EpetraExt Package Browser (Single Doxygen Collection) 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   if (comm.NumProc()==numProc) {
00104     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00105     firstGid = 0;
00106     for (int i=0; i<comm.MyPID(); i++) {
00107       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value
00108       if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00109       firstGid += NumMyElements;
00110     }
00111  
00112     if(fgets(line, lineLength, handle)==0) return(-1); // This PE's NumMyElements value
00113     if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
00114 
00115     for (int i=comm.MyPID()+1; i<numProc; i++) {
00116       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00117     }
00118   }
00119   else {
00120     ierr = 1; // Warning error, different number of processors.
00121 
00122     if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
00123     for (int i=0; i<numProc; i++) {
00124       if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
00125     }
00126 
00127     NumMyElements = NumGlobalElements/comm.NumProc();
00128     firstGid = comm.MyPID()*NumMyElements;
00129     int remainder = NumGlobalElements%comm.NumProc();
00130     if (comm.MyPID()<remainder) NumMyElements++;
00131     int extra = remainder;
00132     if (comm.MyPID()<remainder) extra = comm.MyPID();
00133     firstGid += extra;
00134   }
00135   if(fgets(line, lineLength, handle)==0) return(-1); // Number of rows, columns
00136   if(sscanf(line, "%d %d", &M, &N)==0) return(-1);
00137 
00138   bool doSizes = (N>1);
00139   Epetra_IntSerialDenseVector v1(NumMyElements);
00140   Epetra_IntSerialDenseVector v2(NumMyElements);
00141   for (int i=0; i<firstGid; i++) {
00142     if(fgets(line, lineLength, handle)==0) return(-1); // dump these
00143   }
00144 
00145   if (doSizes) {
00146     for (int i=0; i<NumMyElements; i++) {
00147       if(fgets(line, lineLength, handle)==0) return(-1);
00148       if(sscanf(line, "%d %d", &v1[i], &v2[i])==0) return(-1); // load v1, v2
00149     }
00150   }
00151   else {
00152     for (int i=0; i<NumMyElements; i++) {
00153       if(fgets(line, lineLength, handle)==0) return(-1);
00154       if(sscanf(line, "%d", &v1[i])==0) return(-1); // load v1
00155       v2[i] = MinElementSize; // Fill with constant size
00156     }
00157   }
00158   if (fclose(handle)) return(-1);
00159 
00160   comm.Barrier();
00161 
00162   if (MinElementSize==1 && MaxElementSize==1)
00163     map = new Epetra_Map(-1, NumMyElements, v1.Values(), IndexBase, comm);
00164   else
00165     map = new Epetra_BlockMap(-1, NumMyElements, v1.Values(), v2.Values(), IndexBase, comm);
00166   return(0);
00167 }
00168 
00169 int MatrixMarketFileToRowMap(const char* filename,
00170                              const Epetra_Comm& comm,
00171                              Epetra_BlockMap*& rowmap)
00172 {
00173   FILE* infile = fopen(filename, "r");
00174   MM_typecode matcode;
00175 
00176   int err = mm_read_banner(infile, &matcode);
00177   if (err != 0) return(err);
00178 
00179   if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
00180       !mm_is_real(matcode)   || !mm_is_general(matcode)) {
00181     return(-1);
00182   }
00183 
00184   int numrows, numcols;
00185   err = mm_read_mtx_array_size(infile, &numrows, &numcols);
00186   if (err != 0) return(err);
00187 
00188   fclose(infile);
00189 
00190   rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
00191   return(0);
00192 }
00193 
00194 int MatrixMarketFileToBlockMaps(const char* filename,
00195                                 const Epetra_Comm& comm,
00196                                 Epetra_BlockMap*& rowmap,
00197                                 Epetra_BlockMap*& colmap,
00198                                 Epetra_BlockMap*& rangemap,
00199                                 Epetra_BlockMap*& domainmap)
00200 {
00201   FILE* infile = fopen(filename, "r");
00202   if (infile == NULL) {
00203     return(-1);
00204   }
00205 
00206   MM_typecode matcode;
00207 
00208   int err = mm_read_banner(infile, &matcode);
00209   if (err != 0) return(err);
00210 
00211   if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
00212       !mm_is_real(matcode)   || !mm_is_general(matcode)) {
00213     return(-1);
00214   }
00215 
00216   int numrows, numcols, nnz;
00217   err = mm_read_mtx_crd_size(infile, &numrows, &numcols, &nnz);
00218   if (err != 0) return(err);
00219 
00220   //for this case, we'll assume that the row-map is the same as
00221   //the range-map.
00222   //create row-map and range-map with linear distributions.
00223 
00224   rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
00225   rangemap = new Epetra_BlockMap(numrows, 1, 0, comm);
00226 
00227   int I, J;
00228   double val, imag;
00229 
00230   int num_map_cols = 0, insertPoint, foundOffset;
00231   int allocLen = numcols;
00232   int* map_cols = new int[allocLen];
00233 
00234   //read through all matrix data and construct a list of the column-
00235   //indices that occur in rows that are local to this processor.
00236  
00237   for(int i=0; i<nnz; ++i) {
00238     err = mm_read_mtx_crd_entry(infile, &I, &J, &val,
00239                                 &imag, matcode);
00240 
00241     if (err == 0) {
00242       --I;
00243       --J;
00244       if (rowmap->MyGID(I)) {
00245         foundOffset = Epetra_Util_binary_search(J, map_cols, num_map_cols,
00246                                                 insertPoint);
00247         if (foundOffset < 0) {
00248           Epetra_Util_insert(J, insertPoint, map_cols,
00249                              num_map_cols, allocLen);
00250         }
00251       }
00252     } 
00253   }
00254 
00255   //create colmap with the list of columns associated with rows that are
00256   //local to this processor.
00257   colmap = new Epetra_Map(-1, num_map_cols, map_cols, 0, comm);
00258 
00259   //create domainmap which has a linear distribution
00260   domainmap = new Epetra_BlockMap(numcols, 1, 0, comm);
00261 
00262   delete [] map_cols;
00263 
00264   return(0);
00265 }
00266 
00267 } // namespace EpetraExt
00268 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines