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