EpetraExt_BlockMapOut.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_BlockMapOut.h"
00029 #include "EpetraExt_mmio.h"
00030 #include "Epetra_Comm.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 
00037 using namespace EpetraExt;
00038 namespace EpetraExt {
00039 
00040 int BlockMapToMatrixMarketFile( const char *filename, const Epetra_BlockMap & map, 
00041          const char * mapName,
00042          const char *mapDescription, 
00043          bool writeHeader) {
00044   int M = map.NumGlobalElements();
00045   int N = 1;
00046   if (map.MaxElementSize()>1) N = 2; // Non-trivial block map, store element sizes in second column
00047 
00048   FILE * handle = 0;
00049 
00050   if (map.Comm().MyPID()==0) { // Only PE 0 does this section
00051 
00052     handle = fopen(filename,"w");
00053     if (!handle) return(-1);
00054     MM_typecode matcode;
00055     mm_initialize_typecode(&matcode);
00056     mm_set_matrix(&matcode);
00057     mm_set_array(&matcode);
00058     mm_set_integer(&matcode);
00059 
00060     if (writeHeader==true) { // Only write header if requested (true by default)
00061     
00062       if (mm_write_banner(handle, matcode)) return(-1);
00063       
00064       if (mapName!=0) fprintf(handle, "%% \n%% %s\n", mapName);
00065       if (mapDescription!=0) fprintf(handle, "%% %s\n%% \n", mapDescription);
00066 
00067     }
00068   }
00069     
00070   if (writeHeader==true) { // Only write header if requested (true by default)
00071 
00072     // Make an Epetra_IntVector of length numProc such that all elements are on PE 0 and
00073     // the ith element is NumMyElements from the ith PE
00074 
00075     Epetra_Map map1(-1, 1, 0, map.Comm()); // map with one element on each processor
00076     int length = 0;
00077     if (map.Comm().MyPID()==0) length = map.Comm().NumProc();
00078     Epetra_Map map2(-1, length, 0, map.Comm());
00079     Epetra_Import lengthImporter(map2, map1);
00080     Epetra_IntVector v1(map1);
00081     Epetra_IntVector v2(map2);
00082     v1[0] = map.NumMyElements();
00083     if (v2.Import(v1, lengthImporter, Insert)) return(-1);
00084     if (map.Comm().MyPID()==0) { 
00085       fprintf(handle, "%%Format Version:\n");
00086       //int version = 1; // We may change the format scheme at a later date.
00087       fprintf(handle, "%% %d \n", map.Comm().NumProc());
00088      fprintf(handle, "%%NumProc: Number of processors:\n");
00089       fprintf(handle, "%% %d \n", map.Comm().NumProc());
00090       fprintf(handle, "%%MaxElementSize: Maximum element size:\n");
00091       fprintf(handle, "%% %d \n", map.MaxElementSize());
00092       fprintf(handle, "%%MinElementSize: Minimum element size:\n");
00093       fprintf(handle, "%% %d \n", map.MinElementSize());
00094       fprintf(handle, "%%IndexBase: Index base of map:\n");
00095       fprintf(handle, "%% %d \n", map.IndexBase());
00096       fprintf(handle, "%%NumGlobalElements: Total number of GIDs in map:\n");
00097       fprintf(handle, "%% %d \n", map.NumGlobalElements());
00098       fprintf(handle, "%%NumMyElements: BlockMap lengths per processor:\n");
00099       for ( int i=0; i< v2.MyLength(); i++) fprintf(handle, "%% %d\n", v2[i]);
00100       
00101       if (mm_write_mtx_array_size(handle, M, N)) return(-1);
00102     }
00103   }
00104   if (BlockMapToHandle(handle, map)) return(-1); // Everybody calls this routine
00105   
00106   if (map.Comm().MyPID()==0) // Only PE 0 opened a file
00107     if (fclose(handle)) return(-1);
00108   return(0);
00109 }
00110 
00111 int BlockMapToHandle(FILE * handle, const Epetra_BlockMap & map) {
00112 
00113   const Epetra_Comm & comm = map.Comm();
00114   int numProc = comm.NumProc();
00115   bool doSizes = !map.ConstantElementSize();
00116 
00117   if (numProc==1) {
00118     int * myElements = map.MyGlobalElements();
00119     int * elementSizeList = 0;
00120     if (doSizes) elementSizeList = map.ElementSizeList();
00121     return(writeBlockMap(handle, map.NumGlobalElements(), myElements, elementSizeList, doSizes));
00122   }
00123 
00124   int numRows = map.NumMyElements();
00125   
00126   Epetra_Map allGidsMap(-1, numRows, 0,comm);
00127   
00128   Epetra_IntVector allGids(allGidsMap);
00129   for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
00130   
00131   Epetra_IntVector allSizes(allGidsMap);
00132   for (int i=0; i<numRows; i++) allSizes[i] = map.ElementSize(i);
00133   
00134   // Now construct a Map on PE 0 by strip-mining the rows of the input matrix map.
00135   int numChunks = numProc;
00136   int stripSize = allGids.GlobalLength()/numChunks;
00137   int remainder = allGids.GlobalLength()%numChunks;
00138   int curStart = 0;
00139   int curStripSize = 0;
00140   Epetra_IntSerialDenseVector importGidList;
00141   Epetra_IntSerialDenseVector importSizeList;
00142   if (comm.MyPID()==0) {
00143     importGidList.Size(stripSize+1); // Set size of vector to max needed
00144     if (doSizes) importSizeList.Size(stripSize+1); // Set size of vector to max needed
00145   }
00146   for (int i=0; i<numChunks; i++) {
00147     if (comm.MyPID()==0) { // Only PE 0 does this part
00148       curStripSize = stripSize;
00149       if (i<remainder) curStripSize++; // handle leftovers
00150       for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
00151       curStart += curStripSize;
00152     }
00153     // The following import map will be non-trivial only on PE 0.
00154     Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
00155     Epetra_Import gidImporter(importGidMap, allGidsMap);
00156     
00157     Epetra_IntVector importGids(importGidMap);
00158     if (importGids.Import(allGids, gidImporter, Insert)) return(-1); 
00159     Epetra_IntVector importSizes(importGidMap);
00160     if (doSizes) if (importSizes.Import(allSizes, gidImporter, Insert)) return(-1); 
00161     
00162     // importGids (and importSizes, if non-trivial block map)
00163     // now have a list of GIDs (and sizes, respectively) for the current strip of map.
00164     
00165     int * myElements = importGids.Values();
00166     int * elementSizeList = 0;
00167     if (doSizes) elementSizeList = importSizes.Values();
00168     // Finally we are ready to write this strip of the map to file
00169     writeBlockMap(handle, importGids.MyLength(), myElements, elementSizeList, doSizes);
00170   }
00171   return(0);
00172 }
00173 int writeBlockMap(FILE * handle, int length, const int * v1, const int * v2, bool doSizes) {
00174 
00175   for (int i=0; i<length; i++) {
00176     fprintf(handle, "%d", v1[i]);
00177     if (doSizes) fprintf(handle, " %d", v2[i]);
00178     fprintf(handle, "\n");
00179   }
00180   return(0);
00181 }
00182 } // namespace EpetraExt

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