EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
43#include "Epetra_Comm.h"
44#include "Epetra_CrsMatrix.h"
45#include "Epetra_Map.h"
46#include "Epetra_IntVector.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
77using namespace EpetraExt;
78namespace EpetraExt {
79
80template<typename int_type>
81static void sort_three(int_type* list, int_type *parlista, double *parlistb,
82 int start, int end);
83
85#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
86int MatrixMarketFileToCrsMatrix(const char *filename,
87 const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
88{
90 return(0);
91}
92
94int MatrixMarketFileToCrsMatrix(const char *filename,
95 const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
96{
98 0, 0, 0, 0, transpose));
99 return(0);
100}
102
103int MatrixMarketFileToCrsMatrix(const char *filename,
104 const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
105 const bool verbose)
106{
108 0, 0, 0, 0,
109 transpose, verbose));
110 return(0);
111}
112
114int 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{
120 rowMap.Comm(), A,
121 &rowMap, 0,
122 &rangeMap, &domainMap,
123 transpose, verbose));
124 return(0);
125}
126
128int MatrixMarketFileToCrsMatrix(const char *filename,
129 const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
130 const bool verbose)
131{
133 rowMap.Comm(), A,
134 &rowMap, 0, 0, 0,
135 transpose, verbose));
136 return(0);
137}
138
140int 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{
146 rowMap.Comm(), A,
147 &rowMap, &colMap, 0, 0,
148 transpose, verbose));
149 return(0);
150}
151
153int 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{
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
169int MatrixMarketFileToCrsMatrix64(const char *filename,
170 const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
171{
173 return(0);
174}
175
177int MatrixMarketFileToCrsMatrix64(const char *filename,
178 const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
179{
181 0, 0, 0, 0, transpose));
182 return(0);
183}
185
186int MatrixMarketFileToCrsMatrix64(const char *filename,
187 const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
188 const bool verbose)
189{
191 0, 0, 0, 0,
192 transpose, verbose));
193 return(0);
194}
195
197int 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{
203 rowMap.Comm(), A,
204 &rowMap, 0,
205 &rangeMap, &domainMap,
206 transpose, verbose));
207 return(0);
208}
209
211int MatrixMarketFileToCrsMatrix64(const char *filename,
212 const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
213 const bool verbose)
214{
216 rowMap.Comm(), A,
217 &rowMap, 0, 0, 0,
218 transpose, verbose));
219 return(0);
220}
221
223int 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{
229 rowMap.Comm(), A,
230 &rowMap, &colMap, 0, 0,
231 transpose, verbose));
232 return(0);
233}
234
236int 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{
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
253int 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 {
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
533int 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 {
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
823template<typename int_type>
825 int_type *list, int_type *parlista, double *parlistb,
826 int start, int end, int *equal, int *larger)
827{
828int i;
829int_type key, itmp;
830double 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
863template<typename int_type>
864static void sort_three(int_type* list, int_type *parlista, double *parlistb,
865 int start, int end)
866{
867int 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
879int 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
933int 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
987int 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
1053int 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
#define EPETRA_CHK_ERR(a)
Copy
long long IndexBase64() const
int IndexBase() const
bool UniqueGIDs() const
const Epetra_Comm & Comm() const
const Epetra_Map & RowMap() const
int FillComplete(bool OptimizeDataStorage=true)
const Epetra_Comm & Comm() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
static void SetTracebackMode(int TracebackModeValue)
void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
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 HypreFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix)
int MatrixMarketFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
int MatlabFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
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.
static void sort_three(int_type *list, int_type *parlista, double *parlistb, int start, int end)
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 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)
int MatrixMarketFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
const int N
static void quickpart_list_inc_int(int_type *list, int_type *parlista, double *parlistb, int start, int end, int *equal, int *larger)