Tpetra Matrix/Vector Services Version of the Day
Tpetra_MatrixIO.cpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //          Tpetra: Templated Linear Algebra Services Package
00006 //                 Copyright (2008) Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 // @HEADER
00042 */
00043 
00044 #include "Tpetra_MatrixIO.hpp"
00045 
00046 #include <cstdio>
00047 #include <cstdlib>
00048 #include <cstring>
00049 #include <functional>
00050 #include <algorithm>
00051 #include <iterator>
00052 #include <exception>
00053 #include <string>
00054 #include <cctype>
00055 
00056 bool Tpetra::Utils::parseIfmt(Teuchos::ArrayRCP<char> fmt, int &perline, int &width) {
00057   TEUCHOS_TEST_FOR_EXCEPT(fmt.size() != 0 && fmt[fmt.size()-1] != '\0');
00058   // parses integers n and d out of (nId)
00059   bool error = true;
00060   std::transform(fmt.begin(), fmt.end(), fmt, static_cast < int(*)(int) > (std::toupper));
00061   if (std::sscanf(fmt.getRawPtr(),"(%dI%d)",&perline,&width) == 2) {
00062     error = false;
00063   }
00064   return error;
00065 }
00066 
00067 bool Tpetra::Utils::parseRfmt(Teuchos::ArrayRCP<char> fmt, int &perline, int &width, int &prec, char &valformat) {
00068   TEUCHOS_TEST_FOR_EXCEPT(fmt.size() != 0 && fmt[fmt.size()-1] != '\0');
00069   std::transform(fmt.begin(), fmt.end(), fmt, static_cast < int(*)(int) > (std::toupper));
00070   // find the first left paren '(' and the last right paren ')'
00071   Teuchos::ArrayRCP<char>::iterator firstLeftParen = std::find( fmt.begin(),  fmt.end(), '(');
00072   Teuchos::ArrayRCP<char>::iterator lastRightParen = std::find(std::reverse_iterator<Teuchos::ArrayRCP<char>::iterator>(fmt.end()), 
00073                                                                std::reverse_iterator<Teuchos::ArrayRCP<char>::iterator>(fmt.begin()), 
00074                                                                ')').base()-1;
00075   // select the substring between the parens, including them
00076   // if neither was found, set the string to empty
00077   if (firstLeftParen == fmt.end() || lastRightParen == fmt.begin()) {
00078     fmt.resize(0 + 1);
00079     fmt[0] = '\0';
00080   }
00081   else {
00082     fmt += (firstLeftParen - fmt.begin());
00083     size_t newLen = lastRightParen - firstLeftParen + 1;
00084     fmt.resize(newLen + 1);
00085     fmt[newLen] = '\0';
00086   }
00087   if (std::find(fmt.begin(),fmt.end(),'P') != fmt.end()) {
00088     // not supported
00089     return true;
00090   }
00091   bool error = true;
00092   if (std::sscanf(fmt.getRawPtr(),"(%d%c%d.%d)",&perline,&valformat,&width,&prec) == 4) {
00093     if (valformat == 'E' || valformat == 'D' || valformat == 'F') {
00094       error = false;
00095     }
00096   } 
00097   return error;
00098 }
00099 
00100 
00101 void Tpetra::Utils::readHBHeader(std::ifstream &fin, Teuchos::ArrayRCP<char> &Title, Teuchos::ArrayRCP<char> &Key, Teuchos::ArrayRCP<char> &Type, 
00102                            int &Nrow, int &Ncol, int &Nnzero, int &Nrhs,
00103                            Teuchos::ArrayRCP<char> &Ptrfmt, Teuchos::ArrayRCP<char> &Indfmt, Teuchos::ArrayRCP<char> &Valfmt, Teuchos::ArrayRCP<char> &Rhsfmt, 
00104                            int &Ptrcrd, int &Indcrd, int &Valcrd, int &Rhscrd, Teuchos::ArrayRCP<char> &Rhstype) {
00105   int Totcrd, Neltvl, Nrhsix;
00106   const int MAXLINE = 81;
00107   char line[MAXLINE];
00108   //
00109   Title.resize(72 + 1);  std::fill(Title.begin(),  Title.end(),  '\0');
00110   Key.resize(8 + 1);     std::fill(Key.begin(),    Key.end(),    '\0');
00111   Type.resize(3 + 1);    std::fill(Type.begin(),   Type.end(),   '\0');
00112   Ptrfmt.resize(16 + 1); std::fill(Ptrfmt.begin(), Ptrfmt.end(), '\0');
00113   Indfmt.resize(16 + 1); std::fill(Indfmt.begin(), Indfmt.end(), '\0');
00114   Valfmt.resize(20 + 1); std::fill(Valfmt.begin(), Valfmt.end(), '\0');
00115   Rhsfmt.resize(20 + 1); std::fill(Rhsfmt.begin(), Rhsfmt.end(), '\0');
00116   //
00117   const std::string errStr("Tpetra::Utils::readHBHeader(): Improperly formatted H/B file: ");
00118   /*  First line:   (A72,A8) */
00119   fin.getline(line,MAXLINE);
00120   TEUCHOS_TEST_FOR_EXCEPTION( std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00121   (void)std::sscanf(line, "%72c%8[^\n]", Title.getRawPtr(), Key.getRawPtr());
00122   /*  Second line:  (5I14) or (4I14) */
00123   fin.getline(line,MAXLINE);
00124   TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00125   if ( std::sscanf(line,"%14d%14d%14d%14d%14d",&Totcrd,&Ptrcrd,&Indcrd,&Valcrd,&Rhscrd) != 5 ) {
00126     Rhscrd = 0;
00127     TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%14d%14d%14d%14d",&Totcrd,&Ptrcrd,&Indcrd,&Valcrd) != 4, std::runtime_error, errStr << "error reading pointers (line 2)");
00128   }
00129   /*  Third line:   (A3, 11X, 4I14) */
00130   fin.getline(line,MAXLINE);
00131   TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00132   TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line, "%3c%14i%14i%14i%14i", Type.getRawPtr(),&Nrow,&Ncol,&Nnzero,&Neltvl) != 5 , std::runtime_error, errStr << "error reading matrix meta-data (line 3)");
00133   std::transform(Type.begin(), Type.end(), Type.begin(), static_cast < int(*)(int) > (std::toupper));
00134   /*  Fourth line:  */
00135   fin.getline(line,MAXLINE);
00136   TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00137   if (Rhscrd != 0) {
00138     TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%16c%16c%20c%20c",Ptrfmt.getRawPtr(),Indfmt.getRawPtr(),Valfmt.getRawPtr(),Rhsfmt.getRawPtr()) != 4, std::runtime_error, errStr << "error reading formats (line 4)");
00139   }
00140   else {
00141     TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%16c%16c%20c",Ptrfmt.getRawPtr(),Indfmt.getRawPtr(),Valfmt.getRawPtr()) != 3,                        std::runtime_error, errStr << "error reading formats (line 4)");
00142   }
00143   /*  (Optional) Fifth line: */
00144   if (Rhscrd != 0 ) { 
00145     Rhstype.resize(3 + 1,'\0');
00146     fin.getline(line,MAXLINE);
00147     TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00148     TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%3c%14d%14d", Rhstype.getRawPtr(), &Nrhs, &Nrhsix) != 3, std::runtime_error, errStr << "error reading right-hand-side meta-data (line 5)");
00149   }
00150 }
00151 
00152 
00153 void Tpetra::Utils::readHBInfo(const std::string &filename, int &M, int &N, int &nz, Teuchos::ArrayRCP<char> &Type, int &Nrhs) {
00154   std::ifstream fin;
00155   int Ptrcrd, Indcrd, Valcrd, Rhscrd; 
00156   Teuchos::ArrayRCP<char> Title, Key, Rhstype, Ptrfmt, Indfmt, Valfmt, Rhsfmt;
00157   try {
00158     fin.open(filename.c_str(),std::ifstream::in);
00159     Tpetra::Utils::readHBHeader(fin, Title, Key, Type, M, N, nz, Nrhs,
00160                                 Ptrfmt, Indfmt, Valfmt, Rhsfmt, 
00161                                 Ptrcrd, Indcrd, Valcrd, Rhscrd, Rhstype);
00162     fin.close();
00163   }
00164   catch (std::exception &e) {
00165     TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, 
00166         "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
00167         << e.what() << std::endl);
00168   }
00169 }
00170 
00171 
00172 void Tpetra::Utils::readHBMatDouble(const std::string &filename, int &numRows, int &numCols, int &numNZ, std::string &type, Teuchos::ArrayRCP<int> &colPtrs, Teuchos::ArrayRCP<int> &rowInds, Teuchos::ArrayRCP<double> &vals) {
00173   // NOTE: if complex, allocate enough space for 2*NNZ and interleave real and imaginary parts (real,imag)
00174   //       if pattern, don't touch parameter vals; do not allocate space space for values
00175   try {
00176     std::ifstream fin;
00177     int ptrCrd, indCrd, valCrd;
00178     Teuchos::ArrayRCP<char> Title, Key, Ptrfmt, Indfmt, Valfmt;
00179     const int MAXSIZE = 81;
00180     char lineBuf[MAXSIZE];
00181     // nitty gritty
00182     int ptrsPerLine, ptrWidth, indsPerLine, indWidth, valsPerLine, valWidth, valPrec;
00183     char valFlag;
00184     // 
00185     fin.open(filename.c_str(),std::ifstream::in);
00186     {
00187       // we don't care about RHS-related stuff, so declare those vars in an expiring scope
00188       int Nrhs, rhsCrd;
00189       Teuchos::ArrayRCP<char> Rhstype, Rhsfmt;
00190       Teuchos::ArrayRCP<char> TypeArray;
00191       Tpetra::Utils::readHBHeader(fin, Title, Key, TypeArray, numRows, numCols, numNZ, Nrhs,
00192                                   Ptrfmt, Indfmt, Valfmt, Rhsfmt, 
00193                                   ptrCrd, indCrd, valCrd, rhsCrd, Rhstype);
00194       if (TypeArray.size() > 0) {
00195         type.resize(TypeArray.size()-1);
00196         std::copy(TypeArray.begin(), TypeArray.end(), type.begin());
00197       }
00198     }
00199     const std::string errStr("Tpetra::Utils::readHBHeader(): Improperly formatted H/B file.");
00200     const bool readPatternOnly = (type[0] == 'P' || type[0] == 'p');
00201     const bool readComplex     = (type[0] == 'C' || type[0] == 'c');
00202     /*  Parse the array input formats from Line 3 of HB file  */
00203     TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Ptrfmt,ptrsPerLine,ptrWidth) == true, std::runtime_error,
00204         "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
00205     TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Indfmt,indsPerLine,indWidth) == true, std::runtime_error,
00206         "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
00207     if (readPatternOnly == false) {
00208       TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseRfmt(Valfmt,valsPerLine,valWidth,valPrec,valFlag) == true, std::runtime_error,
00209           "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
00210     }
00211     // special case this: the reason is that the number of colPtrs read is numCols+1, which is non-zero even if numCols == 0
00212     // furthermore, if numCols == 0, there is nothing of interest to read
00213     if (numCols == 0) return;
00214     // allocate space for column pointers, row indices, and values
00215     // if the file is empty, do not touch these ARCPs
00216     colPtrs = Teuchos::arcp<int>(numCols+1);
00217     if (numNZ > 0) {
00218       rowInds = Teuchos::arcp<int>(numNZ);  
00219       if (readPatternOnly == false) {
00220         if (readComplex) {
00221           vals = Teuchos::arcp<double>(2*numNZ); 
00222         }
00223         else {
00224           vals = Teuchos::arcp<double>(numNZ); 
00225         }
00226       }
00227     }
00228     /* Read column pointer array:
00229        Specifically, read ptrCrd number of lines, and on each line, read ptrsPerLine number of integers, each of width ptrWidth
00230        Store these in colPtrs */
00231     {
00232       int colPtrsRead = 0;
00233       char NullSub = '\0';
00234       for (int lno=0; lno < ptrCrd; ++lno) {
00235         fin.getline(lineBuf, MAXSIZE);
00236         TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
00237         char *linePtr = lineBuf;
00238         for (int ptr=0; ptr < ptrsPerLine; ++ptr) {
00239           if (colPtrsRead == numCols + 1) break;
00240           int cptr;
00241           // terminate the string at the end of the current ptr block, saving the character in that location
00242           std::swap(NullSub,linePtr[ptrWidth]);
00243           // read the ptr
00244           std::sscanf(linePtr, "%d", &cptr);
00245           // put the saved character back, and put the '\0' back into NullSub for use again
00246           std::swap(NullSub,linePtr[ptrWidth]);
00247           linePtr += ptrWidth;
00248           colPtrs[colPtrsRead++] = cptr;
00249         }
00250       }
00251       TEUCHOS_TEST_FOR_EXCEPT(colPtrsRead != numCols + 1);
00252     }
00253     /* Read row index array:
00254        Specifically, read indCrd number of lines, and on each line, read indsPerLine number of integers, each of width indWidth
00255        Store these in rowInds */
00256     {
00257       char NullSub = '\0';
00258       int indicesRead = 0;
00259       for (int lno=0; lno < indCrd; ++lno) {
00260         fin.getline(lineBuf, MAXSIZE);
00261         TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
00262         char *linePtr = lineBuf;
00263         for (int indcntr=0; indcntr < indsPerLine; ++indcntr) {
00264           if (indicesRead == numNZ) break;
00265           int ind;
00266           // terminate the string at the end of the current ind block, saving the character in that location
00267           std::swap(NullSub,linePtr[indWidth]);
00268           // read the ind
00269           std::sscanf(linePtr, "%d", &ind);
00270           // put the saved character back, and put the '\0' back into NullSub for use again
00271           std::swap(NullSub,linePtr[indWidth]);
00272           linePtr += indWidth;
00273           rowInds[indicesRead++] = ind;
00274         }
00275       }
00276       TEUCHOS_TEST_FOR_EXCEPT(indicesRead != numNZ);
00277     }
00278     /* Read array of values:
00279        Specifically, read valCrd number of lines, and on each line, read valsPerLine number of real values, each of width/precision valWidth/valPrec
00280        Store these in vals
00281        If readComplex, then read twice as many non-zeros, and interleave real,imag into vals */
00282     if (readPatternOnly == false) {
00283       int totalNumVals;
00284       if (readComplex) {
00285         totalNumVals = 2*numNZ;
00286       }
00287       else {
00288         totalNumVals = numNZ;
00289       }
00290       char NullSub = '\0';
00291       int valsRead = 0;
00292       for (int lno=0; lno < valCrd; ++lno) {
00293         fin.getline(lineBuf, MAXSIZE);
00294         TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
00295         // if valFlag == 'D', then we need to convert [dD] in fp vals into [eE] that scanf can parse
00296         if (valFlag == 'D') std::replace_if(lineBuf, lineBuf+MAXSIZE, std::bind2nd(std::equal_to<char>(), 'D'), 'E'); 
00297         char *linePtr = lineBuf;
00298         for (int valcntr=0; valcntr < valsPerLine; ++valcntr) {
00299           if (valsRead == totalNumVals) break;
00300           double val;
00301           // terminate the string at the end of the current val block, saving the character in that location
00302           std::swap(NullSub,linePtr[valWidth]);
00303           // read the val
00304           std::sscanf(linePtr, "%le", &val);
00305           // put the saved character back, and put the '\0' back into NullSub for use again
00306           std::swap(NullSub,linePtr[valWidth]);
00307           linePtr += valWidth;
00308           vals[valsRead++] = val;
00309         }
00310       }
00311       TEUCHOS_TEST_FOR_EXCEPT(valsRead != totalNumVals);
00312     }
00313     fin.close();
00314   }
00315   catch (std::exception &e) {
00316     TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, 
00317         "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
00318         << e.what() << std::endl);
00319   }
00320 }
00321 
00322 #ifdef HAVE_TPETRA_EXPLICIT_INSTANTIATION
00323 
00324 #include "Tpetra_ETIHelperMacros.h"
00325 #include "Tpetra_MatrixIO_def.hpp"
00326 
00327 namespace Tpetra {
00328   namespace Utils {
00329 
00330     TPETRA_ETI_MANGLING_TYPEDEFS()
00331 
00332     TPETRA_INSTANTIATE_SLGN(TPETRA_MATRIXIO_INSTANT)
00333 
00334   } // namespace Tpetra::Utils
00335 } // namespace Tpetra
00336 
00337 #endif // HAVE_TPETRA_EXPLICIT_INSTANTIATION
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines