EpetraExt  Development
EpetraExt_CrsMatrixIn.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_CrsMatrixIn.h"
43 #include "Epetra_Comm.h"
44 #include "Epetra_CrsMatrix.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_IntVector.h"
47 #include "Epetra_IntSerialDenseVector.h"
48 #include "Epetra_Import.h"
49 #include "Epetra_Time.h"
50 #include "Epetra_Util.h"
51 
52 #if defined(__PGI)
53 #include <sstream>
54 #endif
55 
58 // Functions to read a MatrixMarket file and load it into a matrix.
59 // Adapted from
60 // Trilinos/packages/epetraext/src/inout/EpetraExt_CrsMatrixIn.cpp
61 // Modified by Jon Berry and Karen Devine to make matrix reallocations
62 // more efficient; rather than insert each non-zero separately, we
63 // collect rows of non-zeros for insertion.
64 // Modified by Karen Devine and Steve Plimpton to prevent all processors
65 // from reading the same file at the same time (ouch!); chunks of the file
66 // are read and broadcast by processor zero; each processor then extracts
67 // its nonzeros from the chunk, sorts them by row, and inserts them into
68 // the matrix.
69 // The variable "chunk_read" can be changed to modify the size of the chunks
70 // read from the file. Larger values of chunk_read lead to faster execution
71 // and greater memory use.
74 
75 using namespace EpetraExt;
76 namespace EpetraExt {
77 
78 template<typename int_type>
79 static void sort_three(int_type* list, int_type *parlista, double *parlistb,
80  int start, int end);
81 
83 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
84 int MatrixMarketFileToCrsMatrix(const char *filename,
85  const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
86 {
87  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A));
88  return(0);
89 }
90 
92 int MatrixMarketFileToCrsMatrix(const char *filename,
93  const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
94 {
95  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
96  0, 0, 0, 0, transpose));
97  return(0);
98 }
100 
101 int MatrixMarketFileToCrsMatrix(const char *filename,
102  const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
103  const bool verbose)
104 {
105  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
106  0, 0, 0, 0,
107  transpose, verbose));
108  return(0);
109 }
110 
112 int MatrixMarketFileToCrsMatrix(const char *filename,
113  const Epetra_Map & rowMap, const Epetra_Map& rangeMap,
114  const Epetra_Map& domainMap, Epetra_CrsMatrix * & A, const bool transpose,
115  const bool verbose)
116 {
117  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
118  rowMap.Comm(), A,
119  &rowMap, 0,
120  &rangeMap, &domainMap,
121  transpose, verbose));
122  return(0);
123 }
124 
126 int MatrixMarketFileToCrsMatrix(const char *filename,
127  const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
128  const bool verbose)
129 {
130  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
131  rowMap.Comm(), A,
132  &rowMap, 0, 0, 0,
133  transpose, verbose));
134  return(0);
135 }
136 
138 int MatrixMarketFileToCrsMatrix(const char *filename,
139  const Epetra_Map & rowMap, const Epetra_Map & colMap,
140  Epetra_CrsMatrix * & A, const bool transpose,
141  const bool verbose)
142 {
143  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
144  rowMap.Comm(), A,
145  &rowMap, &colMap, 0, 0,
146  transpose, verbose));
147  return(0);
148 }
149 
151 int MatrixMarketFileToCrsMatrix(const char *filename,
152  const Epetra_Map & rowMap, const Epetra_Map & colMap,
153  const Epetra_Map& rangeMap, const Epetra_Map& domainMap,
154  Epetra_CrsMatrix * & A, const bool transpose,
155  const bool verbose)
156 {
157  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
158  rowMap.Comm(), A,
159  &rowMap, &colMap,
160  &rangeMap, &domainMap,
161  transpose, verbose));
162  return(0);
163 }
164 #endif
165 
166 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
167 int MatrixMarketFileToCrsMatrix64(const char *filename,
168  const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
169 {
170  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A));
171  return(0);
172 }
173 
175 int MatrixMarketFileToCrsMatrix64(const char *filename,
176  const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
177 {
178  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A,
179  0, 0, 0, 0, transpose));
180  return(0);
181 }
183 
184 int MatrixMarketFileToCrsMatrix64(const char *filename,
185  const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
186  const bool verbose)
187 {
188  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A,
189  0, 0, 0, 0,
190  transpose, verbose));
191  return(0);
192 }
193 
195 int MatrixMarketFileToCrsMatrix64(const char *filename,
196  const Epetra_Map & rowMap, const Epetra_Map& rangeMap,
197  const Epetra_Map& domainMap, Epetra_CrsMatrix * & A, const bool transpose,
198  const bool verbose)
199 {
200  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
201  rowMap.Comm(), A,
202  &rowMap, 0,
203  &rangeMap, &domainMap,
204  transpose, verbose));
205  return(0);
206 }
207 
209 int MatrixMarketFileToCrsMatrix64(const char *filename,
210  const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
211  const bool verbose)
212 {
213  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
214  rowMap.Comm(), A,
215  &rowMap, 0, 0, 0,
216  transpose, verbose));
217  return(0);
218 }
219 
221 int MatrixMarketFileToCrsMatrix64(const char *filename,
222  const Epetra_Map & rowMap, const Epetra_Map & colMap,
223  Epetra_CrsMatrix * & A, const bool transpose,
224  const bool verbose)
225 {
226  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
227  rowMap.Comm(), A,
228  &rowMap, &colMap, 0, 0,
229  transpose, verbose));
230  return(0);
231 }
232 
234 int MatrixMarketFileToCrsMatrix64(const char *filename,
235  const Epetra_Map & rowMap, const Epetra_Map & colMap,
236  const Epetra_Map& rangeMap, const Epetra_Map& domainMap,
237  Epetra_CrsMatrix * & A, const bool transpose,
238  const bool verbose)
239 {
240  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
241  rowMap.Comm(), A,
242  &rowMap, &colMap,
243  &rangeMap, &domainMap,
244  transpose, verbose));
245  return(0);
246 }
247 #endif
248 
250 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
251 int MatrixMarketFileToCrsMatrixHandle(const char *filename,
252  const Epetra_Comm & comm,
253  Epetra_CrsMatrix * & A,
254  const Epetra_Map * rowMap,
255  const Epetra_Map * colMap,
256  const Epetra_Map * rangeMap,
257  const Epetra_Map * domainMap,
258  const bool transpose,
259  const bool verbose
260 )
261 {
262  const int chunk_read = 500000; // Modify this variable to change the
263  // size of the chunks read from the file.
264  const int headerlineLength = 257;
265  const int lineLength = 81;
266  const int tokenLength = 35;
267  char line[headerlineLength];
268  char token1[tokenLength];
269  char token2[tokenLength];
270  char token3[tokenLength];
271  char token4[tokenLength];
272  char token5[tokenLength];
273  int M, N, NZ; // Matrix dimensions
274  int me = comm.MyPID();
275 
276  Epetra_Time timer(comm);
277 
278  // Make sure domain and range maps are either both defined or undefined
279  if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
280  EPETRA_CHK_ERR(-3);
281  }
282 
283  // check maps to see if domain and range are 1-to-1
284 
285  if (domainMap!=0) {
286  if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
287  if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
288  }
289  else {
290  // If domain and range are not specified,
291  // rowMap becomes both and must be 1-to-1 if specified
292  if (rowMap!=0) {
293  if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
294  }
295  }
296 
297  FILE * handle = 0;
298  if (me == 0) {
299  if (verbose) std::cout << "Reading MatrixMarket file " << filename << std::endl;
300  handle = fopen(filename,"r"); // Open file
301  if (handle == 0)
302  EPETRA_CHK_ERR(-1); // file not found
303 
304  // Check first line, which should be
305  // %%MatrixMarket matrix coordinate real general
306  if(fgets(line, headerlineLength, handle)==0) {
307  if (handle!=0) fclose(handle);
308  EPETRA_CHK_ERR(-1);
309  }
310  if(sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
311  if (handle!=0) fclose(handle);
312  EPETRA_CHK_ERR(-1);
313  }
314  if (strcmp(token1, "%%MatrixMarket") ||
315  strcmp(token2, "matrix") ||
316  strcmp(token3, "coordinate") ||
317  strcmp(token4, "real") ||
318  strcmp(token5, "general")) {
319  if (handle!=0) fclose(handle);
320  EPETRA_CHK_ERR(-1);
321  }
322 
323  // Next, strip off header lines (which start with "%")
324  do {
325  if(fgets(line, headerlineLength, handle)==0) {
326  if (handle!=0) fclose(handle);
327  EPETRA_CHK_ERR(-1);
328  }
329  } while (line[0] == '%');
330 
331  // Next get problem dimensions: M, N, NZ
332  if(sscanf(line, "%d %d %d", &M, &N, &NZ)==0) {
333  if (handle!=0) fclose(handle);
334  EPETRA_CHK_ERR(-1);
335  }
336  }
337  comm.Broadcast(&M, 1, 0);
338  comm.Broadcast(&N, 1, 0);
339  comm.Broadcast(&NZ, 1, 0);
340 
341  // Now create matrix using user maps if provided.
342 
343 
344  // Now read in chunks of triplets and broadcast them to other processors.
345  char *buffer = new char[chunk_read*lineLength];
346  int nchunk;
347  int nmillion = 0;
348  int nread = 0;
349  int rlen;
350 
351  // Storage for this processor's nonzeros.
352  const int localblock = 100000;
353  int localsize = NZ / comm.NumProc() + localblock;
354  int *iv = (int *) malloc(localsize * sizeof(int));
355  int *jv = (int *) malloc(localsize * sizeof(int));
356  double *vv = (double *) malloc(localsize * sizeof(double));
357  int lnz = 0; // Number of non-zeros on this processor.
358 
359  if (!iv || !jv || !vv)
360  EPETRA_CHK_ERR(-1);
361 
362  Epetra_Map *rowMap1;
363  bool allocatedHere=false;
364  if (rowMap != 0)
365  rowMap1 = const_cast<Epetra_Map *>(rowMap);
366  else {
367  rowMap1 = new Epetra_Map(M, 0, comm);
368  allocatedHere = true;
369  }
370  int ioffset = rowMap1->IndexBase()-1;
371  int joffset = (colMap != 0 ? colMap->IndexBase()-1 : ioffset);
372 
373  int rowmajor = 1; // Assume non-zeros are listed in row-major order,
374  int prevrow = -1; // but test to detect otherwise. If non-zeros
375  // are row-major, we can avoid the sort.
376 
377  while (nread < NZ) {
378  if (NZ-nread > chunk_read) nchunk = chunk_read;
379  else nchunk = NZ - nread;
380 
381  if (me == 0) {
382  char *eof;
383  rlen = 0;
384  for (int i = 0; i < nchunk; i++) {
385  eof = fgets(&buffer[rlen],lineLength,handle);
386  if (eof == NULL) {
387  fprintf(stderr, "%s", "Unexpected end of matrix file.");
388  EPETRA_CHK_ERR(-1);
389  }
390  rlen += strlen(&buffer[rlen]);
391  }
392  buffer[rlen++]='\n';
393  }
394  comm.Broadcast(&rlen, 1, 0);
395  comm.Broadcast(buffer, rlen, 0);
396 
397  buffer[rlen++] = '\0';
398  nread += nchunk;
399 
400  // Process the received data, saving non-zeros belonging on this proc.
401  char *lineptr = buffer;
402  for (rlen = 0; rlen < nchunk; rlen++) {
403  char *next = strchr(lineptr,'\n');
404  int I = atoi(strtok(lineptr," \t\n")) + ioffset;
405  int J = atoi(strtok(NULL," \t\n")) + joffset;
406  double V = atof(strtok(NULL," \t\n"));
407  lineptr = next + 1;
408  if (transpose) {
409  // swap I and J indices.
410  int tmp = I;
411  I = J;
412  J = tmp;
413  }
414  if (rowMap1->MyGID(I)) {
415  // This processor keeps this non-zero.
416  if (lnz >= localsize) {
417  // Need more memory to store nonzeros.
418  localsize += localblock;
419  iv = (int *) realloc(iv, localsize * sizeof(int));
420  jv = (int *) realloc(jv, localsize * sizeof(int));
421  vv = (double *) realloc(vv, localsize * sizeof(double));
422  }
423  iv[lnz] = I;
424  jv[lnz] = J;
425  vv[lnz] = V;
426  lnz++;
427  if (I < prevrow) rowmajor = 0;
428  prevrow = I;
429  }
430  }
431 
432  // Status check
433  if (nread / 1000000 > nmillion) {
434  nmillion++;
435  if (verbose && me == 0) std::cout << nmillion << "M ";
436  }
437  }
438 
439  delete [] buffer;
440 
441  if (!rowmajor) {
442  // Sort into row-major order (by iv) so can insert entire rows at once.
443  // Reorder jv and vv to parallel iv.
444  if (verbose && me == 0) std::cout << std::endl << " Sorting local nonzeros" << std::endl;
445  sort_three(iv, jv, vv, 0, lnz-1);
446  }
447 
448  // Compute number of entries per local row for use in constructor;
449  // saves reallocs in FillComplete.
450 
451  // Now construct the matrix.
452  // Count number of entries in each row so can use StaticProfile=2.
453  if (verbose && me == 0) std::cout << std::endl << " Constructing the matrix" << std::endl;
454  int numRows = rowMap1->NumMyElements();
455  int *numNonzerosPerRow = new int[numRows];
456  for (int i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
457  for (int i = 0; i < lnz; i++)
458  numNonzerosPerRow[rowMap1->LID(iv[i])]++;
459 
460  if (rowMap!=0 && colMap !=0)
461  A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
462  else if (rowMap!=0) {
463  // Construct with StaticProfile=true since we know numNonzerosPerRow.
464  // Less memory will be needed in FillComplete.
465  A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
466  }
467  else {
468  // Construct with StaticProfile=true since we know numNonzerosPerRow.
469  // Less memory will be needed in FillComplete.
470  A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
471  }
472  A->SetTracebackMode(1);
473 
474  // Rows are inserted in ascending global number, and the insertion uses numNonzerosPerRow.
475  // Therefore numNonzerosPerRow must be permuted, as it was created in ascending local order.
476  int *gidList = new int[numRows];
477  for (int i = 0; i < numRows; i++) gidList[i] = rowMap1->GID(i);
478  Epetra_Util::Sort(true,numRows,gidList,0,0,1,&numNonzerosPerRow);
479  delete [] gidList;
480 
481  // Insert the global values into the matrix row-by-row.
482  if (verbose && me == 0) std::cout << " Inserting global values" << std::endl;
483  {
484  int i = 0;
485  for (int sum = 0; i < numRows; i++) {
486  if (numNonzerosPerRow[i]) {
487  int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i],
488  &vv[sum], &jv[sum]);
489  if (ierr<0) EPETRA_CHK_ERR(ierr);
490  sum += numNonzerosPerRow[i];
491  }
492  }
493  }
494 
495  delete [] numNonzerosPerRow;
496  free(iv);
497  free(jv);
498  free(vv);
499 
500  if (verbose && me == 0) std::cout << " Completing matrix fill" << std::endl;
501  if (rangeMap != 0 && domainMap != 0) {
502  EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
503  }
504  else if (M!=N) {
505  Epetra_Map newDomainMap(N,rowMap1->IndexBase(), comm);
506  EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
507  }
508  else {
509  EPETRA_CHK_ERR(A->FillComplete());
510  }
511 
512  if (allocatedHere) delete rowMap1;
513 
514  if (handle!=0) fclose(handle);
515  double dt = timer.ElapsedTime();
516  if (verbose && me == 0) std::cout << "File Read time (secs): " << dt << std::endl;
517  return(0);
518 }
519 #endif
520 
521 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
522 int MatrixMarketFileToCrsMatrixHandle64(const char *filename,
523  const Epetra_Comm & comm,
524  Epetra_CrsMatrix * & A,
525  const Epetra_Map * rowMap,
526  const Epetra_Map * colMap,
527  const Epetra_Map * rangeMap,
528  const Epetra_Map * domainMap,
529  const bool transpose,
530  const bool verbose
531 )
532 {
533  const int chunk_read = 500000; // Modify this variable to change the
534  // size of the chunks read from the file.
535  const int headerlineLength = 257;
536  const int lineLength = 81;
537  const int tokenLength = 35;
538  char line[headerlineLength];
539  char token1[tokenLength];
540  char token2[tokenLength];
541  char token3[tokenLength];
542  char token4[tokenLength];
543  char token5[tokenLength];
544  long long M, N, NZ; // Matrix dimensions
545  int me = comm.MyPID();
546 
547  Epetra_Time timer(comm);
548 
549  // Make sure domain and range maps are either both defined or undefined
550  if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
551  EPETRA_CHK_ERR(-3);
552  }
553 
554  // check maps to see if domain and range are 1-to-1
555 
556  if (domainMap!=0) {
557  if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
558  if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
559  }
560  else {
561  // If domain and range are not specified,
562  // rowMap becomes both and must be 1-to-1 if specified
563  if (rowMap!=0) {
564  if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
565  }
566  }
567 
568  FILE * handle = 0;
569  if (me == 0) {
570  if (verbose) std::cout << "Reading MatrixMarket file " << filename << std::endl;
571  handle = fopen(filename,"r"); // Open file
572  if (handle == 0)
573  EPETRA_CHK_ERR(-1); // file not found
574 
575  // Check first line, which should be
576  // %%MatrixMarket matrix coordinate real general
577  if(fgets(line, headerlineLength, handle)==0) {
578  if (handle!=0) fclose(handle);
579  EPETRA_CHK_ERR(-1);
580  }
581  if(sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
582  if (handle!=0) fclose(handle);
583  EPETRA_CHK_ERR(-1);
584  }
585  if (strcmp(token1, "%%MatrixMarket") ||
586  strcmp(token2, "matrix") ||
587  strcmp(token3, "coordinate") ||
588  strcmp(token4, "real") ||
589  strcmp(token5, "general")) {
590  if (handle!=0) fclose(handle);
591  EPETRA_CHK_ERR(-1);
592  }
593 
594  // Next, strip off header lines (which start with "%")
595  do {
596  if(fgets(line, headerlineLength, handle)==0) {
597  if (handle!=0) fclose(handle);
598  EPETRA_CHK_ERR(-1);
599  }
600  } while (line[0] == '%');
601 
602  // Next get problem dimensions: M, N, NZ
603  if(sscanf(line, "%lld %lld %lld", &M, &N, &NZ)==0) {
604  if (handle!=0) fclose(handle);
605  EPETRA_CHK_ERR(-1);
606  }
607  }
608  comm.Broadcast(&M, 1, 0);
609  comm.Broadcast(&N, 1, 0);
610  comm.Broadcast(&NZ, 1, 0);
611 
612  // Now create matrix using user maps if provided.
613 
614 
615  // Now read in chunks of triplets and broadcast them to other processors.
616  char *buffer = new char[chunk_read*lineLength];
617  long long nchunk;
618  int nmillion = 0;
619  long long nread = 0;
620  int rlen;
621 
622  // Storage for this processor's nonzeros.
623  const int localblock = 100000;
624  int localsize = (int) (NZ / comm.NumProc()) + localblock;
625  long long *iv = (long long *) malloc(localsize * sizeof(long long));
626  long long *jv = (long long *) malloc(localsize * sizeof(long long));
627  double *vv = (double *) malloc(localsize * sizeof(double));
628  int lnz = 0; // Number of non-zeros on this processor.
629 
630  if (!iv || !jv || !vv)
631  EPETRA_CHK_ERR(-1);
632 
633  Epetra_Map *rowMap1;
634  bool allocatedHere=false;
635  if (rowMap != 0)
636  rowMap1 = const_cast<Epetra_Map *>(rowMap);
637  else {
638  rowMap1 = new Epetra_Map(M, 0, comm);
639  allocatedHere = true;
640  }
641  long long ioffset = rowMap1->IndexBase64()-1;
642  long long joffset = (colMap != 0 ? colMap->IndexBase64()-1 : ioffset);
643 
644  int rowmajor = 1; // Assume non-zeros are listed in row-major order,
645  long long prevrow = -1; // but test to detect otherwise. If non-zeros
646  // are row-major, we can avoid the sort.
647 
648  while (nread < NZ) {
649  if (NZ-nread > chunk_read) nchunk = chunk_read;
650  else nchunk = NZ - nread;
651 
652  if (me == 0) {
653  char *eof;
654  rlen = 0;
655  for (int i = 0; i < nchunk; i++) {
656  eof = fgets(&buffer[rlen],lineLength,handle);
657  if (eof == NULL) {
658  fprintf(stderr, "%s", "Unexpected end of matrix file.");
659  EPETRA_CHK_ERR(-1);
660  }
661  rlen += strlen(&buffer[rlen]);
662  }
663  buffer[rlen++]='\n';
664  }
665  comm.Broadcast(&rlen, 1, 0);
666  comm.Broadcast(buffer, rlen, 0);
667 
668  buffer[rlen++] = '\0';
669  nread += nchunk;
670 
671  // Process the received data, saving non-zeros belonging on this proc.
672  char *lineptr = buffer;
673  for (rlen = 0; rlen < nchunk; rlen++) {
674  char *next = strchr(lineptr,'\n');
675  char *endp;
676  const int base = 10;
677 #if defined(_MSC_VER)
678  long long I = _strtoi64(strtok(lineptr," \t\n"), &endp, base) + ioffset;
679  long long J = _strtoi64(strtok(NULL," \t\n"), &endp, base) + joffset;
680 #else
681 #if defined(__PGI)
682  long long I, J;
683  std::istringstream ssI(strtok(lineptr," \t\n"));
684  ssI >> I; I += ioffset;
685  std::istringstream ssJ(strtok(NULL," \t\n"));
686  ssJ >> J; J += joffset;
687 #else
688  long long I = strtoll(strtok(lineptr," \t\n"), &endp, base) + ioffset;
689  long long J = strtoll(strtok(NULL," \t\n"), &endp, base) + joffset;
690 #endif
691 #endif
692  double V = atof(strtok(NULL," \t\n"));
693  lineptr = next + 1;
694  if (transpose) {
695  // swap I and J indices.
696  long long tmp = I;
697  I = J;
698  J = tmp;
699  }
700  if (rowMap1->MyGID(I)) {
701  // This processor keeps this non-zero.
702  if (lnz >= localsize) {
703  // Need more memory to store nonzeros.
704  localsize += localblock;
705  iv = (long long *) realloc(iv, localsize * sizeof(long long));
706  jv = (long long *) realloc(jv, localsize * sizeof(long long));
707  vv = (double *) realloc(vv, localsize * sizeof(double));
708  }
709  iv[lnz] = I;
710  jv[lnz] = J;
711  vv[lnz] = V;
712  lnz++;
713  if (I < prevrow) rowmajor = 0;
714  prevrow = I;
715  }
716  }
717 
718  // Status check
719  if (nread / 1000000 > nmillion) {
720  nmillion++;
721  if (verbose && me == 0) std::cout << nmillion << "M ";
722  }
723  }
724 
725  delete [] buffer;
726 
727  if (!rowmajor) {
728  // Sort into row-major order (by iv) so can insert entire rows at once.
729  // Reorder jv and vv to parallel iv.
730  if (verbose && me == 0) std::cout << std::endl << " Sorting local nonzeros" << std::endl;
731  sort_three(iv, jv, vv, 0, lnz-1);
732  }
733 
734  // Compute number of entries per local row for use in constructor;
735  // saves reallocs in FillComplete.
736 
737  // Now construct the matrix.
738  // Count number of entries in each row so can use StaticProfile=2.
739  if (verbose && me == 0) std::cout << std::endl << " Constructing the matrix" << std::endl;
740  int numRows = rowMap1->NumMyElements();
741  int *numNonzerosPerRow = new int[numRows];
742  for (int i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
743  for (int i = 0; i < lnz; i++)
744  numNonzerosPerRow[rowMap1->LID(iv[i])]++;
745 
746  if (rowMap!=0 && colMap !=0)
747  A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
748  else if (rowMap!=0) {
749  // Construct with StaticProfile=true since we know numNonzerosPerRow.
750  // Less memory will be needed in FillComplete.
751  A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
752  }
753  else {
754  // Construct with StaticProfile=true since we know numNonzerosPerRow.
755  // Less memory will be needed in FillComplete.
756  A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
757  }
758  A->SetTracebackMode(1);
759 
760  // Rows are inserted in ascending global number, and the insertion uses numNonzerosPerRow.
761  // Therefore numNonzerosPerRow must be permuted, as it was created in ascending local order.
762  long long *gidList = new long long[numRows];
763  for (int i = 0; i < numRows; i++) gidList[i] = rowMap1->GID64(i);
764  Epetra_Util::Sort(true,numRows,gidList,0,0,1,&numNonzerosPerRow,0,0);
765  delete [] gidList;
766 
767  // Insert the global values into the matrix row-by-row.
768  if (verbose && me == 0) std::cout << " Inserting global values" << std::endl;
769  {
770  int i = 0;
771  for (int sum = 0; i < numRows; i++) {
772  if (numNonzerosPerRow[i]) {
773  int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i],
774  &vv[sum], &jv[sum]);
775  if (ierr<0) EPETRA_CHK_ERR(ierr);
776  sum += numNonzerosPerRow[i];
777  }
778  }
779  }
780 
781  delete [] numNonzerosPerRow;
782  free(iv);
783  free(jv);
784  free(vv);
785 
786  if (verbose && me == 0) std::cout << " Completing matrix fill" << std::endl;
787  if (rangeMap != 0 && domainMap != 0) {
788  EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
789  }
790  else if (M!=N) {
791  Epetra_Map newDomainMap(N,rowMap1->IndexBase64(), comm);
792  EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
793  }
794  else {
795  EPETRA_CHK_ERR(A->FillComplete());
796  }
797 
798  if (allocatedHere) delete rowMap1;
799 
800  if (handle!=0) fclose(handle);
801  double dt = timer.ElapsedTime();
802  if (verbose && me == 0) std::cout << "File Read time (secs): " << dt << std::endl;
803  return(0);
804 }
805 #endif
806 
808 // Sorting values in array list in increasing order. Criteria is int.
809 // Also rearrange values in arrays parlista and partlistb
810 // to match the new order of list.
811 
812 template<typename int_type>
814  int_type *list, int_type *parlista, double *parlistb,
815  int start, int end, int *equal, int *larger)
816 {
817 int i;
818 int_type key, itmp;
819 double dtmp;
820 
821  key = list ? list[(end+start)/2] : 1;
822 
823  *equal = *larger = start;
824  for (i = start; i <= end; i++)
825  if (list[i] < key) {
826  itmp = parlista[i];
827  parlista[i] = parlista[*larger];
828  parlista[(*larger)] = parlista[*equal];
829  parlista[(*equal)] = itmp;
830  dtmp = parlistb[i];
831  parlistb[i] = parlistb[*larger];
832  parlistb[(*larger)] = parlistb[*equal];
833  parlistb[(*equal)] = dtmp;
834  itmp = list[i];
835  list[i] = list[*larger];
836  list[(*larger)++] = list[*equal];
837  list[(*equal)++] = itmp;
838  }
839  else if (list[i] == key) {
840  itmp = parlista[i];
841  parlista[i] = parlista[*larger];
842  parlista[(*larger)] = itmp;
843  dtmp = parlistb[i];
844  parlistb[i] = parlistb[*larger];
845  parlistb[(*larger)] = dtmp;
846  list[i] = list[*larger];
847  list[(*larger)++] = key;
848  }
849 }
850 
852 template<typename int_type>
853 static void sort_three(int_type* list, int_type *parlista, double *parlistb,
854  int start, int end)
855 {
856 int equal, larger;
857 
858  if (start < end) {
859  quickpart_list_inc_int(list, parlista, parlistb, start, end,
860  &equal, &larger);
861  sort_three(list, parlista, parlistb, start, equal-1);
862  sort_three(list, parlista, parlistb, larger, end);
863  }
864 }
865 
867 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
868 int MatlabFileToCrsMatrix(const char *filename,
869  const Epetra_Comm & comm,
870  Epetra_CrsMatrix * & A)
871 {
872  const int lineLength = 1025;
873  char line[lineLength];
874  int I, J;
875  double V;
876 
877  FILE * handle = 0;
878 
879  handle = fopen(filename,"r"); // Open file
880  if (handle == 0)
881  EPETRA_CHK_ERR(-1); // file not found
882 
883  int numGlobalRows = 0;
884  int numGlobalCols = 0;
885  while(fgets(line, lineLength, handle)!=0) {
886  if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
887  if (I>numGlobalRows) numGlobalRows = I;
888  if (J>numGlobalCols) numGlobalCols = J;
889  }
890 
891  if (handle!=0) fclose(handle);
892  Epetra_Map rangeMap(numGlobalRows, 0, comm);
893  Epetra_Map domainMap(numGlobalCols, 0, comm);
894  A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
895 
896  // Now read in each triplet and store to the local portion of the matrix if the row is owned.
897  const Epetra_Map & rowMap1 = A->RowMap();
898 
899  handle = 0;
900 
901  handle = fopen(filename,"r"); // Open file
902  if (handle == 0)
903  EPETRA_CHK_ERR(-1); // file not found
904 
905  while (fgets(line, lineLength, handle)!=0) {
906  if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
907  I--; J--; // Convert to Zero based
908  if (rowMap1.MyGID(I)) {
909  int ierr = A->InsertGlobalValues(I, 1, &V, &J);
910  if (ierr<0) EPETRA_CHK_ERR(ierr);
911  }
912  }
913 
914  EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
915 
916  if (handle!=0) fclose(handle);
917  return(0);
918 }
919 #endif
920 
921 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
922 int MatlabFileToCrsMatrix64(const char *filename,
923  const Epetra_Comm & comm,
924  Epetra_CrsMatrix * & A)
925 {
926  const int lineLength = 1025;
927  char line[lineLength];
928  long long I, J;
929  double V;
930 
931  FILE * handle = 0;
932 
933  handle = fopen(filename,"r"); // Open file
934  if (handle == 0)
935  EPETRA_CHK_ERR(-1); // file not found
936 
937  long long numGlobalRows = 0;
938  long long numGlobalCols = 0;
939  while(fgets(line, lineLength, handle)!=0) {
940  if(sscanf(line, "%lld %lld %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
941  if (I>numGlobalRows) numGlobalRows = I;
942  if (J>numGlobalCols) numGlobalCols = J;
943  }
944 
945  if (handle!=0) fclose(handle);
946  Epetra_Map rangeMap(numGlobalRows, 0, comm);
947  Epetra_Map domainMap(numGlobalCols, 0, comm);
948  A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
949 
950  // Now read in each triplet and store to the local portion of the matrix if the row is owned.
951  const Epetra_Map & rowMap1 = A->RowMap();
952 
953  handle = 0;
954 
955  handle = fopen(filename,"r"); // Open file
956  if (handle == 0)
957  EPETRA_CHK_ERR(-1); // file not found
958 
959  while (fgets(line, lineLength, handle)!=0) {
960  if(sscanf(line, "%lld %lld %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
961  I--; J--; // Convert to Zero based
962  if (rowMap1.MyGID(I)) {
963  int ierr = A->InsertGlobalValues(I, 1, &V, &J);
964  if (ierr<0) EPETRA_CHK_ERR(ierr);
965  }
966  }
967 
968  EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
969 
970  if (handle!=0) fclose(handle);
971  return(0);
972 }
973 #endif
974 
975 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
976 int HypreFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix){
977  int MyPID = comm.MyPID();
978  // This double will be in the format we want for the extension besides the leading zero
979  double filePID = (double)MyPID/(double)100000;
980  std::ostringstream stream;
981  // Using setprecision() puts it in the std::string
982  stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
983  // Then just ignore the first character
984  std::string fileName(filename);
985  fileName += stream.str().substr(1,7);
986  // Open the file
987  std::ifstream file(fileName.c_str());
988  std::string line;
989  if(file.is_open()){
990  std::getline(file, line);
991  int ilower, iupper;
992  std::istringstream istream(line);
993  // The first line of the file has the beginning and ending rows
994  istream >> ilower;
995  istream >> iupper;
996  // Using those we can create a row map
997  Epetra_Map RowMap(-1, iupper-ilower+1, 0, comm);
998  Matrix = new Epetra_CrsMatrix(Copy, RowMap, 0);
999  int currRow = -1;
1000  int counter = 0;
1001  std::vector<int> indices;
1002  std::vector<double> values;
1003  while(!file.eof()){
1004  std::getline(file, line);
1005  std::istringstream lineStr(line);
1006  int row, col;
1007  double val;
1008  lineStr >> row;
1009  lineStr >> col;
1010  lineStr >> val;
1011  if(currRow == -1) currRow = row; // First line
1012  if(row == currRow){
1013  // add to the vector
1014  counter = counter + 1;
1015  indices.push_back(col);
1016  values.push_back(val);
1017  } else {
1018  Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
1019  indices.clear();
1020  values.clear();
1021  counter = 0;
1022  currRow = row;
1023  // make a new vector
1024  indices.push_back(col);
1025  values.push_back(val);
1026  counter = counter + 1;
1027  }
1028  }
1029  Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
1030  Matrix->Comm().Barrier();
1031  Matrix->FillComplete();
1032  file.close();
1033  return 0;
1034  } else {
1035  std::cout << "\nERROR:\nCouldn't open " << fileName << ".\n";
1036  return -1;
1037  }
1038 }
1039 #endif
1040 
1041 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1042 int HypreFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix){
1043  int MyPID = comm.MyPID();
1044  // This double will be in the format we want for the extension besides the leading zero
1045  double filePID = (double)MyPID/(double)100000;
1046  std::ostringstream stream;
1047  // Using setprecision() puts it in the std::string
1048  stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
1049  // Then just ignore the first character
1050  std::string fileName(filename);
1051  fileName += stream.str().substr(1,7);
1052  // Open the file
1053  std::ifstream file(fileName.c_str());
1054  std::string line;
1055  if(file.is_open()){
1056  std::getline(file, line);
1057  int ilower, iupper;
1058  std::istringstream istream(line);
1059  // The first line of the file has the beginning and ending rows
1060  istream >> ilower;
1061  istream >> iupper;
1062  // Using those we can create a row map
1063  Epetra_Map RowMap(-1LL, iupper-ilower+1, 0LL, comm);
1064  Matrix = new Epetra_CrsMatrix(Copy, RowMap, 0);
1065  long long currRow = -1;
1066  int counter = 0;
1067  std::vector<long long> indices;
1068  std::vector<double> values;
1069  while(!file.eof()){
1070  std::getline(file, line);
1071  std::istringstream lineStr(line);
1072  long long row, col;
1073  double val;
1074  lineStr >> row;
1075  lineStr >> col;
1076  lineStr >> val;
1077  if(currRow == -1) currRow = row; // First line
1078  if(row == currRow){
1079  // add to the vector
1080  counter = counter + 1;
1081  indices.push_back(col);
1082  values.push_back(val);
1083  } else {
1084  Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
1085  indices.clear();
1086  values.clear();
1087  counter = 0;
1088  currRow = row;
1089  // make a new vector
1090  indices.push_back(col);
1091  values.push_back(val);
1092  counter = counter + 1;
1093  }
1094  }
1095  Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
1096  Matrix->Comm().Barrier();
1097  Matrix->FillComplete();
1098  file.close();
1099  return 0;
1100  } else {
1101  std::cout << "\nERROR:\nCouldn't open " << fileName << ".\n";
1102  return -1;
1103  }
1104 }
1105 #endif
1106 
1107 } // namespace EpetraExt
1108 
int MatrixMarketFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
int HypreFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
static void quickpart_list_inc_int(int_type *list, int_type *parlista, double *parlistb, int start, int end, int *equal, int *larger)
int MatrixMarketFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
int MatlabFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
int MatrixMarketFileToCrsMatrixHandle(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A, const Epetra_Map *rowMap, const Epetra_Map *colMap, const Epetra_Map *rangeMap, const Epetra_Map *domainMap, const bool transpose, const bool verbose)
int MatlabFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
Constructs an Epetra_CrsMatrix object from a Matlab format file, distributes rows evenly across proce...
int HypreFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix)
Constructs an Epetra_CrsMatrix object from a Hypre Matrix Print command, the row map is specified...
int MatrixMarketFileToCrsMatrixHandle64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A, const Epetra_Map *rowMap, const Epetra_Map *colMap, const Epetra_Map *rangeMap, const Epetra_Map *domainMap, const bool transpose, const bool verbose)
static void sort_three(int_type *list, int_type *parlista, double *parlistb, int start, int end)