Tpetra Matrix/Vector Services Version of the Day
Tpetra_MatrixIO.cpp
00001 #include "Tpetra_MatrixIO.hpp"
00002 
00003 #include <cstdio>
00004 #include <cstdlib>
00005 #include <cstring>
00006 #include <functional>
00007 #include <algorithm>
00008 #include <iterator>
00009 #include <exception>
00010 #include <string>
00011 #include <cctype>
00012 
00013 bool Tpetra::Utils::parseIfmt(Teuchos::ArrayRCP<char> fmt, int &perline, int &width) {
00014   TEST_FOR_EXCEPT(fmt.size() != 0 && fmt[fmt.size()-1] != '\0');
00015   // parses integers n and d out of (nId)
00016   bool error = true;
00017   std::transform(fmt.begin(), fmt.end(), fmt, static_cast < int(*)(int) > (std::toupper));
00018   if (std::sscanf(fmt.getRawPtr(),"(%dI%d)",&perline,&width) == 2) {
00019     error = false;
00020   }
00021   return error;
00022 }
00023 
00024 bool Tpetra::Utils::parseRfmt(Teuchos::ArrayRCP<char> fmt, int &perline, int &width, int &prec, char &valformat) {
00025   TEST_FOR_EXCEPT(fmt.size() != 0 && fmt[fmt.size()-1] != '\0');
00026   std::transform(fmt.begin(), fmt.end(), fmt, static_cast < int(*)(int) > (std::toupper));
00027   // find the first left paren '(' and the last right paren ')'
00028   Teuchos::ArrayRCP<char>::iterator firstLeftParen = std::find( fmt.begin(),  fmt.end(), '(');
00029   Teuchos::ArrayRCP<char>::iterator lastRightParen = std::find(std::reverse_iterator<Teuchos::ArrayRCP<char>::iterator>(fmt.end()), 
00030                                                                std::reverse_iterator<Teuchos::ArrayRCP<char>::iterator>(fmt.begin()), 
00031                                                                ')').base()-1;
00032   // select the substring between the parens, including them
00033   // if neither was found, set the string to empty
00034   if (firstLeftParen == fmt.end() || lastRightParen == fmt.begin()) {
00035     fmt.resize(0 + 1);
00036     fmt[0] = '\0';
00037   }
00038   else {
00039     fmt += (firstLeftParen - fmt.begin());
00040     size_t newLen = lastRightParen - firstLeftParen + 1;
00041     fmt.resize(newLen + 1);
00042     fmt[newLen] = '\0';
00043   }
00044   if (std::find(fmt.begin(),fmt.end(),'P') != fmt.end()) {
00045     // not supported
00046     return true;
00047   }
00048   bool error = true;
00049   if (std::sscanf(fmt.getRawPtr(),"(%d%c%d.%d)",&perline,&valformat,&width,&prec) == 4) {
00050     if (valformat == 'E' || valformat == 'D' || valformat == 'F') {
00051       error = false;
00052     }
00053   } 
00054   return error;
00055 }
00056 
00057 
00058 void Tpetra::Utils::readHBHeader(std::ifstream &fin, Teuchos::ArrayRCP<char> &Title, Teuchos::ArrayRCP<char> &Key, Teuchos::ArrayRCP<char> &Type, 
00059                            int &Nrow, int &Ncol, int &Nnzero, int &Nrhs,
00060                            Teuchos::ArrayRCP<char> &Ptrfmt, Teuchos::ArrayRCP<char> &Indfmt, Teuchos::ArrayRCP<char> &Valfmt, Teuchos::ArrayRCP<char> &Rhsfmt, 
00061                            int &Ptrcrd, int &Indcrd, int &Valcrd, int &Rhscrd, Teuchos::ArrayRCP<char> &Rhstype) {
00062   int Totcrd, Neltvl, Nrhsix;
00063   const int MAXLINE = 81;
00064   char line[MAXLINE];
00065   //
00066   Title.resize(72 + 1);  std::fill(Title.begin(),  Title.end(),  '\0');
00067   Key.resize(8 + 1);     std::fill(Key.begin(),    Key.end(),    '\0');
00068   Type.resize(3 + 1);    std::fill(Type.begin(),   Type.end(),   '\0');
00069   Ptrfmt.resize(16 + 1); std::fill(Ptrfmt.begin(), Ptrfmt.end(), '\0');
00070   Indfmt.resize(16 + 1); std::fill(Indfmt.begin(), Indfmt.end(), '\0');
00071   Valfmt.resize(20 + 1); std::fill(Valfmt.begin(), Valfmt.end(), '\0');
00072   Rhsfmt.resize(20 + 1); std::fill(Rhsfmt.begin(), Rhsfmt.end(), '\0');
00073   //
00074   const std::string errStr("Tpetra::Utils::readHBHeader(): Improperly formatted H/B file: ");
00075   /*  First line:   (A72,A8) */
00076   fin.getline(line,MAXLINE);
00077   TEST_FOR_EXCEPTION( std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00078   (void)std::sscanf(line, "%72c%8[^\n]", Title.getRawPtr(), Key.getRawPtr());
00079   /*  Second line:  (5I14) or (4I14) */
00080   fin.getline(line,MAXLINE);
00081   TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00082   if ( std::sscanf(line,"%14d%14d%14d%14d%14d",&Totcrd,&Ptrcrd,&Indcrd,&Valcrd,&Rhscrd) != 5 ) {
00083     Rhscrd = 0;
00084     TEST_FOR_EXCEPTION(std::sscanf(line,"%14d%14d%14d%14d",&Totcrd,&Ptrcrd,&Indcrd,&Valcrd) != 4, std::runtime_error, errStr << "error reading pointers (line 2)");
00085   }
00086   /*  Third line:   (A3, 11X, 4I14) */
00087   fin.getline(line,MAXLINE);
00088   TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00089   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)");
00090   std::transform(Type.begin(), Type.end(), Type.begin(), static_cast < int(*)(int) > (std::toupper));
00091   /*  Fourth line:  */
00092   fin.getline(line,MAXLINE);
00093   TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00094   if (Rhscrd != 0) {
00095     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)");
00096   }
00097   else {
00098     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)");
00099   }
00100   /*  (Optional) Fifth line: */
00101   if (Rhscrd != 0 ) { 
00102     Rhstype.resize(3 + 1,'\0');
00103     fin.getline(line,MAXLINE);
00104     TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
00105     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)");
00106   }
00107 }
00108 
00109 
00110 void Tpetra::Utils::readHBInfo(const std::string &filename, int &M, int &N, int &nz, Teuchos::ArrayRCP<char> &Type, int &Nrhs) {
00111   std::ifstream fin;
00112   int Ptrcrd, Indcrd, Valcrd, Rhscrd; 
00113   Teuchos::ArrayRCP<char> Title, Key, Rhstype, Ptrfmt, Indfmt, Valfmt, Rhsfmt;
00114   try {
00115     fin.open(filename.c_str(),std::ifstream::in);
00116     Tpetra::Utils::readHBHeader(fin, Title, Key, Type, M, N, nz, Nrhs,
00117                                 Ptrfmt, Indfmt, Valfmt, Rhsfmt, 
00118                                 Ptrcrd, Indcrd, Valcrd, Rhscrd, Rhstype);
00119     fin.close();
00120   }
00121   catch (std::exception &e) {
00122     TEST_FOR_EXCEPTION(true, std::runtime_error, 
00123         "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
00124         << e.what() << std::endl);
00125   }
00126 }
00127 
00128 
00129 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) {
00130   // NOTE: if complex, allocate enough space for 2*NNZ and interleave real and imaginary parts (real,imag)
00131   //       if pattern, don't touch parameter vals; do not allocate space space for values
00132   try {
00133     std::ifstream fin;
00134     int ptrCrd, indCrd, valCrd;
00135     Teuchos::ArrayRCP<char> Title, Key, Ptrfmt, Indfmt, Valfmt;
00136     const int MAXSIZE = 81;
00137     char lineBuf[MAXSIZE];
00138     // nitty gritty
00139     int ptrsPerLine, ptrWidth, indsPerLine, indWidth, valsPerLine, valWidth, valPrec;
00140     char valFlag;
00141     // 
00142     fin.open(filename.c_str(),std::ifstream::in);
00143     {
00144       // we don't care about RHS-related stuff, so declare those vars in an expiring scope
00145       int Nrhs, rhsCrd;
00146       Teuchos::ArrayRCP<char> Rhstype, Rhsfmt;
00147       Teuchos::ArrayRCP<char> TypeArray;
00148       Tpetra::Utils::readHBHeader(fin, Title, Key, TypeArray, numRows, numCols, numNZ, Nrhs,
00149                                   Ptrfmt, Indfmt, Valfmt, Rhsfmt, 
00150                                   ptrCrd, indCrd, valCrd, rhsCrd, Rhstype);
00151       if (TypeArray.size() > 0) {
00152         type.resize(TypeArray.size()-1);
00153         std::copy(TypeArray.begin(), TypeArray.end(), type.begin());
00154       }
00155     }
00156     const std::string errStr("Tpetra::Utils::readHBHeader(): Improperly formatted H/B file.");
00157     const bool readPatternOnly = (type[0] == 'P' || type[0] == 'p');
00158     const bool readComplex     = (type[0] == 'C' || type[0] == 'c');
00159     /*  Parse the array input formats from Line 3 of HB file  */
00160     TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Ptrfmt,ptrsPerLine,ptrWidth) == true, std::runtime_error,
00161         "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
00162     TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Indfmt,indsPerLine,indWidth) == true, std::runtime_error,
00163         "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
00164     if (readPatternOnly == false) {
00165       TEST_FOR_EXCEPTION( Tpetra::Utils::parseRfmt(Valfmt,valsPerLine,valWidth,valPrec,valFlag) == true, std::runtime_error,
00166           "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
00167     }
00168     // special case this: the reason is that the number of colPtrs read is numCols+1, which is non-zero even if numCols == 0
00169     // furthermore, if numCols == 0, there is nothing of interest to read
00170     if (numCols == 0) return;
00171     // allocate space for column pointers, row indices, and values
00172     // if the file is empty, do not touch these ARCPs
00173     colPtrs = Teuchos::arcp<int>(numCols+1);
00174     if (numNZ > 0) {
00175       rowInds = Teuchos::arcp<int>(numNZ);  
00176       if (readPatternOnly == false) {
00177         if (readComplex) {
00178           vals = Teuchos::arcp<double>(2*numNZ); 
00179         }
00180         else {
00181           vals = Teuchos::arcp<double>(numNZ); 
00182         }
00183       }
00184     }
00185     /* Read column pointer array:
00186        Specifically, read ptrCrd number of lines, and on each line, read ptrsPerLine number of integers, each of width ptrWidth
00187        Store these in colPtrs */
00188     {
00189       int colPtrsRead = 0;
00190       char NullSub = '\0';
00191       for (int lno=0; lno < ptrCrd; ++lno) {
00192         fin.getline(lineBuf, MAXSIZE);
00193         TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
00194         char *linePtr = lineBuf;
00195         for (int ptr=0; ptr < ptrsPerLine; ++ptr) {
00196           if (colPtrsRead == numCols + 1) break;
00197           int cptr;
00198           // terminate the string at the end of the current ptr block, saving the character in that location
00199           std::swap(NullSub,linePtr[ptrWidth]);
00200           // read the ptr
00201           std::sscanf(linePtr, "%d", &cptr);
00202           // put the saved character back, and put the '\0' back into NullSub for use again
00203           std::swap(NullSub,linePtr[ptrWidth]);
00204           linePtr += ptrWidth;
00205           colPtrs[colPtrsRead++] = cptr;
00206         }
00207       }
00208       TEST_FOR_EXCEPT(colPtrsRead != numCols + 1);
00209     }
00210     /* Read row index array:
00211        Specifically, read indCrd number of lines, and on each line, read indsPerLine number of integers, each of width indWidth
00212        Store these in rowInds */
00213     {
00214       char NullSub = '\0';
00215       int indicesRead = 0;
00216       for (int lno=0; lno < indCrd; ++lno) {
00217         fin.getline(lineBuf, MAXSIZE);
00218         TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
00219         char *linePtr = lineBuf;
00220         for (int indcntr=0; indcntr < indsPerLine; ++indcntr) {
00221           if (indicesRead == numNZ) break;
00222           int ind;
00223           // terminate the string at the end of the current ind block, saving the character in that location
00224           std::swap(NullSub,linePtr[indWidth]);
00225           // read the ind
00226           std::sscanf(linePtr, "%d", &ind);
00227           // put the saved character back, and put the '\0' back into NullSub for use again
00228           std::swap(NullSub,linePtr[indWidth]);
00229           linePtr += indWidth;
00230           rowInds[indicesRead++] = ind;
00231         }
00232       }
00233       TEST_FOR_EXCEPT(indicesRead != numNZ);
00234     }
00235     /* Read array of values:
00236        Specifically, read valCrd number of lines, and on each line, read valsPerLine number of real values, each of width/precision valWidth/valPrec
00237        Store these in vals
00238        If readComplex, then read twice as many non-zeros, and interleave real,imag into vals */
00239     if (readPatternOnly == false) {
00240       int totalNumVals;
00241       if (readComplex) {
00242         totalNumVals = 2*numNZ;
00243       }
00244       else {
00245         totalNumVals = numNZ;
00246       }
00247       char NullSub = '\0';
00248       int valsRead = 0;
00249       for (int lno=0; lno < valCrd; ++lno) {
00250         fin.getline(lineBuf, MAXSIZE);
00251         TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
00252         // if valFlag == 'D', then we need to convert [dD] in fp vals into [eE] that scanf can parse
00253         if (valFlag == 'D') std::replace_if(lineBuf, lineBuf+MAXSIZE, std::bind2nd(std::equal_to<char>(), 'D'), 'E'); 
00254         char *linePtr = lineBuf;
00255         for (int valcntr=0; valcntr < valsPerLine; ++valcntr) {
00256           if (valsRead == totalNumVals) break;
00257           double val;
00258           // terminate the string at the end of the current val block, saving the character in that location
00259           std::swap(NullSub,linePtr[valWidth]);
00260           // read the val
00261           std::sscanf(linePtr, "%le", &val);
00262           // put the saved character back, and put the '\0' back into NullSub for use again
00263           std::swap(NullSub,linePtr[valWidth]);
00264           linePtr += valWidth;
00265           vals[valsRead++] = val;
00266         }
00267       }
00268       TEST_FOR_EXCEPT(valsRead != totalNumVals);
00269     }
00270     fin.close();
00271   }
00272   catch (std::exception &e) {
00273     TEST_FOR_EXCEPTION(true, std::runtime_error, 
00274         "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
00275         << e.what() << std::endl);
00276   }
00277 }
00278 
00279 #ifdef HAVE_TPETRA_EXPLICIT_INSTANTIATION
00280 
00281 #include "Tpetra_MatrixIO_def.hpp"
00282 
00283 #include <Kokkos_ConfigDefs.hpp>
00284 #include <Kokkos_SerialNode.hpp>
00285 #ifdef HAVE_KOKKOS_TBB
00286 #  include <Kokkos_TBBNode.hpp>
00287 #endif
00288 #ifdef HAVE_KOKKOS_THREADPOOL
00289 #  include <Kokkos_TPINode.hpp>
00290 #endif
00291 #ifdef HAVE_KOKKOS_THRUST
00292 #  include <Kokkos_ThrustGPUNode.hpp>
00293 #endif
00294 
00295 namespace Tpetra {
00296   namespace Utils {
00297 
00298 #if defined(HAVE_TPETRA_INST_FLOAT)
00299   TPETRA_MATRIXIO_INSTANT(float,int,int,Kokkos::SerialNode)
00300 # ifdef HAVE_KOKKOS_TBB
00301     TPETRA_MATRIXIO_INSTANT(float,int,int,Kokkos::TBBNode)
00302 # endif
00303 # ifdef HAVE_KOKKOS_THREADPOOL
00304     TPETRA_MATRIXIO_INSTANT(float,int,int,Kokkos::TPINode)
00305 # endif
00306 # if defined(HAVE_KOKKOS_THRUST) && defined(HAVE_KOKKOS_CUDA_FLOAT)
00307     TPETRA_MATRIXIO_INSTANT(float,int,int,Kokkos::ThrustGPUNode)
00308 # endif
00309 #endif
00310 
00311 #if defined(HAVE_TPETRA_INST_DOUBLE)
00312   TPETRA_MATRIXIO_INSTANT(double,int,int,Kokkos::SerialNode)
00313 # ifdef HAVE_KOKKOS_TBB
00314     TPETRA_MATRIXIO_INSTANT(double,int,int,Kokkos::TBBNode)
00315 # endif
00316 # ifdef HAVE_KOKKOS_THREADPOOL
00317     TPETRA_MATRIXIO_INSTANT(double,int,int,Kokkos::TPINode)
00318 # endif
00319 # if defined(HAVE_KOKKOS_THRUST) && defined(HAVE_KOKKOS_CUDA_DOUBLE)
00320     TPETRA_MATRIXIO_INSTANT(double,int,int,Kokkos::ThrustGPUNode)
00321 # endif
00322 #endif
00323 
00324 } // namespace Tpetra::Utils
00325 } // namespace Tpetra
00326 
00327 #endif // HAVE_TPETRA_EXPLICIT_INSTANTIATION
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines