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_Util.h"
00036 
00037 using namespace EpetraExt;
00038 namespace EpetraExt {
00039 
00040 int MatrixMarketFileToCrsMatrix( const char *filename, const Epetra_Comm & comm, Epetra_CrsMatrix * & A) {
00041     
00042   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A));
00043   return(0);
00044 }
00045   
00046 int MatrixMarketFileToCrsMatrix(const char *filename, const Epetra_Map & rowMap,
00047           const Epetra_Map& rangeMap, const Epetra_Map& domainMap, Epetra_CrsMatrix * & A) {
00048   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, rowMap.Comm(), A, &rowMap, 0, &rangeMap, &domainMap));
00049   return(0);
00050 }
00051 
00052 int MatrixMarketFileToCrsMatrix( const char *filename, const Epetra_Map & rowMap, Epetra_CrsMatrix * & A) {
00053 
00054   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, rowMap.Comm(), A, &rowMap));
00055   return(0);
00056 }
00057 
00058 int MatrixMarketFileToCrsMatrix( const char *filename, const Epetra_Map & rowMap, const Epetra_Map & colMap, Epetra_CrsMatrix * & A) {
00059 
00060   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, rowMap.Comm(), A, &rowMap, &colMap));
00061   return(0);
00062 }
00063 
00064 int MatrixMarketFileToCrsMatrix(const char *filename, const Epetra_Map & rowMap, const Epetra_Map & colMap,
00065                                 const Epetra_Map& rangeMap, const Epetra_Map& domainMap, Epetra_CrsMatrix * & A) {
00066   EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, rowMap.Comm(), A, &rowMap, &colMap, &rangeMap, &domainMap));
00067   return(0);
00068 }
00069 
00070 int MatrixMarketFileToCrsMatrixHandle(const char *filename,
00071               const Epetra_Comm & comm,
00072                                       Epetra_CrsMatrix * & A,
00073               const Epetra_Map * rowMap,
00074               const Epetra_Map * colMap,
00075                                       const Epetra_Map * rangeMap,
00076                                       const Epetra_Map * domainMap)
00077 {
00078 
00079   const int lineLength = 1025;
00080   const int tokenLength = 35;
00081   char line[lineLength];
00082   char token1[tokenLength];
00083   char token2[tokenLength];
00084   char token3[tokenLength];
00085   char token4[tokenLength];
00086   char token5[tokenLength];
00087   int M, N, NZ;
00088 
00089   // Make sure domain and range maps are either both defined or undefined 
00090   if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {EPETRA_CHK_ERR(-3);}
00091 
00092   // check maps to see if domain and range are 1-to-1
00093 
00094   if (domainMap!=0) {
00095     if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00096     if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00097   }
00098   else {
00099     // If domain and range are not specified, rowMap becomes both and must be 1-to-1 if specified
00100     if (rowMap!=0) {
00101       if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00102     }
00103   }
00104       
00105 
00106   FILE * handle = 0;
00107 
00108   handle = fopen(filename,"r");  // Open file
00109   if (handle == 0)
00110     EPETRA_CHK_ERR(-1); // file not found
00111 
00112   // Check first line, which should be "%%MatrixMarket matrix coordinate real general" (without quotes)
00113   if(fgets(line, lineLength, handle)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00114   if(sscanf(line, "%s %s %s %s %s", token1, token2, token3, token4, token5 )==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00115   if (strcmp(token1, "%%MatrixMarket") ||
00116       strcmp(token2, "matrix") ||
00117       strcmp(token3, "coordinate") ||
00118       strcmp(token4, "real") ||
00119       strcmp(token5, "general")) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00120 
00121   // Next, strip off header lines (which start with "%")
00122   do {
00123     if(fgets(line, lineLength, handle)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00124   } while (line[0] == '%');
00125 
00126   // Next get problem dimensions: M, N, NZ
00127   if(sscanf(line, "%d %d %d", &M, &N, &NZ)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00128 
00129   // Now create matrix using user maps if provided.
00130 
00131   if (rowMap!=0 && colMap !=0)
00132     A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, 0);
00133   else if (rowMap!=0)
00134     A = new Epetra_CrsMatrix(Copy, *rowMap, 0);
00135   else {
00136     Epetra_Map newRowMap(M,0, comm);
00137     A = new Epetra_CrsMatrix(Copy, newRowMap, 0);
00138   }
00139   // Now read in each triplet and store to the local portion of the matrix if the row is owned.
00140   int I, J;
00141   double V;
00142   const Epetra_Map & rowMap1 = A->RowMap();
00143   const Epetra_Map & colMap1 = A->ColMap();
00144   int ioffset = rowMap1.IndexBase()-1;
00145   int joffset = colMap1.IndexBase()-1;
00146   
00147   for (int i=0; i<NZ; i++) {
00148     if(fgets(line, lineLength, handle)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);};
00149     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00150     I+=ioffset; J+=joffset; // Convert to Zero based
00151     if (rowMap1.MyGID(I)) {
00152       int ierr = A->InsertGlobalValues(I, 1, &V, &J);
00153       if (ierr<0) EPETRA_CHK_ERR(ierr);
00154     }
00155   }
00156     
00157   if (rangeMap != 0 && domainMap != 0) {
00158     A->FillComplete(*domainMap, *rangeMap);
00159   }
00160   else if (M!=N) {
00161     Epetra_Map newDomainMap(N,rowMap1.IndexBase(), comm);
00162     A->FillComplete(newDomainMap, rowMap1);
00163   }
00164   else {
00165     A->FillComplete();
00166   }
00167   
00168   if (handle!=0) fclose(handle);
00169   return(0);
00170 }
00171 
00172 int MatlabFileToCrsMatrix(const char *filename,
00173         const Epetra_Comm & comm,
00174         Epetra_CrsMatrix * & A)
00175 {
00176   const int lineLength = 1025;
00177   char line[lineLength];
00178   int I, J;
00179   double V;
00180 
00181   FILE * handle = 0;
00182 
00183   handle = fopen(filename,"r");  // Open file
00184   if (handle == 0)
00185     EPETRA_CHK_ERR(-1); // file not found
00186 
00187   int numGlobalRows = 0;
00188   int numGlobalCols = 0;
00189   while(fgets(line, lineLength, handle)!=0) {
00190     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00191     if (I>numGlobalRows) numGlobalRows = I;
00192     if (J>numGlobalCols) numGlobalCols = J;
00193   }
00194 
00195   if (handle!=0) fclose(handle);
00196   Epetra_Map rangeMap(numGlobalRows, 0, comm);
00197   Epetra_Map domainMap(numGlobalCols, 0, comm);
00198   A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
00199 
00200   // Now read in each triplet and store to the local portion of the matrix if the row is owned.
00201   const Epetra_Map & rowMap1 = A->RowMap();
00202   
00203   handle = 0;
00204 
00205   handle = fopen(filename,"r");  // Open file
00206   if (handle == 0)
00207     EPETRA_CHK_ERR(-1); // file not found
00208 
00209   while (fgets(line, lineLength, handle)!=0) {
00210     if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00211     I--; J--; // Convert to Zero based
00212     if (rowMap1.MyGID(I)) {
00213       int ierr = A->InsertGlobalValues(I, 1, &V, &J);
00214       if (ierr<0) EPETRA_CHK_ERR(ierr);
00215     }
00216   }
00217     
00218   EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
00219 
00220   if (handle!=0) fclose(handle);
00221   return(0);
00222 }
00223 } // namespace EpetraExt
00224 

Generated on Tue Oct 20 12:45:29 2009 for EpetraExt by doxygen 1.4.7