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 
00171   Epetra_Time timer(comm);
00172 
00173   // Make sure domain and range maps are either both defined or undefined 
00174   if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
00175     EPETRA_CHK_ERR(-3);
00176   }
00177 
00178   // check maps to see if domain and range are 1-to-1
00179 
00180   if (domainMap!=0) {
00181     if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00182     if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00183   }
00184   else {
00185     // If domain and range are not specified, 
00186     // rowMap becomes both and must be 1-to-1 if specified
00187     if (rowMap!=0) {
00188       if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00189     }
00190   }
00191       
00192   FILE * handle = 0;
00193   if (me == 0) {
00194     if (verbose) cout << "Reading MatrixMarket file " << filename << endl;
00195     handle = fopen(filename,"r");  // Open file
00196     if (handle == 0)
00197       EPETRA_CHK_ERR(-1); // file not found
00198 
00199     // Check first line, which should be 
00200     // %%MatrixMarket matrix coordinate real general
00201     if(fgets(line, headerlineLength, handle)==0) {
00202       if (handle!=0) fclose(handle); 
00203       EPETRA_CHK_ERR(-1);
00204     }
00205     if(sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
00206       if (handle!=0) fclose(handle); 
00207       EPETRA_CHK_ERR(-1);
00208     }
00209     if (strcmp(token1, "%%MatrixMarket") ||
00210         strcmp(token2, "matrix") ||
00211         strcmp(token3, "coordinate") ||
00212         strcmp(token4, "real") ||
00213         strcmp(token5, "general")) {
00214       if (handle!=0) fclose(handle); 
00215       EPETRA_CHK_ERR(-1);
00216     }
00217 
00218     // Next, strip off header lines (which start with "%")
00219     do {
00220       if(fgets(line, headerlineLength, handle)==0) {
00221         if (handle!=0) fclose(handle); 
00222         EPETRA_CHK_ERR(-1);
00223       }
00224     } while (line[0] == '%');
00225 
00226     // Next get problem dimensions: M, N, NZ
00227     if(sscanf(line, "%d %d %d", &M, &N, &NZ)==0) {
00228       if (handle!=0) fclose(handle); 
00229       EPETRA_CHK_ERR(-1);
00230     }
00231   }
00232   comm.Broadcast(&M, 1, 0);
00233   comm.Broadcast(&N, 1, 0);
00234   comm.Broadcast(&NZ, 1, 0);
00235 
00236   // Now create matrix using user maps if provided.
00237 
00238 
00239   // Now read in chunks of triplets and broadcast them to other processors.
00240   char *buffer = new char[chunk_read*lineLength];
00241   int nchunk; 
00242   int nmillion = 0;
00243   int nread = 0;
00244   int rlen;
00245 
00246   // Storage for this processor's nonzeros.
00247   const int localblock = 100000;
00248   int localsize = NZ / comm.NumProc() + localblock;
00249   int *iv = (int *) malloc(localsize * sizeof(int));
00250   int *jv = (int *) malloc(localsize * sizeof(int));
00251   double *vv = (double *) malloc(localsize * sizeof(double));
00252   int lnz = 0;   //  Number of non-zeros on this processor.
00253 
00254   if (!iv || !jv || !vv) 
00255     EPETRA_CHK_ERR(-1);
00256 
00257   Epetra_Map *rowMap1;
00258   bool allocatedHere=false;
00259   if (rowMap != 0) 
00260     rowMap1 = const_cast<Epetra_Map *>(rowMap);
00261   else {
00262     rowMap1 = new Epetra_Map(M, 0, comm);
00263     allocatedHere = true;
00264   }
00265   int ioffset = rowMap1->IndexBase()-1;
00266   int joffset = (colMap != 0 ? colMap->IndexBase()-1 : ioffset);
00267 
00268   int rowmajor = 1;  // Assume non-zeros are listed in row-major order, 
00269   int prevrow = -1;  // but test to detect otherwise.  If non-zeros
00270                      // are row-major, we can avoid the sort.
00271 
00272   while (nread < NZ) {
00273     if (NZ-nread > chunk_read) nchunk = chunk_read;
00274     else nchunk = NZ - nread;
00275 
00276     if (me == 0) {
00277       char *eof;
00278       rlen = 0;
00279       for (int i = 0; i < nchunk; i++) {
00280         eof = fgets(&buffer[rlen],lineLength,handle);
00281         if (eof == NULL) {
00282           fprintf(stderr, "%s", "Unexpected end of matrix file.");
00283           EPETRA_CHK_ERR(-1);
00284         }
00285         rlen += strlen(&buffer[rlen]);
00286       }
00287       buffer[rlen++]='\n';
00288     }
00289     comm.Broadcast(&rlen, 1, 0);
00290     comm.Broadcast(buffer, rlen, 0);
00291 
00292     buffer[rlen++] = '\0';
00293     nread += nchunk;
00294 
00295     // Process the received data, saving non-zeros belonging on this proc.
00296     char *lineptr = buffer;
00297     for (rlen = 0; rlen < nchunk; rlen++) {
00298       char *next = strchr(lineptr,'\n');
00299       int I = atoi(strtok(lineptr," \t\n")) + ioffset;
00300       int J = atoi(strtok(NULL," \t\n")) + joffset;
00301       double V = atof(strtok(NULL," \t\n"));
00302       lineptr = next + 1;
00303       if (transpose) {
00304         // swap I and J indices.
00305         int tmp = I;
00306         I = J;
00307         J = tmp;
00308       }
00309       if (rowMap1->MyGID(I)) {
00310         //  This processor keeps this non-zero.
00311         if (lnz >= localsize) {  
00312           // Need more memory to store nonzeros.
00313           localsize += localblock;
00314           iv = (int *) realloc(iv, localsize * sizeof(int));
00315           jv = (int *) realloc(jv, localsize * sizeof(int));
00316           vv = (double *) realloc(vv, localsize * sizeof(double));
00317         }
00318         iv[lnz] = I;
00319         jv[lnz] = J;
00320         vv[lnz] = V;
00321         lnz++;
00322         if (I < prevrow) rowmajor = 0;
00323         prevrow = I;
00324       }
00325     }
00326 
00327     // Status check
00328     if (nread / 1000000 > nmillion) {
00329       nmillion++;
00330       if (verbose && me == 0) cout << nmillion << "M ";
00331     }
00332   }
00333 
00334   delete [] buffer;
00335 
00336   if (!rowmajor) {
00337     // Sort into row-major order (by iv) so can insert entire rows at once.
00338     // Reorder jv and vv to parallel iv.
00339     if (verbose && me == 0) cout << endl << "   Sorting local nonzeros" << endl;
00340     sort_three(iv, jv, vv, 0, lnz-1);
00341   }
00342 
00343   //  Compute number of entries per local row for use in constructor;
00344   //  saves reallocs in FillComplete.
00345 
00346   //  Now construct the matrix.
00347   //  Count number of entries in each row so can use StaticProfile=2.
00348   if (verbose && me == 0) cout << endl << "   Constructing the matrix" << endl;
00349   int numRows = rowMap1->NumMyElements();
00350   int *numNonzerosPerRow = new int[numRows];
00351   for (i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
00352   for (i = 0; i < lnz; i++) 
00353     numNonzerosPerRow[rowMap1->LID(iv[i])]++;
00354 
00355   if (rowMap!=0 && colMap !=0) 
00356     A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
00357   else if (rowMap!=0) {
00358     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00359     // Less memory will be needed in FillComplete.
00360     A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
00361   }
00362   else {
00363     // Construct with StaticProfile=true since we know numNonzerosPerRow.
00364     // Less memory will be needed in FillComplete.
00365     A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
00366   }
00367   A->SetTracebackMode(1);
00368 
00369   // Rows are inserted in ascending global number, and the insertion uses numNonzerosPerRow.
00370   // Therefore numNonzerosPerRow must be permuted, as it was created in ascending local order.
00371   int *gidList = new int[numRows];
00372   for (i = 0; i < numRows; i++) gidList[i] = rowMap1->GID(i);
00373   Epetra_Util::Sort(true,numRows,gidList,0,0,1,&numNonzerosPerRow);
00374   delete [] gidList;
00375 
00376   //  Insert the global values into the matrix row-by-row.
00377   if (verbose && me == 0) cout << "   Inserting global values" << endl;
00378   for (int sum = 0, i = 0; i < numRows; i++) {
00379     if (numNonzerosPerRow[i]) {
00380       int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i], 
00381                                        &vv[sum], &jv[sum]);
00382       if (ierr<0) EPETRA_CHK_ERR(ierr);
00383       sum += numNonzerosPerRow[i];
00384     }
00385   }
00386 
00387   delete [] numNonzerosPerRow;
00388   free(iv);
00389   free(jv);
00390   free(vv);
00391     
00392   if (verbose && me == 0) cout << "   Completing matrix fill" << endl;
00393   if (rangeMap != 0 && domainMap != 0) {
00394     EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
00395   }
00396   else if (M!=N) {
00397     Epetra_Map newDomainMap(N,rowMap1->IndexBase(), comm);
00398     EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
00399   }
00400   else {
00401     EPETRA_CHK_ERR(A->FillComplete());
00402   }
00403 
00404   if (allocatedHere) delete rowMap1;
00405   
00406   if (handle!=0) fclose(handle);
00407   double dt = timer.ElapsedTime();
00408   if (verbose && me == 0) cout << "File Read time (secs):  " << dt << endl;
00409   return(0);
00410 }
00411 
00413 // Sorting values in array list in increasing order. Criteria is int. 
00414 // Also rearrange values in arrays parlista and partlistb
00415 // to match the new order of list. 
00416 
00417 static void quickpart_list_inc_int (
00418   int *list, int *parlista, double *parlistb,
00419   int start, int end, int *equal, int *larger)
00420 {
00421 int i, key, itmp;
00422 double dtmp;
00423 
00424   key = list ? list[(end+start)/2] : 1;
00425 
00426   *equal = *larger = start;
00427   for (i = start; i <= end; i++)
00428     if (list[i] < key) {
00429       itmp                = parlista[i];
00430       parlista[i]         = parlista[*larger];
00431       parlista[(*larger)] = parlista[*equal];
00432       parlista[(*equal)]  = itmp;
00433       dtmp                = parlistb[i];
00434       parlistb[i]         = parlistb[*larger];
00435       parlistb[(*larger)] = parlistb[*equal];
00436       parlistb[(*equal)]  = dtmp;
00437       itmp                = list[i];
00438       list[i]             = list[*larger];
00439       list[(*larger)++]   = list[*equal];
00440       list[(*equal)++]    = itmp;
00441     }
00442     else if (list[i] == key) {
00443       itmp                = parlista[i];
00444       parlista[i]         = parlista[*larger];
00445       parlista[(*larger)] = itmp;
00446       dtmp                = parlistb[i];
00447       parlistb[i]         = parlistb[*larger];
00448       parlistb[(*larger)] = dtmp;
00449       list[i]             = list[*larger];
00450       list[(*larger)++]   = key;
00451     }
00452 }
00453 
00455 static void sort_three(int* list, int *parlista, double *parlistb, 
00456                        int start, int end)
00457 {
00458 int  equal, larger;
00459 
00460   if (start < end) {
00461     quickpart_list_inc_int(list, parlista, parlistb, start, end, 
00462                            &equal, &larger);
00463     sort_three(list, parlista, parlistb, start,  equal-1);
00464     sort_three(list, parlista, parlistb, larger, end);
00465   }
00466 }
00467 
00469 int MatlabFileToCrsMatrix(const char *filename,
00470         const Epetra_Comm & comm,
00471         Epetra_CrsMatrix * & A)
00472 {
00473   const int lineLength = 1025;
00474   char line[lineLength];
00475   int I, J;
00476   double V;
00477 
00478   FILE * handle = 0;
00479 
00480   handle = fopen(filename,"r");  // Open file
00481   if (handle == 0)
00482     EPETRA_CHK_ERR(-1); // file not found
00483 
00484   int numGlobalRows = 0;
00485   int numGlobalCols = 0;
00486   while(fgets(line, lineLength, handle)!=0) {
00487     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00488     if (I>numGlobalRows) numGlobalRows = I;
00489     if (J>numGlobalCols) numGlobalCols = J;
00490   }
00491 
00492   if (handle!=0) fclose(handle);
00493   Epetra_Map rangeMap(numGlobalRows, 0, comm);
00494   Epetra_Map domainMap(numGlobalCols, 0, comm);
00495   A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
00496 
00497   // Now read in each triplet and store to the local portion of the matrix if the row is owned.
00498   const Epetra_Map & rowMap1 = A->RowMap();
00499   
00500   handle = 0;
00501 
00502   handle = fopen(filename,"r");  // Open file
00503   if (handle == 0)
00504     EPETRA_CHK_ERR(-1); // file not found
00505 
00506   while (fgets(line, lineLength, handle)!=0) {
00507     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00508     I--; J--; // Convert to Zero based
00509     if (rowMap1.MyGID(I)) {
00510       int ierr = A->InsertGlobalValues(I, 1, &V, &J);
00511       if (ierr<0) EPETRA_CHK_ERR(ierr);
00512     }
00513   }
00514     
00515   EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
00516 
00517   if (handle!=0) fclose(handle);
00518   return(0);
00519 }
00520 
00521 int HypreFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix){
00522   int MyPID = comm.MyPID();
00523   // This double will be in the format we want for the extension besides the leading zero
00524   double filePID = (double)MyPID/(double)100000;
00525   std::ostringstream stream;
00526   // Using setprecision() puts it in the string
00527   stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
00528   // Then just ignore the first character
00529   std::string fileName(filename);
00530   fileName += stream.str().substr(1,7);
00531   // Open the file
00532   std::ifstream file(fileName.c_str());
00533   string line;
00534   if(file.is_open()){
00535     std::getline(file, line);
00536     int ilower, iupper;
00537     std::istringstream istream(line);
00538     // The first line of the file has the beginning and ending rows
00539     istream >> ilower;
00540     istream >> iupper;
00541     // Using those we can create a row map
00542     Epetra_Map RowMap(-1, iupper-ilower+1, 0, comm);
00543     Matrix = new Epetra_CrsMatrix(Copy, RowMap, 0);
00544     int currRow = -1;
00545     int counter = 0;
00546     std::vector<int> indices;
00547     std::vector<double> values;
00548     while(!file.eof()){
00549       std::getline(file, line);
00550       std::istringstream lineStr(line);
00551       int row, col;
00552       double val;
00553       lineStr >> row;
00554       lineStr >> col;
00555       lineStr >> val;
00556       if(currRow == -1) currRow = row; // First line
00557       if(row == currRow){
00558         // add to the vector
00559         counter = counter + 1;
00560         indices.push_back(col);
00561         values.push_back(val);
00562       } else {
00563         Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
00564         indices.clear();
00565         values.clear();
00566         counter = 0;
00567         currRow = row;
00568         // make a new vector
00569         indices.push_back(col);
00570         values.push_back(val);
00571         counter = counter + 1;
00572       }
00573     }
00574     Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
00575     Matrix->Comm().Barrier();
00576     Matrix->FillComplete();
00577     file.close();
00578     return 0;
00579   } else {
00580     cout << "\nERROR:\nCouldn't open " << fileName << ".\n";
00581     return -1;
00582   }
00583 }
00584 
00585 } // namespace EpetraExt
00586 

Generated on Tue Jul 13 09:23:05 2010 for EpetraExt by  doxygen 1.4.7