EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
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"
53#include "Epetra_Import.h"
54
55using namespace EpetraExt;
56namespace EpetraExt {
57
58int 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
129template<typename int_type>
130int 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
194int 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
210int 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
220int 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
#define mm_set_array(typecode)
char MM_typecode[4]
#define mm_initialize_typecode(typecode)
#define mm_set_integer(typecode)
#define mm_set_matrix(typecode)
Insert
long long IndexBase64() const
int MaxElementSize() const
int * ElementSizeList() const
long long GID64(int LID) const
int ElementSize() const
long long NumGlobalElements64() const
bool GlobalIndicesInt() const
int MinElementSize() const
const Epetra_Comm & Comm() const
bool ConstantElementSize() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
virtual int NumProc() const=0
virtual int MyPID() const=0
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Size(int Length_in)
int MyLength() const
int * Values() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int mm_write_mtx_array_size(FILE *f, long long M, long long N)
int mm_write_banner(FILE *f, MM_typecode matcode)
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.
int TBlockMapToHandle(FILE *handle, const Epetra_BlockMap &map)
const int N
int BlockMapToHandle(FILE *handle, const Epetra_BlockMap &map)
int writeBlockMap(FILE *handle, long long length, const int *v1, const int *v2, bool doSizes)