EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_BlockMapIn.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 "Epetra_Comm.h"
44#include "Epetra_Util.h"
45#include "Epetra_BlockMap.h"
46#include "Epetra_Map.h"
47#include "Epetra_IntVector.h"
48#include "Epetra_IntSerialDenseVector.h"
49#include "Epetra_Import.h"
50#include "EpetraExt_mmio.h"
51
52using namespace EpetraExt;
53namespace EpetraExt {
54
55#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
56int MatrixMarketFileToMap( const char *filename, const Epetra_Comm & comm, Epetra_Map * & map) {
57
58 Epetra_BlockMap * bmap;
59 if (MatrixMarketFileToBlockMap(filename, comm, bmap)) return(-1);
60 map = dynamic_cast<Epetra_Map *>(bmap);
61 return(0);
62}
63#endif
64
65#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
66int MatrixMarketFileToMap64( const char *filename, const Epetra_Comm & comm, Epetra_Map * & map) {
67
68 Epetra_BlockMap * bmap;
69 if (MatrixMarketFileToBlockMap64(filename, comm, bmap)) return(-1);
70 map = dynamic_cast<Epetra_Map *>(bmap);
71 return(0);
72}
73#endif
74
75#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
76int MatrixMarketFileToBlockMap( const char *filename, const Epetra_Comm & comm, Epetra_BlockMap * & map) {
77
78 const int lineLength = 1025;
79 char line[lineLength];
80 char token[lineLength];
81 int M, N, numProc, MaxElementSize, MinElementSize, NumMyElements, IndexBase, NumGlobalElements, firstGid;
82
83 FILE * handle = 0;
84
85 bool inHeader = true;
86
87 handle = fopen(filename,"r");
88 if (handle == 0)
89 EPETRA_CHK_ERR(-1); // file not found
90
91 while (inHeader) {
92 if(fgets(line, lineLength, handle)==0) return(-1);
93 if(sscanf(line, "%s", token)==0) return(-1);
94 if (!strcmp(token, "%NumProc:")) inHeader = false;
95 }
96
97 if(fgets(line, lineLength, handle)==0) return(-1); // numProc value
98 if(sscanf(line, "%s %d", token, &numProc)==0) return(-1);
99
100 if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize header line
101 if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize value
102 if(sscanf(line, "%s %d", token, &MaxElementSize)==0) return(-1);
103
104 if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize header line
105 if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize value
106 if(sscanf(line, "%s %d", token, &MinElementSize)==0) return(-1);
107
108 if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase header line
109 if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase value
110 if(sscanf(line, "%s %d", token, &IndexBase)==0) return(-1);
111
112 if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements header line
113 if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements value
114 if(sscanf(line, "%s %d", token, &NumGlobalElements)==0) return(-1);
115
116 int ierr = 0;
117 (void) ierr; // mfh 13 Jan 2013: Forestall compiler warning for unused var.
118 if (comm.NumProc()==numProc) {
119 if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
120 firstGid = 0;
121 for (int i=0; i<comm.MyPID(); i++) {
122 if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value
123 if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
124 firstGid += NumMyElements;
125 }
126
127 if(fgets(line, lineLength, handle)==0) return(-1); // This PE's NumMyElements value
128 if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
129
130 for (int i=comm.MyPID()+1; i<numProc; i++) {
131 if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
132 }
133 }
134 else {
135 ierr = 1; // Warning error, different number of processors.
136
137 if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
138 for (int i=0; i<numProc; i++) {
139 if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
140 }
141
142 NumMyElements = NumGlobalElements/comm.NumProc();
143 firstGid = comm.MyPID()*NumMyElements;
144 int remainder = NumGlobalElements%comm.NumProc();
145 if (comm.MyPID()<remainder) NumMyElements++;
146 int extra = remainder;
147 if (comm.MyPID()<remainder) extra = comm.MyPID();
148 firstGid += extra;
149 }
150 if(fgets(line, lineLength, handle)==0) return(-1); // Number of rows, columns
151 if(sscanf(line, "%d %d", &M, &N)==0) return(-1);
152
153 bool doSizes = (N>1);
154 Epetra_IntSerialDenseVector v1(NumMyElements);
155 Epetra_IntSerialDenseVector v2(NumMyElements);
156 for (int i=0; i<firstGid; i++) {
157 if(fgets(line, lineLength, handle)==0) return(-1); // dump these
158 }
159
160 if (doSizes) {
161 for (int i=0; i<NumMyElements; i++) {
162 if(fgets(line, lineLength, handle)==0) return(-1);
163 if(sscanf(line, "%d %d", &v1[i], &v2[i])==0) return(-1); // load v1, v2
164 }
165 }
166 else {
167 for (int i=0; i<NumMyElements; i++) {
168 if(fgets(line, lineLength, handle)==0) return(-1);
169 if(sscanf(line, "%d", &v1[i])==0) return(-1); // load v1
170 v2[i] = MinElementSize; // Fill with constant size
171 }
172 }
173 if (fclose(handle)) return(-1);
174
175 comm.Barrier();
176
177 if (MinElementSize==1 && MaxElementSize==1)
178 map = new Epetra_Map(-1, NumMyElements, v1.Values(), IndexBase, comm);
179 else
180 map = new Epetra_BlockMap(-1, NumMyElements, v1.Values(), v2.Values(), IndexBase, comm);
181 return(0);
182}
183#endif
184
185#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
186int MatrixMarketFileToBlockMap64( const char *filename, const Epetra_Comm & comm, Epetra_BlockMap * & map) {
187
188 const int lineLength = 1025;
189 char line[lineLength];
190 char token[lineLength];
191 int numProc, MaxElementSize, MinElementSize, NumMyElements;
192 long long M, N, IndexBase, NumGlobalElements, firstGid;
193
194 FILE * handle = 0;
195
196 bool inHeader = true;
197
198 handle = fopen(filename,"r");
199 if (handle == 0)
200 EPETRA_CHK_ERR(-1); // file not found
201
202 while (inHeader) {
203 if(fgets(line, lineLength, handle)==0) return(-1);
204 if(sscanf(line, "%s", token)==0) return(-1);
205 if (!strcmp(token, "%NumProc:")) inHeader = false;
206 }
207
208 if(fgets(line, lineLength, handle)==0) return(-1); // numProc value
209 if(sscanf(line, "%s %d", token, &numProc)==0) return(-1);
210
211 if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize header line
212 if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize value
213 if(sscanf(line, "%s %d", token, &MaxElementSize)==0) return(-1);
214
215 if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize header line
216 if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize value
217 if(sscanf(line, "%s %d", token, &MinElementSize)==0) return(-1);
218
219 if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase header line
220 if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase value
221 if(sscanf(line, "%s %lld", token, &IndexBase)==0) return(-1);
222
223 if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements header line
224 if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements value
225 if(sscanf(line, "%s %lld", token, &NumGlobalElements)==0) return(-1);
226
227 int ierr = 0;
228 (void) ierr; // mfh 13 Jan 2013: Forestall compiler warning for unused var.
229 if (comm.NumProc()==numProc) {
230 if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
231 firstGid = 0;
232 for (int i=0; i<comm.MyPID(); i++) {
233 if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value
234 if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
235 firstGid += NumMyElements;
236 }
237
238 if(fgets(line, lineLength, handle)==0) return(-1); // This PE's NumMyElements value
239 if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
240
241 for (int i=comm.MyPID()+1; i<numProc; i++) {
242 if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
243 }
244 }
245 else {
246 ierr = 1; // Warning error, different number of processors.
247
248 if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
249 for (int i=0; i<numProc; i++) {
250 if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
251 }
252
253 NumMyElements = (int) (NumGlobalElements/comm.NumProc());
254 firstGid = ((long long) comm.MyPID())*NumMyElements;
255 int remainder = (int) (NumGlobalElements%comm.NumProc());
256 if (comm.MyPID()<remainder) NumMyElements++;
257 int extra = remainder;
258 if (comm.MyPID()<remainder) extra = comm.MyPID();
259 firstGid += extra;
260 }
261 if(fgets(line, lineLength, handle)==0) return(-1); // Number of rows, columns
262 if(sscanf(line, "%lld %lld", &M, &N)==0) return(-1);
263
264 bool doSizes = (N>1);
265 Epetra_LongLongSerialDenseVector v1(NumMyElements);
266 Epetra_IntSerialDenseVector v2(NumMyElements);
267 for (long long i=0; i<firstGid; i++) {
268 if(fgets(line, lineLength, handle)==0) return(-1); // dump these
269 }
270
271 if (doSizes) {
272 for (int i=0; i<NumMyElements; i++) {
273 if(fgets(line, lineLength, handle)==0) return(-1);
274 if(sscanf(line, "%lld %d", &v1[i], &v2[i])==0) return(-1); // load v1, v2
275 }
276 }
277 else {
278 for (int i=0; i<NumMyElements; i++) {
279 if(fgets(line, lineLength, handle)==0) return(-1);
280 if(sscanf(line, "%lld", &v1[i])==0) return(-1); // load v1
281 v2[i] = MinElementSize; // Fill with constant size
282 }
283 }
284 if (fclose(handle)) return(-1);
285
286 comm.Barrier();
287
288 if (MinElementSize==1 && MaxElementSize==1)
289 map = new Epetra_Map(-1LL, NumMyElements, v1.Values(), IndexBase, comm);
290 else
291 map = new Epetra_BlockMap(-1LL, NumMyElements, v1.Values(), v2.Values(), IndexBase, comm);
292 return(0);
293}
294#endif
295
296#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
297int MatrixMarketFileToRowMap(const char* filename,
298 const Epetra_Comm& comm,
299 Epetra_BlockMap*& rowmap)
300{
301 FILE* infile = fopen(filename, "r");
302 MM_typecode matcode;
303
304 int err = mm_read_banner(infile, &matcode);
305 if (err != 0) return(err);
306
307 if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
308 !mm_is_real(matcode) || !mm_is_general(matcode)) {
309 return(-1);
310 }
311
312 int numrows, numcols;
313 err = mm_read_mtx_array_size(infile, &numrows, &numcols);
314 if (err != 0) return(err);
315
316 fclose(infile);
317
318 rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
319 return(0);
320}
321#endif
322
323#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
324int MatrixMarketFileToBlockMaps(const char* filename,
325 const Epetra_Comm& comm,
326 Epetra_BlockMap*& rowmap,
327 Epetra_BlockMap*& colmap,
328 Epetra_BlockMap*& rangemap,
329 Epetra_BlockMap*& domainmap)
330{
331 FILE* infile = fopen(filename, "r");
332 if (infile == NULL) {
333 return(-1);
334 }
335
336 MM_typecode matcode;
337
338 int err = mm_read_banner(infile, &matcode);
339 if (err != 0) return(err);
340
341 if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
342 !mm_is_real(matcode) || !mm_is_general(matcode)) {
343 return(-1);
344 }
345
346 int numrows, numcols, nnz;
347 err = mm_read_mtx_crd_size(infile, &numrows, &numcols, &nnz);
348 if (err != 0) return(err);
349
350 //for this case, we'll assume that the row-map is the same as
351 //the range-map.
352 //create row-map and range-map with linear distributions.
353
354 rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
355 rangemap = new Epetra_BlockMap(numrows, 1, 0, comm);
356
357 int I, J;
358 double val, imag;
359
360 int num_map_cols = 0, insertPoint, foundOffset;
361 int allocLen = numcols;
362 int* map_cols = new int[allocLen];
363
364 //read through all matrix data and construct a list of the column-
365 //indices that occur in rows that are local to this processor.
366
367 for(int i=0; i<nnz; ++i) {
368 err = mm_read_mtx_crd_entry(infile, &I, &J, &val,
369 &imag, matcode);
370
371 if (err == 0) {
372 --I;
373 --J;
374 if (rowmap->MyGID(I)) {
375 foundOffset = Epetra_Util_binary_search(J, map_cols, num_map_cols,
376 insertPoint);
377 if (foundOffset < 0) {
378 Epetra_Util_insert(J, insertPoint, map_cols,
379 num_map_cols, allocLen);
380 }
381 }
382 }
383 }
384
385 //create colmap with the list of columns associated with rows that are
386 //local to this processor.
387 colmap = new Epetra_Map(-1, num_map_cols, map_cols, 0, comm);
388
389 //create domainmap which has a linear distribution
390 domainmap = new Epetra_BlockMap(numcols, 1, 0, comm);
391
392 delete [] map_cols;
393
394 return(0);
395}
396#endif
397
398#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
399int MatrixMarketFileToBlockMaps64(const char* filename,
400 const Epetra_Comm& comm,
401 Epetra_BlockMap*& rowmap,
402 Epetra_BlockMap*& colmap,
403 Epetra_BlockMap*& rangemap,
404 Epetra_BlockMap*& domainmap)
405{
406 FILE* infile = fopen(filename, "r");
407 if (infile == NULL) {
408 return(-1);
409 }
410
411 MM_typecode matcode;
412
413 int err = mm_read_banner(infile, &matcode);
414 if (err != 0) return(err);
415
416 if (!mm_is_matrix(matcode) || !mm_is_coordinate(matcode) ||
417 !mm_is_real(matcode) || !mm_is_general(matcode)) {
418 return(-1);
419 }
420
421 long long numrows, numcols, nnz;
422 err = mm_read_mtx_crd_size(infile, &numrows, &numcols, &nnz);
423 if (err != 0) return(err);
424
425 //for this case, we'll assume that the row-map is the same as
426 //the range-map.
427 //create row-map and range-map with linear distributions.
428
429 rowmap = new Epetra_BlockMap(numrows, 1, 0, comm);
430 rangemap = new Epetra_BlockMap(numrows, 1, 0, comm);
431
432 long long I, J;
433 double val, imag;
434
435 int num_map_cols = 0, insertPoint, foundOffset;
436 int allocLen = 0;
437 if(numcols > (long long) std::numeric_limits<int>::max())
438 allocLen = std::numeric_limits<int>::max();
439 else
440 allocLen = (int) numcols;
441
442 long long* map_cols = new long long[allocLen];
443
444 //read through all matrix data and construct a list of the column-
445 //indices that occur in rows that are local to this processor.
446
447 for(int i=0; i<nnz; ++i) {
448 err = mm_read_mtx_crd_entry(infile, &I, &J, &val,
449 &imag, matcode);
450
451 if (err == 0) {
452 --I;
453 --J;
454 if (rowmap->MyGID(I)) {
455 foundOffset = Epetra_Util_binary_search(J, map_cols, num_map_cols,
456 insertPoint);
457 if (foundOffset < 0) {
458 Epetra_Util_insert(J, insertPoint, map_cols,
459 num_map_cols, allocLen);
460 }
461 }
462 }
463 }
464
465 //create colmap with the list of columns associated with rows that are
466 //local to this processor.
467 colmap = new Epetra_Map(-1LL, num_map_cols, map_cols, 0, comm);
468
469 //create domainmap which has a linear distribution
470 domainmap = new Epetra_BlockMap(numcols, 1, 0, comm);
471
472 delete [] map_cols;
473
474 return(0);
475}
476#endif
477
478} // namespace EpetraExt
479
char MM_typecode[4]
#define mm_is_matrix(typecode)
#define mm_is_real(typecode)
#define mm_is_coordinate(typecode)
#define mm_is_general(typecode)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int MatrixMarketFileToBlockMaps(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&rowmap, Epetra_BlockMap *&colmap, Epetra_BlockMap *&rangemap, Epetra_BlockMap *&domainmap)
Constructs row,col,range and domain maps from a matrix-market matrix file.
int MatrixMarketFileToMap64(const char *filename, const Epetra_Comm &comm, Epetra_Map *&map)
int MatrixMarketFileToBlockMaps64(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&rowmap, Epetra_BlockMap *&colmap, Epetra_BlockMap *&rangemap, Epetra_BlockMap *&domainmap)
Constructs row,col,range and domain maps from a matrix-market matrix file.
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
int MatrixMarketFileToRowMap(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&rowmap)
int mm_read_banner(FILE *f, MM_typecode *matcode)
int MatrixMarketFileToMap(const char *filename, const Epetra_Comm &comm, Epetra_Map *&map)
Constructs an Epetra_BlockMap object from a Matrix Market format file.
int MatrixMarketFileToBlockMap(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&map)
Constructs an Epetra_BlockMap object from a Matrix Market format file.
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *imag, MM_typecode matcode)
int MatrixMarketFileToBlockMap64(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&map)