EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_BlockMapOut.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 #include "Epetra_ConfigDefs.h"
42 #include "EpetraExt_BlockMapOut.h"
43 #include "EpetraExt_mmio.h"
44 #include "Epetra_Comm.h"
45 #include "Epetra_BlockMap.h"
46 #include "Epetra_Map.h"
47 #include "Epetra_IntVector.h"
48 #include "Epetra_LongLongVector.h"
49 #include "Epetra_GIDTypeVector.h"
50 #include "Epetra_IntSerialDenseVector.h"
51 #include "Epetra_LongLongSerialDenseVector.h"
52 #include "Epetra_GIDTypeSerialDenseVector.h"
53 #include "Epetra_Import.h"
54 
55 using namespace EpetraExt;
56 namespace EpetraExt {
57 
58 int BlockMapToMatrixMarketFile( const char *filename, const Epetra_BlockMap & map,
59  const char * mapName,
60  const char *mapDescription,
61  bool writeHeader) {
62  long long M = map.NumGlobalElements64();
63  int N = 1;
64  if (map.MaxElementSize()>1) N = 2; // Non-trivial block map, store element sizes in second column
65 
66  FILE * handle = 0;
67 
68  if (map.Comm().MyPID()==0) { // Only PE 0 does this section
69 
70  handle = fopen(filename,"w");
71  if (!handle) return(-1);
72  MM_typecode matcode;
73  mm_initialize_typecode(&matcode);
74  mm_set_matrix(&matcode);
75  mm_set_array(&matcode);
76  mm_set_integer(&matcode);
77 
78  if (writeHeader==true) { // Only write header if requested (true by default)
79 
80  if (mm_write_banner(handle, matcode)) return(-1);
81 
82  if (mapName!=0) fprintf(handle, "%% \n%% %s\n", mapName);
83  if (mapDescription!=0) fprintf(handle, "%% %s\n%% \n", mapDescription);
84 
85  }
86  }
87 
88  if (writeHeader==true) { // Only write header if requested (true by default)
89 
90  // Make an Epetra_IntVector of length numProc such that all elements are on PE 0 and
91  // the ith element is NumMyElements from the ith PE
92 
93  Epetra_Map map1(-1, 1, 0, map.Comm()); // map with one element on each processor
94  int length = 0;
95  if (map.Comm().MyPID()==0) length = map.Comm().NumProc();
96  Epetra_Map map2(-1, length, 0, map.Comm());
97  Epetra_Import lengthImporter(map2, map1);
98  Epetra_IntVector v1(map1);
99  Epetra_IntVector v2(map2);
100  v1[0] = map.NumMyElements();
101  if (v2.Import(v1, lengthImporter, Insert)) return(-1);
102  if (map.Comm().MyPID()==0) {
103  fprintf(handle, "%s", "%Format Version:\n");
104  //int version = 1; // We may change the format scheme at a later date.
105  fprintf(handle, "%% %d \n", map.Comm().NumProc());
106  fprintf(handle, "%s", "%NumProc: Number of processors:\n");
107  fprintf(handle, "%% %d \n", map.Comm().NumProc());
108  fprintf(handle, "%s", "%MaxElementSize: Maximum element size:\n");
109  fprintf(handle, "%% %d \n", map.MaxElementSize());
110  fprintf(handle, "%s", "%MinElementSize: Minimum element size:\n");
111  fprintf(handle, "%% %d \n", map.MinElementSize());
112  fprintf(handle, "%s", "%IndexBase: Index base of map:\n");
113  fprintf(handle, "%% %lld \n", map.IndexBase64());
114  fprintf(handle, "%s", "%NumGlobalElements: Total number of GIDs in map:\n");
115  fprintf(handle, "%% %lld \n", map.NumGlobalElements64());
116  fprintf(handle, "%s", "%NumMyElements: BlockMap lengths per processor:\n");
117  for ( int i=0; i< v2.MyLength(); i++) fprintf(handle, "%% %d\n", v2[i]);
118 
119  if (mm_write_mtx_array_size(handle, M, N)) return(-1);
120  }
121  }
122  if (BlockMapToHandle(handle, map)) return(-1); // Everybody calls this routine
123 
124  if (map.Comm().MyPID()==0) // Only PE 0 opened a file
125  if (fclose(handle)) return(-1);
126  return(0);
127 }
128 
129 template<typename int_type>
130 int TBlockMapToHandle(FILE * handle, const Epetra_BlockMap & map) {
131 
132  const Epetra_Comm & comm = map.Comm();
133  int numProc = comm.NumProc();
134  bool doSizes = !map.ConstantElementSize();
135 
136  if (numProc==1) {
137  int_type * myElements = 0;
138  map.MyGlobalElementsPtr(myElements);
139  int * elementSizeList = 0;
140  if (doSizes) elementSizeList = map.ElementSizeList();
141  return(writeBlockMap(handle, map.NumGlobalElements64(), myElements, elementSizeList, doSizes));
142  }
143 
144  int numRows = map.NumMyElements();
145 
146  Epetra_Map allGidsMap((int_type) -1, numRows, (int_type) 0,comm);
147 
148  typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
149  for (int i=0; i<numRows; i++) allGids[i] = (int_type) map.GID64(i);
150 
151  Epetra_IntVector allSizes(allGidsMap);
152  for (int i=0; i<numRows; i++) allSizes[i] = map.ElementSize(i);
153 
154  // Now construct a Map on PE 0 by strip-mining the rows of the input matrix map.
155  int numChunks = numProc;
156  int stripSize = allGids.GlobalLength64()/numChunks;
157  int remainder = allGids.GlobalLength64()%numChunks;
158  int curStart = 0;
159  int curStripSize = 0;
161  Epetra_IntSerialDenseVector importSizeList;
162  if (comm.MyPID()==0) {
163  importGidList.Size(stripSize+1); // Set size of vector to max needed
164  if (doSizes) importSizeList.Size(stripSize+1); // Set size of vector to max needed
165  }
166  for (int i=0; i<numChunks; i++) {
167  if (comm.MyPID()==0) { // Only PE 0 does this part
168  curStripSize = stripSize;
169  if (i<remainder) curStripSize++; // handle leftovers
170  for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
171  curStart += curStripSize;
172  }
173  // The following import map will be non-trivial only on PE 0.
174  Epetra_Map importGidMap((int_type) -1, curStripSize, importGidList.Values(), 0, comm);
175  Epetra_Import gidImporter(importGidMap, allGidsMap);
176 
177  typename Epetra_GIDTypeVector<int_type>::impl importGids(importGidMap);
178  if (importGids.Import(allGids, gidImporter, Insert)) return(-1);
179  Epetra_IntVector importSizes(importGidMap);
180  if (doSizes) if (importSizes.Import(allSizes, gidImporter, Insert)) return(-1);
181 
182  // importGids (and importSizes, if non-trivial block map)
183  // now have a list of GIDs (and sizes, respectively) for the current strip of map.
184 
185  int_type * myElements = importGids.Values();
186  int * elementSizeList = 0;
187  if (doSizes) elementSizeList = importSizes.Values();
188  // Finally we are ready to write this strip of the map to file
189  writeBlockMap(handle, importGids.MyLength(), myElements, elementSizeList, doSizes);
190  }
191  return(0);
192 }
193 
194 int BlockMapToHandle(FILE * handle, const Epetra_BlockMap & map) {
195 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
196  if(map.GlobalIndicesInt()) {
197  return TBlockMapToHandle<int>(handle, map);
198  }
199  else
200 #endif
201 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
202  if(map.GlobalIndicesLongLong()) {
203  return TBlockMapToHandle<long long>(handle, map);
204  }
205  else
206 #endif
207  throw "EpetraExt::BlockMapToHandle: GlobalIndices type unknown";
208 }
209 
210 int writeBlockMap(FILE * handle, long long length, const int * v1, const int * v2, bool doSizes) {
211 
212  for (long long i=0; i<length; i++) {
213  fprintf(handle, "%d", v1[i]);
214  if (doSizes) fprintf(handle, " %d", v2[i]);
215  fprintf(handle, "%s", "\n");
216  }
217  return(0);
218 }
219 
220 int writeBlockMap(FILE * handle, long long length, const long long * v1, const int * v2, bool doSizes) {
221 
222  for (long long i=0; i<length; i++) {
223  fprintf(handle, "%lld", v1[i]);
224  if (doSizes) fprintf(handle, " %d", v2[i]);
225  fprintf(handle, "%s", "\n");
226  }
227  return(0);
228 }
229 
230 } // namespace EpetraExt
231 
int writeBlockMap(FILE *handle, long long length, const int *v1, const int *v2, bool doSizes)
int ElementSize() const
bool GlobalIndicesLongLong() const
int BlockMapToMatrixMarketFile(const char *filename, const Epetra_BlockMap &map, const char *mapName, const char *mapDescription, bool writeHeader)
Writes an Epetra_BlockMap or Epetra_Map object to a Matrix Market format file.
#define mm_set_matrix(typecode)
int Size(int Length_in)
bool ConstantElementSize() const
int mm_write_banner(FILE *f, MM_typecode matcode)
#define mm_set_integer(typecode)
bool GlobalIndicesInt() const
int * Values() const
virtual int MyPID() const =0
int NumMyElements() const
int * ElementSizeList() const
char MM_typecode[4]
int BlockMapToHandle(FILE *handle, const Epetra_BlockMap &map)
const Epetra_Comm & Comm() const
virtual int NumProc() const =0
int TBlockMapToHandle(FILE *handle, const Epetra_BlockMap &map)
int MaxElementSize() const
#define mm_set_array(typecode)
int mm_write_mtx_array_size(FILE *f, long long M, long long N)
#define mm_initialize_typecode(typecode)
int MyLength() const
int MinElementSize() const