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

Generated on Wed May 12 21:24:45 2010 for EpetraExt by  doxygen 1.4.7