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