EpetraExt Development
EpetraExt_BlockMapIn.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_BlockMapIn.h"
00043 #include "Epetra_Comm.h"
00044 #include "Epetra_Util.h"
00045 #include "Epetra_BlockMap.h"
00046 #include "Epetra_Map.h"
00047 #include "Epetra_IntVector.h"
00048 #include "Epetra_IntSerialDenseVector.h"
00049 #include "Epetra_Import.h"
00050 #include "EpetraExt_mmio.h"
00051 
00052 using namespace EpetraExt;
00053 namespace EpetraExt {
00054 
00055 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00056 int MatrixMarketFileToMap( const char *filename, const Epetra_Comm & comm, Epetra_Map * & map) {
00057 
00058   Epetra_BlockMap * bmap;
00059   if (MatrixMarketFileToBlockMap(filename, comm, bmap)) return(-1);
00060   map = dynamic_cast<Epetra_Map *>(bmap);
00061   return(0);
00062 }
00063 #endif
00064 
00065 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00066 int MatrixMarketFileToMap64( const char *filename, const Epetra_Comm & comm, Epetra_Map * & map) {
00067 
00068   Epetra_BlockMap * bmap;
00069   if (MatrixMarketFileToBlockMap64(filename, comm, bmap)) return(-1);
00070   map = dynamic_cast<Epetra_Map *>(bmap);
00071   return(0);
00072 }
00073 #endif
00074 
00075 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00076 int MatrixMarketFileToBlockMap( const char *filename, const Epetra_Comm & comm, Epetra_BlockMap * & map) {
00077 
00078   const int lineLength = 1025;
00079   char line[lineLength];
00080   char token[lineLength];
00081   int M, N, numProc, MaxElementSize, MinElementSize, NumMyElements, IndexBase, NumGlobalElements, firstGid;
00082 
00083   FILE * handle = 0;
00084 
00085   bool inHeader = true;
00086 
00087   handle = fopen(filename,"r");
00088   if (handle == 0)
00089     EPETRA_CHK_ERR(-1); // file not found
00090 
00091   while (inHeader) {
00092     if(fgets(line, lineLength, handle)==0) return(-1);
00093     if(sscanf(line, "%s", token)==0) return(-1);
00094     if (!strcmp(token, "%NumProc:")) inHeader = false;
00095   }
00096 
00097   if(fgets(line, lineLength, handle)==0) return(-1); // numProc value
00098   if(sscanf(line, "%s %d", token, &numProc)==0) return(-1);
00099 
00100   if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize header line
00101   if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize value
00102   if(sscanf(line, "%s %d", token, &MaxElementSize)==0) return(-1);
00103 
00104   if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize header line
00105   if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize value
00106   if(sscanf(line, "%s %d", token, &MinElementSize)==0) return(-1);
00107 
00108   if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase header line
00109   if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase value
00110   if(sscanf(line, "%s %d", token, &IndexBase)==0) return(-1);
00111 
00112   if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements header line
00113   if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements value
00114   if(sscanf(line, "%s %d", token, &NumGlobalElements)==0) return(-1);
00115 
00116   int ierr = 0;
00117   (void) ierr; // mfh 13 Jan 2013: Forestall compiler warning for unused var.
00118   if (comm.NumProc()==numProc) {
00119     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00120     firstGid = 0;
00121     for (int i=0; i<comm.MyPID(); i++) {
00122       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value
00123       if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00124       firstGid += NumMyElements;
00125     }
00126  
00127     if(fgets(line, lineLength, handle)==0) return(-1); // This PE's NumMyElements value
00128     if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00129 
00130     for (int i=comm.MyPID()+1; i<numProc; i++) {
00131       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00132     }
00133   }
00134   else {
00135     ierr = 1; // Warning error, different number of processors.
00136 
00137     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00138     for (int i=0; i<numProc; i++) {
00139       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00140     }
00141 
00142     NumMyElements = NumGlobalElements/comm.NumProc();
00143     firstGid = comm.MyPID()*NumMyElements;
00144     int remainder = NumGlobalElements%comm.NumProc();
00145     if (comm.MyPID()<remainder) NumMyElements++;
00146     int extra = remainder;
00147     if (comm.MyPID()<remainder) extra = comm.MyPID();
00148     firstGid += extra;
00149   }
00150   if(fgets(line, lineLength, handle)==0) return(-1); // Number of rows, columns
00151   if(sscanf(line, "%d %d", &M, &N)==0) return(-1);
00152 
00153   bool doSizes = (N>1);
00154   Epetra_IntSerialDenseVector v1(NumMyElements);
00155   Epetra_IntSerialDenseVector v2(NumMyElements);
00156   for (int i=0; i<firstGid; i++) {
00157     if(fgets(line, lineLength, handle)==0) return(-1); // dump these
00158   }
00159 
00160   if (doSizes) {
00161     for (int i=0; i<NumMyElements; i++) {
00162       if(fgets(line, lineLength, handle)==0) return(-1);
00163       if(sscanf(line, "%d %d", &v1[i], &v2[i])==0) return(-1); // load v1, v2
00164     }
00165   }
00166   else {
00167     for (int i=0; i<NumMyElements; i++) {
00168       if(fgets(line, lineLength, handle)==0) return(-1);
00169       if(sscanf(line, "%d", &v1[i])==0) return(-1); // load v1
00170       v2[i] = MinElementSize; // Fill with constant size
00171     }
00172   }
00173   if (fclose(handle)) return(-1);
00174 
00175   comm.Barrier();
00176 
00177   if (MinElementSize==1 && MaxElementSize==1)
00178     map = new Epetra_Map(-1, NumMyElements, v1.Values(), IndexBase, comm);
00179   else
00180     map = new Epetra_BlockMap(-1, NumMyElements, v1.Values(), v2.Values(), IndexBase, comm);
00181   return(0);
00182 }
00183 #endif
00184 
00185 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00186 int MatrixMarketFileToBlockMap64( const char *filename, const Epetra_Comm & comm, Epetra_BlockMap * & map) {
00187 
00188   const int lineLength = 1025;
00189   char line[lineLength];
00190   char token[lineLength];
00191   int numProc, MaxElementSize, MinElementSize, NumMyElements;
00192   long long M, N, IndexBase, NumGlobalElements, firstGid;
00193 
00194   FILE * handle = 0;
00195 
00196   bool inHeader = true;
00197 
00198   handle = fopen(filename,"r");
00199   if (handle == 0)
00200     EPETRA_CHK_ERR(-1); // file not found
00201 
00202   while (inHeader) {
00203     if(fgets(line, lineLength, handle)==0) return(-1);
00204     if(sscanf(line, "%s", token)==0) return(-1);
00205     if (!strcmp(token, "%NumProc:")) inHeader = false;
00206   }
00207 
00208   if(fgets(line, lineLength, handle)==0) return(-1); // numProc value
00209   if(sscanf(line, "%s %d", token, &numProc)==0) return(-1);
00210 
00211   if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize header line
00212   if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize value
00213   if(sscanf(line, "%s %d", token, &MaxElementSize)==0) return(-1);
00214 
00215   if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize header line
00216   if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize value
00217   if(sscanf(line, "%s %d", token, &MinElementSize)==0) return(-1);
00218 
00219   if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase header line
00220   if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase value
00221   if(sscanf(line, "%s %lld", token, &IndexBase)==0) return(-1);
00222 
00223   if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements header line
00224   if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements value
00225   if(sscanf(line, "%s %lld", token, &NumGlobalElements)==0) return(-1);
00226 
00227   int ierr = 0;
00228   (void) ierr; // mfh 13 Jan 2013: Forestall compiler warning for unused var.
00229   if (comm.NumProc()==numProc) {
00230     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00231     firstGid = 0;
00232     for (int i=0; i<comm.MyPID(); i++) {
00233       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value
00234       if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00235       firstGid += NumMyElements;
00236     }
00237  
00238     if(fgets(line, lineLength, handle)==0) return(-1); // This PE's NumMyElements value
00239     if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00240 
00241     for (int i=comm.MyPID()+1; i<numProc; i++) {
00242       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00243     }
00244   }
00245   else {
00246     ierr = 1; // Warning error, different number of processors.
00247 
00248     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00249     for (int i=0; i<numProc; i++) {
00250       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00251     }
00252 
00253     NumMyElements = (int) (NumGlobalElements/comm.NumProc());
00254     firstGid = ((long long) comm.MyPID())*NumMyElements;
00255     int remainder = (int) (NumGlobalElements%comm.NumProc());
00256     if (comm.MyPID()<remainder) NumMyElements++;
00257     int extra = remainder;
00258     if (comm.MyPID()<remainder) extra = comm.MyPID();
00259     firstGid += extra;
00260   }
00261   if(fgets(line, lineLength, handle)==0) return(-1); // Number of rows, columns
00262   if(sscanf(line, "%lld %lld", &M, &N)==0) return(-1);
00263 
00264   bool doSizes = (N>1);
00265   Epetra_LongLongSerialDenseVector v1(NumMyElements);
00266   Epetra_IntSerialDenseVector v2(NumMyElements);
00267   for (long long i=0; i<firstGid; i++) {
00268     if(fgets(line, lineLength, handle)==0) return(-1); // dump these
00269   }
00270 
00271   if (doSizes) {
00272     for (int i=0; i<NumMyElements; i++) {
00273       if(fgets(line, lineLength, handle)==0) return(-1);
00274       if(sscanf(line, "%lld %d", &v1[i], &v2[i])==0) return(-1); // load v1, v2
00275     }
00276   }
00277   else {
00278     for (int i=0; i<NumMyElements; i++) {
00279       if(fgets(line, lineLength, handle)==0) return(-1);
00280       if(sscanf(line, "%lld", &v1[i])==0) return(-1); // load v1
00281       v2[i] = MinElementSize; // Fill with constant size
00282     }
00283   }
00284   if (fclose(handle)) return(-1);
00285 
00286   comm.Barrier();
00287 
00288   if (MinElementSize==1 && MaxElementSize==1)
00289     map = new Epetra_Map(-1LL, NumMyElements, v1.Values(), IndexBase, comm);
00290   else
00291     map = new Epetra_BlockMap(-1LL, NumMyElements, v1.Values(), v2.Values(), IndexBase, comm);
00292   return(0);
00293 }
00294 #endif
00295 
00296 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00297 int MatrixMarketFileToRowMap(const char* filename,
00298                              const Epetra_Comm& comm,
00299                              Epetra_BlockMap*& rowmap)
00300 {
00301   FILE* infile = fopen(filename, "r");
00302   MM_typecode matcode;
00303 
00304   int err = mm_read_banner(infile, &matcode);
00305   if (err != 0) return(err);
00306 
00307   if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
00308       !mm_is_real(matcode)   || !mm_is_general(matcode)) {
00309     return(-1);
00310   }
00311 
00312   int numrows, numcols;
00313   err = mm_read_mtx_array_size(infile, &numrows, &numcols);
00314   if (err != 0) return(err);
00315 
00316   fclose(infile);
00317 
00318   rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
00319   return(0);
00320 }
00321 #endif
00322 
00323 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00324 int MatrixMarketFileToBlockMaps(const char* filename,
00325                                 const Epetra_Comm& comm,
00326                                 Epetra_BlockMap*& rowmap,
00327                                 Epetra_BlockMap*& colmap,
00328                                 Epetra_BlockMap*& rangemap,
00329                                 Epetra_BlockMap*& domainmap)
00330 {
00331   FILE* infile = fopen(filename, "r");
00332   if (infile == NULL) {
00333     return(-1);
00334   }
00335 
00336   MM_typecode matcode;
00337 
00338   int err = mm_read_banner(infile, &matcode);
00339   if (err != 0) return(err);
00340 
00341   if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
00342       !mm_is_real(matcode)   || !mm_is_general(matcode)) {
00343     return(-1);
00344   }
00345 
00346   int numrows, numcols, nnz;
00347   err = mm_read_mtx_crd_size(infile, &numrows, &numcols, &nnz);
00348   if (err != 0) return(err);
00349 
00350   //for this case, we'll assume that the row-map is the same as
00351   //the range-map.
00352   //create row-map and range-map with linear distributions.
00353 
00354   rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
00355   rangemap = new Epetra_BlockMap(numrows, 1, 0, comm);
00356 
00357   int I, J;
00358   double val, imag;
00359 
00360   int num_map_cols = 0, insertPoint, foundOffset;
00361   int allocLen = numcols;
00362   int* map_cols = new int[allocLen];
00363 
00364   //read through all matrix data and construct a list of the column-
00365   //indices that occur in rows that are local to this processor.
00366  
00367   for(int i=0; i<nnz; ++i) {
00368     err = mm_read_mtx_crd_entry(infile, &I, &J, &val,
00369                                 &imag, matcode);
00370 
00371     if (err == 0) {
00372       --I;
00373       --J;
00374       if (rowmap->MyGID(I)) {
00375         foundOffset = Epetra_Util_binary_search(J, map_cols, num_map_cols,
00376                                                 insertPoint);
00377         if (foundOffset < 0) {
00378           Epetra_Util_insert(J, insertPoint, map_cols,
00379                              num_map_cols, allocLen);
00380         }
00381       }
00382     } 
00383   }
00384 
00385   //create colmap with the list of columns associated with rows that are
00386   //local to this processor.
00387   colmap = new Epetra_Map(-1, num_map_cols, map_cols, 0, comm);
00388 
00389   //create domainmap which has a linear distribution
00390   domainmap = new Epetra_BlockMap(numcols, 1, 0, comm);
00391 
00392   delete [] map_cols;
00393 
00394   return(0);
00395 }
00396 #endif
00397 
00398 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00399 int MatrixMarketFileToBlockMaps64(const char* filename,
00400                                 const Epetra_Comm& comm,
00401                                 Epetra_BlockMap*& rowmap,
00402                                 Epetra_BlockMap*& colmap,
00403                                 Epetra_BlockMap*& rangemap,
00404                                 Epetra_BlockMap*& domainmap)
00405 {
00406   FILE* infile = fopen(filename, "r");
00407   if (infile == NULL) {
00408     return(-1);
00409   }
00410 
00411   MM_typecode matcode;
00412 
00413   int err = mm_read_banner(infile, &matcode);
00414   if (err != 0) return(err);
00415 
00416   if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
00417       !mm_is_real(matcode)   || !mm_is_general(matcode)) {
00418     return(-1);
00419   }
00420 
00421   long long numrows, numcols, nnz;
00422   err = mm_read_mtx_crd_size(infile, &numrows, &numcols, &nnz);
00423   if (err != 0) return(err);
00424 
00425   //for this case, we'll assume that the row-map is the same as
00426   //the range-map.
00427   //create row-map and range-map with linear distributions.
00428 
00429   rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
00430   rangemap = new Epetra_BlockMap(numrows, 1, 0, comm);
00431 
00432   long long I, J;
00433   double val, imag;
00434 
00435   int num_map_cols = 0, insertPoint, foundOffset;
00436   int allocLen = 0;
00437   if(numcols > (long long) std::numeric_limits<int>::max())
00438     allocLen = std::numeric_limits<int>::max();
00439   else
00440     allocLen = (int) numcols;
00441 
00442   long long* map_cols = new long long[allocLen];
00443 
00444   //read through all matrix data and construct a list of the column-
00445   //indices that occur in rows that are local to this processor.
00446  
00447   for(int i=0; i<nnz; ++i) {
00448     err = mm_read_mtx_crd_entry(infile, &I, &J, &val,
00449                                 &imag, matcode);
00450 
00451     if (err == 0) {
00452       --I;
00453       --J;
00454       if (rowmap->MyGID(I)) {
00455         foundOffset = Epetra_Util_binary_search(J, map_cols, num_map_cols,
00456                                                 insertPoint);
00457         if (foundOffset < 0) {
00458           Epetra_Util_insert(J, insertPoint, map_cols,
00459                              num_map_cols, allocLen);
00460         }
00461       }
00462     } 
00463   }
00464 
00465   //create colmap with the list of columns associated with rows that are
00466   //local to this processor.
00467   colmap = new Epetra_Map(-1LL, num_map_cols, map_cols, 0, comm);
00468 
00469   //create domainmap which has a linear distribution
00470   domainmap = new Epetra_BlockMap(numcols, 1, 0, comm);
00471 
00472   delete [] map_cols;
00473 
00474   return(0);
00475 }
00476 #endif
00477 
00478 } // namespace EpetraExt
00479 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines