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 "EpetraExt_BlockMapIn.h"
00042 #include "Epetra_Comm.h"
00043 #include "Epetra_Util.h"
00044 #include "Epetra_BlockMap.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_IntVector.h"
00047 #include "Epetra_IntSerialDenseVector.h"
00048 #include "Epetra_Import.h"
00049 #include "EpetraExt_mmio.h"
00050 
00051 using namespace EpetraExt;
00052 namespace EpetraExt {
00053 
00054 int MatrixMarketFileToMap( const char *filename, const Epetra_Comm & comm, Epetra_Map * & map) {
00055 
00056   Epetra_BlockMap * bmap;
00057   if (MatrixMarketFileToBlockMap(filename, comm, bmap)) return(-1);
00058   map = dynamic_cast<Epetra_Map *>(bmap);
00059   return(0);
00060 }
00061 
00062 int MatrixMarketFileToBlockMap( const char *filename, const Epetra_Comm & comm, Epetra_BlockMap * & map) {
00063 
00064   const int lineLength = 1025;
00065   char line[lineLength];
00066   char token[lineLength];
00067   int M, N, numProc, MaxElementSize, MinElementSize, NumMyElements, IndexBase, NumGlobalElements, firstGid;
00068 
00069   FILE * handle = 0;
00070 
00071   bool inHeader = true;
00072 
00073   handle = fopen(filename,"r");
00074   if (handle == 0)
00075     EPETRA_CHK_ERR(-1); // file not found
00076 
00077   while (inHeader) {
00078     if(fgets(line, lineLength, handle)==0) return(-1);
00079     if(sscanf(line, "%s", token)==0) return(-1);
00080     if (!strcmp(token, "%NumProc:")) inHeader = false;
00081   }
00082 
00083   if(fgets(line, lineLength, handle)==0) return(-1); // numProc value
00084   if(sscanf(line, "%s %d", token, &numProc)==0) return(-1);
00085 
00086   if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize header line
00087   if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize value
00088   if(sscanf(line, "%s %d", token, &MaxElementSize)==0) return(-1);
00089 
00090   if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize header line
00091   if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize value
00092   if(sscanf(line, "%s %d", token, &MinElementSize)==0) return(-1);
00093 
00094   if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase header line
00095   if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase value
00096   if(sscanf(line, "%s %d", token, &IndexBase)==0) return(-1);
00097 
00098   if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements header line
00099   if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements value
00100   if(sscanf(line, "%s %d", token, &NumGlobalElements)==0) return(-1);
00101 
00102   int ierr = 0;
00103   (void) ierr; // mfh 13 Jan 2013: Forestall compiler warning for unused var.
00104   if (comm.NumProc()==numProc) {
00105     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00106     firstGid = 0;
00107     for (int i=0; i<comm.MyPID(); i++) {
00108       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value
00109       if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00110       firstGid += NumMyElements;
00111     }
00112  
00113     if(fgets(line, lineLength, handle)==0) return(-1); // This PE's NumMyElements value
00114     if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00115 
00116     for (int i=comm.MyPID()+1; i<numProc; i++) {
00117       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00118     }
00119   }
00120   else {
00121     ierr = 1; // Warning error, different number of processors.
00122 
00123     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00124     for (int i=0; i<numProc; i++) {
00125       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00126     }
00127 
00128     NumMyElements = NumGlobalElements/comm.NumProc();
00129     firstGid = comm.MyPID()*NumMyElements;
00130     int remainder = NumGlobalElements%comm.NumProc();
00131     if (comm.MyPID()<remainder) NumMyElements++;
00132     int extra = remainder;
00133     if (comm.MyPID()<remainder) extra = comm.MyPID();
00134     firstGid += extra;
00135   }
00136   if(fgets(line, lineLength, handle)==0) return(-1); // Number of rows, columns
00137   if(sscanf(line, "%d %d", &M, &N)==0) return(-1);
00138 
00139   bool doSizes = (N>1);
00140   Epetra_IntSerialDenseVector v1(NumMyElements);
00141   Epetra_IntSerialDenseVector v2(NumMyElements);
00142   for (int i=0; i<firstGid; i++) {
00143     if(fgets(line, lineLength, handle)==0) return(-1); // dump these
00144   }
00145 
00146   if (doSizes) {
00147     for (int i=0; i<NumMyElements; i++) {
00148       if(fgets(line, lineLength, handle)==0) return(-1);
00149       if(sscanf(line, "%d %d", &v1[i], &v2[i])==0) return(-1); // load v1, v2
00150     }
00151   }
00152   else {
00153     for (int i=0; i<NumMyElements; i++) {
00154       if(fgets(line, lineLength, handle)==0) return(-1);
00155       if(sscanf(line, "%d", &v1[i])==0) return(-1); // load v1
00156       v2[i] = MinElementSize; // Fill with constant size
00157     }
00158   }
00159   if (fclose(handle)) return(-1);
00160 
00161   comm.Barrier();
00162 
00163   if (MinElementSize==1 && MaxElementSize==1)
00164     map = new Epetra_Map(-1, NumMyElements, v1.Values(), IndexBase, comm);
00165   else
00166     map = new Epetra_BlockMap(-1, NumMyElements, v1.Values(), v2.Values(), IndexBase, comm);
00167   return(0);
00168 }
00169 
00170 int MatrixMarketFileToRowMap(const char* filename,
00171                              const Epetra_Comm& comm,
00172                              Epetra_BlockMap*& rowmap)
00173 {
00174   FILE* infile = fopen(filename, "r");
00175   MM_typecode matcode;
00176 
00177   int err = mm_read_banner(infile, &matcode);
00178   if (err != 0) return(err);
00179 
00180   if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
00181       !mm_is_real(matcode)   || !mm_is_general(matcode)) {
00182     return(-1);
00183   }
00184 
00185   int numrows, numcols;
00186   err = mm_read_mtx_array_size(infile, &numrows, &numcols);
00187   if (err != 0) return(err);
00188 
00189   fclose(infile);
00190 
00191   rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
00192   return(0);
00193 }
00194 
00195 int MatrixMarketFileToBlockMaps(const char* filename,
00196                                 const Epetra_Comm& comm,
00197                                 Epetra_BlockMap*& rowmap,
00198                                 Epetra_BlockMap*& colmap,
00199                                 Epetra_BlockMap*& rangemap,
00200                                 Epetra_BlockMap*& domainmap)
00201 {
00202   FILE* infile = fopen(filename, "r");
00203   if (infile == NULL) {
00204     return(-1);
00205   }
00206 
00207   MM_typecode matcode;
00208 
00209   int err = mm_read_banner(infile, &matcode);
00210   if (err != 0) return(err);
00211 
00212   if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
00213       !mm_is_real(matcode)   || !mm_is_general(matcode)) {
00214     return(-1);
00215   }
00216 
00217   int numrows, numcols, nnz;
00218   err = mm_read_mtx_crd_size(infile, &numrows, &numcols, &nnz);
00219   if (err != 0) return(err);
00220 
00221   //for this case, we'll assume that the row-map is the same as
00222   //the range-map.
00223   //create row-map and range-map with linear distributions.
00224 
00225   rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
00226   rangemap = new Epetra_BlockMap(numrows, 1, 0, comm);
00227 
00228   int I, J;
00229   double val, imag;
00230 
00231   int num_map_cols = 0, insertPoint, foundOffset;
00232   int allocLen = numcols;
00233   int* map_cols = new int[allocLen];
00234 
00235   //read through all matrix data and construct a list of the column-
00236   //indices that occur in rows that are local to this processor.
00237  
00238   for(int i=0; i<nnz; ++i) {
00239     err = mm_read_mtx_crd_entry(infile, &I, &J, &val,
00240                                 &imag, matcode);
00241 
00242     if (err == 0) {
00243       --I;
00244       --J;
00245       if (rowmap->MyGID(I)) {
00246         foundOffset = Epetra_Util_binary_search(J, map_cols, num_map_cols,
00247                                                 insertPoint);
00248         if (foundOffset < 0) {
00249           Epetra_Util_insert(J, insertPoint, map_cols,
00250                              num_map_cols, allocLen);
00251         }
00252       }
00253     } 
00254   }
00255 
00256   //create colmap with the list of columns associated with rows that are
00257   //local to this processor.
00258   colmap = new Epetra_Map(-1, num_map_cols, map_cols, 0, comm);
00259 
00260   //create domainmap which has a linear distribution
00261   domainmap = new Epetra_BlockMap(numcols, 1, 0, comm);
00262 
00263   delete [] map_cols;
00264 
00265   return(0);
00266 }
00267 
00268 } // namespace EpetraExt
00269 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines