EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
42 #include "EpetraExt_BlockMapIn.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"
49 #include "Epetra_Import.h"
50 #include "EpetraExt_mmio.h"
51 
52 using namespace EpetraExt;
53 namespace EpetraExt {
54 
55 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
56 int 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
66 int 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
76 int 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
186 int 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
297 int 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
324 int 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
399 int 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 
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
#define EPETRA_CHK_ERR(a)
virtual void Barrier() const =0
int MatrixMarketFileToMap64(const char *filename, const Epetra_Comm &comm, Epetra_Map *&map)
virtual int MyPID() const =0
int mm_read_banner(FILE *f, MM_typecode *matcode)
const int N
#define mm_is_coordinate(typecode)
int MatrixMarketFileToMap(const char *filename, const Epetra_Comm &comm, Epetra_Map *&map)
Constructs an Epetra_BlockMap object from a Matrix Market format file.
int MatrixMarketFileToRowMap(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&rowmap)
char MM_typecode[4]
bool MyGID(int GID_in) const
#define mm_is_real(typecode)
int Epetra_Util_binary_search(T item, const T *list, int len, int &insertPoint)
virtual int NumProc() const =0
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *imag, MM_typecode matcode)
#define mm_is_matrix(typecode)
int MatrixMarketFileToBlockMap64(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&map)
#define mm_is_general(typecode)
int MatrixMarketFileToBlockMap(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&map)
Constructs an Epetra_BlockMap object from a Matrix Market format file.
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 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 Epetra_Util_insert(T item, int offset, T *&list, int &usedLength, int &allocatedLength, int allocChunkSize=32)