EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_BlockMapOut.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_BlockMapOut.h"
00043 #include "EpetraExt_mmio.h"
00044 #include "Epetra_Comm.h"
00045 #include "Epetra_BlockMap.h"
00046 #include "Epetra_Map.h"
00047 #include "Epetra_IntVector.h"
00048 #include "Epetra_LongLongVector.h"
00049 #include "Epetra_GIDTypeVector.h"
00050 #include "Epetra_IntSerialDenseVector.h"
00051 #include "Epetra_LongLongSerialDenseVector.h"
00052 #include "Epetra_GIDTypeSerialDenseVector.h"
00053 #include "Epetra_Import.h"
00054 
00055 using namespace EpetraExt;
00056 namespace EpetraExt {
00057 
00058 int BlockMapToMatrixMarketFile( const char *filename, const Epetra_BlockMap & map, 
00059          const char * mapName,
00060          const char *mapDescription, 
00061          bool writeHeader) {
00062   long long M = map.NumGlobalElements64();
00063   int N = 1;
00064   if (map.MaxElementSize()>1) N = 2; // Non-trivial block map, store element sizes in second column
00065 
00066   FILE * handle = 0;
00067 
00068   if (map.Comm().MyPID()==0) { // Only PE 0 does this section
00069 
00070     handle = fopen(filename,"w");
00071     if (!handle) return(-1);
00072     MM_typecode matcode;
00073     mm_initialize_typecode(&matcode);
00074     mm_set_matrix(&matcode);
00075     mm_set_array(&matcode);
00076     mm_set_integer(&matcode);
00077 
00078     if (writeHeader==true) { // Only write header if requested (true by default)
00079     
00080       if (mm_write_banner(handle, matcode)) return(-1);
00081       
00082       if (mapName!=0) fprintf(handle, "%% \n%% %s\n", mapName);
00083       if (mapDescription!=0) fprintf(handle, "%% %s\n%% \n", mapDescription);
00084 
00085     }
00086   }
00087     
00088   if (writeHeader==true) { // Only write header if requested (true by default)
00089 
00090     // Make an Epetra_IntVector of length numProc such that all elements are on PE 0 and
00091     // the ith element is NumMyElements from the ith PE
00092 
00093     Epetra_Map map1(-1, 1, 0, map.Comm()); // map with one element on each processor
00094     int length = 0;
00095     if (map.Comm().MyPID()==0) length = map.Comm().NumProc();
00096     Epetra_Map map2(-1, length, 0, map.Comm());
00097     Epetra_Import lengthImporter(map2, map1);
00098     Epetra_IntVector v1(map1);
00099     Epetra_IntVector v2(map2);
00100     v1[0] = map.NumMyElements();
00101     if (v2.Import(v1, lengthImporter, Insert)) return(-1);
00102     if (map.Comm().MyPID()==0) { 
00103       fprintf(handle, "%s", "%Format Version:\n");
00104       //int version = 1; // We may change the format scheme at a later date.
00105       fprintf(handle, "%% %d \n", map.Comm().NumProc());
00106       fprintf(handle, "%s", "%NumProc: Number of processors:\n");
00107       fprintf(handle, "%% %d \n", map.Comm().NumProc());
00108       fprintf(handle, "%s", "%MaxElementSize: Maximum element size:\n");
00109       fprintf(handle, "%% %d \n", map.MaxElementSize());
00110       fprintf(handle, "%s", "%MinElementSize: Minimum element size:\n");
00111       fprintf(handle, "%% %d \n", map.MinElementSize());
00112       fprintf(handle, "%s", "%IndexBase: Index base of map:\n");
00113       fprintf(handle, "%% %lld \n", map.IndexBase64());
00114       fprintf(handle, "%s", "%NumGlobalElements: Total number of GIDs in map:\n");
00115       fprintf(handle, "%% %lld \n", map.NumGlobalElements64());
00116       fprintf(handle, "%s", "%NumMyElements: BlockMap lengths per processor:\n");
00117       for ( int i=0; i< v2.MyLength(); i++) fprintf(handle, "%% %d\n", v2[i]);
00118       
00119       if (mm_write_mtx_array_size(handle, M, N)) return(-1);
00120     }
00121   }
00122   if (BlockMapToHandle(handle, map)) return(-1); // Everybody calls this routine
00123   
00124   if (map.Comm().MyPID()==0) // Only PE 0 opened a file
00125     if (fclose(handle)) return(-1);
00126   return(0);
00127 }
00128 
00129 template<typename int_type>
00130 int TBlockMapToHandle(FILE * handle, const Epetra_BlockMap & map) {
00131 
00132   const Epetra_Comm & comm = map.Comm();
00133   int numProc = comm.NumProc();
00134   bool doSizes = !map.ConstantElementSize();
00135 
00136   if (numProc==1) {
00137     int_type * myElements = 0;
00138   map.MyGlobalElementsPtr(myElements);
00139     int * elementSizeList = 0;
00140     if (doSizes) elementSizeList = map.ElementSizeList();
00141     return(writeBlockMap(handle, map.NumGlobalElements64(), myElements, elementSizeList, doSizes));
00142   }
00143 
00144   int numRows = map.NumMyElements();
00145   
00146   Epetra_Map allGidsMap((int_type) -1, numRows, (int_type) 0,comm);
00147   
00148   typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
00149   for (int i=0; i<numRows; i++) allGids[i] = (int_type) map.GID64(i);
00150   
00151   Epetra_IntVector allSizes(allGidsMap);
00152   for (int i=0; i<numRows; i++) allSizes[i] = map.ElementSize(i);
00153   
00154   // Now construct a Map on PE 0 by strip-mining the rows of the input matrix map.
00155   int numChunks = numProc;
00156   int stripSize = allGids.GlobalLength64()/numChunks;
00157   int remainder = allGids.GlobalLength64()%numChunks;
00158   int curStart = 0;
00159   int curStripSize = 0;
00160   typename Epetra_GIDTypeSerialDenseVector<int_type>::impl importGidList;
00161   Epetra_IntSerialDenseVector importSizeList;
00162   if (comm.MyPID()==0) {
00163     importGidList.Size(stripSize+1); // Set size of vector to max needed
00164     if (doSizes) importSizeList.Size(stripSize+1); // Set size of vector to max needed
00165   }
00166   for (int i=0; i<numChunks; i++) {
00167     if (comm.MyPID()==0) { // Only PE 0 does this part
00168       curStripSize = stripSize;
00169       if (i<remainder) curStripSize++; // handle leftovers
00170       for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
00171       curStart += curStripSize;
00172     }
00173     // The following import map will be non-trivial only on PE 0.
00174     Epetra_Map importGidMap((int_type) -1, curStripSize, importGidList.Values(), 0, comm);
00175     Epetra_Import gidImporter(importGidMap, allGidsMap);
00176     
00177     typename Epetra_GIDTypeVector<int_type>::impl importGids(importGidMap);
00178     if (importGids.Import(allGids, gidImporter, Insert)) return(-1); 
00179     Epetra_IntVector importSizes(importGidMap);
00180     if (doSizes) if (importSizes.Import(allSizes, gidImporter, Insert)) return(-1); 
00181     
00182     // importGids (and importSizes, if non-trivial block map)
00183     // now have a list of GIDs (and sizes, respectively) for the current strip of map.
00184     
00185     int_type * myElements = importGids.Values();
00186     int * elementSizeList = 0;
00187     if (doSizes) elementSizeList = importSizes.Values();
00188     // Finally we are ready to write this strip of the map to file
00189     writeBlockMap(handle, importGids.MyLength(), myElements, elementSizeList, doSizes);
00190   }
00191   return(0);
00192 }
00193 
00194 int BlockMapToHandle(FILE * handle, const Epetra_BlockMap & map) {
00195 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00196   if(map.GlobalIndicesInt()) {
00197     return TBlockMapToHandle<int>(handle, map);
00198   }
00199   else
00200 #endif
00201 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00202   if(map.GlobalIndicesLongLong()) {
00203     return TBlockMapToHandle<long long>(handle, map);
00204   }
00205   else
00206 #endif
00207     throw "EpetraExt::BlockMapToHandle: GlobalIndices type unknown";
00208 }
00209 
00210 int writeBlockMap(FILE * handle, long long length, const int * v1, const int * v2, bool doSizes) {
00211 
00212   for (long long i=0; i<length; i++) {
00213     fprintf(handle, "%d", v1[i]);
00214     if (doSizes) fprintf(handle, " %d", v2[i]);
00215     fprintf(handle, "%s", "\n");
00216   }
00217   return(0);
00218 }
00219 
00220 int writeBlockMap(FILE * handle, long long length, const long long * v1, const int * v2, bool doSizes) {
00221 
00222   for (long long i=0; i<length; i++) {
00223     fprintf(handle, "%lld", v1[i]);
00224     if (doSizes) fprintf(handle, " %d", v2[i]);
00225     fprintf(handle, "%s", "\n");
00226   }
00227   return(0);
00228 }
00229 
00230 } // namespace EpetraExt
00231 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines