Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
MatrixMarket_TpetraNew.hpp
1
2#ifndef __MATRIXMARKET_TPETRA_NEW_HPP
3#define __MATRIXMARKET_TPETRA_NEW_HPP
4
6// Functions to read a MatrixMarket file and load it into a matrix.
7// Adapted from
8// Trilinos/packages/epetraext/src/inout/EpetraExt_CrsMatrixIn.cpp
9// Modified by Jon Berry and Karen Devine to make matrix reallocations
10// more efficient; rather than insert each non-zero separately, we
11// collect rows of non-zeros for insertion.
12// Modified by Karen Devine and Steve Plimpton to prevent all processors
13// from reading the same file at the same time (ouch!); chunks of the file
14// are read and broadcast by processor zero; each processor then extracts
15// its nonzeros from the chunk, sorts them by row, and inserts them into
16// the matrix.
17// The variable "chunkSize" can be changed to modify the size of the chunks
18// read from the file. Larger values of chunkSize lead to faster execution
19// and greater memory use.
21
22private:
23
24template <typename global_ordinal_type, typename scalar_type>
25static
26Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> >
27buildDistribution(
28 std::string &distribution, // string indicating whether to use 1D, 2D or
29 // file-based matrix distribution
30 size_t nRow, // Number of global matrix rows
31 size_t nCol, // Number of global matrix rows
32 const Teuchos::ParameterList &params, // Parameters to the file reading
33 const Teuchos::RCP<const Teuchos::Comm<int> > &comm
34 // communicator to be used in maps
35)
36{
37 // Function to build the sets of GIDs for 1D or 2D matrix distribution.
38 // Output is a Distribution object.
39
40 int me = comm->getRank();
41
42 using basedist_t = Distribution<global_ordinal_type,scalar_type>;
43 basedist_t *retval;
44
45 bool randomize = false; // Randomly permute rows and columns before storing
46 {
47 const Teuchos::ParameterEntry *pe = params.getEntryPtr("randomize");
48 if (pe != NULL)
49 randomize = pe->getValue<bool>(&randomize);
50 }
51
52 std::string partitionFile = ""; // File to reassign rows to parts
53 {
54 const Teuchos::ParameterEntry *pe = params.getEntryPtr("partitionFile");
55 if (pe != NULL)
56 partitionFile = pe->getValue<std::string>(&partitionFile);
57 }
58
59 if (distribution == "2D") { // Generate 2D distribution
60 if (partitionFile != "") {
61 // Generate 2D distribution with vector partition file
62 TEUCHOS_TEST_FOR_EXCEPTION(randomize, std::logic_error,
63 "Randomization not available for 2DVec distributions.");
64 Distribution2DVec<global_ordinal_type,scalar_type> *dist =
65 new Distribution2DVec<global_ordinal_type, scalar_type>(
66 nRow, comm, params, partitionFile);
67 retval = dynamic_cast<basedist_t *>(dist);
68 }
69 else if (randomize) {
70 // Randomly assign rows/columns to processors
71 Distribution2DRandom<global_ordinal_type,scalar_type> *dist =
72 new Distribution2DRandom<global_ordinal_type,scalar_type>(
73 nRow, comm, params);
74 retval = dynamic_cast<basedist_t *>(dist);
75 }
76 else {
77 // The vector will be distributed linearly, with extra vector entries
78 // spread among processors to maintain balance in rows and columns.
79 Distribution2DLinear<global_ordinal_type,scalar_type> *dist =
80 new Distribution2DLinear<global_ordinal_type,scalar_type>(
81 nRow, comm, params);
82 retval = dynamic_cast<basedist_t *>(dist);
83 }
84 }
85 else if (distribution == "1D") {
86 if (partitionFile != "") {
87 // Generate 1D row-based distribution with vector partition file
88 Distribution1DVec<global_ordinal_type,scalar_type> *dist =
89 new Distribution1DVec<global_ordinal_type,scalar_type>(
90 nRow, comm, params, partitionFile);
91 retval = dynamic_cast<basedist_t*>(dist);
92 }
93 else if (randomize) {
94 // Randomly assign rows to processors.
95 Distribution1DRandom<global_ordinal_type,scalar_type> *dist =
96 new Distribution1DRandom<global_ordinal_type,scalar_type>(
97 nRow, comm, params);
98 retval = dynamic_cast<basedist_t *>(dist);
99 }
100 else {
101 // Linear map similar to Trilinos default.
102 Distribution1DLinear<global_ordinal_type,scalar_type> *dist =
103 new Distribution1DLinear<global_ordinal_type,scalar_type>(
104 nRow, comm, params);
105 retval = dynamic_cast<basedist_t *>(dist);
106 }
107 }
108 else if (distribution == "LowerTriangularBlock") {
109 if (randomize && me == 0) {
110 std::cout << "WARNING: Randomization not available for "
111 << "LowerTriangularBlock distributions." << std::endl;
112 }
113
114 DistributionLowerTriangularBlock<global_ordinal_type,scalar_type> *dist =
115 new DistributionLowerTriangularBlock<global_ordinal_type,scalar_type>(
116 nRow, comm, params);
117 retval = dynamic_cast<basedist_t*>(dist);
118 }
119 else if (distribution == "MMFile") {
120 // Generate user-defined distribution from Matrix-Market file
121 if (randomize && me == 0) {
122 std::cout << "WARNING: Randomization not available for MMFile "
123 << "distributions." << std::endl;
124 }
125 DistributionMMFile<global_ordinal_type,scalar_type> *dist =
126 new DistributionMMFile<global_ordinal_type,scalar_type>(
127 nRow, comm, params);
128 retval = dynamic_cast<basedist_t*>(dist);
129 }
130
131 else {
132 std::cout << "ERROR: Invalid distribution method " << distribution
133 << std::endl;
134 exit(-1);
135 }
136 return Teuchos::rcp<basedist_t>(retval);
137}
138
139static
140void
141readMatrixMarket(
142 const std::string &filename, // MatrixMarket file to read
143 const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
144 const Teuchos::ParameterList &params,
145 size_t &nRow,
146 size_t &nCol,
147 typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t &localNZ,
148 Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist
149)
150{
151 bool useTimers = false; // should we time various parts of the reader?
152 {
153 const Teuchos::ParameterEntry *pe = params.getEntryPtr("useTimers");
154 if (pe != NULL)
155 useTimers = pe->getValue<bool>(&useTimers);
156 }
157
158 Teuchos::RCP<Teuchos::TimeMonitor> timer = Teuchos::null;
159 if (useTimers)
160 timer = rcp(new Teuchos::TimeMonitor(
161 *Teuchos::TimeMonitor::getNewTimer("RMM parameterSetup")));
162
163 int me = comm->getRank();
164
165 bool verbose = false; // Print status as reading
166 {
167 const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
168 if (pe != NULL)
169 verbose = pe->getValue<bool>(&verbose);
170 }
171
172 size_t chunkSize = 500000; // Number of lines to read / broadcast at once
173 {
174 const Teuchos::ParameterEntry *pe = params.getEntryPtr("chunkSize");
175 if (pe != NULL)
176 chunkSize = pe->getValue<size_t>(&chunkSize);
177 }
178
179 bool symmetrize = false; // Symmetrize the matrix
180 {
181 const Teuchos::ParameterEntry *pe = params.getEntryPtr("symmetrize");
182 if (pe != NULL)
183 symmetrize = pe->getValue<bool>(&symmetrize);
184 }
185
186 bool transpose = false; // Store the transpose
187 {
188 const Teuchos::ParameterEntry *pe = params.getEntryPtr("transpose");
189 if (pe != NULL)
190 transpose = pe->getValue<bool>(&transpose);
191 }
192
193 std::string diagonal = ""; // Are diagonal entries required (so add them)
194 // or ignored (so delete them) or neither?
195 // Default is neither.
196 {
197 const Teuchos::ParameterEntry *pe = params.getEntryPtr("diagonal");
198 if (pe != NULL)
199 diagonal = pe->getValue<std::string>(&diagonal);
200 }
201 bool ignoreDiagonal = (diagonal == "exclude");
202 bool requireDiagonal = (diagonal == "require");
203
204 std::string distribution = "1D"; // Default distribution is 1D row-based
205 {
206 const Teuchos::ParameterEntry *pe = params.getEntryPtr("distribution");
207 if (pe != NULL)
208 distribution = pe->getValue<std::string>(&distribution);
209 }
210
211 if (useTimers) {
212 timer = Teuchos::null;
213 timer = rcp(new Teuchos::TimeMonitor(
214 *Teuchos::TimeMonitor::getNewTimer("RMM header")));
215 }
216
217 if (verbose && me == 0)
218 std::cout << "Reading matrix market file... " << filename << std::endl;
219
220 FILE *fp = NULL;
221 size_t dim[3] = {0,0,0}; // nRow, nCol, nNz as read from MatrixMarket
222 MM_typecode mmcode;
223
224 if (me == 0) {
225
226 fp = fopen(filename.c_str(), "r");
227
228 if (fp == NULL) {
229 std::cout << "Error: cannot open file " << filename << std::endl;
230 }
231 else {
232 // Read MatrixMarket banner to get type of matrix
233 if (mm_read_banner(fp, &mmcode) != 0) {
234 std::cout << "Error: bad MatrixMarket banner " << std::endl;
235 }
236 else {
237 // Error check for file types that this function can read
238 if (!mm_is_valid(mmcode) ||
239 mm_is_dense(mmcode) || mm_is_array(mmcode) ||
240 mm_is_complex(mmcode) ||
241 mm_is_skew(mmcode) || mm_is_hermitian(mmcode)) {
242 std::cout << "Error: matrix type is not supported by this reader\n "
243 << "This reader supports sparse, coordinate, "
244 << "real, pattern, integer, symmetric, general"
245 << std::endl;
246 }
247 else {
248 // Read nRow, nCol, nNz from MatrixMarket file
249 if (mm_read_mtx_crd_size(fp, &dim[0], &dim[1], &dim[2]) != 0) {
250 std::cout << "Error: bad matrix dimensions " << std::endl;
251 dim[0] = dim[1] = dim[2] = 0;
252 }
253 }
254 }
255 }
256 }
257
258 // Broadcast relevant info
259 // Bad input if dim[0] or dim[1] still is zero after broadcast;
260 // all procs throw together
261 Teuchos::broadcast<int, size_t>(*comm, 0, 3, dim);
262 if (dim[0] == 0 || dim[1] == 0) {
263 throw std::runtime_error("Error: bad matrix header information");
264 }
265 Teuchos::broadcast<int, char>(*comm, 0, sizeof(MM_typecode), mmcode);
266
267 nRow = dim[0];
268 nCol = dim[1];
269
270 TEUCHOS_TEST_FOR_EXCEPTION(nRow != nCol, std::logic_error,
271 "This overload of readSparseFile requires nRow == nCol "
272 << "(nRow = " << nRow << ", nCol = " << nCol << "); "
273 << "For now, use a different overload of readSparseFile until this "
274 << "overload is fixed to read rectangular matrices. "
275 << "See Trilinos github issues #7045 and #8472.");
276
277 size_t nNz = dim[2];
278 bool patternInput = mm_is_pattern(mmcode);
279 bool symmetricInput = mm_is_symmetric(mmcode);
280 if (symmetricInput) symmetrize = true;
281
282 if (verbose && me == 0)
283 std::cout << "Matrix market file... "
284 << "\n symmetrize = " << symmetrize
285 << "\n transpose = " << transpose
286 << "\n change diagonal = " << diagonal
287 << "\n distribution = " << distribution
288 << std::endl;
289
290 if (useTimers) {
291 timer = Teuchos::null;
292 timer = rcp(new Teuchos::TimeMonitor(
293 *Teuchos::TimeMonitor::getNewTimer("RMM distribution")));
294 }
295
296 // Create distribution based on nRow, nCol, npRow, npCol
297 dist = buildDistribution<global_ordinal_type,scalar_type>(distribution,
298 nRow, nCol, params,
299 comm);
300 if (useTimers) {
301 timer = Teuchos::null;
302 timer = rcp(new Teuchos::TimeMonitor(
303 *Teuchos::TimeMonitor::getNewTimer("RMM readChunks")));
304 }
305
306 std::set<global_ordinal_type> diagset;
307 // If diagonal == require, this set keeps track of
308 // which diagonal entries already existing so we can
309 // add those that don't
310
311 using nzindex_t =
312 typename Distribution<global_ordinal_type,scalar_type>::NZindex_t;
313
314 // Chunk information and buffers
315 const int maxLineLength = 81;
316 char *buffer = new char[chunkSize*maxLineLength+1];
317 size_t nChunk;
318 size_t nMillion = 0;
319 size_t nRead = 0;
320 size_t rlen;
321
322 auto timerRead = Teuchos::TimeMonitor::getNewTimer("RMM readNonzeros");
323 auto timerSelect = Teuchos::TimeMonitor::getNewTimer("RMM selectNonzeros");
324 // Read chunks until the entire file is read
325 Teuchos::RCP<Teuchos::TimeMonitor> innerTimer = Teuchos::null;
326 while (nRead < nNz) {
327
328 innerTimer = rcp(new Teuchos::TimeMonitor(*timerRead));
329
330 if (nNz-nRead > chunkSize) nChunk = chunkSize;
331 else nChunk = (nNz - nRead);
332
333 // Processor 0 reads a chunk
334 if (me == 0) {
335 char *eof;
336 rlen = 0;
337 for (size_t i = 0; i < nChunk; i++) {
338 eof = fgets(&buffer[rlen],maxLineLength,fp);
339 if (eof == NULL) {
340 std::cout << "Unexpected end of matrix file." << std::endl;
341 std::cout.flush();
342 exit(-1);
343 }
344 rlen += strlen(&buffer[rlen]);
345 }
346 buffer[rlen++]='\n';
347 }
348
349 // Processor 0 broadcasts the chunk
350 Teuchos::broadcast<int, size_t>(*comm, 0, 1, &rlen);
351 Teuchos::broadcast<int, char>(*comm, 0, rlen, buffer);
352
353 buffer[rlen++] = '\0';
354 nRead += nChunk;
355
356 innerTimer = Teuchos::null;
357 innerTimer = rcp(new Teuchos::TimeMonitor(*timerSelect));
358
359 // All processors check the received data, saving nonzeros belonging to them
360 char *lineptr = buffer;
361 for (rlen = 0; rlen < nChunk; rlen++) {
362
363 char *next = strchr(lineptr,'\n');
364 global_ordinal_type I = atol(strtok(lineptr," \t\n"))
365 - 1; // since MatrixMarket files are one-based
366 global_ordinal_type J = atol(strtok(NULL," \t\n"))
367 - 1; // since MatrixMarket files are one-based
368
369 scalar_type V = (patternInput ? -1. : atof(strtok(NULL," \t\n")));
370 lineptr = next + 1;
371
372 // Special processing of nonzero
373 if ((I == J) && ignoreDiagonal) continue;
374
375 if (transpose) std::swap<global_ordinal_type>(I,J);
376
377 // Add nonzero (I,J) to the map if it should be on this processor
378 // Some file-based distributions have processor assignment stored as
379 // the non-zero's value, so pass the value to Mine.
380 if (dist->Mine(I,J,int(V))) {
381 nzindex_t idx = std::make_pair(I,J);
382 localNZ[idx] = V;
383 if (requireDiagonal && (I == J)) diagset.insert(I);
384 }
385
386 // If symmetrizing, add (J,I) to the map if it should be on this processor
387 // Some file-based distributions have processor assignment stored as
388 // the non-zero's value, so pass the value to Mine.
389 if (symmetrize && (I != J) && dist->Mine(J,I,int(V))) {
390 // Add entry (J, I) if need to symmetrize
391 // This processor keeps this non-zero.
392 nzindex_t idx = std::make_pair(J,I);
393 localNZ[idx] = V;
394 }
395 }
396
397 // Status check
398 if (verbose) {
399 if (nRead / 1000000 > nMillion) {
400 nMillion++;
401 if (me == 0) std::cout << nMillion << "M ";
402 }
403 }
404
405 innerTimer = Teuchos::null;
406 }
407
408 if (verbose && me == 0)
409 std::cout << std::endl << nRead << " nonzeros read " << std::endl;
410
411 if (fp != NULL) fclose(fp);
412 delete [] buffer;
413
414 if (useTimers) {
415 timer = Teuchos::null;
416 timer = rcp(new Teuchos::TimeMonitor(
417 *Teuchos::TimeMonitor::getNewTimer("RMM diagonal")));
418 }
419
420 // Add diagonal entries if they are required.
421 // Check to make sure they are all here; add them if they are not.
422 if (requireDiagonal) {
423 if (dist->DistType() == MMFile) {
424 // All diagonal entries should be present in the file; we cannot create
425 // them for file-based data distributions.
426 // Do an error check to ensure they are all there.
427 size_t localss = diagset.size();
428 size_t globalss;
429 Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, 1,
430 &localss, &globalss);
431 TEUCHOS_TEST_FOR_EXCEPTION(globalss != nRow, std::logic_error,
432 "File-based nonzero distribution requires all diagonal "
433 << "entries to be present in the file. \n"
434 << "Number of diagonal entries in file = " << globalss << "\n"
435 << "Number of matrix rows = " << nRow << "\n");
436 }
437 else {
438 for (global_ordinal_type i = 0;
439 i < static_cast<global_ordinal_type>(nRow); i++)
440 {
441 if (dist->Mine(i,i)) {
442 if (diagset.find(i) == diagset.end()) {
443 nzindex_t idx = std::make_pair(i,i);
444 localNZ[idx] = 0;
445 }
446 }
447 }
448 }
449 }
450 // Done with diagset; free its memory
451 std::set<global_ordinal_type>().swap(diagset);
452
453 if (useTimers)
454 timer = Teuchos::null;
455}
456
459// (This format is used by the upcoming readBinary function) //
461//
462// File format:
463// #vertices #edges src_1 dst_1 src_2 dst_2 ... src_#edges dst_#edges
464//
465// Types:
466// #edges: unsigned long long
467// everything else: unsigned int
468//
469// Base of indexing:
470// 1
471//
473//
474// A code example that converts MatrixMarket format into this format:
475//
476// typedef unsigned int ord_type;
477// typedef unsigned long long size_type;
478//
479// std::string line;
480// std::ifstream in(infilename);
481// std::ofstream out(outfilename, std::ios::out | std::ios::binary);
482//
483// ord_type nv, dummy, edge[2];
484// size_type ne;
485//
486// do
487// std::getline (in, line);
488// while(line[0] == '%');
489//
490// std::stringstream ss(line);
491// ss >> nv >> dummy >> ne;
492// out.write((char *)&nv, sizeof(ord_type));
493// out.write((char *)&ne, sizeof(size_type));
494//
495// while(std::getline(in, line)) {
496// std::stringstream ss2(line);
497// ss2 >> edge[0] >> edge[1];
498// out.write((char *)edge, sizeof(ord_type)*2);
499//
500// }
501// in.close();
502// out.close();
503//
505
506static
507void
508readBinary(
509 const std::string &filename, // MatrixMarket file to read
510 const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
511 const Teuchos::ParameterList &params,
512 size_t &nRow,
513 size_t &nCol,
514 typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t &localNZ,
515 Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist
516)
517{
518
519 int me = comm->getRank();
520
521 bool verbose = false; // Print status as reading
522 {
523 const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
524 if (pe != NULL)
525 verbose = pe->getValue<bool>(&verbose);
526 }
527
528 size_t chunkSize = 500000; // Number of lines to read / broadcast at once
529 {
530 const Teuchos::ParameterEntry *pe = params.getEntryPtr("chunkSize");
531 if (pe != NULL)
532 chunkSize = pe->getValue<size_t>(&chunkSize);
533 }
534
535 bool symmetrize = false; // Symmetrize the matrix
536 {
537 const Teuchos::ParameterEntry *pe = params.getEntryPtr("symmetrize");
538 if (pe != NULL)
539 symmetrize = pe->getValue<bool>(&symmetrize);
540 }
541
542 bool transpose = false; // Store the transpose
543 {
544 const Teuchos::ParameterEntry *pe = params.getEntryPtr("transpose");
545 if (pe != NULL)
546 transpose = pe->getValue<bool>(&transpose);
547 }
548
549 std::string diagonal = ""; // Are diagonal entries required (so add them)
550 // or ignored (so delete them) or neither?
551 // Default is neither.
552 {
553 const Teuchos::ParameterEntry *pe = params.getEntryPtr("diagonal");
554 if (pe != NULL)
555 diagonal = pe->getValue<std::string>(&diagonal);
556 }
557 bool ignoreDiagonal = (diagonal == "exclude");
558 bool requireDiagonal = (diagonal == "require");
559
560 std::string distribution = "1D"; // Default distribution is 1D row-based
561 {
562 const Teuchos::ParameterEntry *pe = params.getEntryPtr("distribution");
563 if (pe != NULL)
564 distribution = pe->getValue<std::string>(&distribution);
565 }
566
567 if (verbose && me == 0)
568 std::cout << "Reading binary file... " << filename << std::endl;
569
570 FILE *fp = NULL;
571 size_t dim[3] = {0,0,0}; // Expected to read nRow and nNz, nCol = nRow
572
573 if (me == 0) {
574
575 fp = fopen(filename.c_str(), "rb");
576
577 if (fp == NULL) {
578 std::cout << "Error: cannot open file " << filename << std::endl;
579 }
580 else {
581 // The header in the binary file contains only numVertices and numEdges
582 unsigned int nv = 0;
583 unsigned long long ne = 0;
584 if (fread(&nv, sizeof(unsigned int), 1, fp) != 1)
585 throw std::runtime_error("Error: bad number of read values.");
586 if (fread(&ne, sizeof(unsigned long long), 1, fp) != 1)
587 throw std::runtime_error("Error: bad number of read values.");
588 dim[0] = nv;
589 dim[1] = dim[0];
590 dim[2] = ne;
591 }
592 }
593
594 // Broadcast relevant info
595 // Bad input if dim[0] or dim[1] still is zero after broadcast;
596 // all procs throw together
597 Teuchos::broadcast<int, size_t>(*comm, 0, 3, dim);
598 if (dim[0] == 0 || dim[1] == 0) {
599 throw std::runtime_error("Error: bad matrix header information");
600 }
601
602 nRow = dim[0];
603 nCol = dim[1];
604 size_t nNz = dim[2];
605
606 if (verbose && me == 0)
607 std::cout << "Binary file... "
608 << "\n symmetrize = " << symmetrize
609 << "\n transpose = " << transpose
610 << "\n change diagonal = " << diagonal
611 << "\n distribution = " << distribution
612 << std::endl;
613
614 // Create distribution based on nRow, nCol, npRow, npCol
615 dist = buildDistribution<global_ordinal_type,scalar_type>(distribution,
616 nRow, nCol, params,
617 comm);
618
619 std::set<global_ordinal_type> diagset;
620 // If diagonal == require, this set keeps track of
621 // which diagonal entries already existing so we can
622 // add those that don't
623
624 using nzindex_t =
625 typename Distribution<global_ordinal_type,scalar_type>::NZindex_t;
626
627 // Chunk information and buffers
628 unsigned int *buffer = new unsigned int[chunkSize*2];
629 size_t nChunk;
630 size_t nMillion = 0;
631 size_t nRead = 0;
632 size_t rlen;
633 const scalar_type ONE = Teuchos::ScalarTraits<scalar_type>::one();
634
635
636 // Read chunks until the entire file is read
637 while (nRead < nNz) {
638 if (nNz-nRead > chunkSize) nChunk = chunkSize;
639 else nChunk = (nNz - nRead);
640
641 // Processor 0 reads a chunk
642 if (me == 0) {
643 size_t ret = 0;
644 rlen = 0;
645 for (size_t i = 0; i < nChunk; i++) {
646 ret = fread(&buffer[rlen], sizeof(unsigned int), 2, fp);
647 if (ret == 0) {
648 std::cout << "Unexpected end of matrix file." << std::endl;
649 std::cout.flush();
650 exit(-1);
651 }
652 rlen += 2;
653 }
654 }
655
656 // Processor 0 broadcasts the chunk
657 Teuchos::broadcast<int, size_t>(*comm, 0, 1, &rlen);
658 Teuchos::broadcast<int, unsigned int>(*comm, 0, rlen, buffer);
659
660 nRead += nChunk;
661
662 // All processors check the received data, saving nonzeros belonging to them
663 for (rlen = 0; rlen < nChunk; rlen++) {
664
665 global_ordinal_type I = buffer[2*rlen]-1;
666 global_ordinal_type J = buffer[2*rlen+1]-1;
667
668 // Special processing of nonzero
669 if ((I == J) && ignoreDiagonal) continue;
670
671 if (transpose) std::swap<global_ordinal_type>(I,J);
672
673 // Add nonzero (I,J) to the map if it should be on this processor
674 // Some file-based distributions have processor assignment stored as
675 // the non-zero's value, so pass the value to Mine.
676 if (dist->Mine(I,J,ONE)) {
677 nzindex_t idx = std::make_pair(I,J);
678 localNZ[idx] = ONE; // For now, the input binary format does not
679 // support numeric values, so we insert one.
680 if (requireDiagonal && (I == J)) diagset.insert(I);
681 }
682
683 // If symmetrizing, add (J,I) to the map if it should be on this processor
684 // Some file-based distributions have processor assignment stored as
685 // the non-zero's value, so pass the value to Mine.
686 if (symmetrize && (I != J) && dist->Mine(J,I,ONE)) {
687 // Add entry (J, I) if need to symmetrize
688 // This processor keeps this non-zero.
689 nzindex_t idx = std::make_pair(J,I);
690 localNZ[idx] = ONE;
691 }
692 }
693
694 // Status check
695 if (verbose) {
696 if (nRead / 1000000 > nMillion) {
697 nMillion++;
698 if (me == 0) std::cout << nMillion << "M ";
699 }
700 }
701 }
702
703 if (verbose && me == 0)
704 std::cout << std::endl << nRead << " nonzeros read " << std::endl;
705
706 if (fp != NULL) fclose(fp);
707 delete [] buffer;
708
709 // Add diagonal entries if they are required.
710 // Check to make sure they are all here; add them if they are not.
711 if (requireDiagonal) {
712 if (dist->DistType() == MMFile) {
713 // All diagonal entries should be present in the file; we cannot create
714 // them for file-based data distributions.
715 // Do an error check to ensure they are all there.
716 size_t localss = diagset.size();
717 size_t globalss;
718 Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, 1,
719 &localss, &globalss);
720 TEUCHOS_TEST_FOR_EXCEPTION(globalss != nRow, std::logic_error,
721 "File-based nonzero distribution requires all diagonal "
722 << "entries to be present in the file. \n"
723 << "Number of diagonal entries in file = " << globalss << "\n"
724 << "Number of matrix rows = " << nRow << "\n");
725 }
726 else {
727 for (global_ordinal_type i = 0;
728 i < static_cast<global_ordinal_type>(nRow); i++)
729 {
730 if (dist->Mine(i,i)) {
731 if (diagset.find(i) == diagset.end()) {
732 nzindex_t idx = std::make_pair(i,i);
733 localNZ[idx] = 0;
734 }
735 }
736 }
737 }
738 }
739 // Done with diagset; free its memory
740 std::set<global_ordinal_type>().swap(diagset);
741
742}
743
746// (This format is used by the upcoming readPerProcessBinary function) //
748//
749//
750// FILE FORMAT:
751// globalNumRows globalNumCols localNumNnzs row_1 col_1 ... row_n col_n,
752// where n = localNumNnzs
753//
754//
755// ASSUMPTIONS:
756// - The nonzeros should be sorted by their row indices within each file.
757// - Nonzeros have global row and column indices.
758// - The user specifies the base <filename>, where rank i reads <filename>.i.cooBin.
759//
760//
761// TYPES:
762// localNumNnzs: unsigned long long
763// everything else: unsigned int
764//
765//
766// BASE OF INDEXING: 1
767//
768//
770static
771void
772readPerProcessBinary(
773 const std::string &filename,
774 const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
775 const Teuchos::ParameterList &params,
776 size_t &nRow,
777 size_t &nCol,
778 typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t &localNZ,
779 Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist,
780 unsigned int* &buffer,
781 size_t &nNz
782)
783{
784
785 int me = comm->getRank();
786
787 bool verbose = false; // Print status as reading
788 {
789 const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
790 if (pe != NULL)
791 verbose = pe->getValue<bool>(&verbose);
792 }
793
794 std::string distribution = "1D"; // Default distribution is 1D row-based
795 {
796 const Teuchos::ParameterEntry *pe = params.getEntryPtr("distribution");
797 if (pe != NULL)
798 distribution = pe->getValue<std::string>(&distribution);
799 }
800
801 if (verbose && me == 0)
802 std::cout << "Reading per-process binary files... " << filename << std::endl;
803
804
805 std::string rankFileName = filename + "." + std::to_string(me) + ".cooBin";
806 FILE *fp = NULL;
807
808 fp = fopen(rankFileName.c_str(), "rb");
809 if (fp == NULL) {
810 std::cout << "Error: cannot open file " << filename << std::endl;
811 throw std::runtime_error("Error: non-existing input file: " + rankFileName);
812 }
813
814 // The header in each per-process file: globalNumRows globalNumCols localNumNonzeros
815 unsigned int globalNumRows = 0, globalNumCols = 0;
816 unsigned long long localNumNonzeros = 0;
817 if (fread(&globalNumRows, sizeof(unsigned int), 1, fp) != 1)
818 throw std::runtime_error("Error: bad number of read values.");
819 if (fread(&globalNumCols, sizeof(unsigned int), 1, fp) != 1)
820 throw std::runtime_error("Error: bad number of read values.");
821 if (fread(&localNumNonzeros, sizeof(unsigned long long), 1, fp) != 1)
822 throw std::runtime_error("Error: bad number of read values.");
823
824 nRow = static_cast<size_t>(globalNumRows);
825 nCol = static_cast<size_t>(globalNumCols);
826 nNz = static_cast<size_t>(localNumNonzeros);
827
828 // Fill the simple buffer array instead of a std::map
829 // S. Acer: With large graphs, we can't afford std::map
830 buffer = new unsigned int[nNz*2];
831
832 if(nNz > 0) {
833 size_t ret = fread(buffer, sizeof(unsigned int), 2*nNz, fp);
834 if (ret == 0) {
835 std::cout << "Unexpected end of matrix file: " << rankFileName << std::endl;
836 std::cout.flush();
837 delete [] buffer;
838 exit(-1);
839 }
840 }
841 if (fp != NULL) fclose(fp);
842
843 // This barrier is not necessary but useful for debugging
844 comm->barrier();
845 if(verbose && me == 0)
846 std::cout << "All ranks finished reading their nonzeros from their individual files\n";
847
848 // Create distribution based on nRow, nCol, npRow, npCol
849 dist = buildDistribution<global_ordinal_type,scalar_type>(distribution,
850 nRow, nCol, params,
851 comm);
852}
853
854public:
855
856// This is the default interface.
857static Teuchos::RCP<sparse_matrix_type>
858readSparseFile(
859 const std::string &filename, // MatrixMarket file to read
860 const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
861 const Teuchos::ParameterList &params
862)
863{
864 Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > dist;
865 return readSparseFile(filename, comm, params, dist);
866}
867
868// This version has the Distribution object as an output parameter.
869// S. Acer needs the distribution object to get the chunk cuts from
870// LowerTriangularBlock distribution.
871static Teuchos::RCP<sparse_matrix_type>
872readSparseFile(
873 const std::string &filename, // MatrixMarket file to read
874 const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
875 const Teuchos::ParameterList &params,
876 Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist
877)
878{
879 bool useTimers = false; // should we time various parts of the reader?
880 {
881 const Teuchos::ParameterEntry *pe = params.getEntryPtr("useTimers");
882 if (pe != NULL)
883 useTimers = pe->getValue<bool>(&useTimers);
884 }
885
886 Teuchos::RCP<Teuchos::TimeMonitor> timer = Teuchos::null;
887 if (useTimers)
888 timer = rcp(new Teuchos::TimeMonitor(
889 *Teuchos::TimeMonitor::getNewTimer("RSF parameterSetup")));
890
891 int me = comm->getRank();
892
893 // Check parameters to determine how to process the matrix while reading
894 // TODO: Add validators for the parameters
895
896 bool verbose = false; // Print status as reading
897 {
898 const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
899 if (pe != NULL)
900 verbose = pe->getValue<bool>(&verbose);
901 }
902
903 bool callFillComplete = true; // should we fillComplete the new CrsMatrix?
904 {
905 const Teuchos::ParameterEntry *pe = params.getEntryPtr("callFillComplete");
906 if (pe != NULL)
907 callFillComplete = pe->getValue<bool>(&callFillComplete);
908 }
909
910 // Don't want to use MatrixMarket's coordinate reader, because don't want
911 // entire matrix on one processor.
912 // Instead, Proc 0 reads nonzeros in chunks and broadcasts chunks to all
913 // processors.
914 // All processors insert nonzeros they own into a std::map
915
916 // Storage for this processor's nonzeros.
917 using localNZmap_t =
918 typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t;
919 localNZmap_t localNZ;
920
921 bool binary = false; // should we read a binary file?
922 {
923 const Teuchos::ParameterEntry *pe = params.getEntryPtr("binary");
924 if (pe != NULL)
925 binary = pe->getValue<bool>(&binary);
926 }
927
928 bool readPerProcess = false; // should we read a separate file per process?
929 {
930 const Teuchos::ParameterEntry *pe = params.getEntryPtr("readPerProcess");
931 if (pe != NULL)
932 readPerProcess = pe->getValue<bool>(&readPerProcess);
933 }
934
935 if (useTimers) {
936 const char *timername = (binary?"RSF readBinary":"RSF readMatrixMarket");
937 timer = Teuchos::null;
938 timer = rcp(new Teuchos::TimeMonitor(
939 *Teuchos::TimeMonitor::getNewTimer(timername)));
940 }
941
942 // Read nonzeros from the given file(s)
943 size_t nRow = 0, nCol = 0;
944 unsigned int *buffer=0; size_t nNz = 0;
945 if(binary){
946 if(readPerProcess)
947 readPerProcessBinary(filename, comm, params, nRow, nCol, localNZ, dist, buffer, nNz);
948 else
949 readBinary(filename, comm, params, nRow, nCol, localNZ, dist);
950 }
951 else
952 readMatrixMarket(filename, comm, params, nRow, nCol, localNZ, dist);
953
954 if(readPerProcess == false){
955
956 // Redistribute nonzeros as needed to satisfy the Distribution
957 // For most Distributions, this is a no-op
958 if (useTimers) {
959 timer = Teuchos::null;
960 timer = rcp(new Teuchos::TimeMonitor(
961 *Teuchos::TimeMonitor::getNewTimer("RSF redistribute")));
962 }
963
964 dist->Redistribute(localNZ);
965 }
966
967 if (useTimers) {
968 timer = Teuchos::null;
969 timer = rcp(new Teuchos::TimeMonitor(
970 *Teuchos::TimeMonitor::getNewTimer("RSF nonzerosConstruction")));
971 }
972
973 // Now construct the matrix.
974 // Count number of entries in each row for best use of StaticProfile
975 if (verbose && me == 0)
976 std::cout << std::endl << "Constructing the matrix" << std::endl;
977
978 Teuchos::Array<global_ordinal_type> rowIdx;
979 Teuchos::Array<size_t> nnzPerRow;
980 Teuchos::Array<global_ordinal_type> colIdx;
981 Teuchos::Array<scalar_type> val;
982 Teuchos::Array<size_t> offsets;
983
984 if(readPerProcess) {
985 global_ordinal_type prevI = std::numeric_limits<global_ordinal_type>::max();
986 for (size_t it = 0; it < nNz; it++){
987 global_ordinal_type I = buffer[2*it]-1;
988 global_ordinal_type J = buffer[2*it+1]-1;
989 if (prevI != I) {
990 prevI = I;
991 rowIdx.push_back(I);
992 nnzPerRow.push_back(0);
993 }
994 nnzPerRow.back()++;
995 colIdx.push_back(J);
996 }
997 delete [] buffer;
998 }
999 else {
1000 // Exploit fact that map has entries sorted by I, then J
1001 global_ordinal_type prevI = std::numeric_limits<global_ordinal_type>::max();
1002 for (auto it = localNZ.begin(); it != localNZ.end(); it++) {
1003 global_ordinal_type I = it->first.first;
1004 global_ordinal_type J = it->first.second;
1005 scalar_type V = it->second;
1006 if (prevI != I) {
1007 prevI = I;
1008 rowIdx.push_back(I);
1009 nnzPerRow.push_back(0);
1010 }
1011 nnzPerRow.back()++;
1012 colIdx.push_back(J);
1013 val.push_back(V);
1014 }
1015
1016 // Done with localNZ map; free it.
1017 localNZmap_t().swap(localNZ);
1018 }
1019
1020 // Compute prefix sum in offsets array
1021 offsets.resize(rowIdx.size() + 1);
1022 offsets[0] = 0;
1023 for (size_t row = 0; row < static_cast<size_t>(rowIdx.size()); row++)
1024 offsets[row+1] = offsets[row] + nnzPerRow[row];
1025
1026 if (useTimers) {
1027 timer = Teuchos::null;
1028 timer = rcp(new Teuchos::TimeMonitor(
1029 *Teuchos::TimeMonitor::getNewTimer("RSF insertNonzeros")));
1030 }
1031
1032 // Create a new RowMap with only rows having non-zero entries.
1033 size_t dummy = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
1034 Teuchos::RCP<const Tpetra::Map<> > rowMap =
1035 Teuchos::rcp(new Tpetra::Map<>(dummy, rowIdx(), 0, comm));
1036
1037 Teuchos::RCP<sparse_matrix_type> A =
1038 rcp(new sparse_matrix_type(rowMap, nnzPerRow()));
1039
1040 // Insert the global values into the matrix row-by-row.
1041 if (verbose && me == 0)
1042 std::cout << "Inserting global values" << std::endl;
1043
1044 if(readPerProcess){
1045 const scalar_type ONE = Teuchos::ScalarTraits<scalar_type>::one();
1046 for (int i = 0; i < rowIdx.size(); i++) {
1047 size_t nnz = nnzPerRow[i];
1048 size_t off = offsets[i];
1049 val.assign(nnz, ONE);
1050 // ReadPerProcess routine does not read any numeric values from the file,
1051 // So we insert dummy values here.
1052 A->insertGlobalValues(rowIdx[i], colIdx(off, nnz), val());
1053 }
1054 }
1055 else {
1056 for (int i = 0; i < rowIdx.size(); i++) {
1057 size_t nnz = nnzPerRow[i];
1058 size_t off = offsets[i];
1059 A->insertGlobalValues(rowIdx[i], colIdx(off, nnz), val(off, nnz));
1060 }
1061 }
1062
1063 // free memory that is no longer needed
1064 Teuchos::Array<size_t>().swap(nnzPerRow);
1065 Teuchos::Array<size_t>().swap(offsets);
1066 Teuchos::Array<global_ordinal_type>().swap(rowIdx);
1067 Teuchos::Array<global_ordinal_type>().swap(colIdx);
1068 Teuchos::Array<scalar_type>().swap(val);
1069
1070 if (useTimers)
1071 timer = Teuchos::null;
1072
1073 if (callFillComplete) {
1074
1075 if (verbose && me == 0)
1076 std::cout << "Building vector maps" << std::endl;
1077
1078 if (useTimers) {
1079 timer = Teuchos::null;
1080 timer = rcp(new Teuchos::TimeMonitor(
1081 *Teuchos::TimeMonitor::getNewTimer("RSF buildVectorMaps")));
1082 }
1083
1084 // Build domain map that corresponds to distribution
1085 Teuchos::Array<global_ordinal_type> vectorSet;
1086 for (global_ordinal_type i = 0;
1087 i < static_cast<global_ordinal_type>(nCol); i++)
1088 if (dist->VecMine(i)) vectorSet.push_back(i);
1089
1090 dummy = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
1091 Teuchos::RCP<const Tpetra::Map<> > domainMap =
1092 Teuchos::rcp(new Tpetra::Map<>(dummy, vectorSet(), 0, comm));
1093
1094 Teuchos::Array<global_ordinal_type>().swap(vectorSet);
1095
1096 // Build range map that corresponds to distribution
1097 for (global_ordinal_type i = 0;
1098 i < static_cast<global_ordinal_type>(nRow); i++)
1099 if (dist->VecMine(i)) vectorSet.push_back(i);
1100
1101 dummy = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
1102 Teuchos::RCP<const Tpetra::Map<> > rangeMap =
1103 Teuchos::rcp(new Tpetra::Map<>(dummy, vectorSet(), 0, comm));
1104
1105 Teuchos::Array<global_ordinal_type>().swap(vectorSet);
1106
1107 // FillComplete the matrix
1108 if (useTimers) {
1109 timer = Teuchos::null;
1110 timer = rcp(new Teuchos::TimeMonitor(
1111 *Teuchos::TimeMonitor::getNewTimer("RSF fillComplete")));
1112 }
1113
1114 if (verbose && me == 0)
1115 std::cout << "Calling fillComplete" << std::endl;
1116
1117 A->fillComplete(domainMap, rangeMap);
1118
1119 if (useTimers)
1120 timer = Teuchos::null;
1121
1122 if (verbose) {
1123 std::cout << "\nRank " << A->getComm()->getRank()
1124 << " nRows " << A->getLocalNumRows()
1125 << " minRow " << A->getRowMap()->getMinGlobalIndex()
1126 << " maxRow " << A->getRowMap()->getMaxGlobalIndex()
1127 << "\nRank " << A->getComm()->getRank()
1128 << " nCols " << A->getLocalNumCols()
1129 << " minCol " << A->getColMap()->getMinGlobalIndex()
1130 << " maxCol " << A->getColMap()->getMaxGlobalIndex()
1131 << "\nRank " << A->getComm()->getRank()
1132 << " nDomain " << A->getDomainMap()->getLocalNumElements()
1133 << " minDomain " << A->getDomainMap()->getMinGlobalIndex()
1134 << " maxDomain " << A->getDomainMap()->getMaxGlobalIndex()
1135 << "\nRank " << A->getComm()->getRank()
1136 << " nRange " << A->getRangeMap()->getLocalNumElements()
1137 << " minRange " << A->getRangeMap()->getMinGlobalIndex()
1138 << " maxRange " << A->getRangeMap()->getMaxGlobalIndex()
1139 << "\nRank " << A->getComm()->getRank()
1140 << " nEntries " << A->getLocalNumEntries()
1141 << std::endl;
1142 }
1143 }
1144
1145 return A;
1146}
1147
1148
1149
1150#endif
1151
Struct that holds views of the contents of a CrsMatrix.