EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_CrsMatrixIn.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 #include "Epetra_ConfigDefs.h"
00042 #include "EpetraExt_CrsMatrixIn.h"
00043 #include "Epetra_Comm.h"
00044 #include "Epetra_CrsMatrix.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_IntVector.h"
00047 #include "Epetra_IntSerialDenseVector.h"
00048 #include "Epetra_Import.h"
00049 #include "Epetra_Time.h"
00050 #include "Epetra_Util.h"
00051 
00052 #if defined(__PGI)
00053 #include <sstream>
00054 #endif
00055 
00058 // Functions to read a MatrixMarket file and load it into a matrix.
00059 // Adapted from 
00060 //    Trilinos/packages/epetraext/src/inout/EpetraExt_CrsMatrixIn.cpp
00061 // Modified by Jon Berry and Karen Devine to make matrix reallocations 
00062 // more efficient; rather than insert each non-zero separately, we 
00063 // collect rows of non-zeros for insertion.
00064 // Modified by Karen Devine and Steve Plimpton to prevent all processors 
00065 // from reading the same file at the same time (ouch!); chunks of the file 
00066 // are read and broadcast by processor zero; each processor then extracts 
00067 // its nonzeros from the chunk, sorts them by row, and inserts them into 
00068 // the matrix.
00069 // The variable "chunk_read" can be changed to modify the size of the chunks
00070 // read from the file.  Larger values of chunk_read lead to faster execution
00071 // and greater memory use.
00074 
00075 using namespace EpetraExt;
00076 namespace EpetraExt {
00077 
00078 template<typename int_type>
00079 static void sort_three(int_type* list, int_type *parlista, double *parlistb, 
00080                        int start, int end);
00081 
00083 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00084 int MatrixMarketFileToCrsMatrix(const char *filename, 
00085     const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
00086 {
00087   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A));
00088   return(0);
00089 }
00090 
00092 int MatrixMarketFileToCrsMatrix(const char *filename, 
00093     const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
00094 {
00095   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
00096                                                    0, 0, 0, 0, transpose));
00097   return(0);
00098 }
00100   
00101 int MatrixMarketFileToCrsMatrix(const char *filename, 
00102     const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
00103     const bool verbose) 
00104 {
00105   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
00106                                                    0, 0, 0, 0,
00107                                                    transpose, verbose));
00108   return(0);
00109 }
00110   
00112 int MatrixMarketFileToCrsMatrix(const char *filename, 
00113     const Epetra_Map & rowMap, const Epetra_Map& rangeMap, 
00114     const Epetra_Map& domainMap, Epetra_CrsMatrix * & A, const bool transpose,
00115     const bool verbose) 
00116 {
00117   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00118                                                    rowMap.Comm(), A,
00119                                                    &rowMap, 0, 
00120                                                    &rangeMap, &domainMap,
00121                                                    transpose, verbose));
00122   return(0);
00123 }
00124 
00126 int MatrixMarketFileToCrsMatrix(const char *filename, 
00127     const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
00128     const bool verbose) 
00129 {
00130   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00131                                                    rowMap.Comm(), A,
00132                                                    &rowMap, 0, 0, 0, 
00133                                                    transpose, verbose));
00134   return(0);
00135 }
00136 
00138 int MatrixMarketFileToCrsMatrix(const char *filename, 
00139     const Epetra_Map & rowMap, const Epetra_Map & colMap, 
00140     Epetra_CrsMatrix * & A, const bool transpose,
00141     const bool verbose) 
00142 {
00143   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00144                                                    rowMap.Comm(), A,
00145                                                    &rowMap, &colMap, 0, 0,
00146                                                    transpose, verbose));
00147   return(0);
00148 }
00149 
00151 int MatrixMarketFileToCrsMatrix(const char *filename, 
00152     const Epetra_Map & rowMap, const Epetra_Map & colMap,
00153     const Epetra_Map& rangeMap, const Epetra_Map& domainMap, 
00154     Epetra_CrsMatrix * & A, const bool transpose,
00155     const bool verbose) 
00156 {
00157   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00158                                                    rowMap.Comm(), A, 
00159                                                    &rowMap, &colMap, 
00160                                                    &rangeMap, &domainMap,
00161                                                    transpose, verbose));
00162   return(0);
00163 }
00164 #endif
00165 
00166 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00167 int MatrixMarketFileToCrsMatrix64(const char *filename, 
00168     const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
00169 {
00170   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A));
00171   return(0);
00172 }
00173 
00175 int MatrixMarketFileToCrsMatrix64(const char *filename, 
00176     const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
00177 {
00178   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A,
00179                                                    0, 0, 0, 0, transpose));
00180   return(0);
00181 }
00183   
00184 int MatrixMarketFileToCrsMatrix64(const char *filename, 
00185     const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
00186     const bool verbose) 
00187 {
00188   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A,
00189                                                    0, 0, 0, 0,
00190                                                    transpose, verbose));
00191   return(0);
00192 }
00193   
00195 int MatrixMarketFileToCrsMatrix64(const char *filename, 
00196     const Epetra_Map & rowMap, const Epetra_Map& rangeMap, 
00197     const Epetra_Map& domainMap, Epetra_CrsMatrix * & A, const bool transpose,
00198     const bool verbose) 
00199 {
00200   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, 
00201                                                    rowMap.Comm(), A,
00202                                                    &rowMap, 0, 
00203                                                    &rangeMap, &domainMap,
00204                                                    transpose, verbose));
00205   return(0);
00206 }
00207 
00209 int MatrixMarketFileToCrsMatrix64(const char *filename, 
00210     const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
00211     const bool verbose) 
00212 {
00213   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, 
00214                                                    rowMap.Comm(), A,
00215                                                    &rowMap, 0, 0, 0, 
00216                                                    transpose, verbose));
00217   return(0);
00218 }
00219 
00221 int MatrixMarketFileToCrsMatrix64(const char *filename, 
00222     const Epetra_Map & rowMap, const Epetra_Map & colMap, 
00223     Epetra_CrsMatrix * & A, const bool transpose,
00224     const bool verbose) 
00225 {
00226   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, 
00227                                                    rowMap.Comm(), A,
00228                                                    &rowMap, &colMap, 0, 0,
00229                                                    transpose, verbose));
00230   return(0);
00231 }
00232 
00234 int MatrixMarketFileToCrsMatrix64(const char *filename, 
00235     const Epetra_Map & rowMap, const Epetra_Map & colMap,
00236     const Epetra_Map& rangeMap, const Epetra_Map& domainMap, 
00237     Epetra_CrsMatrix * & A, const bool transpose,
00238     const bool verbose) 
00239 {
00240   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, 
00241                                                    rowMap.Comm(), A, 
00242                                                    &rowMap, &colMap, 
00243                                                    &rangeMap, &domainMap,
00244                                                    transpose, verbose));
00245   return(0);
00246 }
00247 #endif
00248 
00250 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00251 int MatrixMarketFileToCrsMatrixHandle(const char *filename,
00252   const Epetra_Comm & comm,
00253   Epetra_CrsMatrix * & A,
00254   const Epetra_Map * rowMap,
00255   const Epetra_Map * colMap,
00256   const Epetra_Map * rangeMap,
00257   const Epetra_Map * domainMap,
00258   const bool transpose,
00259   const bool verbose
00260 )
00261 {
00262   const int chunk_read = 500000;  //  Modify this variable to change the
00263                                   //  size of the chunks read from the file.
00264   const int headerlineLength = 257;
00265   const int lineLength = 81;
00266   const int tokenLength = 35;
00267   char line[headerlineLength];
00268   char token1[tokenLength];
00269   char token2[tokenLength];
00270   char token3[tokenLength];
00271   char token4[tokenLength];
00272   char token5[tokenLength];
00273   int M, N, NZ;      // Matrix dimensions
00274   int i;
00275   int me = comm.MyPID();
00276 
00277   Epetra_Time timer(comm);
00278 
00279   // Make sure domain and range maps are either both defined or undefined 
00280   if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
00281     EPETRA_CHK_ERR(-3);
00282   }
00283 
00284   // check maps to see if domain and range are 1-to-1
00285 
00286   if (domainMap!=0) {
00287     if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00288     if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00289   }
00290   else {
00291     // If domain and range are not specified, 
00292     // rowMap becomes both and must be 1-to-1 if specified
00293     if (rowMap!=0) {
00294       if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00295     }
00296   }
00297       
00298   FILE * handle = 0;
00299   if (me == 0) {
00300     if (verbose) std::cout << "Reading MatrixMarket file " << filename << std::endl;
00301     handle = fopen(filename,"r");  // Open file
00302     if (handle == 0)
00303       EPETRA_CHK_ERR(-1); // file not found
00304 
00305     // Check first line, which should be 
00306     // %%MatrixMarket matrix coordinate real general
00307     if(fgets(line, headerlineLength, handle)==0) {
00308       if (handle!=0) fclose(handle); 
00309       EPETRA_CHK_ERR(-1);
00310     }
00311     if(sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
00312       if (handle!=0) fclose(handle); 
00313       EPETRA_CHK_ERR(-1);
00314     }
00315     if (strcmp(token1, "%%MatrixMarket") ||
00316         strcmp(token2, "matrix") ||
00317         strcmp(token3, "coordinate") ||
00318         strcmp(token4, "real") ||
00319         strcmp(token5, "general")) {
00320       if (handle!=0) fclose(handle); 
00321       EPETRA_CHK_ERR(-1);
00322     }
00323 
00324     // Next, strip off header lines (which start with "%")
00325     do {
00326       if(fgets(line, headerlineLength, handle)==0) {
00327         if (handle!=0) fclose(handle); 
00328         EPETRA_CHK_ERR(-1);
00329       }
00330     } while (line[0] == '%');
00331 
00332     // Next get problem dimensions: M, N, NZ
00333     if(sscanf(line, "%d %d %d", &M, &N, &NZ)==0) {
00334       if (handle!=0) fclose(handle); 
00335       EPETRA_CHK_ERR(-1);
00336     }
00337   }
00338   comm.Broadcast(&M, 1, 0);
00339   comm.Broadcast(&N, 1, 0);
00340   comm.Broadcast(&NZ, 1, 0);
00341 
00342   // Now create matrix using user maps if provided.
00343 
00344 
00345   // Now read in chunks of triplets and broadcast them to other processors.
00346   char *buffer = new char[chunk_read*lineLength];
00347   int nchunk; 
00348   int nmillion = 0;
00349   int nread = 0;
00350   int rlen;
00351 
00352   // Storage for this processor's nonzeros.
00353   const int localblock = 100000;
00354   int localsize = NZ / comm.NumProc() + localblock;
00355   int *iv = (int *) malloc(localsize * sizeof(int));
00356   int *jv = (int *) malloc(localsize * sizeof(int));
00357   double *vv = (double *) malloc(localsize * sizeof(double));
00358   int lnz = 0;   //  Number of non-zeros on this processor.
00359 
00360   if (!iv || !jv || !vv) 
00361     EPETRA_CHK_ERR(-1);
00362 
00363   Epetra_Map *rowMap1;
00364   bool allocatedHere=false;
00365   if (rowMap != 0) 
00366     rowMap1 = const_cast<Epetra_Map *>(rowMap);
00367   else {
00368     rowMap1 = new Epetra_Map(M, 0, comm);
00369     allocatedHere = true;
00370   }
00371   int ioffset = rowMap1->IndexBase()-1;
00372   int joffset = (colMap != 0 ? colMap->IndexBase()-1 : ioffset);
00373 
00374   int rowmajor = 1;  // Assume non-zeros are listed in row-major order, 
00375   int prevrow = -1;  // but test to detect otherwise.  If non-zeros
00376                      // are row-major, we can avoid the sort.
00377 
00378   while (nread < NZ) {
00379     if (NZ-nread > chunk_read) nchunk = chunk_read;
00380     else nchunk = NZ - nread;
00381 
00382     if (me == 0) {
00383       char *eof;
00384       rlen = 0;
00385       for (int i = 0; i < nchunk; i++) {
00386         eof = fgets(&buffer[rlen],lineLength,handle);
00387         if (eof == NULL) {
00388           fprintf(stderr, "%s", "Unexpected end of matrix file.");
00389           EPETRA_CHK_ERR(-1);
00390         }
00391         rlen += strlen(&buffer[rlen]);
00392       }
00393       buffer[rlen++]='\n';
00394     }
00395     comm.Broadcast(&rlen, 1, 0);
00396     comm.Broadcast(buffer, rlen, 0);
00397 
00398     buffer[rlen++] = '\0';
00399     nread += nchunk;
00400 
00401     // Process the received data, saving non-zeros belonging on this proc.
00402     char *lineptr = buffer;
00403     for (rlen = 0; rlen < nchunk; rlen++) {
00404       char *next = strchr(lineptr,'\n');
00405       int I = atoi(strtok(lineptr," \t\n")) + ioffset;
00406       int J = atoi(strtok(NULL," \t\n")) + joffset;
00407       double V = atof(strtok(NULL," \t\n"));
00408       lineptr = next + 1;
00409       if (transpose) {
00410         // swap I and J indices.
00411         int tmp = I;
00412         I = J;
00413         J = tmp;
00414       }
00415       if (rowMap1->MyGID(I)) {
00416         //  This processor keeps this non-zero.
00417         if (lnz >= localsize) {  
00418           // Need more memory to store nonzeros.
00419           localsize += localblock;
00420           iv = (int *) realloc(iv, localsize * sizeof(int));
00421           jv = (int *) realloc(jv, localsize * sizeof(int));
00422           vv = (double *) realloc(vv, localsize * sizeof(double));
00423         }
00424         iv[lnz] = I;
00425         jv[lnz] = J;
00426         vv[lnz] = V;
00427         lnz++;
00428         if (I < prevrow) rowmajor = 0;
00429         prevrow = I;
00430       }
00431     }
00432 
00433     // Status check
00434     if (nread / 1000000 > nmillion) {
00435       nmillion++;
00436       if (verbose && me == 0) std::cout << nmillion << "M ";
00437     }
00438   }
00439 
00440   delete [] buffer;
00441 
00442   if (!rowmajor) {
00443     // Sort into row-major order (by iv) so can insert entire rows at once.
00444     // Reorder jv and vv to parallel iv.
00445     if (verbose && me == 0) std::cout << std::endl << "   Sorting local nonzeros" << std::endl;
00446     sort_three(iv, jv, vv, 0, lnz-1);
00447   }
00448 
00449   //  Compute number of entries per local row for use in constructor;
00450   //  saves reallocs in FillComplete.
00451 
00452   //  Now construct the matrix.
00453   //  Count number of entries in each row so can use StaticProfile=2.
00454   if (verbose && me == 0) std::cout << std::endl << "   Constructing the matrix" << std::endl;
00455   int numRows = rowMap1->NumMyElements();
00456   int *numNonzerosPerRow = new int[numRows];
00457   for (i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
00458   for (i = 0; i < lnz; i++) 
00459     numNonzerosPerRow[rowMap1->LID(iv[i])]++;
00460 
00461   if (rowMap!=0 && colMap !=0) 
00462     A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
00463   else if (rowMap!=0) {
00464     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00465     // Less memory will be needed in FillComplete.
00466     A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
00467   }
00468   else {
00469     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00470     // Less memory will be needed in FillComplete.
00471     A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
00472   }
00473   A->SetTracebackMode(1);
00474 
00475   // Rows are inserted in ascending global number, and the insertion uses numNonzerosPerRow.
00476   // Therefore numNonzerosPerRow must be permuted, as it was created in ascending local order.
00477   int *gidList = new int[numRows];
00478   for (i = 0; i < numRows; i++) gidList[i] = rowMap1->GID(i);
00479   Epetra_Util::Sort(true,numRows,gidList,0,0,1,&numNonzerosPerRow);
00480   delete [] gidList;
00481 
00482   //  Insert the global values into the matrix row-by-row.
00483   if (verbose && me == 0) std::cout << "   Inserting global values" << std::endl;
00484   for (int sum = 0, i = 0; i < numRows; i++) {
00485     if (numNonzerosPerRow[i]) {
00486       int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i], 
00487                                        &vv[sum], &jv[sum]);
00488       if (ierr<0) EPETRA_CHK_ERR(ierr);
00489       sum += numNonzerosPerRow[i];
00490     }
00491   }
00492 
00493   delete [] numNonzerosPerRow;
00494   free(iv);
00495   free(jv);
00496   free(vv);
00497     
00498   if (verbose && me == 0) std::cout << "   Completing matrix fill" << std::endl;
00499   if (rangeMap != 0 && domainMap != 0) {
00500     EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
00501   }
00502   else if (M!=N) {
00503     Epetra_Map newDomainMap(N,rowMap1->IndexBase(), comm);
00504     EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
00505   }
00506   else {
00507     EPETRA_CHK_ERR(A->FillComplete());
00508   }
00509 
00510   if (allocatedHere) delete rowMap1;
00511   
00512   if (handle!=0) fclose(handle);
00513   double dt = timer.ElapsedTime();
00514   if (verbose && me == 0) std::cout << "File Read time (secs):  " << dt << std::endl;
00515   return(0);
00516 }
00517 #endif
00518 
00519 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00520 int MatrixMarketFileToCrsMatrixHandle64(const char *filename,
00521   const Epetra_Comm & comm,
00522   Epetra_CrsMatrix * & A,
00523   const Epetra_Map * rowMap,
00524   const Epetra_Map * colMap,
00525   const Epetra_Map * rangeMap,
00526   const Epetra_Map * domainMap,
00527   const bool transpose,
00528   const bool verbose
00529 )
00530 {
00531   const int chunk_read = 500000;  //  Modify this variable to change the
00532                                   //  size of the chunks read from the file.
00533   const int headerlineLength = 257;
00534   const int lineLength = 81;
00535   const int tokenLength = 35;
00536   char line[headerlineLength];
00537   char token1[tokenLength];
00538   char token2[tokenLength];
00539   char token3[tokenLength];
00540   char token4[tokenLength];
00541   char token5[tokenLength];
00542   long long M, N, NZ;      // Matrix dimensions
00543   int i;
00544   int me = comm.MyPID();
00545 
00546   Epetra_Time timer(comm);
00547 
00548   // Make sure domain and range maps are either both defined or undefined 
00549   if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
00550     EPETRA_CHK_ERR(-3);
00551   }
00552 
00553   // check maps to see if domain and range are 1-to-1
00554 
00555   if (domainMap!=0) {
00556     if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00557     if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00558   }
00559   else {
00560     // If domain and range are not specified, 
00561     // rowMap becomes both and must be 1-to-1 if specified
00562     if (rowMap!=0) {
00563       if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00564     }
00565   }
00566       
00567   FILE * handle = 0;
00568   if (me == 0) {
00569     if (verbose) std::cout << "Reading MatrixMarket file " << filename << std::endl;
00570     handle = fopen(filename,"r");  // Open file
00571     if (handle == 0)
00572       EPETRA_CHK_ERR(-1); // file not found
00573 
00574     // Check first line, which should be 
00575     // %%MatrixMarket matrix coordinate real general
00576     if(fgets(line, headerlineLength, handle)==0) {
00577       if (handle!=0) fclose(handle); 
00578       EPETRA_CHK_ERR(-1);
00579     }
00580     if(sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
00581       if (handle!=0) fclose(handle); 
00582       EPETRA_CHK_ERR(-1);
00583     }
00584     if (strcmp(token1, "%%MatrixMarket") ||
00585         strcmp(token2, "matrix") ||
00586         strcmp(token3, "coordinate") ||
00587         strcmp(token4, "real") ||
00588         strcmp(token5, "general")) {
00589       if (handle!=0) fclose(handle); 
00590       EPETRA_CHK_ERR(-1);
00591     }
00592 
00593     // Next, strip off header lines (which start with "%")
00594     do {
00595       if(fgets(line, headerlineLength, handle)==0) {
00596         if (handle!=0) fclose(handle); 
00597         EPETRA_CHK_ERR(-1);
00598       }
00599     } while (line[0] == '%');
00600 
00601     // Next get problem dimensions: M, N, NZ
00602     if(sscanf(line, "%lld %lld %lld", &M, &N, &NZ)==0) {
00603       if (handle!=0) fclose(handle); 
00604       EPETRA_CHK_ERR(-1);
00605     }
00606   }
00607   comm.Broadcast(&M, 1, 0);
00608   comm.Broadcast(&N, 1, 0);
00609   comm.Broadcast(&NZ, 1, 0);
00610 
00611   // Now create matrix using user maps if provided.
00612 
00613 
00614   // Now read in chunks of triplets and broadcast them to other processors.
00615   char *buffer = new char[chunk_read*lineLength];
00616   long long nchunk; 
00617   int nmillion = 0;
00618   long long nread = 0;
00619   int rlen;
00620 
00621   // Storage for this processor's nonzeros.
00622   const int localblock = 100000;
00623   int localsize = (int) (NZ / comm.NumProc()) + localblock;
00624   long long *iv = (long long *) malloc(localsize * sizeof(long long));
00625   long long *jv = (long long *) malloc(localsize * sizeof(long long));
00626   double *vv = (double *) malloc(localsize * sizeof(double));
00627   int lnz = 0;   //  Number of non-zeros on this processor.
00628 
00629   if (!iv || !jv || !vv) 
00630     EPETRA_CHK_ERR(-1);
00631 
00632   Epetra_Map *rowMap1;
00633   bool allocatedHere=false;
00634   if (rowMap != 0) 
00635     rowMap1 = const_cast<Epetra_Map *>(rowMap);
00636   else {
00637     rowMap1 = new Epetra_Map(M, 0, comm);
00638     allocatedHere = true;
00639   }
00640   long long ioffset = rowMap1->IndexBase64()-1;
00641   long long joffset = (colMap != 0 ? colMap->IndexBase64()-1 : ioffset);
00642 
00643   int rowmajor = 1;  // Assume non-zeros are listed in row-major order, 
00644   long long prevrow = -1;  // but test to detect otherwise.  If non-zeros
00645                      // are row-major, we can avoid the sort.
00646 
00647   while (nread < NZ) {
00648     if (NZ-nread > chunk_read) nchunk = chunk_read;
00649     else nchunk = NZ - nread;
00650 
00651     if (me == 0) {
00652       char *eof;
00653       rlen = 0;
00654       for (int i = 0; i < nchunk; i++) {
00655         eof = fgets(&buffer[rlen],lineLength,handle);
00656         if (eof == NULL) {
00657           fprintf(stderr, "%s", "Unexpected end of matrix file.");
00658           EPETRA_CHK_ERR(-1);
00659         }
00660         rlen += strlen(&buffer[rlen]);
00661       }
00662       buffer[rlen++]='\n';
00663     }
00664     comm.Broadcast(&rlen, 1, 0);
00665     comm.Broadcast(buffer, rlen, 0);
00666 
00667     buffer[rlen++] = '\0';
00668     nread += nchunk;
00669 
00670     // Process the received data, saving non-zeros belonging on this proc.
00671     char *lineptr = buffer;
00672     for (rlen = 0; rlen < nchunk; rlen++) {
00673       char *next = strchr(lineptr,'\n');
00674       char *endp;
00675       const int base = 10;
00676 #if defined(_MSC_VER)
00677       long long I = _strtoi64(strtok(lineptr," \t\n"), &endp, base) + ioffset;
00678       long long J = _strtoi64(strtok(NULL," \t\n"), &endp, base) + joffset;
00679 #else
00680 #if defined(__PGI)
00681       long long I, J;
00682       std::istringstream ssI(strtok(lineptr," \t\n"));
00683       ssI >> I; I += ioffset;
00684       std::istringstream ssJ(strtok(NULL," \t\n"));
00685       ssJ >> J; J += joffset;
00686 #else
00687       long long I = strtoll(strtok(lineptr," \t\n"), &endp, base) + ioffset;
00688       long long J = strtoll(strtok(NULL," \t\n"), &endp, base) + joffset;
00689 #endif
00690 #endif
00691       double V = atof(strtok(NULL," \t\n"));
00692       lineptr = next + 1;
00693       if (transpose) {
00694         // swap I and J indices.
00695         long long tmp = I;
00696         I = J;
00697         J = tmp;
00698       }
00699       if (rowMap1->MyGID(I)) {
00700         //  This processor keeps this non-zero.
00701         if (lnz >= localsize) {  
00702           // Need more memory to store nonzeros.
00703           localsize += localblock;
00704           iv = (long long *) realloc(iv, localsize * sizeof(long long));
00705           jv = (long long *) realloc(jv, localsize * sizeof(long long));
00706           vv = (double *) realloc(vv, localsize * sizeof(double));
00707         }
00708         iv[lnz] = I;
00709         jv[lnz] = J;
00710         vv[lnz] = V;
00711         lnz++;
00712         if (I < prevrow) rowmajor = 0;
00713         prevrow = I;
00714       }
00715     }
00716 
00717     // Status check
00718     if (nread / 1000000 > nmillion) {
00719       nmillion++;
00720       if (verbose && me == 0) std::cout << nmillion << "M ";
00721     }
00722   }
00723 
00724   delete [] buffer;
00725 
00726   if (!rowmajor) {
00727     // Sort into row-major order (by iv) so can insert entire rows at once.
00728     // Reorder jv and vv to parallel iv.
00729     if (verbose && me == 0) std::cout << std::endl << "   Sorting local nonzeros" << std::endl;
00730     sort_three(iv, jv, vv, 0, lnz-1);
00731   }
00732 
00733   //  Compute number of entries per local row for use in constructor;
00734   //  saves reallocs in FillComplete.
00735 
00736   //  Now construct the matrix.
00737   //  Count number of entries in each row so can use StaticProfile=2.
00738   if (verbose && me == 0) std::cout << std::endl << "   Constructing the matrix" << std::endl;
00739   int numRows = rowMap1->NumMyElements();
00740   int *numNonzerosPerRow = new int[numRows];
00741   for (i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
00742   for (i = 0; i < lnz; i++) 
00743     numNonzerosPerRow[rowMap1->LID(iv[i])]++;
00744 
00745   if (rowMap!=0 && colMap !=0) 
00746     A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
00747   else if (rowMap!=0) {
00748     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00749     // Less memory will be needed in FillComplete.
00750     A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
00751   }
00752   else {
00753     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00754     // Less memory will be needed in FillComplete.
00755     A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
00756   }
00757   A->SetTracebackMode(1);
00758 
00759   // Rows are inserted in ascending global number, and the insertion uses numNonzerosPerRow.
00760   // Therefore numNonzerosPerRow must be permuted, as it was created in ascending local order.
00761   long long *gidList = new long long[numRows];
00762   for (i = 0; i < numRows; i++) gidList[i] = rowMap1->GID64(i);
00763   Epetra_Util::Sort(true,numRows,gidList,0,0,1,&numNonzerosPerRow,0,0);
00764   delete [] gidList;
00765 
00766   //  Insert the global values into the matrix row-by-row.
00767   if (verbose && me == 0) std::cout << "   Inserting global values" << std::endl;
00768   for (int sum = 0, i = 0; i < numRows; i++) {
00769     if (numNonzerosPerRow[i]) {
00770       int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i], 
00771                                        &vv[sum], &jv[sum]);
00772       if (ierr<0) EPETRA_CHK_ERR(ierr);
00773       sum += numNonzerosPerRow[i];
00774     }
00775   }
00776 
00777   delete [] numNonzerosPerRow;
00778   free(iv);
00779   free(jv);
00780   free(vv);
00781     
00782   if (verbose && me == 0) std::cout << "   Completing matrix fill" << std::endl;
00783   if (rangeMap != 0 && domainMap != 0) {
00784     EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
00785   }
00786   else if (M!=N) {
00787     Epetra_Map newDomainMap(N,rowMap1->IndexBase64(), comm);
00788     EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
00789   }
00790   else {
00791     EPETRA_CHK_ERR(A->FillComplete());
00792   }
00793 
00794   if (allocatedHere) delete rowMap1;
00795   
00796   if (handle!=0) fclose(handle);
00797   double dt = timer.ElapsedTime();
00798   if (verbose && me == 0) std::cout << "File Read time (secs):  " << dt << std::endl;
00799   return(0);
00800 }
00801 #endif
00802 
00804 // Sorting values in array list in increasing order. Criteria is int. 
00805 // Also rearrange values in arrays parlista and partlistb
00806 // to match the new order of list. 
00807 
00808 template<typename int_type>
00809 static void quickpart_list_inc_int (
00810   int_type *list, int_type *parlista, double *parlistb,
00811   int start, int end, int *equal, int *larger)
00812 {
00813 int i;
00814 int_type key, itmp;
00815 double dtmp;
00816 
00817   key = list ? list[(end+start)/2] : 1;
00818 
00819   *equal = *larger = start;
00820   for (i = start; i <= end; i++)
00821     if (list[i] < key) {
00822       itmp                = parlista[i];
00823       parlista[i]         = parlista[*larger];
00824       parlista[(*larger)] = parlista[*equal];
00825       parlista[(*equal)]  = itmp;
00826       dtmp                = parlistb[i];
00827       parlistb[i]         = parlistb[*larger];
00828       parlistb[(*larger)] = parlistb[*equal];
00829       parlistb[(*equal)]  = dtmp;
00830       itmp                = list[i];
00831       list[i]             = list[*larger];
00832       list[(*larger)++]   = list[*equal];
00833       list[(*equal)++]    = itmp;
00834     }
00835     else if (list[i] == key) {
00836       itmp                = parlista[i];
00837       parlista[i]         = parlista[*larger];
00838       parlista[(*larger)] = itmp;
00839       dtmp                = parlistb[i];
00840       parlistb[i]         = parlistb[*larger];
00841       parlistb[(*larger)] = dtmp;
00842       list[i]             = list[*larger];
00843       list[(*larger)++]   = key;
00844     }
00845 }
00846 
00848 template<typename int_type>
00849 static void sort_three(int_type* list, int_type *parlista, double *parlistb, 
00850                        int start, int end)
00851 {
00852 int  equal, larger;
00853 
00854   if (start < end) {
00855     quickpart_list_inc_int(list, parlista, parlistb, start, end, 
00856                            &equal, &larger);
00857     sort_three(list, parlista, parlistb, start,  equal-1);
00858     sort_three(list, parlista, parlistb, larger, end);
00859   }
00860 }
00861 
00863 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00864 int MatlabFileToCrsMatrix(const char *filename,
00865         const Epetra_Comm & comm,
00866         Epetra_CrsMatrix * & A)
00867 {
00868   const int lineLength = 1025;
00869   char line[lineLength];
00870   int I, J;
00871   double V;
00872 
00873   FILE * handle = 0;
00874 
00875   handle = fopen(filename,"r");  // Open file
00876   if (handle == 0)
00877     EPETRA_CHK_ERR(-1); // file not found
00878 
00879   int numGlobalRows = 0;
00880   int numGlobalCols = 0;
00881   while(fgets(line, lineLength, handle)!=0) {
00882     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00883     if (I>numGlobalRows) numGlobalRows = I;
00884     if (J>numGlobalCols) numGlobalCols = J;
00885   }
00886 
00887   if (handle!=0) fclose(handle);
00888   Epetra_Map rangeMap(numGlobalRows, 0, comm);
00889   Epetra_Map domainMap(numGlobalCols, 0, comm);
00890   A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
00891 
00892   // Now read in each triplet and store to the local portion of the matrix if the row is owned.
00893   const Epetra_Map & rowMap1 = A->RowMap();
00894   
00895   handle = 0;
00896 
00897   handle = fopen(filename,"r");  // Open file
00898   if (handle == 0)
00899     EPETRA_CHK_ERR(-1); // file not found
00900 
00901   while (fgets(line, lineLength, handle)!=0) {
00902     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00903     I--; J--; // Convert to Zero based
00904     if (rowMap1.MyGID(I)) {
00905       int ierr = A->InsertGlobalValues(I, 1, &V, &J);
00906       if (ierr<0) EPETRA_CHK_ERR(ierr);
00907     }
00908   }
00909     
00910   EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
00911 
00912   if (handle!=0) fclose(handle);
00913   return(0);
00914 }
00915 #endif
00916 
00917 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00918 int MatlabFileToCrsMatrix64(const char *filename,
00919         const Epetra_Comm & comm,
00920         Epetra_CrsMatrix * & A)
00921 {
00922   const int lineLength = 1025;
00923   char line[lineLength];
00924   long long I, J;
00925   double V;
00926 
00927   FILE * handle = 0;
00928 
00929   handle = fopen(filename,"r");  // Open file
00930   if (handle == 0)
00931     EPETRA_CHK_ERR(-1); // file not found
00932 
00933   long long numGlobalRows = 0;
00934   long long numGlobalCols = 0;
00935   while(fgets(line, lineLength, handle)!=0) {
00936     if(sscanf(line, "%lld %lld %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00937     if (I>numGlobalRows) numGlobalRows = I;
00938     if (J>numGlobalCols) numGlobalCols = J;
00939   }
00940 
00941   if (handle!=0) fclose(handle);
00942   Epetra_Map rangeMap(numGlobalRows, 0, comm);
00943   Epetra_Map domainMap(numGlobalCols, 0, comm);
00944   A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
00945 
00946   // Now read in each triplet and store to the local portion of the matrix if the row is owned.
00947   const Epetra_Map & rowMap1 = A->RowMap();
00948   
00949   handle = 0;
00950 
00951   handle = fopen(filename,"r");  // Open file
00952   if (handle == 0)
00953     EPETRA_CHK_ERR(-1); // file not found
00954 
00955   while (fgets(line, lineLength, handle)!=0) {
00956     if(sscanf(line, "%lld %lld %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00957     I--; J--; // Convert to Zero based
00958     if (rowMap1.MyGID(I)) {
00959       int ierr = A->InsertGlobalValues(I, 1, &V, &J);
00960       if (ierr<0) EPETRA_CHK_ERR(ierr);
00961     }
00962   }
00963     
00964   EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
00965 
00966   if (handle!=0) fclose(handle);
00967   return(0);
00968 }
00969 #endif
00970 
00971 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00972 int HypreFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix){
00973   int MyPID = comm.MyPID();
00974   // This double will be in the format we want for the extension besides the leading zero
00975   double filePID = (double)MyPID/(double)100000;
00976   std::ostringstream stream;
00977   // Using setprecision() puts it in the std::string
00978   stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
00979   // Then just ignore the first character
00980   std::string fileName(filename);
00981   fileName += stream.str().substr(1,7);
00982   // Open the file
00983   std::ifstream file(fileName.c_str());
00984   std::string line;
00985   if(file.is_open()){
00986     std::getline(file, line);
00987     int ilower, iupper;
00988     std::istringstream istream(line);
00989     // The first line of the file has the beginning and ending rows
00990     istream >> ilower;
00991     istream >> iupper;
00992     // Using those we can create a row map
00993     Epetra_Map RowMap(-1, iupper-ilower+1, 0, comm);
00994     Matrix = new Epetra_CrsMatrix(Copy, RowMap, 0);
00995     int currRow = -1;
00996     int counter = 0;
00997     std::vector<int> indices;
00998     std::vector<double> values;
00999     while(!file.eof()){
01000       std::getline(file, line);
01001       std::istringstream lineStr(line);
01002       int row, col;
01003       double val;
01004       lineStr >> row;
01005       lineStr >> col;
01006       lineStr >> val;
01007       if(currRow == -1) currRow = row; // First line
01008       if(row == currRow){
01009         // add to the vector
01010         counter = counter + 1;
01011         indices.push_back(col);
01012         values.push_back(val);
01013       } else {
01014         Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
01015         indices.clear();
01016         values.clear();
01017         counter = 0;
01018         currRow = row;
01019         // make a new vector
01020         indices.push_back(col);
01021         values.push_back(val);
01022         counter = counter + 1;
01023       }
01024     }
01025     Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
01026     Matrix->Comm().Barrier();
01027     Matrix->FillComplete();
01028     file.close();
01029     return 0;
01030   } else {
01031     std::cout << "\nERROR:\nCouldn't open " << fileName << ".\n";
01032     return -1;
01033   }
01034 }
01035 #endif
01036 
01037 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01038 int HypreFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix){
01039   int MyPID = comm.MyPID();
01040   // This double will be in the format we want for the extension besides the leading zero
01041   double filePID = (double)MyPID/(double)100000;
01042   std::ostringstream stream;
01043   // Using setprecision() puts it in the std::string
01044   stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
01045   // Then just ignore the first character
01046   std::string fileName(filename);
01047   fileName += stream.str().substr(1,7);
01048   // Open the file
01049   std::ifstream file(fileName.c_str());
01050   std::string line;
01051   if(file.is_open()){
01052     std::getline(file, line);
01053     int ilower, iupper;
01054     std::istringstream istream(line);
01055     // The first line of the file has the beginning and ending rows
01056     istream >> ilower;
01057     istream >> iupper;
01058     // Using those we can create a row map
01059     Epetra_Map RowMap(-1LL, iupper-ilower+1, 0LL, comm);
01060     Matrix = new Epetra_CrsMatrix(Copy, RowMap, 0);
01061     long long currRow = -1;
01062     int counter = 0;
01063     std::vector<long long> indices;
01064     std::vector<double> values;
01065     while(!file.eof()){
01066       std::getline(file, line);
01067       std::istringstream lineStr(line);
01068       long long row, col;
01069       double val;
01070       lineStr >> row;
01071       lineStr >> col;
01072       lineStr >> val;
01073       if(currRow == -1) currRow = row; // First line
01074       if(row == currRow){
01075         // add to the vector
01076         counter = counter + 1;
01077         indices.push_back(col);
01078         values.push_back(val);
01079       } else {
01080         Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
01081         indices.clear();
01082         values.clear();
01083         counter = 0;
01084         currRow = row;
01085         // make a new vector
01086         indices.push_back(col);
01087         values.push_back(val);
01088         counter = counter + 1;
01089       }
01090     }
01091     Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
01092     Matrix->Comm().Barrier();
01093     Matrix->FillComplete();
01094     file.close();
01095     return 0;
01096   } else {
01097     std::cout << "\nERROR:\nCouldn't open " << fileName << ".\n";
01098     return -1;
01099   }
01100 }
01101 #endif
01102 
01103 } // namespace EpetraExt
01104 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines