EpetraExt_CrsMatrixIn.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 #include "EpetraExt_CrsMatrixIn.h"
00029 #include "Epetra_Comm.h"
00030 #include "Epetra_CrsMatrix.h"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_IntVector.h"
00033 #include "Epetra_IntSerialDenseVector.h"
00034 #include "Epetra_Import.h"
00035 #include "Epetra_Time.h"
00036 #include "Epetra_Util.h"
00037 
00040 // Functions to read a MatrixMarket file and load it into a matrix.
00041 // Adapted from 
00042 //    Trilinos/packages/epetraext/src/inout/EpetraExt_CrsMatrixIn.cpp
00043 // Modified by Jon Berry and Karen Devine to make matrix reallocations 
00044 // more efficient; rather than insert each non-zero separately, we 
00045 // collect rows of non-zeros for insertion.
00046 // Modified by Karen Devine and Steve Plimpton to prevent all processors 
00047 // from reading the same file at the same time (ouch!); chunks of the file 
00048 // are read and broadcast by processor zero; each processor then extracts 
00049 // its nonzeros from the chunk, sorts them by row, and inserts them into 
00050 // the matrix.
00051 // The variable "chunk_read" can be changed to modify the size of the chunks
00052 // read from the file.  Larger values of chunk_read lead to faster execution
00053 // and greater memory use.
00056 
00057 using namespace EpetraExt;
00058 namespace EpetraExt {
00059 
00060 static void sort_three(int* list, int *parlista, double *parlistb, 
00061                        int start, int end);
00063 int MatrixMarketFileToCrsMatrix(const char *filename, 
00064     const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
00065 {
00066   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A));
00067   return(0);
00068 }
00069 
00071 int MatrixMarketFileToCrsMatrix(const char *filename, 
00072     const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
00073 {
00074   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
00075                                                    0, 0, 0, 0, transpose));
00076   return(0);
00077 }
00079   
00080 int MatrixMarketFileToCrsMatrix(const char *filename, 
00081     const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
00082     const bool verbose) 
00083 {
00084   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
00085                                                    0, 0, 0, 0,
00086                                                    transpose, verbose));
00087   return(0);
00088 }
00089   
00091 int MatrixMarketFileToCrsMatrix(const char *filename, 
00092     const Epetra_Map & rowMap, const Epetra_Map& rangeMap, 
00093     const Epetra_Map& domainMap, Epetra_CrsMatrix * & A, const bool transpose,
00094     const bool verbose) 
00095 {
00096   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00097                                                    rowMap.Comm(), A,
00098                                                    &rowMap, 0, 
00099                                                    &rangeMap, &domainMap,
00100                                                    transpose, verbose));
00101   return(0);
00102 }
00103 
00105 int MatrixMarketFileToCrsMatrix(const char *filename, 
00106     const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
00107     const bool verbose) 
00108 {
00109   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00110                                                    rowMap.Comm(), A,
00111                                                    &rowMap, 0, 0, 0, 
00112                                                    transpose, verbose));
00113   return(0);
00114 }
00115 
00117 int MatrixMarketFileToCrsMatrix(const char *filename, 
00118     const Epetra_Map & rowMap, const Epetra_Map & colMap, 
00119     Epetra_CrsMatrix * & A, const bool transpose,
00120     const bool verbose) 
00121 {
00122   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00123                                                    rowMap.Comm(), A,
00124                                                    &rowMap, &colMap, 0, 0,
00125                                                    transpose, verbose));
00126   return(0);
00127 }
00128 
00130 int MatrixMarketFileToCrsMatrix(const char *filename, 
00131     const Epetra_Map & rowMap, const Epetra_Map & colMap,
00132     const Epetra_Map& rangeMap, const Epetra_Map& domainMap, 
00133     Epetra_CrsMatrix * & A, const bool transpose,
00134     const bool verbose) 
00135 {
00136   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, 
00137                                                    rowMap.Comm(), A, 
00138                                                    &rowMap, &colMap, 
00139                                                    &rangeMap, &domainMap,
00140                                                    transpose, verbose));
00141   return(0);
00142 }
00143 
00145 int MatrixMarketFileToCrsMatrixHandle(const char *filename,
00146   const Epetra_Comm & comm,
00147   Epetra_CrsMatrix * & A,
00148   const Epetra_Map * rowMap,
00149   const Epetra_Map * colMap,
00150   const Epetra_Map * rangeMap,
00151   const Epetra_Map * domainMap,
00152   const bool transpose,
00153   const bool verbose
00154 )
00155 {
00156   const int chunk_read = 500000;  //  Modify this variable to change the
00157                                   //  size of the chunks read from the file.
00158   const int headerlineLength = 257;
00159   const int lineLength = 81;
00160   const int tokenLength = 35;
00161   char line[lineLength];
00162   char token1[tokenLength];
00163   char token2[tokenLength];
00164   char token3[tokenLength];
00165   char token4[tokenLength];
00166   char token5[tokenLength];
00167   int M, N, NZ;      // Matrix dimensions
00168   int i;
00169   int me = comm.MyPID();
00170   int nproc = comm.NumProc();
00171 
00172   Epetra_Time timer(comm);
00173 
00174   // Make sure domain and range maps are either both defined or undefined 
00175   if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
00176     EPETRA_CHK_ERR(-3);
00177   }
00178 
00179   // check maps to see if domain and range are 1-to-1
00180 
00181   if (domainMap!=0) {
00182     if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00183     if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00184   }
00185   else {
00186     // If domain and range are not specified, 
00187     // rowMap becomes both and must be 1-to-1 if specified
00188     if (rowMap!=0) {
00189       if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00190     }
00191   }
00192       
00193   FILE * handle = 0;
00194   if (me == 0) {
00195     if (verbose) cout << "Reading MatrixMarket file " << filename << endl;
00196     handle = fopen(filename,"r");  // Open file
00197     if (handle == 0)
00198       EPETRA_CHK_ERR(-1); // file not found
00199 
00200     // Check first line, which should be 
00201     // %%MatrixMarket matrix coordinate real general
00202     if(fgets(line, headerlineLength, handle)==0) {
00203       if (handle!=0) fclose(handle); 
00204       EPETRA_CHK_ERR(-1);
00205     }
00206     if(sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
00207       if (handle!=0) fclose(handle); 
00208       EPETRA_CHK_ERR(-1);
00209     }
00210     if (strcmp(token1, "%%MatrixMarket") ||
00211         strcmp(token2, "matrix") ||
00212         strcmp(token3, "coordinate") ||
00213         strcmp(token4, "real") ||
00214         strcmp(token5, "general")) {
00215       if (handle!=0) fclose(handle); 
00216       EPETRA_CHK_ERR(-1);
00217     }
00218 
00219     // Next, strip off header lines (which start with "%")
00220     do {
00221       if(fgets(line, headerlineLength, handle)==0) {
00222         if (handle!=0) fclose(handle); 
00223         EPETRA_CHK_ERR(-1);
00224       }
00225     } while (line[0] == '%');
00226 
00227     // Next get problem dimensions: M, N, NZ
00228     if(sscanf(line, "%d %d %d", &M, &N, &NZ)==0) {
00229       if (handle!=0) fclose(handle); 
00230       EPETRA_CHK_ERR(-1);
00231     }
00232   }
00233   comm.Broadcast(&M, 1, 0);
00234   comm.Broadcast(&N, 1, 0);
00235   comm.Broadcast(&NZ, 1, 0);
00236 
00237   // Now create matrix using user maps if provided.
00238 
00239 
00240   // Now read in chunks of triplets and broadcast them to other processors.
00241   char *buffer = new char[chunk_read*lineLength];
00242   int nchunk; 
00243   int nmillion = 0;
00244   int nread = 0;
00245   int rlen;
00246 
00247   // Storage for this processor's nonzeros.
00248   const int localblock = 100000;
00249   int localsize = NZ / comm.NumProc() + localblock;
00250   int *iv = (int *) malloc(localsize * sizeof(int));
00251   int *jv = (int *) malloc(localsize * sizeof(int));
00252   double *vv = (double *) malloc(localsize * sizeof(double));
00253   int lnz = 0;   //  Number of non-zeros on this processor.
00254 
00255   if (!iv || !jv || !vv) 
00256     EPETRA_CHK_ERR(-1);
00257 
00258   Epetra_Map *rowMap1;
00259   if (rowMap != 0) 
00260     rowMap1 = const_cast<Epetra_Map *>(rowMap);
00261   else
00262     rowMap1 = new Epetra_Map(M, 0, comm);
00263   int ioffset = rowMap1->IndexBase()-1;
00264   int joffset = (colMap != 0 ? colMap->IndexBase()-1 : ioffset);
00265 
00266   int rowmajor = 1;  // Assume non-zeros are listed in row-major order, 
00267   int prevrow = -1;  // but test to detect otherwise.  If non-zeros
00268                      // are row-major, we can avoid the sort.
00269 
00270   while (nread < NZ) {
00271     if (NZ-nread > chunk_read) nchunk = chunk_read;
00272     else nchunk = NZ - nread;
00273 
00274     if (me == 0) {
00275       char *eof;
00276       rlen = 0;
00277       for (int i = 0; i < nchunk; i++) {
00278         eof = fgets(&buffer[rlen],lineLength,handle);
00279         if (eof == NULL) {
00280           fprintf(stderr, "Unexpected end of matrix file.");
00281           EPETRA_CHK_ERR(-1);
00282         }
00283         rlen += strlen(&buffer[rlen]);
00284       }
00285       buffer[rlen++]='\n';
00286     }
00287     comm.Broadcast(&rlen, 1, 0);
00288     comm.Broadcast(buffer, rlen, 0);
00289 
00290     buffer[rlen++] = '\0';
00291     nread += nchunk;
00292 
00293     // Process the received data, saving non-zeros belonging on this proc.
00294     char *lineptr = buffer;
00295     for (rlen = 0; rlen < nchunk; rlen++) {
00296       char *next = strchr(lineptr,'\n');
00297       int I = atoi(strtok(lineptr," \t\n")) + ioffset;
00298       int J = atoi(strtok(NULL," \t\n")) + joffset;
00299       double V = atof(strtok(NULL," \t\n"));
00300       lineptr = next + 1;
00301       if (transpose) {
00302         // swap I and J indices.
00303         int tmp = I;
00304         I = J;
00305         J = tmp;
00306       }
00307       if (rowMap1->MyGID(I)) {
00308         //  This processor keeps this non-zero.
00309         if (lnz >= localsize) {  
00310           // Need more memory to store nonzeros.
00311           localsize += localblock;
00312           iv = (int *) realloc(iv, localsize * sizeof(int));
00313           jv = (int *) realloc(jv, localsize * sizeof(int));
00314           vv = (double *) realloc(vv, localsize * sizeof(double));
00315         }
00316         iv[lnz] = I;
00317         jv[lnz] = J;
00318         vv[lnz] = V;
00319         lnz++;
00320         if (I < prevrow) rowmajor = 0;
00321         prevrow = I;
00322       }
00323     }
00324 
00325     // Status check
00326     if (nread / 1000000 > nmillion) {
00327       nmillion++;
00328       if (verbose && me == 0) cout << nmillion << "M ";
00329     }
00330   }
00331 
00332   delete [] buffer;
00333 
00334   if (!rowmajor) {
00335     // Sort into row-major order (by iv) so can insert entire rows at once.
00336     // Reorder jv and vv to parallel iv.
00337     if (verbose && me == 0) cout << endl << "   Sorting local nonzeros" << endl;
00338     sort_three(iv, jv, vv, 0, lnz-1);
00339   }
00340 
00341   //  Compute number of entries per local row for use in constructor;
00342   //  saves reallocs in FillComplete.
00343 
00344   //  Now construct the matrix.
00345   //  Count number of entries in each row so can use StaticProfile=2.
00346   if (verbose && me == 0) cout << endl << "   Constructing the matrix" << endl;
00347   int numRows = rowMap1->NumMyElements();
00348   int *numNonzerosPerRow = new int[numRows];
00349   for (i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
00350   for (i = 0; i < lnz; i++) 
00351     numNonzerosPerRow[rowMap1->LID(iv[i])]++;
00352 
00353   if (rowMap!=0 && colMap !=0) 
00354     A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
00355   else if (rowMap!=0) {
00356     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00357     // Less memory will be needed in FillComplete.
00358     A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
00359   }
00360   else {
00361     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00362     // Less memory will be needed in FillComplete.
00363     A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
00364   }
00365   A->SetTracebackMode(2);
00366 
00367   //  Insert the global values into the matrix row-by-row.
00368   if (verbose && me == 0) cout << "   Inserting global values" << endl;
00369   for (int sum = 0, i = 0; i < numRows; i++) {
00370     if (numNonzerosPerRow[i]) {
00371       int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i], 
00372                                        &vv[sum], &jv[sum]);
00373       if (ierr<0) EPETRA_CHK_ERR(ierr);
00374       sum += numNonzerosPerRow[i];
00375     }
00376   }
00377 
00378   delete [] numNonzerosPerRow;
00379   free(iv); iv=NULL;
00380   free(jv); jv=NULL;
00381   free(vv); vv=NULL;
00382     
00383   if (verbose && me == 0) cout << "   Completing matrix fill" << endl;
00384   if (rangeMap != 0 && domainMap != 0) {
00385     EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
00386   }
00387   else if (M!=N) {
00388     Epetra_Map newDomainMap(N,rowMap1->IndexBase(), comm);
00389     EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
00390   }
00391   else {
00392     EPETRA_CHK_ERR(A->FillComplete());
00393   }
00394   
00395   if (handle!=0) fclose(handle);
00396   double dt = timer.ElapsedTime();
00397   if (verbose && me == 0) cout << "File Read time (secs):  " << dt << endl;
00398   return(0);
00399 }
00400 
00402 // Sorting values in array list in increasing order. Criteria is int. 
00403 // Also rearrange values in arrays parlista and partlistb
00404 // to match the new order of list. 
00405 
00406 static void quickpart_list_inc_int (
00407   int *list, int *parlista, double *parlistb,
00408   int start, int end, int *equal, int *larger)
00409 {
00410 int i, key, itmp;
00411 double dtmp;
00412 
00413   key = list ? list[(end+start)/2] : 1;
00414 
00415   *equal = *larger = start;
00416   for (i = start; i <= end; i++)
00417     if (list[i] < key) {
00418       itmp                = parlista[i];
00419       parlista[i]         = parlista[*larger];
00420       parlista[(*larger)] = parlista[*equal];
00421       parlista[(*equal)]  = itmp;
00422       dtmp                = parlistb[i];
00423       parlistb[i]         = parlistb[*larger];
00424       parlistb[(*larger)] = parlistb[*equal];
00425       parlistb[(*equal)]  = dtmp;
00426       itmp                = list[i];
00427       list[i]             = list[*larger];
00428       list[(*larger)++]   = list[*equal];
00429       list[(*equal)++]    = itmp;
00430     }
00431     else if (list[i] == key) {
00432       itmp                = parlista[i];
00433       parlista[i]         = parlista[*larger];
00434       parlista[(*larger)] = itmp;
00435       dtmp                = parlistb[i];
00436       parlistb[i]         = parlistb[*larger];
00437       parlistb[(*larger)] = dtmp;
00438       list[i]             = list[*larger];
00439       list[(*larger)++]   = key;
00440     }
00441 }
00442 
00444 static void sort_three(int* list, int *parlista, double *parlistb, 
00445                        int start, int end)
00446 {
00447 int  equal, larger;
00448 
00449   if (start < end) {
00450     quickpart_list_inc_int(list, parlista, parlistb, start, end, 
00451                            &equal, &larger);
00452     sort_three(list, parlista, parlistb, start,  equal-1);
00453     sort_three(list, parlista, parlistb, larger, end);
00454   }
00455 }
00456 
00458 int MatlabFileToCrsMatrix(const char *filename,
00459         const Epetra_Comm & comm,
00460         Epetra_CrsMatrix * & A)
00461 {
00462   const int lineLength = 1025;
00463   char line[lineLength];
00464   int I, J;
00465   double V;
00466 
00467   FILE * handle = 0;
00468 
00469   handle = fopen(filename,"r");  // Open file
00470   if (handle == 0)
00471     EPETRA_CHK_ERR(-1); // file not found
00472 
00473   int numGlobalRows = 0;
00474   int numGlobalCols = 0;
00475   while(fgets(line, lineLength, handle)!=0) {
00476     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00477     if (I>numGlobalRows) numGlobalRows = I;
00478     if (J>numGlobalCols) numGlobalCols = J;
00479   }
00480 
00481   if (handle!=0) fclose(handle);
00482   Epetra_Map rangeMap(numGlobalRows, 0, comm);
00483   Epetra_Map domainMap(numGlobalCols, 0, comm);
00484   A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
00485 
00486   // Now read in each triplet and store to the local portion of the matrix if the row is owned.
00487   const Epetra_Map & rowMap1 = A->RowMap();
00488   
00489   handle = 0;
00490 
00491   handle = fopen(filename,"r");  // Open file
00492   if (handle == 0)
00493     EPETRA_CHK_ERR(-1); // file not found
00494 
00495   while (fgets(line, lineLength, handle)!=0) {
00496     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00497     I--; J--; // Convert to Zero based
00498     if (rowMap1.MyGID(I)) {
00499       int ierr = A->InsertGlobalValues(I, 1, &V, &J);
00500       if (ierr<0) EPETRA_CHK_ERR(ierr);
00501     }
00502   }
00503     
00504   EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
00505 
00506   if (handle!=0) fclose(handle);
00507   return(0);
00508 }
00509 } // namespace EpetraExt
00510 

Generated on Wed May 12 21:40:37 2010 for EpetraExt by  doxygen 1.4.7