143 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
158 typedef typename SparseMatrixType::global_ordinal_type
194 typedef Teuchos::ArrayRCP<int>::size_type size_type;
206 static Teuchos::RCP<const map_type>
212 pComm, GloballyDistributed));
242 static Teuchos::RCP<const map_type>
243 makeRowMap(
const Teuchos::RCP<const map_type>&
pRowMap,
251 pComm, GloballyDistributed));
254 std::invalid_argument,
255 "The specified row map is not distributed, "
256 "but the given communicator includes more than one process (in "
258 <<
pComm->getSize() <<
" processes).");
260 "The specified row Map's communicator (pRowMap->getComm()) "
261 "differs from the given separately supplied communicator pComm.");
280 static Teuchos::RCP<const map_type>
281 makeDomainMap(
const Teuchos::RCP<const map_type>& pRangeMap,
289 if (numRows == numCols) {
292 return createUniformContigMapWithNode<LO, GO, NT>(numCols,
293 pRangeMap->getComm());
370 distribute(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
371 Teuchos::ArrayRCP<size_t>& myRowPtr,
372 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
373 Teuchos::ArrayRCP<scalar_type>& myValues,
374 const Teuchos::RCP<const map_type>& pRowMap,
375 Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
376 Teuchos::ArrayRCP<size_t>& rowPtr,
377 Teuchos::ArrayRCP<global_ordinal_type>& colInd,
378 Teuchos::ArrayRCP<scalar_type>& values,
379 const bool debug =
false) {
383 using Teuchos::ArrayRCP;
384 using Teuchos::ArrayView;
387 using Teuchos::CommRequest;
390 using Teuchos::receive;
393 const bool extraDebug =
false;
395 const int numProcs = pComm->getSize();
396 const int myRank = pComm->getRank();
397 const int rootRank = 0;
404 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
405 const size_type myNumRows = myRows.size();
406 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<size_t>(myNumRows) !=
407 pRowMap->getLocalNumElements(),
409 "pRowMap->getLocalElementList().size() = "
411 <<
" != pRowMap->getLocalNumElements() = "
412 << pRowMap->getLocalNumElements() <<
". "
413 "Please report this bug to the Tpetra developers.");
414 TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
416 "On Proc 0: numEntriesPerRow.size() = "
417 << numEntriesPerRow.size()
418 <<
" != pRowMap->getLocalElementList().size() = "
419 << myNumRows <<
". Please report this bug to the "
420 "Tpetra developers.");
424 myNumEntriesPerRow = arcp<size_t>(myNumRows);
426 if (myRank != rootRank) {
430 send(*pComm, myNumRows, rootRank);
431 if (myNumRows != 0) {
435 send(*pComm,
static_cast<int>(myNumRows),
436 myRows.getRawPtr(), rootRank);
446 receive(*pComm, rootRank,
447 static_cast<int>(myNumRows),
448 myNumEntriesPerRow.getRawPtr());
453 std::accumulate(myNumEntriesPerRow.begin(),
454 myNumEntriesPerRow.end(), 0);
460 myColInd = arcp<GO>(myNumEntries);
461 myValues = arcp<scalar_type>(myNumEntries);
462 if (myNumEntries > 0) {
465 receive(*pComm, rootRank,
466 static_cast<int>(myNumEntries),
467 myColInd.getRawPtr());
468 receive(*pComm, rootRank,
469 static_cast<int>(myNumEntries),
470 myValues.getRawPtr());
476 cerr <<
"-- Proc 0: Copying my data from global arrays" << endl;
480 for (size_type k = 0; k < myNumRows; ++k) {
481 const GO myCurRow = myRows[k];
483 myNumEntriesPerRow[k] = numEntriesInThisRow;
485 if (extraDebug && debug) {
486 cerr <<
"Proc " << pRowMap->getComm()->getRank()
487 <<
": myNumEntriesPerRow[0.." << (myNumRows - 1) <<
"] = [";
488 for (size_type k = 0; k < myNumRows; ++k) {
489 cerr << myNumEntriesPerRow[k];
490 if (k < myNumRows - 1) {
498 std::accumulate(myNumEntriesPerRow.begin(),
499 myNumEntriesPerRow.end(), 0);
501 cerr <<
"-- Proc 0: I own " << myNumRows <<
" rows and "
502 << myNumEntries <<
" entries" << endl;
504 myColInd = arcp<GO>(myNumEntries);
505 myValues = arcp<scalar_type>(myNumEntries);
513 for (size_type k = 0; k < myNumRows;
514 myCurPos += myNumEntriesPerRow[k], ++k) {
516 const GO myRow = myRows[k];
517 const size_t curPos = rowPtr[myRow];
520 if (curNumEntries > 0) {
521 ArrayView<GO> colIndView = colInd(curPos, curNumEntries);
522 ArrayView<GO> myColIndView = myColInd(myCurPos, curNumEntries);
523 std::copy(colIndView.begin(), colIndView.end(),
524 myColIndView.begin());
526 ArrayView<scalar_type> valuesView =
527 values(curPos, curNumEntries);
528 ArrayView<scalar_type> myValuesView =
529 myValues(myCurPos, curNumEntries);
530 std::copy(valuesView.begin(), valuesView.end(),
531 myValuesView.begin());
536 for (
int p = 1; p < numProcs; ++p) {
538 cerr <<
"-- Proc 0: Processing proc " << p << endl;
541 size_type theirNumRows = 0;
546 receive(*pComm, p, &theirNumRows);
548 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
549 << theirNumRows <<
" rows" << endl;
551 if (theirNumRows != 0) {
556 ArrayRCP<GO> theirRows = arcp<GO>(theirNumRows);
557 receive(*pComm, p, as<int>(theirNumRows),
558 theirRows.getRawPtr());
567 const global_size_t numRows = pRowMap->getGlobalNumElements();
568 const GO indexBase = pRowMap->getIndexBase();
569 bool theirRowsValid =
true;
570 for (size_type k = 0; k < theirNumRows; ++k) {
571 if (theirRows[k] < indexBase ||
572 as<global_size_t>(theirRows[k] - indexBase) >= numRows) {
573 theirRowsValid =
false;
576 if (!theirRowsValid) {
577 TEUCHOS_TEST_FOR_EXCEPTION(
578 !theirRowsValid, std::logic_error,
579 "Proc " << p <<
" has at least one invalid row index. "
580 "Here are all of them: "
581 << Teuchos::toString(theirRows()) <<
". Valid row index "
582 "range (zero-based): [0, "
583 << (numRows - 1) <<
"].");
598 ArrayRCP<size_t> theirNumEntriesPerRow;
599 theirNumEntriesPerRow = arcp<size_t>(theirNumRows);
600 for (size_type k = 0; k < theirNumRows; ++k) {
601 theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
608 send(*pComm,
static_cast<int>(theirNumRows),
609 theirNumEntriesPerRow.getRawPtr(), p);
613 std::accumulate(theirNumEntriesPerRow.begin(),
614 theirNumEntriesPerRow.end(), 0);
617 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
618 << theirNumEntries <<
" entries" << endl;
623 if (theirNumEntries == 0) {
632 ArrayRCP<GO> theirColInd(theirNumEntries);
633 ArrayRCP<scalar_type> theirValues(theirNumEntries);
641 for (size_type k = 0; k < theirNumRows;
642 theirCurPos += theirNumEntriesPerRow[k], k++) {
644 const GO theirRow = theirRows[k];
650 if (curNumEntries > 0) {
651 ArrayView<GO> colIndView =
652 colInd(curPos, curNumEntries);
653 ArrayView<GO> theirColIndView =
654 theirColInd(theirCurPos, curNumEntries);
655 std::copy(colIndView.begin(), colIndView.end(),
656 theirColIndView.begin());
658 ArrayView<scalar_type> valuesView =
659 values(curPos, curNumEntries);
660 ArrayView<scalar_type> theirValuesView =
661 theirValues(theirCurPos, curNumEntries);
662 std::copy(valuesView.begin(), valuesView.end(),
663 theirValuesView.begin());
670 send(*pComm,
static_cast<int>(theirNumEntries),
671 theirColInd.getRawPtr(), p);
672 send(*pComm,
static_cast<int>(theirNumEntries),
673 theirValues.getRawPtr(), p);
676 cerr <<
"-- Proc 0: Finished with proc " << p << endl;
684 numEntriesPerRow = null;
689 if (debug && myRank == 0) {
690 cerr <<
"-- Proc 0: About to fill in myRowPtr" << endl;
698 myRowPtr = arcp<size_t>(myNumRows + 1);
700 for (size_type k = 1; k < myNumRows + 1; ++k) {
701 myRowPtr[k] = myRowPtr[k - 1] + myNumEntriesPerRow[k - 1];
703 if (extraDebug && debug) {
704 cerr <<
"Proc " << Teuchos::rank(*(pRowMap->getComm()))
705 <<
": myRowPtr[0.." << myNumRows <<
"] = [";
706 for (size_type k = 0; k < myNumRows + 1; ++k) {
716 if (debug && myRank == 0) {
717 cerr <<
"-- Proc 0: Done with distribute" << endl;
734 static Teuchos::RCP<sparse_matrix_type>
735 makeMatrix(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
736 Teuchos::ArrayRCP<size_t>& myRowPtr,
737 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
738 Teuchos::ArrayRCP<scalar_type>& myValues,
739 const Teuchos::RCP<const map_type>& pRowMap,
740 const Teuchos::RCP<const map_type>& pRangeMap,
741 const Teuchos::RCP<const map_type>& pDomainMap,
742 const bool callFillComplete =
true) {
745 using Teuchos::ArrayView;
757 TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
758 "makeMatrix: myRowPtr array is null. "
759 "Please report this bug to the Tpetra developers.");
760 TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
761 "makeMatrix: domain map is null. "
762 "Please report this bug to the Tpetra developers.");
763 TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
764 "makeMatrix: range map is null. "
765 "Please report this bug to the Tpetra developers.");
766 TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
767 "makeMatrix: row map is null. "
768 "Please report this bug to the Tpetra developers.");
772 RCP<sparse_matrix_type> A =
777 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
778 const size_type myNumRows = myRows.size();
781 const GO indexBase = pRowMap->getIndexBase();
782 for (size_type i = 0; i < myNumRows; ++i) {
783 const size_type myCurPos = myRowPtr[i];
785 ArrayView<GO> curColInd = myColInd.view(myCurPos, curNumEntries);
786 ArrayView<scalar_type> curValues = myValues.view(myCurPos, curNumEntries);
789 for (size_type k = 0; k < curNumEntries; ++k) {
790 curColInd[k] += indexBase;
793 if (curNumEntries > 0) {
794 A->insertGlobalValues(myRows[i], curColInd, curValues);
801 myNumEntriesPerRow = null;
806 if (callFillComplete) {
807 A->fillComplete(pDomainMap, pRangeMap);
817 static Teuchos::RCP<sparse_matrix_type>
818 makeMatrix(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
819 Teuchos::ArrayRCP<size_t>& myRowPtr,
820 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
821 Teuchos::ArrayRCP<scalar_type>& myValues,
822 const Teuchos::RCP<const map_type>& pRowMap,
823 const Teuchos::RCP<const map_type>& pRangeMap,
824 const Teuchos::RCP<const map_type>& pDomainMap,
825 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
826 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams) {
829 using Teuchos::ArrayView;
841 TEUCHOS_TEST_FOR_EXCEPTION(
842 myRowPtr.is_null(), std::logic_error,
843 "makeMatrix: myRowPtr array is null. "
844 "Please report this bug to the Tpetra developers.");
845 TEUCHOS_TEST_FOR_EXCEPTION(
846 pDomainMap.is_null(), std::logic_error,
847 "makeMatrix: domain map is null. "
848 "Please report this bug to the Tpetra developers.");
849 TEUCHOS_TEST_FOR_EXCEPTION(
850 pRangeMap.is_null(), std::logic_error,
851 "makeMatrix: range map is null. "
852 "Please report this bug to the Tpetra developers.");
853 TEUCHOS_TEST_FOR_EXCEPTION(
854 pRowMap.is_null(), std::logic_error,
855 "makeMatrix: row map is null. "
856 "Please report this bug to the Tpetra developers.");
860 RCP<sparse_matrix_type> A =
866 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
867 const size_type myNumRows = myRows.size();
870 const GO indexBase = pRowMap->getIndexBase();
871 for (size_type i = 0; i < myNumRows; ++i) {
872 const size_type myCurPos = myRowPtr[i];
874 ArrayView<GO> curColInd = myColInd.view(myCurPos, curNumEntries);
875 ArrayView<scalar_type> curValues = myValues.view(myCurPos, curNumEntries);
878 for (size_type k = 0; k < curNumEntries; ++k) {
879 curColInd[k] += indexBase;
881 if (curNumEntries > 0) {
882 A->insertGlobalValues(myRows[i], curColInd, curValues);
889 myNumEntriesPerRow = null;
894 A->fillComplete(pDomainMap, pRangeMap, fillCompleteParams);
902 static Teuchos::RCP<sparse_matrix_type>
903 makeMatrix(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
904 Teuchos::ArrayRCP<size_t>& myRowPtr,
905 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
906 Teuchos::ArrayRCP<scalar_type>& myValues,
907 const Teuchos::RCP<const map_type>& rowMap,
908 Teuchos::RCP<const map_type>& colMap,
909 const Teuchos::RCP<const map_type>& domainMap,
910 const Teuchos::RCP<const map_type>& rangeMap,
911 const bool callFillComplete =
true) {
912 using Teuchos::ArrayView;
921 RCP<sparse_matrix_type> A;
922 if (colMap.is_null()) {
930 ArrayView<const GO> myRows = rowMap->getLocalElementList();
931 const size_type myNumRows = myRows.size();
934 const GO indexBase = rowMap->getIndexBase();
935 for (size_type i = 0; i < myNumRows; ++i) {
936 const size_type myCurPos = myRowPtr[i];
937 const size_type curNumEntries = as<size_type>(myNumEntriesPerRow[i]);
938 ArrayView<GO> curColInd = myColInd.view(myCurPos, curNumEntries);
939 ArrayView<scalar_type> curValues = myValues.view(myCurPos, curNumEntries);
942 for (size_type k = 0; k < curNumEntries; ++k) {
943 curColInd[k] += indexBase;
945 if (curNumEntries > 0) {
946 A->insertGlobalValues(myRows[i], curColInd, curValues);
953 myNumEntriesPerRow = null;
958 if (callFillComplete) {
959 A->fillComplete(domainMap, rangeMap);
960 if (colMap.is_null()) {
961 colMap = A->getColMap();
984 static Teuchos::RCP<const Teuchos::MatrixMarket::Banner>
985 readBanner(std::istream& in,
987 const bool tolerant =
false,
989 const bool isGraph =
false) {
994 using Teuchos::MatrixMarket::Banner;
995 typedef Teuchos::ScalarTraits<scalar_type> STS;
1001 const bool readFailed = !getline(in, line);
1002 TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1003 "Failed to get Matrix Market banner line from input.");
1010 pBanner = rcp(
new Banner(line, tolerant));
1011 }
catch (std::exception& e) {
1012 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1013 "Matrix Market banner line contains syntax error(s): "
1017 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() !=
"matrix",
1018 std::invalid_argument,
1019 "The Matrix Market file does not contain "
1020 "matrix data. Its Banner (first) line says that its object type is \""
1021 << pBanner->matrixType() <<
"\", rather than the required \"matrix\".");
1025 TEUCHOS_TEST_FOR_EXCEPTION(
1026 !STS::isComplex && pBanner->dataType() ==
"complex",
1027 std::invalid_argument,
1028 "The Matrix Market file contains complex-valued data, but you are "
1029 "trying to read it into a matrix containing entries of the real-"
1030 "valued Scalar type \""
1031 << Teuchos::TypeNameTraits<scalar_type>::name() <<
"\".");
1032 TEUCHOS_TEST_FOR_EXCEPTION(
1034 pBanner->dataType() !=
"real" &&
1035 pBanner->dataType() !=
"complex" &&
1036 pBanner->dataType() !=
"integer",
1037 std::invalid_argument,
1038 "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
1039 "Matrix Market file may not contain a \"pattern\" matrix. A "
1040 "pattern matrix is really just a graph with no weights. It "
1041 "should be stored in a CrsGraph, not a CrsMatrix.");
1043 TEUCHOS_TEST_FOR_EXCEPTION(
1045 pBanner->dataType() !=
"pattern",
1046 std::invalid_argument,
1047 "When reading Matrix Market data into a Tpetra::CrsGraph, the "
1048 "Matrix Market file must contain a \"pattern\" matrix.");
1075 static Teuchos::Tuple<global_ordinal_type, 3>
1076 readCoordDims(std::istream& in,
1078 const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1080 const bool tolerant =
false,
1081 const bool =
false) {
1082 using Teuchos::Tuple;
1083 using Teuchos::MatrixMarket::readCoordinateDimensions;
1088 Tuple<global_ordinal_type, 3> dims;
1094 bool success =
false;
1095 if (pComm->getRank() == 0) {
1096 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() !=
"coordinate",
1097 std::invalid_argument,
1098 "The Tpetra::CrsMatrix Matrix Market reader "
1099 "only accepts \"coordinate\" (sparse) matrix data.");
1103 success = readCoordinateDimensions(in, numRows, numCols,
1104 numNonzeros, lineNumber,
1110 dims[2] = numNonzeros;
1118 int the_success = success ? 1 : 0;
1119 Teuchos::broadcast(*pComm, 0, &the_success);
1120 success = (the_success == 1);
1125 Teuchos::broadcast(*pComm, 0, dims);
1132 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1133 "Error reading Matrix Market sparse matrix: failed to read "
1134 "coordinate matrix dimensions.");
1149 typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type>> adder_type;
1151 typedef Teuchos::MatrixMarket::SymmetrizingGraphAdder<Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type>> graph_adder_type;
1178 static Teuchos::RCP<adder_type>
1179 makeAdder(
const Teuchos::RCP<
const Teuchos::Comm<int>>& pComm,
1180 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1181 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1182 const bool tolerant =
false,
1183 const bool debug =
false) {
1184 if (pComm->getRank() == 0) {
1185 typedef Teuchos::MatrixMarket::Raw::Adder<
scalar_type,
1188 Teuchos::RCP<raw_adder_type> pRaw =
1189 Teuchos::rcp(
new raw_adder_type(dims[0], dims[1], dims[2],
1191 return Teuchos::rcp(
new adder_type(pRaw, pBanner->symmType()));
1193 return Teuchos::null;
1222 static Teuchos::RCP<graph_adder_type>
1223 makeGraphAdder(
const Teuchos::RCP<
const Teuchos::Comm<int>>& pComm,
1224 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1225 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1226 const bool tolerant =
false,
1227 const bool debug =
false) {
1228 if (pComm->getRank() == 0) {
1229 typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
1230 Teuchos::RCP<raw_adder_type> pRaw =
1231 Teuchos::rcp(
new raw_adder_type(dims[0], dims[1], dims[2],
1233 return Teuchos::rcp(
new graph_adder_type(pRaw, pBanner->symmType()));
1235 return Teuchos::null;
1240 static Teuchos::RCP<sparse_graph_type>
1241 readSparseGraphHelper(std::istream& in,
1242 const Teuchos::RCP<
const Teuchos::Comm<int>>& pComm,
1243 const Teuchos::RCP<const map_type>& rowMap,
1244 Teuchos::RCP<const map_type>& colMap,
1245 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1246 const bool tolerant,
1252 using Teuchos::Tuple;
1253 using Teuchos::MatrixMarket::Banner;
1255 const int myRank = pComm->getRank();
1256 const int rootRank = 0;
1261 size_t lineNumber = 1;
1263 if (debug && myRank == rootRank) {
1264 cerr <<
"Matrix Market reader: readGraph:" << endl
1265 <<
"-- Reading banner line" << endl;
1274 RCP<const Banner> pBanner;
1280 int bannerIsCorrect = 1;
1281 std::ostringstream errMsg;
1283 if (myRank == rootRank) {
1286 pBanner = readBanner(in, lineNumber, tolerant, debug,
true);
1287 }
catch (std::exception& e) {
1288 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
1289 "threw an exception: "
1291 bannerIsCorrect = 0;
1294 if (bannerIsCorrect) {
1299 if (!tolerant && pBanner->matrixType() !=
"coordinate") {
1300 bannerIsCorrect = 0;
1301 errMsg <<
"The Matrix Market input file must contain a "
1302 "\"coordinate\"-format sparse graph in order to create a "
1303 "Tpetra::CrsGraph object from it, but the file's matrix "
1305 << pBanner->matrixType() <<
"\" instead.";
1310 if (tolerant && pBanner->matrixType() ==
"array") {
1311 bannerIsCorrect = 0;
1312 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
1313 "format sparse graph in order to create a Tpetra::CrsGraph "
1314 "object from it, but the file's matrix type is \"array\" "
1315 "instead. That probably means the file contains dense matrix "
1322 broadcast(*pComm, rootRank, ptr(&bannerIsCorrect));
1329 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1330 std::invalid_argument, errMsg.str());
1332 if (debug && myRank == rootRank) {
1333 cerr <<
"-- Reading dimensions line" << endl;
1341 Tuple<global_ordinal_type, 3> dims =
1342 readCoordDims(in, lineNumber, pBanner, pComm, tolerant, debug);
1344 if (debug && myRank == rootRank) {
1345 cerr <<
"-- Making Adder for collecting graph data" << endl;
1352 RCP<graph_adder_type> pAdder =
1353 makeGraphAdder(pComm, pBanner, dims, tolerant, debug);
1355 if (debug && myRank == rootRank) {
1356 cerr <<
"-- Reading graph data" << endl;
1366 int readSuccess = 1;
1367 std::ostringstream errMsg;
1368 if (myRank == rootRank) {
1371 typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
1374 reader_type reader(pAdder);
1377 std::pair<bool, std::vector<size_t>> results =
1378 reader.read(in, lineNumber, tolerant, debug);
1379 readSuccess = results.first ? 1 : 0;
1380 }
catch (std::exception& e) {
1385 broadcast(*pComm, rootRank, ptr(&readSuccess));
1394 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1395 "Failed to read the Matrix Market sparse graph file: "
1399 if (debug && myRank == rootRank) {
1400 cerr <<
"-- Successfully read the Matrix Market data" << endl;
1413 if (debug && myRank == rootRank) {
1414 cerr <<
"-- Tolerant mode: rebroadcasting graph dimensions"
1416 <<
"----- Dimensions before: "
1417 << dims[0] <<
" x " << dims[1]
1421 Tuple<global_ordinal_type, 2> updatedDims;
1422 if (myRank == rootRank) {
1429 std::max(dims[0], pAdder->getAdder()->numRows());
1430 updatedDims[1] = pAdder->getAdder()->numCols();
1432 broadcast(*pComm, rootRank, updatedDims);
1433 dims[0] = updatedDims[0];
1434 dims[1] = updatedDims[1];
1435 if (debug && myRank == rootRank) {
1436 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
1449 if (myRank == rootRank) {
1456 if (dims[0] < pAdder->getAdder()->numRows()) {
1460 broadcast(*pComm, 0, ptr(&dimsMatch));
1461 if (dimsMatch == 0) {
1468 Tuple<global_ordinal_type, 2> addersDims;
1469 if (myRank == rootRank) {
1470 addersDims[0] = pAdder->getAdder()->numRows();
1471 addersDims[1] = pAdder->getAdder()->numCols();
1473 broadcast(*pComm, 0, addersDims);
1474 TEUCHOS_TEST_FOR_EXCEPTION(
1475 dimsMatch == 0, std::runtime_error,
1476 "The graph metadata says that the graph is " << dims[0] <<
" x "
1477 << dims[1] <<
", but the actual data says that the graph is "
1478 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
1479 "data includes more rows than reported in the metadata. This "
1480 "is not allowed when parsing in strict mode. Parse the graph in "
1481 "tolerant mode to ignore the metadata when it disagrees with the "
1487 RCP<map_type> proc0Map;
1489 if (Teuchos::is_null(rowMap)) {
1492 indexBase = rowMap->getIndexBase();
1494 if (myRank == rootRank) {
1495 proc0Map = rcp(
new map_type(dims[0], dims[0], indexBase, pComm));
1497 proc0Map = rcp(
new map_type(dims[0], 0, indexBase, pComm));
1501 std::map<global_ordinal_type, size_t> numEntriesPerRow_map;
1502 if (myRank == rootRank) {
1503 const auto& entries = pAdder()->getAdder()->getEntries();
1508 for (
const auto& entry : entries) {
1510 ++numEntriesPerRow_map[gblRow];
1514 Teuchos::Array<size_t> numEntriesPerRow(proc0Map->getLocalNumElements());
1515 for (
const auto& ent : numEntriesPerRow_map) {
1517 numEntriesPerRow[lclRow] = ent.second;
1524 std::map<global_ordinal_type, size_t> empty_map;
1525 std::swap(numEntriesPerRow_map, empty_map);
1528 RCP<sparse_graph_type> proc0Graph =
1530 constructorParams));
1531 if (myRank == rootRank) {
1532 typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
1535 const std::vector<element_type>& entries =
1536 pAdder->getAdder()->getEntries();
1539 for (
size_t curPos = 0; curPos < entries.size(); curPos++) {
1540 const element_type& curEntry = entries[curPos];
1543 Teuchos::ArrayView<const global_ordinal_type> colView(&curCol, 1);
1544 proc0Graph->insertGlobalIndices(curRow, colView);
1547 proc0Graph->fillComplete();
1549 RCP<sparse_graph_type> distGraph;
1550 if (Teuchos::is_null(rowMap)) {
1552 RCP<map_type> distMap =
1553 rcp(
new map_type(dims[0], 0, pComm, GloballyDistributed));
1559 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1560 import_type importer(proc0Map, distMap);
1563 distGraph->doImport(*proc0Graph, importer,
INSERT);
1568 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1569 import_type importer(proc0Map, rowMap);
1572 distGraph->doImport(*proc0Graph, importer,
INSERT);
1602 static Teuchos::RCP<sparse_graph_type>
1607 const bool debug =
false) {
1646 static Teuchos::RCP<sparse_graph_type>
1648 const Teuchos::RCP<
const Teuchos::Comm<int>>&
pComm,
1652 const bool debug =
false) {
1700 static Teuchos::RCP<sparse_graph_type>
1702 const Teuchos::RCP<const map_type>& rowMap,
1703 Teuchos::RCP<const map_type>& colMap,
1704 const Teuchos::RCP<const map_type>& domainMap,
1705 const Teuchos::RCP<const map_type>&
rangeMap,
1708 const bool debug =
false) {
1710 "Input rowMap must be nonnull.");
1712 if (comm.is_null()) {
1715 return Teuchos::null;
1749 static Teuchos::RCP<sparse_graph_type>
1751 const Teuchos::RCP<
const Teuchos::Comm<int>>&
pComm,
1754 const bool debug =
false) {
1759 Teuchos::RCP<sparse_graph_type>
graph =
1760 readSparseGraphHelper(
in,
pComm,
1764 graph->fillComplete();
1798 static Teuchos::RCP<sparse_graph_type>
1800 const Teuchos::RCP<
const Teuchos::Comm<int>>&
pComm,
1804 const bool debug =
false) {
1807 Teuchos::RCP<sparse_graph_type>
graph =
1808 readSparseGraphHelper(
in,
pComm,
1855 static Teuchos::RCP<sparse_graph_type>
1857 const Teuchos::RCP<const map_type>& rowMap,
1858 Teuchos::RCP<const map_type>& colMap,
1859 const Teuchos::RCP<const map_type>& domainMap,
1860 const Teuchos::RCP<const map_type>&
rangeMap,
1863 const bool debug =
false) {
1864 Teuchos::RCP<sparse_graph_type>
graph =
1865 readSparseGraphHelper(
in, rowMap->getComm(),
1866 rowMap, colMap, Teuchos::null,
tolerant,
1874#include "MatrixMarket_TpetraNew.hpp"
1899 static Teuchos::RCP<sparse_matrix_type>
1904 const bool debug =
false) {
1941 static Teuchos::RCP<sparse_matrix_type>
1947 const bool debug =
false) {
1991 static Teuchos::RCP<sparse_matrix_type>
1993 const Teuchos::RCP<const map_type>& rowMap,
1994 Teuchos::RCP<const map_type>& colMap,
1995 const Teuchos::RCP<const map_type>& domainMap,
1996 const Teuchos::RCP<const map_type>&
rangeMap,
1999 const bool debug =
false) {
2001 rowMap.is_null(), std::invalid_argument,
2002 "Row Map must be nonnull.");
2037 static Teuchos::RCP<sparse_matrix_type>
2039 const Teuchos::RCP<
const Teuchos::Comm<int>>&
pComm,
2042 const bool debug =
false) {
2045 using Teuchos::arcp;
2046 using Teuchos::ArrayRCP;
2047 using Teuchos::broadcast;
2048 using Teuchos::null;
2051 using Teuchos::REDUCE_MAX;
2052 using Teuchos::reduceAll;
2053 using Teuchos::Tuple;
2054 using Teuchos::MatrixMarket::Banner;
2055 typedef Teuchos::ScalarTraits<scalar_type> STS;
2067 cerr <<
"Matrix Market reader: readSparse:" <<
endl
2068 <<
"-- Reading banner line" <<
endl;
2084 std::ostringstream
errMsg;
2090 }
catch (std::exception&
e) {
2091 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2092 "threw an exception: "
2104 errMsg <<
"The Matrix Market input file must contain a "
2105 "\"coordinate\"-format sparse matrix in order to create a "
2106 "Tpetra::CrsMatrix object from it, but the file's matrix "
2108 <<
pBanner->matrixType() <<
"\" instead.";
2115 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2116 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2117 "object from it, but the file's matrix type is \"array\" "
2118 "instead. That probably means the file contains dense matrix "
2133 std::invalid_argument,
errMsg.str());
2136 cerr <<
"-- Reading dimensions line" <<
endl;
2148 cerr <<
"-- Making Adder for collecting matrix data" <<
endl;
2157 cerr <<
"-- Reading matrix data" <<
endl;
2168 std::ostringstream
errMsg;
2172 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2178 std::pair<bool, std::vector<size_t>>
results =
2181 }
catch (std::exception&
e) {
2196 "Failed to read the Matrix Market sparse matrix file: "
2201 cerr <<
"-- Successfully read the Matrix Market data" <<
endl;
2215 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2217 <<
"----- Dimensions before: "
2230 std::max(
dims[0],
pAdder->getAdder()->numRows());
2237 cerr <<
"----- Dimensions after: " <<
dims[0] <<
" x "
2257 if (
dims[0] <
pAdder->getAdder()->numRows()) {
2277 "The matrix metadata says that the matrix is " <<
dims[0] <<
" x "
2278 <<
dims[1] <<
", but the actual data says that the matrix is "
2280 "data includes more rows than reported in the metadata. This "
2281 "is not allowed when parsing in strict mode. Parse the matrix in "
2282 "tolerant mode to ignore the metadata when it disagrees with the "
2288 cerr <<
"-- Converting matrix data into CSR format on Proc 0" <<
endl;
2309 std::ostringstream
errMsg;
2313 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2327 pAdder->getAdder()->merge();
2330 const std::vector<element_type>& entries =
2331 pAdder->getAdder()->getEntries();
2334 const size_t numEntries = (
size_t)entries.size();
2337 cerr <<
"----- Proc 0: Matrix has numRows=" <<
numRows
2338 <<
" rows and numEntries=" << numEntries
2339 <<
" entries." <<
endl;
2362 "Row indices are out of order, even though they are supposed "
2363 "to be sorted. curRow = "
2364 <<
curRow <<
", prvRow = "
2365 <<
prvRow <<
", at curPos = " <<
curPos <<
". Please report "
2366 "this bug to the Tpetra developers.");
2382 catch (std::exception&
e) {
2384 errMsg <<
"Failed to merge sparse matrix entries and convert to "
2394 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
2397 cerr <<
"(only showing first and last few entries) ";
2402 for (size_type
k = 0;
k < 2; ++
k) {
2410 for (size_type
k = 0;
k <
numRows - 1; ++
k) {
2418 cerr <<
"----- Proc 0: rowPtr ";
2420 cerr <<
"(only showing first and last few entries) ";
2424 for (size_type
k = 0;
k < 2; ++
k) {
2448 cerr <<
"-- Making range, domain, and row maps" <<
endl;
2461 cerr <<
"-- Distributing the matrix data" <<
endl;
2484 cerr <<
"-- Inserting matrix entries on each processor";
2486 cerr <<
" and calling fillComplete()";
2508 "Reader::makeMatrix() returned a null pointer on at least one "
2509 "process. Please report this bug to the Tpetra developers.");
2512 "Reader::makeMatrix() returned a null pointer. "
2513 "Please report this bug to the Tpetra developers.");
2534 cerr <<
"-- Matrix is "
2536 <<
" with " <<
pMatrix->getGlobalNumEntries()
2537 <<
" entries, and index base "
2543 cerr <<
"-- Proc " <<
p <<
" owns "
2544 <<
pMatrix->getLocalNumCols() <<
" columns, and "
2545 <<
pMatrix->getLocalNumEntries() <<
" entries." <<
endl;
2553 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
2587 static Teuchos::RCP<sparse_matrix_type>
2589 const Teuchos::RCP<
const Teuchos::Comm<int>>&
pComm,
2593 const bool debug =
false) {
2596 using Teuchos::arcp;
2597 using Teuchos::ArrayRCP;
2598 using Teuchos::broadcast;
2599 using Teuchos::null;
2602 using Teuchos::reduceAll;
2603 using Teuchos::Tuple;
2604 using Teuchos::MatrixMarket::Banner;
2605 typedef Teuchos::ScalarTraits<scalar_type> STS;
2617 cerr <<
"Matrix Market reader: readSparse:" <<
endl
2618 <<
"-- Reading banner line" <<
endl;
2634 std::ostringstream
errMsg;
2640 }
catch (std::exception&
e) {
2641 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2642 "threw an exception: "
2654 errMsg <<
"The Matrix Market input file must contain a "
2655 "\"coordinate\"-format sparse matrix in order to create a "
2656 "Tpetra::CrsMatrix object from it, but the file's matrix "
2658 <<
pBanner->matrixType() <<
"\" instead.";
2665 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2666 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2667 "object from it, but the file's matrix type is \"array\" "
2668 "instead. That probably means the file contains dense matrix "
2683 std::invalid_argument,
errMsg.str());
2686 cerr <<
"-- Reading dimensions line" <<
endl;
2698 cerr <<
"-- Making Adder for collecting matrix data" <<
endl;
2707 cerr <<
"-- Reading matrix data" <<
endl;
2718 std::ostringstream
errMsg;
2722 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2728 std::pair<bool, std::vector<size_t>>
results =
2731 }
catch (std::exception&
e) {
2746 "Failed to read the Matrix Market sparse matrix file: "
2751 cerr <<
"-- Successfully read the Matrix Market data" <<
endl;
2765 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2767 <<
"----- Dimensions before: "
2780 std::max(
dims[0],
pAdder->getAdder()->numRows());
2787 cerr <<
"----- Dimensions after: " <<
dims[0] <<
" x "
2807 if (
dims[0] <
pAdder->getAdder()->numRows()) {
2827 "The matrix metadata says that the matrix is " <<
dims[0] <<
" x "
2828 <<
dims[1] <<
", but the actual data says that the matrix is "
2830 "data includes more rows than reported in the metadata. This "
2831 "is not allowed when parsing in strict mode. Parse the matrix in "
2832 "tolerant mode to ignore the metadata when it disagrees with the "
2838 cerr <<
"-- Converting matrix data into CSR format on Proc 0" <<
endl;
2859 std::ostringstream
errMsg;
2863 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2877 pAdder->getAdder()->merge();
2880 const std::vector<element_type>& entries =
2881 pAdder->getAdder()->getEntries();
2884 const size_t numEntries = (
size_t)entries.size();
2887 cerr <<
"----- Proc 0: Matrix has numRows=" <<
numRows
2888 <<
" rows and numEntries=" << numEntries
2889 <<
" entries." <<
endl;
2912 "Row indices are out of order, even though they are supposed "
2913 "to be sorted. curRow = "
2914 <<
curRow <<
", prvRow = "
2915 <<
prvRow <<
", at curPos = " <<
curPos <<
". Please report "
2916 "this bug to the Tpetra developers.");
2932 catch (std::exception&
e) {
2934 errMsg <<
"Failed to merge sparse matrix entries and convert to "
2944 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
2947 cerr <<
"(only showing first and last few entries) ";
2952 for (size_type
k = 0;
k < 2; ++
k) {
2960 for (size_type
k = 0;
k <
numRows - 1; ++
k) {
2968 cerr <<
"----- Proc 0: rowPtr ";
2970 cerr <<
"(only showing first and last few entries) ";
2974 for (size_type
k = 0;
k < 2; ++
k) {
2998 cerr <<
"-- Making range, domain, and row maps" <<
endl;
3011 cerr <<
"-- Distributing the matrix data" <<
endl;
3034 cerr <<
"-- Inserting matrix entries on each process "
3035 "and calling fillComplete()"
3045 Teuchos::RCP<sparse_matrix_type>
pMatrix =
3057 "Reader::makeMatrix() returned a null pointer on at least one "
3058 "process. Please report this bug to the Tpetra developers.");
3061 "Reader::makeMatrix() returned a null pointer. "
3062 "Please report this bug to the Tpetra developers.");
3078 cerr <<
"-- Matrix is "
3080 <<
" with " <<
pMatrix->getGlobalNumEntries()
3081 <<
" entries, and index base "
3087 cerr <<
"-- Proc " <<
p <<
" owns "
3088 <<
pMatrix->getLocalNumCols() <<
" columns, and "
3089 <<
pMatrix->getLocalNumEntries() <<
" entries." <<
endl;
3096 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3142 static Teuchos::RCP<sparse_matrix_type>
3144 const Teuchos::RCP<const map_type>& rowMap,
3145 Teuchos::RCP<const map_type>& colMap,
3146 const Teuchos::RCP<const map_type>& domainMap,
3147 const Teuchos::RCP<const map_type>&
rangeMap,
3150 const bool debug =
false) {
3153 using Teuchos::arcp;
3154 using Teuchos::ArrayRCP;
3155 using Teuchos::ArrayView;
3157 using Teuchos::broadcast;
3158 using Teuchos::Comm;
3159 using Teuchos::null;
3162 using Teuchos::reduceAll;
3163 using Teuchos::Tuple;
3164 using Teuchos::MatrixMarket::Banner;
3165 typedef Teuchos::ScalarTraits<scalar_type> STS;
3176 rowMap.is_null(), std::invalid_argument,
3177 "Row Map must be nonnull.");
3179 rangeMap.is_null(), std::invalid_argument,
3180 "Range Map must be nonnull.");
3182 domainMap.is_null(), std::invalid_argument,
3183 "Domain Map must be nonnull.");
3185 rowMap->getComm().getRawPtr() !=
pComm.getRawPtr(),
3186 std::invalid_argument,
3187 "The specified row Map's communicator (rowMap->getComm())"
3188 "differs from the given separately supplied communicator pComm.");
3190 domainMap->getComm().getRawPtr() !=
pComm.getRawPtr(),
3191 std::invalid_argument,
3192 "The specified domain Map's communicator (domainMap->getComm())"
3193 "differs from the given separately supplied communicator pComm.");
3196 std::invalid_argument,
3197 "The specified range Map's communicator (rangeMap->getComm())"
3198 "differs from the given separately supplied communicator pComm.");
3206 cerr <<
"Matrix Market reader: readSparse:" <<
endl
3207 <<
"-- Reading banner line" <<
endl;
3223 std::ostringstream
errMsg;
3229 }
catch (std::exception&
e) {
3230 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
3231 "threw an exception: "
3243 errMsg <<
"The Matrix Market input file must contain a "
3244 "\"coordinate\"-format sparse matrix in order to create a "
3245 "Tpetra::CrsMatrix object from it, but the file's matrix "
3247 <<
pBanner->matrixType() <<
"\" instead.";
3254 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
3255 "format sparse matrix in order to create a Tpetra::CrsMatrix "
3256 "object from it, but the file's matrix type is \"array\" "
3257 "instead. That probably means the file contains dense matrix "
3272 std::invalid_argument,
errMsg.str());
3275 cerr <<
"-- Reading dimensions line" <<
endl;
3287 cerr <<
"-- Making Adder for collecting matrix data" <<
endl;
3298 cerr <<
"-- Reading matrix data" <<
endl;
3309 std::ostringstream
errMsg;
3313 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3319 std::pair<bool, std::vector<size_t>>
results =
3322 }
catch (std::exception&
e) {
3337 "Failed to read the Matrix Market sparse matrix file: "
3342 cerr <<
"-- Successfully read the Matrix Market data" <<
endl;
3356 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
3358 <<
"----- Dimensions before: "
3371 std::max(
dims[0],
pAdder->getAdder()->numRows());
3378 cerr <<
"----- Dimensions after: " <<
dims[0] <<
" x "
3398 if (
dims[0] <
pAdder->getAdder()->numRows()) {
3418 "The matrix metadata says that the matrix is " <<
dims[0] <<
" x "
3419 <<
dims[1] <<
", but the actual data says that the matrix is "
3421 "data includes more rows than reported in the metadata. This "
3422 "is not allowed when parsing in strict mode. Parse the matrix in "
3423 "tolerant mode to ignore the metadata when it disagrees with the "
3429 cerr <<
"-- Converting matrix data into CSR format on Proc 0" <<
endl;
3450 std::ostringstream
errMsg;
3454 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3468 pAdder->getAdder()->merge();
3471 const std::vector<element_type>& entries =
3472 pAdder->getAdder()->getEntries();
3475 const size_t numEntries = (
size_t)entries.size();
3478 cerr <<
"----- Proc 0: Matrix has numRows=" <<
numRows
3479 <<
" rows and numEntries=" << numEntries
3480 <<
" entries." <<
endl;
3503 "Row indices are out of order, even though they are supposed "
3504 "to be sorted. curRow = "
3505 <<
curRow <<
", prvRow = "
3506 <<
prvRow <<
", at curPos = " <<
curPos <<
". Please report "
3507 "this bug to the Tpetra developers.");
3521 catch (std::exception&
e) {
3523 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3533 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3536 cerr <<
"(only showing first and last few entries) ";
3541 for (size_type
k = 0;
k < 2; ++
k) {
3549 for (size_type
k = 0;
k <
numRows - 1; ++
k) {
3557 cerr <<
"----- Proc 0: rowPtr ";
3559 cerr <<
"(only showing first and last few entries) ";
3563 for (size_type
k = 0;
k < 2; ++
k) {
3577 cerr <<
"----- Proc 0: colInd = [";
3596 cerr <<
"-- Verifying Maps" <<
endl;
3600 std::invalid_argument,
3601 "The range Map has " <<
rangeMap->getMaxAllGlobalIndex()
3602 <<
" max entry, but the matrix has a global number of rows " <<
dims[0]
3606 std::invalid_argument,
3607 "The domain Map has " << domainMap->getMaxAllGlobalIndex()
3608 <<
" max entry, but the matrix has a global number of columns "
3632 cerr <<
"-- Proc 0: Filling gather matrix" <<
endl;
3692 if (colMap.is_null()) {
3718 cerr <<
"-- Matrix is "
3720 <<
" with " <<
A->getGlobalNumEntries()
3721 <<
" entries, and index base "
3722 <<
A->getIndexBase() <<
"." <<
endl;
3727 cerr <<
"-- Proc " <<
p <<
" owns "
3728 <<
A->getLocalNumCols() <<
" columns, and "
3729 <<
A->getLocalNumEntries() <<
" entries." <<
endl;
3737 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3772 static Teuchos::RCP<multivector_type>
3775 Teuchos::RCP<const map_type>&
map,
3777 const bool debug =
false,
3778 const bool binary =
false) {
3797 std::ios_base::openmode
mode = std::ios_base::in) {
3805 if (comm->getRank() == 0) {
3816 "Could not open input file '" +
filename +
"' on root rank 0.");
3850 static Teuchos::RCP<vector_type>
3853 Teuchos::RCP<const map_type>&
map,
3855 const bool debug =
false) {
3929 static Teuchos::RCP<multivector_type>
3932 Teuchos::RCP<const map_type>&
map,
3934 const bool debug =
false,
3935 const bool binary =
false) {
3936 Teuchos::RCP<Teuchos::FancyOStream>
err =
3937 Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
3942 static Teuchos::RCP<vector_type>
3945 Teuchos::RCP<const map_type>&
map,
3947 const bool debug =
false) {
3948 Teuchos::RCP<Teuchos::FancyOStream>
err =
3949 Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
3974 static Teuchos::RCP<const map_type>
3978 const bool debug =
false,
3979 const bool binary =
false) {
3986 template <
class MultiVectorScalarType>
3991 readDenseImpl(std::istream&
in,
3993 Teuchos::RCP<const map_type>&
map,
3994 const Teuchos::RCP<Teuchos::FancyOStream>&
err,
3996 const bool debug =
false,
3997 const bool binary =
false) {
3999 using Teuchos::ArrayRCP;
4001 using Teuchos::broadcast;
4002 using Teuchos::outArg;
4004 using Teuchos::Tuple;
4005 using Teuchos::MatrixMarket::Banner;
4006 using Teuchos::MatrixMarket::checkCommentLine;
4011 typedef Teuchos::ScalarTraits<ST> STS;
4012 typedef typename STS::magnitudeType
MT;
4013 typedef Teuchos::ScalarTraits<MT> STM;
4018 const int myRank = comm->getRank();
4020 if (!
err.is_null()) {
4026 if (!err.is_null()) {
4061 size_t lineNumber = 1;
4064 std::ostringstream exMsg;
4065 int localBannerReadSuccess = 1;
4066 int localDimsReadSuccess = 1;
4072 *err << myRank <<
": readDenseImpl: Reading banner line (dense)" << endl;
4078 RCP<const Banner> pBanner;
4080 pBanner = readBanner(in, lineNumber, tolerant, debug);
4081 }
catch (std::exception& e) {
4083 localBannerReadSuccess = 0;
4086 if (localBannerReadSuccess) {
4087 if (pBanner->matrixType() !=
"array") {
4088 exMsg <<
"The Matrix Market file does not contain dense matrix "
4089 "data. Its banner (first) line says that its matrix type is \""
4090 << pBanner->matrixType() <<
"\", rather that the required "
4092 localBannerReadSuccess = 0;
4093 }
else if (pBanner->dataType() ==
"pattern") {
4094 exMsg <<
"The Matrix Market file's banner (first) "
4095 "line claims that the matrix's data type is \"pattern\". This does "
4096 "not make sense for a dense matrix, yet the file reports the matrix "
4097 "as dense. The only valid data types for a dense matrix are "
4098 "\"real\", \"complex\", and \"integer\".";
4099 localBannerReadSuccess = 0;
4103 dims[2] = encodeDataType(pBanner->dataType());
4109 if (localBannerReadSuccess) {
4111 *err << myRank <<
": readDenseImpl: Reading dimensions line (dense)" << endl;
4120 bool commentLine =
true;
4122 while (commentLine) {
4125 if (in.eof() || in.fail()) {
4126 exMsg <<
"Unable to get array dimensions line (at all) from line "
4127 << lineNumber <<
" of input stream. The input stream "
4128 <<
"claims that it is "
4129 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4130 localDimsReadSuccess = 0;
4133 if (getline(in, line)) {
4140 size_t start = 0, size = 0;
4141 commentLine = checkCommentLine(line, start, size, lineNumber, tolerant);
4148 std::istringstream istr(line);
4152 if (istr.eof() || istr.fail()) {
4153 exMsg <<
"Unable to read any data from line " << lineNumber
4154 <<
" of input; the line should contain the matrix dimensions "
4155 <<
"\"<numRows> <numCols>\".";
4156 localDimsReadSuccess = 0;
4161 exMsg <<
"Failed to get number of rows from line "
4162 << lineNumber <<
" of input; the line should contains the "
4163 <<
"matrix dimensions \"<numRows> <numCols>\".";
4164 localDimsReadSuccess = 0;
4166 dims[0] = theNumRows;
4168 exMsg <<
"No more data after number of rows on line "
4169 << lineNumber <<
" of input; the line should contain the "
4170 <<
"matrix dimensions \"<numRows> <numCols>\".";
4171 localDimsReadSuccess = 0;
4176 exMsg <<
"Failed to get number of columns from line "
4177 << lineNumber <<
" of input; the line should contain "
4178 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4179 localDimsReadSuccess = 0;
4181 dims[1] = theNumCols;
4190 in.read(
reinterpret_cast<char*
>(&numRows),
sizeof(numRows));
4191 in.read(
reinterpret_cast<char*
>(&numCols),
sizeof(numCols));
4192 dims[0] = Teuchos::as<GO>(numRows);
4193 dims[1] = Teuchos::as<GO>(numCols);
4194 if ((
typeid(ST) ==
typeid(
double)) || Teuchos::ScalarTraits<ST>::isOrdinal) {
4196 }
else if (Teuchos::ScalarTraits<ST>::isComplex) {
4199 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
4200 "Unrecognized Matrix Market data type. "
4201 "We should never get here. "
4202 "Please report this bug to the Tpetra developers.");
4204 localDimsReadSuccess =
true;
4205 localDimsReadSuccess =
true;
4212 Tuple<GO, 5> bannerDimsReadResult;
4214 bannerDimsReadResult[0] = dims[0];
4215 bannerDimsReadResult[1] = dims[1];
4216 bannerDimsReadResult[2] = dims[2];
4217 bannerDimsReadResult[3] = localBannerReadSuccess;
4218 bannerDimsReadResult[4] = localDimsReadSuccess;
4222 broadcast(*comm, 0, bannerDimsReadResult);
4224 TEUCHOS_TEST_FOR_EXCEPTION(
4225 bannerDimsReadResult[3] == 0, std::runtime_error,
4226 "Failed to read banner line: " << exMsg.str());
4227 TEUCHOS_TEST_FOR_EXCEPTION(
4228 bannerDimsReadResult[4] == 0, std::runtime_error,
4229 "Failed to read matrix dimensions line: " << exMsg.str());
4231 dims[0] = bannerDimsReadResult[0];
4232 dims[1] = bannerDimsReadResult[1];
4233 dims[2] = bannerDimsReadResult[2];
4238 const size_t numCols =
static_cast<size_t>(dims[1]);
4243 RCP<const map_type> proc0Map;
4244 if (map.is_null()) {
4248 map = createUniformContigMapWithNode<LO, GO, NT>(numRows, comm);
4249 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4251 proc0Map = createContigMapWithNode<LO, GO, NT>(numRows, localNumRows,
4254 proc0Map = Details::computeGatherMap<map_type>(map, err, debug);
4258 RCP<MV> X = createMultiVector<ST, LO, GO, NT>(proc0Map, numCols);
4264 int localReadDataSuccess = 1;
4269 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4274 TEUCHOS_TEST_FOR_EXCEPTION(
4275 !X->isConstantStride(), std::logic_error,
4276 "Can't get a 1-D view of the entries of the MultiVector X on "
4277 "Process 0, because the stride between the columns of X is not "
4278 "constant. This shouldn't happen because we just created X and "
4279 "haven't filled it in yet. Please report this bug to the Tpetra "
4286 ArrayRCP<ST> X_view = X->get1dViewNonConst();
4287 TEUCHOS_TEST_FOR_EXCEPTION(
4288 as<global_size_t>(X_view.size()) < numRows * numCols,
4290 "The view of X has size " << X_view.size() <<
" which is not enough to "
4291 "accommodate the expected number of entries numRows*numCols = "
4292 << numRows <<
"*" << numCols <<
" = " << numRows * numCols <<
". "
4293 "Please report this bug to the Tpetra developers.");
4294 const size_t stride = X->getStride();
4301 const bool isComplex = (dims[2] == 1);
4302 size_type count = 0, curRow = 0, curCol = 0;
4305 while (getline(in, line)) {
4309 size_t start = 0, size = 0;
4310 const bool commentLine =
4311 checkCommentLine(line, start, size, lineNumber, tolerant);
4318 if (count >= X_view.size()) {
4322 TEUCHOS_TEST_FOR_EXCEPTION(
4323 count >= X_view.size(),
4325 "The Matrix Market input stream has more data in it than "
4326 "its metadata reported. Current line number is "
4327 << lineNumber <<
".");
4336 const size_t pos = line.substr(start, size).find(
':');
4337 if (pos != std::string::npos) {
4341 std::istringstream istr(line.substr(start, size));
4345 if (istr.eof() || istr.fail()) {
4352 TEUCHOS_TEST_FOR_EXCEPTION(
4353 !tolerant && (istr.eof() || istr.fail()),
4355 "Line " << lineNumber <<
" of the Matrix Market file is "
4356 "empty, or we cannot read from it for some other reason.");
4359 ST val = STS::zero();
4362 MT real = STM::zero(), imag = STM::zero();
4379 TEUCHOS_TEST_FOR_EXCEPTION(
4380 !tolerant && istr.eof(), std::runtime_error,
4381 "Failed to get the real part of a complex-valued matrix "
4383 << lineNumber <<
" of the Matrix Market "
4389 }
else if (istr.eof()) {
4390 TEUCHOS_TEST_FOR_EXCEPTION(
4391 !tolerant && istr.eof(), std::runtime_error,
4392 "Missing imaginary part of a complex-valued matrix entry "
4394 << lineNumber <<
" of the Matrix Market file.");
4400 TEUCHOS_TEST_FOR_EXCEPTION(
4401 !tolerant && istr.fail(), std::runtime_error,
4402 "Failed to get the imaginary part of a complex-valued "
4403 "matrix entry from line "
4404 << lineNumber <<
" of the "
4405 "Matrix Market file.");
4412 TEUCHOS_TEST_FOR_EXCEPTION(
4413 !tolerant && istr.fail(), std::runtime_error,
4414 "Failed to get a real-valued matrix entry from line "
4415 << lineNumber <<
" of the Matrix Market file.");
4418 if (istr.fail() && tolerant) {
4424 TEUCHOS_TEST_FOR_EXCEPTION(
4425 !tolerant && istr.fail(), std::runtime_error,
4426 "Failed to read matrix data from line " << lineNumber
4427 <<
" of the Matrix Market file.");
4430 Teuchos::MatrixMarket::details::assignScalar<ST>(val, real, imag);
4432 curRow = count % numRows;
4433 curCol = count / numRows;
4434 X_view[curRow + curCol * stride] = val;
4439 TEUCHOS_TEST_FOR_EXCEPTION(
4440 !tolerant &&
static_cast<global_size_t>(count) < numRows * numCols,
4442 "The Matrix Market metadata reports that the dense matrix is "
4443 << numRows <<
" x " << numCols <<
", and thus has "
4444 << numRows * numCols <<
" total entries, but we only found " << count
4445 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
4446 }
catch (std::exception& e) {
4448 localReadDataSuccess = 0;
4453 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4458 TEUCHOS_TEST_FOR_EXCEPTION(
4459 !X->isConstantStride(), std::logic_error,
4460 "Can't get a 1-D view of the entries of the MultiVector X on "
4461 "Process 0, because the stride between the columns of X is not "
4462 "constant. This shouldn't happen because we just created X and "
4463 "haven't filled it in yet. Please report this bug to the Tpetra "
4470 auto X_view = X->getLocalViewHost(Access::OverwriteAll);
4472 TEUCHOS_TEST_FOR_EXCEPTION(
4473 as<global_size_t>(X_view.extent(0)) < numRows,
4475 "The view of X has " << X_view.extent(0) <<
" rows which is not enough to "
4476 "accommodate the expected number of entries numRows = "
4478 "Please report this bug to the Tpetra developers.");
4479 TEUCHOS_TEST_FOR_EXCEPTION(
4480 as<global_size_t>(X_view.extent(1)) < numCols,
4482 "The view of X has " << X_view.extent(1) <<
" colums which is not enough to "
4483 "accommodate the expected number of entries numRows = "
4485 "Please report this bug to the Tpetra developers.");
4487 for (
size_t curRow = 0; curRow < numRows; ++curRow) {
4488 for (
size_t curCol = 0; curCol < numCols; ++curCol) {
4489 if (Teuchos::ScalarTraits<ST>::isOrdinal) {
4491 in.read(
reinterpret_cast<char*
>(&val),
sizeof(val));
4492 X_view(curRow, curCol) = val;
4495 in.read(
reinterpret_cast<char*
>(&val),
sizeof(val));
4496 X_view(curRow, curCol) = val;
4504 *err << myRank <<
": readDenseImpl: done reading data" << endl;
4508 int globalReadDataSuccess = localReadDataSuccess;
4509 broadcast(*comm, 0, outArg(globalReadDataSuccess));
4510 TEUCHOS_TEST_FOR_EXCEPTION(
4511 globalReadDataSuccess == 0, std::runtime_error,
4512 "Failed to read the multivector's data: " << exMsg.str());
4517 if (comm->getSize() == 1 && map.is_null()) {
4519 if (!err.is_null()) {
4523 *err << myRank <<
": readDenseImpl: done" << endl;
4525 if (!err.is_null()) {
4532 *err << myRank <<
": readDenseImpl: Creating target MV" << endl;
4536 RCP<MV> Y = createMultiVector<ST, LO, GO, NT>(map, numCols);
4539 *err << myRank <<
": readDenseImpl: Creating Export" << endl;
4545 Export<LO, GO, NT> exporter(proc0Map, map, err);
4548 *err << myRank <<
": readDenseImpl: Exporting" << endl;
4551 Y->doExport(*X, exporter,
INSERT);
4553 if (!err.is_null()) {
4557 *err << myRank <<
": readDenseImpl: done" << endl;
4559 if (!err.is_null()) {
4567 template <
class VectorScalarType>
4572 readVectorImpl(std::istream& in,
4574 Teuchos::RCP<const map_type>& map,
4575 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4576 const bool tolerant =
false,
4577 const bool debug =
false) {
4580 using Teuchos::broadcast;
4581 using Teuchos::outArg;
4583 using Teuchos::Tuple;
4584 using Teuchos::MatrixMarket::Banner;
4585 using Teuchos::MatrixMarket::checkCommentLine;
4586 typedef VectorScalarType ST;
4590 typedef Teuchos::ScalarTraits<ST> STS;
4591 typedef typename STS::magnitudeType MT;
4592 typedef Teuchos::ScalarTraits<MT> STM;
4597 const int myRank = comm->getRank();
4599 if (!err.is_null()) {
4603 *err << myRank <<
": readVectorImpl" << endl;
4605 if (!err.is_null()) {
4639 size_t lineNumber = 1;
4642 std::ostringstream exMsg;
4643 int localBannerReadSuccess = 1;
4644 int localDimsReadSuccess = 1;
4649 *err << myRank <<
": readVectorImpl: Reading banner line (dense)" << endl;
4655 RCP<const Banner> pBanner;
4657 pBanner = readBanner(in, lineNumber, tolerant, debug);
4658 }
catch (std::exception& e) {
4660 localBannerReadSuccess = 0;
4663 if (localBannerReadSuccess) {
4664 if (pBanner->matrixType() !=
"array") {
4665 exMsg <<
"The Matrix Market file does not contain dense matrix "
4666 "data. Its banner (first) line says that its matrix type is \""
4667 << pBanner->matrixType() <<
"\", rather that the required "
4669 localBannerReadSuccess = 0;
4670 }
else if (pBanner->dataType() ==
"pattern") {
4671 exMsg <<
"The Matrix Market file's banner (first) "
4672 "line claims that the matrix's data type is \"pattern\". This does "
4673 "not make sense for a dense matrix, yet the file reports the matrix "
4674 "as dense. The only valid data types for a dense matrix are "
4675 "\"real\", \"complex\", and \"integer\".";
4676 localBannerReadSuccess = 0;
4680 dims[2] = encodeDataType(pBanner->dataType());
4686 if (localBannerReadSuccess) {
4688 *err << myRank <<
": readVectorImpl: Reading dimensions line (dense)" << endl;
4697 bool commentLine =
true;
4699 while (commentLine) {
4702 if (in.eof() || in.fail()) {
4703 exMsg <<
"Unable to get array dimensions line (at all) from line "
4704 << lineNumber <<
" of input stream. The input stream "
4705 <<
"claims that it is "
4706 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4707 localDimsReadSuccess = 0;
4710 if (getline(in, line)) {
4717 size_t start = 0, size = 0;
4718 commentLine = checkCommentLine(line, start, size, lineNumber, tolerant);
4725 std::istringstream istr(line);
4729 if (istr.eof() || istr.fail()) {
4730 exMsg <<
"Unable to read any data from line " << lineNumber
4731 <<
" of input; the line should contain the matrix dimensions "
4732 <<
"\"<numRows> <numCols>\".";
4733 localDimsReadSuccess = 0;
4738 exMsg <<
"Failed to get number of rows from line "
4739 << lineNumber <<
" of input; the line should contains the "
4740 <<
"matrix dimensions \"<numRows> <numCols>\".";
4741 localDimsReadSuccess = 0;
4743 dims[0] = theNumRows;
4745 exMsg <<
"No more data after number of rows on line "
4746 << lineNumber <<
" of input; the line should contain the "
4747 <<
"matrix dimensions \"<numRows> <numCols>\".";
4748 localDimsReadSuccess = 0;
4753 exMsg <<
"Failed to get number of columns from line "
4754 << lineNumber <<
" of input; the line should contain "
4755 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4756 localDimsReadSuccess = 0;
4758 dims[1] = theNumCols;
4768 exMsg <<
"File does not contain a 1D Vector";
4769 localDimsReadSuccess = 0;
4775 Tuple<GO, 5> bannerDimsReadResult;
4777 bannerDimsReadResult[0] = dims[0];
4778 bannerDimsReadResult[1] = dims[1];
4779 bannerDimsReadResult[2] = dims[2];
4780 bannerDimsReadResult[3] = localBannerReadSuccess;
4781 bannerDimsReadResult[4] = localDimsReadSuccess;
4786 broadcast(*comm, 0, bannerDimsReadResult);
4788 TEUCHOS_TEST_FOR_EXCEPTION(
4789 bannerDimsReadResult[3] == 0, std::runtime_error,
4790 "Failed to read banner line: " << exMsg.str());
4791 TEUCHOS_TEST_FOR_EXCEPTION(
4792 bannerDimsReadResult[4] == 0, std::runtime_error,
4793 "Failed to read matrix dimensions line: " << exMsg.str());
4795 dims[0] = bannerDimsReadResult[0];
4796 dims[1] = bannerDimsReadResult[1];
4797 dims[2] = bannerDimsReadResult[2];
4802 const size_t numCols =
static_cast<size_t>(dims[1]);
4807 RCP<const map_type> proc0Map;
4808 if (map.is_null()) {
4812 map = createUniformContigMapWithNode<LO, GO, NT>(numRows, comm);
4813 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4815 proc0Map = createContigMapWithNode<LO, GO, NT>(numRows, localNumRows,
4818 proc0Map = Details::computeGatherMap<map_type>(map, err, debug);
4822 RCP<MV> X = createVector<ST, LO, GO, NT>(proc0Map);
4828 int localReadDataSuccess = 1;
4832 *err << myRank <<
": readVectorImpl: Reading matrix data (dense)"
4837 TEUCHOS_TEST_FOR_EXCEPTION(
4838 !X->isConstantStride(), std::logic_error,
4839 "Can't get a 1-D view of the entries of the MultiVector X on "
4840 "Process 0, because the stride between the columns of X is not "
4841 "constant. This shouldn't happen because we just created X and "
4842 "haven't filled it in yet. Please report this bug to the Tpetra "
4849 Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst();
4850 TEUCHOS_TEST_FOR_EXCEPTION(
4851 as<global_size_t>(X_view.size()) < numRows * numCols,
4853 "The view of X has size " << X_view <<
" which is not enough to "
4854 "accommodate the expected number of entries numRows*numCols = "
4855 << numRows <<
"*" << numCols <<
" = " << numRows * numCols <<
". "
4856 "Please report this bug to the Tpetra developers.");
4857 const size_t stride = X->getStride();
4864 const bool isComplex = (dims[2] == 1);
4865 size_type count = 0, curRow = 0, curCol = 0;
4868 while (getline(in, line)) {
4872 size_t start = 0, size = 0;
4873 const bool commentLine =
4874 checkCommentLine(line, start, size, lineNumber, tolerant);
4881 if (count >= X_view.size()) {
4885 TEUCHOS_TEST_FOR_EXCEPTION(
4886 count >= X_view.size(),
4888 "The Matrix Market input stream has more data in it than "
4889 "its metadata reported. Current line number is "
4890 << lineNumber <<
".");
4899 const size_t pos = line.substr(start, size).find(
':');
4900 if (pos != std::string::npos) {
4904 std::istringstream istr(line.substr(start, size));
4908 if (istr.eof() || istr.fail()) {
4915 TEUCHOS_TEST_FOR_EXCEPTION(
4916 !tolerant && (istr.eof() || istr.fail()),
4918 "Line " << lineNumber <<
" of the Matrix Market file is "
4919 "empty, or we cannot read from it for some other reason.");
4922 ST val = STS::zero();
4925 MT real = STM::zero(), imag = STM::zero();
4942 TEUCHOS_TEST_FOR_EXCEPTION(
4943 !tolerant && istr.eof(), std::runtime_error,
4944 "Failed to get the real part of a complex-valued matrix "
4946 << lineNumber <<
" of the Matrix Market "
4952 }
else if (istr.eof()) {
4953 TEUCHOS_TEST_FOR_EXCEPTION(
4954 !tolerant && istr.eof(), std::runtime_error,
4955 "Missing imaginary part of a complex-valued matrix entry "
4957 << lineNumber <<
" of the Matrix Market file.");
4963 TEUCHOS_TEST_FOR_EXCEPTION(
4964 !tolerant && istr.fail(), std::runtime_error,
4965 "Failed to get the imaginary part of a complex-valued "
4966 "matrix entry from line "
4967 << lineNumber <<
" of the "
4968 "Matrix Market file.");
4975 TEUCHOS_TEST_FOR_EXCEPTION(
4976 !tolerant && istr.fail(), std::runtime_error,
4977 "Failed to get a real-valued matrix entry from line "
4978 << lineNumber <<
" of the Matrix Market file.");
4981 if (istr.fail() && tolerant) {
4987 TEUCHOS_TEST_FOR_EXCEPTION(
4988 !tolerant && istr.fail(), std::runtime_error,
4989 "Failed to read matrix data from line " << lineNumber
4990 <<
" of the Matrix Market file.");
4993 Teuchos::MatrixMarket::details::assignScalar<ST>(val, real, imag);
4995 curRow = count % numRows;
4996 curCol = count / numRows;
4997 X_view[curRow + curCol * stride] = val;
5002 TEUCHOS_TEST_FOR_EXCEPTION(
5003 !tolerant &&
static_cast<global_size_t>(count) < numRows * numCols,
5005 "The Matrix Market metadata reports that the dense matrix is "
5006 << numRows <<
" x " << numCols <<
", and thus has "
5007 << numRows * numCols <<
" total entries, but we only found " << count
5008 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
5009 }
catch (std::exception& e) {
5011 localReadDataSuccess = 0;
5016 *err << myRank <<
": readVectorImpl: done reading data" << endl;
5020 int globalReadDataSuccess = localReadDataSuccess;
5021 broadcast(*comm, 0, outArg(globalReadDataSuccess));
5022 TEUCHOS_TEST_FOR_EXCEPTION(
5023 globalReadDataSuccess == 0, std::runtime_error,
5024 "Failed to read the multivector's data: " << exMsg.str());
5029 if (comm->getSize() == 1 && map.is_null()) {
5031 if (!err.is_null()) {
5035 *err << myRank <<
": readVectorImpl: done" << endl;
5037 if (!err.is_null()) {
5044 *err << myRank <<
": readVectorImpl: Creating target MV" << endl;
5048 RCP<MV> Y = createVector<ST, LO, GO, NT>(map);
5051 *err << myRank <<
": readVectorImpl: Creating Export" << endl;
5057 Export<LO, GO, NT> exporter(proc0Map, map, err);
5060 *err << myRank <<
": readVectorImpl: Exporting" << endl;
5063 Y->doExport(*X, exporter,
INSERT);
5065 if (!err.is_null()) {
5069 *err << myRank <<
": readVectorImpl: done" << endl;
5071 if (!err.is_null()) {
5100 static Teuchos::RCP<const map_type>
5104 const bool debug =
false,
5105 const bool binary =
false) {
5106 Teuchos::RCP<Teuchos::FancyOStream>
err =
5107 Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
5137 static Teuchos::RCP<const map_type>
5140 const Teuchos::RCP<Teuchos::FancyOStream>&
err,
5142 const bool debug =
false,
5143 const bool binary =
false) {
5145 using Teuchos::arcp;
5146 using Teuchos::Array;
5147 using Teuchos::ArrayRCP;
5149 using Teuchos::broadcast;
5150 using Teuchos::Comm;
5151 using Teuchos::CommRequest;
5152 using Teuchos::inOutArg;
5153 using Teuchos::ireceive;
5154 using Teuchos::isend;
5155 using Teuchos::outArg;
5157 using Teuchos::receive;
5158 using Teuchos::REDUCE_MIN;
5159 using Teuchos::reduceAll;
5160 using Teuchos::SerialComm;
5161 using Teuchos::toString;
5162 using Teuchos::wait;
5170 const int numProcs = comm->getSize();
5171 const int myRank = comm->getRank();
5173 if (
err.is_null()) {
5177 std::ostringstream
os;
5181 if (
err.is_null()) {
5208 std::ostringstream
exMsg;
5225 if (
data.is_null()) {
5227 exMsg <<
"readDenseImpl() returned null." <<
endl;
5229 }
catch (std::exception&
e) {
5253 if (
data->getNumVectors() == 2) {
5261 const int pid =
static_cast<int>(
PIDs[0]);
5265 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5266 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5277 const int pid =
static_cast<int>(
PIDs[
k]);
5280 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5281 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5291 }
else if (
data->getNumVectors() == 1) {
5292 if (
data->getGlobalLength() % 2 !=
static_cast<GST>(0)) {
5294 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data has the "
5295 "wrong format (for Map format 2.0). The global number of rows "
5296 "in the MultiVector must be even (divisible by 2)."
5306 const int pid =
static_cast<int>(
theData[1]);
5309 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5310 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5320 const int pid =
static_cast<int>(
theData[2 *
k + 1]);
5323 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5324 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5336 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data must have "
5337 "either 1 column (for the new Map format 2.0) or 2 columns (for "
5338 "the old Map format 1.0).";
5361 "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
5362 "following exception message: "
5391 typename std::map<int, Array<GO>>::iterator
it =
pid2gids.find(0);
5418 typename std::map<int, Array<GO>>::iterator
it =
pid2gids.find(
p);
5436 std::ostringstream
os;
5440 const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid();
5451 }
catch (std::exception&
e) {
5453 errStrm <<
"Process " << comm->getRank() <<
" threw an exception: "
5454 <<
e.what() << std::endl;
5457 errStrm <<
"Process " << comm->getRank() <<
" threw an exception "
5458 "not a subclass of std::exception"
5461 Teuchos::reduceAll<int, int>(*comm, Teuchos::REDUCE_MIN,
5468 if (
err.is_null()) {
5472 std::ostringstream
os;
5476 if (
err.is_null()) {
5494 encodeDataType(
const std::string&
dataType) {
5497 }
else if (
dataType ==
"complex") {
5499 }
else if (dataType ==
"pattern") {
5505 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
5506 "Unrecognized Matrix Market data type \"" << dataType
5507 <<
"\". We should never get here. "
5508 "Please report this bug to the Tpetra developers.");
5542 static Teuchos::RCP<sparse_matrix_type>
5545 const Teuchos::RCP<const map_type>& rowMap,
5546 Teuchos::RCP<const map_type>& colMap,
5547 const Teuchos::RCP<const map_type>& domainMap,
5548 const Teuchos::RCP<const map_type>&
rangeMap,
5552 const bool debug =
false) {
5556 using STS =
typename Teuchos::ScalarTraits<ST>;
5557 using Teuchos::arcp;
5558 using Teuchos::ArrayRCP;
5567 rowMap.is_null(), std::invalid_argument,
5568 "Row Map must be nonnull.");
5569 Teuchos::RCP<const Teuchos::Comm<int>> comm = rowMap->getComm();
5571 "The input row map's communicator (Teuchos::Comm object) is null.");
5573 rangeMap.is_null(), std::invalid_argument,
5574 "Range Map must be nonnull.");
5576 domainMap.is_null(), std::invalid_argument,
5577 "Domain Map must be nonnull.");
5579 domainMap->getComm().getRawPtr() != comm.getRawPtr(),
5580 std::invalid_argument,
5581 "The specified domain Map's communicator (domainMap->getComm())"
5582 "differs from row Map's communicator");
5584 rangeMap->getComm().getRawPtr() != comm.getRawPtr(),
5585 std::invalid_argument,
5586 "The specified range Map's communicator (rangeMap->getComm())"
5587 "differs from row Map's communicator");
5590 const int myRank = comm->getRank();
5591 const int numProc = comm->getSize();
5602 std::ostringstream
errMsg;
5619 using Teuchos::MatrixMarket::Banner;
5624 }
catch (std::exception&
e) {
5625 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
5626 "threw an exception: "
5637 errMsg <<
"The Matrix Market input file must contain a "
5638 "\"coordinate\"-format sparse matrix in order to create a "
5639 "Tpetra::CrsMatrix object from it, but the file's matrix "
5641 <<
pBanner->matrixType() <<
"\" instead.";
5648 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
5649 "format sparse matrix in order to create a Tpetra::CrsMatrix "
5650 "object from it, but the file's matrix type is \"array\" "
5651 "instead. That probably means the file contains dense matrix "
5657 using Teuchos::MatrixMarket::readCoordinateDimensions;
5664 "# rows in file does not match rowmap.");
5666 "# rows in file does not match colmap.");
5669 typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type>
raw_adder_type;
5671 Teuchos::RCP<raw_adder_type>
pRaw =
5676 std::cerr <<
"-- Reading matrix data" << std::endl;
5681 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
5690 }
catch (std::exception&
e) {
5697 typedef Teuchos::MatrixMarket::Raw::Element<scalar_type, global_ordinal_type>
element_type;
5700 pAdder->getAdder()->merge();
5703 const std::vector<element_type>& entries =
pAdder->getAdder()->getEntries();
5706 const size_t numEntries = (
size_t)entries.size();
5709 std::cerr <<
"----- Proc " <<
myRank <<
": Matrix has numRows=" <<
numRows
5710 <<
" rows and numEntries=" << numEntries
5711 <<
" entries." << std::endl;
5728 LO
INVALID = Teuchos::OrdinalTraits<LO>::invalid();
5737 "Current global row " <<
curRow <<
" is invalid.");
5740 "Row indices are out of order, even though they are supposed "
5741 "to be sorted. curRow = "
5744 "this bug to the Tpetra developers.");
5769 if (colMap.is_null()) {
5771 for (
size_t i = 0;
i < rowMap->getLocalNumElements();
i++) {
5772 GO
g_row = rowMap->getGlobalElement(
i);
5780 throw std::runtime_error(
"Reading with a column map is not yet implemented");
5830 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
5898 const bool debug =
false) {
5901 "The input matrix's communicator (Teuchos::Comm object) is null.");
5902 const int myRank = comm->getRank();
5915 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
5918 const bool debug =
false) {
5920 "The input matrix is null.");
5947 const bool debug =
false) {
5954 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
5955 const bool debug =
false) {
5994 const bool debug =
false) {
5997 using Teuchos::ArrayView;
5998 using Teuchos::Comm;
5999 using Teuchos::FancyOStream;
6000 using Teuchos::getFancyOStream;
6001 using Teuchos::null;
6003 using Teuchos::rcpFromRef;
6007 using STS =
typename Teuchos::ScalarTraits<ST>;
6013 Teuchos::SetScientific<ST>
sci(
out);
6018 comm.is_null(), std::invalid_argument,
6019 "The input matrix's communicator (Teuchos::Comm object) is null.");
6020 const int myRank = comm->getRank();
6026 std::ostringstream
os;
6043 const global_size_t numCols = domainMap->getMaxAllGlobalIndex() + 1 - domainMap->getIndexBase();
6045 if (debug &&
myRank == 0) {
6046 std::ostringstream
os;
6047 os <<
"-- Input sparse matrix is:"
6050 << (
matrix.isGloballyIndexed() ?
"Globally" :
"Locally")
6051 <<
" indexed." <<
endl
6052 <<
"---- Its row map has " << rowMap->getGlobalNumElements()
6053 <<
" elements." <<
endl
6054 <<
"---- Its col map has " << colMap->getGlobalNumElements()
6055 <<
" elements." <<
endl;
6062 std::ostringstream
os;
6063 os <<
"-- " <<
myRank <<
": making gatherRowMap" <<
endl;
6067 Details::computeGatherMap(rowMap,
err, debug);
6077 Details::computeGatherMap(colMap,
err, debug);
6088 static_cast<size_t>(0)));
6097 domainMap->getIndexBase(),
6107 cerr <<
"-- Output sparse matrix is:"
6108 <<
"---- " <<
newMatrix->getRangeMap()->getGlobalNumElements()
6110 <<
newMatrix->getDomainMap()->getGlobalNumElements()
6114 << (
newMatrix->isGloballyIndexed() ?
"Globally" :
"Locally")
6115 <<
" indexed." <<
endl
6116 <<
"---- Its row map has "
6117 <<
newMatrix->getRowMap()->getGlobalNumElements()
6118 <<
" elements, with index base "
6120 <<
"---- Its col map has "
6121 <<
newMatrix->getColMap()->getGlobalNumElements()
6122 <<
" elements, with index base "
6124 <<
"---- Element count of output matrix's column Map may differ "
6125 <<
"from that of the input matrix's column Map, if some columns "
6126 <<
"of the matrix contain no entries." <<
endl;
6139 out <<
"%%MatrixMarket matrix coordinate "
6140 << (STS::isComplex ?
"complex" :
"real")
6141 <<
" general" <<
endl;
6158 out <<
newMatrix->getRangeMap()->getMaxAllGlobalIndex() + 1 -
newMatrix->getRangeMap()->getIndexBase() <<
" "
6159 <<
newMatrix->getDomainMap()->getMaxAllGlobalIndex() + 1 -
newMatrix->getDomainMap()->getIndexBase() <<
" "
6182 typename sparse_matrix_type::global_inds_host_view_type
ind;
6183 typename sparse_matrix_type::values_host_view_type
val;
6185 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6191 if (STS::isComplex) {
6200 using OTG = Teuchos::OrdinalTraits<GO>;
6209 "Failed to convert the supposed local row index "
6211 "Please report this bug to the Tpetra developers.");
6212 typename sparse_matrix_type::local_inds_host_view_type
ind;
6213 typename sparse_matrix_type::values_host_view_type
val;
6215 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6221 "On local row " <<
localRowIndex <<
" of the sparse matrix: "
6222 "Failed to convert the supposed local column index "
6223 <<
ind(
ii) <<
" into a global column index. Please report "
6224 "this bug to the Tpetra developers.");
6229 if (STS::isComplex) {
6244 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
6247 const bool debug =
false) {
6249 "The input matrix is null.");
6288 const bool debug =
false) {
6291 using Teuchos::ArrayView;
6292 using Teuchos::Comm;
6293 using Teuchos::FancyOStream;
6294 using Teuchos::getFancyOStream;
6295 using Teuchos::null;
6297 using Teuchos::rcpFromRef;
6305 auto rowMap =
graph.getRowMap();
6306 if (rowMap.is_null()) {
6309 auto comm = rowMap->getComm();
6310 if (comm.is_null()) {
6313 const int myRank = comm->getRank();
6319 std::ostringstream
os;
6331 auto colMap =
graph.getColMap();
6332 auto domainMap =
graph.getDomainMap();
6336 const global_size_t numCols = domainMap->getGlobalNumElements();
6338 if (debug &&
myRank == 0) {
6339 std::ostringstream
os;
6340 os <<
"-- Input sparse graph is:"
6341 <<
"---- " <<
numRows <<
" x " << numCols <<
" with "
6342 <<
graph.getGlobalNumEntries() <<
" entries;" <<
endl
6344 << (
graph.isGloballyIndexed() ?
"Globally" :
"Locally")
6345 <<
" indexed." <<
endl
6346 <<
"---- Its row Map has " << rowMap->getGlobalNumElements()
6347 <<
" elements." <<
endl
6348 <<
"---- Its col Map has " << colMap->getGlobalNumElements()
6349 <<
" elements." <<
endl;
6356 std::ostringstream
os;
6357 os <<
"-- " <<
myRank <<
": making gatherRowMap" <<
endl;
6378 static_cast<size_t>(0));
6387 domainMap->getIndexBase(),
6397 cerr <<
"-- Output sparse graph is:"
6398 <<
"---- " <<
newGraph.getRangeMap()->getGlobalNumElements()
6400 <<
newGraph.getDomainMap()->getGlobalNumElements()
6402 <<
newGraph.getGlobalNumEntries() <<
" entries;" <<
endl
6404 << (
newGraph.isGloballyIndexed() ?
"Globally" :
"Locally")
6405 <<
" indexed." <<
endl
6406 <<
"---- Its row map has "
6407 <<
newGraph.getRowMap()->getGlobalNumElements()
6408 <<
" elements, with index base "
6410 <<
"---- Its col map has "
6411 <<
newGraph.getColMap()->getGlobalNumElements()
6412 <<
" elements, with index base "
6414 <<
"---- Element count of output graph's column Map may differ "
6415 <<
"from that of the input matrix's column Map, if some columns "
6416 <<
"of the matrix contain no entries." <<
endl;
6430 out <<
"%%MatrixMarket matrix coordinate pattern general" <<
endl;
6448 out <<
newGraph.getRangeMap()->getMaxAllGlobalIndex() + 1 -
newGraph.getRangeMap()->getIndexBase() <<
" "
6449 <<
newGraph.getDomainMap()->getMaxAllGlobalIndex() + 1 -
newGraph.getDomainMap()->getIndexBase() <<
" "
6464 if (
newGraph.isGloballyIndexed()) {
6472 typename crs_graph_type::global_inds_host_view_type
ind;
6474 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6484 typedef Teuchos::OrdinalTraits<GO>
OTG;
6493 "to convert the supposed local row index "
6494 <<
localRowIndex <<
" into a global row index. Please report this bug to the "
6495 "Tpetra developers.");
6496 typename crs_graph_type::local_inds_host_view_type
ind;
6498 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6504 "On local row " <<
localRowIndex <<
" of the sparse graph: "
6505 "Failed to convert the supposed local column index "
6506 <<
ind(
ii) <<
" into a global column index. Please report "
6507 "this bug to the Tpetra developers.");
6527 const bool debug =
false) {
6570 const bool debug =
false) {
6571 auto comm =
graph.getComm();
6572 if (comm.is_null()) {
6580 const int myRank = comm->getRank();
6597 const bool debug =
false) {
6611 const Teuchos::RCP<const crs_graph_type>&
pGraph,
6614 const bool debug =
false) {
6629 const Teuchos::RCP<const crs_graph_type>&
pGraph,
6630 const bool debug =
false) {
6659 const bool debug =
false) {
6666 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
6667 const bool debug =
false) {
6704 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6705 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
6724 const Teuchos::RCP<const multivector_type>&
X,
6727 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6728 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
6730 X.is_null(), std::invalid_argument,
6731 "Tpetra::MatrixMarket::"
6732 "writeDenseFile: The input MultiVector X is null.");
6744 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6745 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
6756 const Teuchos::RCP<const multivector_type>&
X,
6757 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6758 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
6760 X.is_null(), std::invalid_argument,
6761 "Tpetra::MatrixMarket::"
6762 "writeDenseFile: The input MultiVector X is null.");
6801 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6802 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
6804 using Teuchos::Comm;
6805 using Teuchos::outArg;
6806 using Teuchos::REDUCE_MAX;
6807 using Teuchos::reduceAll;
6815 const bool debug = !
dbg.is_null();
6818 std::ostringstream
os;
6828 const size_t numVecs =
X.getNumVectors();
6830 writeDenseColumn(
out, *(
X.getVector(
j)),
err,
dbg);
6835 std::ostringstream
os;
6848 static std::ofstream openOutFileOnRankZero(
6851 const std::ios_base::openmode
mode = std::ios_base::out) {
6870 "Could not open output file '" +
filename +
"' on root rank 0.");
6901 writeDenseHeader(std::ostream& out,
6903 const std::string& matrixName,
6904 const std::string& matrixDescription,
6905 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6906 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null) {
6908 using Teuchos::Comm;
6909 using Teuchos::outArg;
6911 using Teuchos::REDUCE_MAX;
6912 using Teuchos::reduceAll;
6913 typedef Teuchos::ScalarTraits<scalar_type> STS;
6914 const char prefix[] =
"Tpetra::MatrixMarket::writeDenseHeader: ";
6924 const bool debug = !dbg.is_null();
6928 std::ostringstream os;
6929 os << myRank <<
": writeDenseHeader" << endl;
6948 std::ostringstream hdr;
6950 std::string dataType;
6951 if (STS::isComplex) {
6952 dataType =
"complex";
6953 }
else if (STS::isOrdinal) {
6954 dataType =
"integer";
6958 hdr <<
"%%MatrixMarket matrix array " << dataType <<
" general"
6963 if (matrixName !=
"") {
6964 printAsComment(hdr, matrixName);
6966 if (matrixDescription !=
"") {
6967 printAsComment(hdr, matrixDescription);
6970 hdr << X.getGlobalLength() <<
" " << X.getNumVectors() << endl;
6974 }
catch (std::exception& e) {
6975 if (!err.is_null()) {
6976 *err << prefix <<
"While writing the Matrix Market header, "
6977 "Process 0 threw an exception: "
6978 << e.what() << endl;
6987 reduceAll<int, int>(*comm, REDUCE_MAX, lclErr, outArg(gblErr));
6988 TEUCHOS_TEST_FOR_EXCEPTION(
6989 gblErr == 1, std::runtime_error, prefix <<
"Some error occurred "
6990 "which prevented this method from completing.");
6994 *dbg << myRank <<
": writeDenseHeader: Done" << endl;
7017 writeDenseColumn(std::ostream& out,
7019 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7020 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null) {
7022 using Teuchos::arcp;
7023 using Teuchos::Array;
7024 using Teuchos::ArrayRCP;
7025 using Teuchos::ArrayView;
7026 using Teuchos::Comm;
7027 using Teuchos::CommRequest;
7028 using Teuchos::ireceive;
7029 using Teuchos::isend;
7030 using Teuchos::outArg;
7032 using Teuchos::REDUCE_MAX;
7033 using Teuchos::reduceAll;
7034 using Teuchos::TypeNameTraits;
7035 using Teuchos::wait;
7036 typedef Teuchos::ScalarTraits<scalar_type> STS;
7038 const Comm<int>& comm = *(X.getMap()->getComm());
7039 const int myRank = comm.getRank();
7040 const int numProcs = comm.getSize();
7047 const bool debug = !dbg.is_null();
7051 std::ostringstream os;
7052 os << myRank <<
": writeDenseColumn" << endl;
7061 Teuchos::SetScientific<scalar_type> sci(out);
7063 const size_t myNumRows = X.getLocalLength();
7064 const size_t numCols = X.getNumVectors();
7067 const int sizeTag = 1337;
7068 const int dataTag = 1338;
7129 RCP<CommRequest<int>> sendReqSize, sendReqData;
7135 Array<ArrayRCP<size_t>> recvSizeBufs(3);
7136 Array<ArrayRCP<scalar_type>> recvDataBufs(3);
7137 Array<RCP<CommRequest<int>>> recvSizeReqs(3);
7138 Array<RCP<CommRequest<int>>> recvDataReqs(3);
7141 ArrayRCP<size_t> sendDataSize(1);
7142 sendDataSize[0] = myNumRows;
7146 std::ostringstream os;
7147 os << myRank <<
": Post receive-size receives from "
7148 "Procs 1 and 2: tag = "
7153 recvSizeBufs[0].resize(1);
7157 (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid();
7158 recvSizeBufs[1].resize(1);
7159 (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid();
7160 recvSizeBufs[2].resize(1);
7161 (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid();
7164 ireceive<int, size_t>(recvSizeBufs[1], 1, sizeTag, comm);
7168 ireceive<int, size_t>(recvSizeBufs[2], 2, sizeTag, comm);
7170 }
else if (myRank == 1 || myRank == 2) {
7172 std::ostringstream os;
7173 os << myRank <<
": Post send-size send: size = "
7174 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7181 sendReqSize = isend<int, size_t>(sendDataSize, 0, sizeTag, comm);
7184 std::ostringstream os;
7185 os << myRank <<
": Not posting my send-size send yet" << endl;
7194 std::ostringstream os;
7195 os << myRank <<
": Pack my entries" << endl;
7198 ArrayRCP<scalar_type> sendDataBuf;
7200 sendDataBuf = arcp<scalar_type>(myNumRows * numCols);
7201 X.get1dCopy(sendDataBuf(), myNumRows);
7202 }
catch (std::exception& e) {
7204 if (!err.is_null()) {
7205 std::ostringstream os;
7206 os <<
"Process " << myRank <<
": Attempt to pack my MultiVector "
7207 "entries threw an exception: "
7208 << e.what() << endl;
7213 std::ostringstream os;
7214 os << myRank <<
": Done packing my entries" << endl;
7223 *dbg << myRank <<
": Post send-data send: tag = " << dataTag
7226 sendReqData = isend<int, scalar_type>(sendDataBuf, 0, dataTag, comm);
7234 std::ostringstream os;
7235 os << myRank <<
": Write my entries" << endl;
7241 const size_t printNumRows = myNumRows;
7242 ArrayView<const scalar_type> printData = sendDataBuf();
7243 const size_t printStride = printNumRows;
7244 if (
static_cast<size_t>(printData.size()) < printStride * numCols) {
7246 if (!err.is_null()) {
7247 std::ostringstream os;
7248 os <<
"Process " << myRank <<
": My MultiVector data's size "
7249 << printData.size() <<
" does not match my local dimensions "
7250 << printStride <<
" x " << numCols <<
"." << endl;
7257 for (
size_t col = 0; col < numCols; ++col) {
7258 for (
size_t row = 0; row < printNumRows; ++row) {
7259 if (STS::isComplex) {
7260 out << STS::real(printData[row + col * printStride]) <<
" "
7261 << STS::imag(printData[row + col * printStride]) << endl;
7263 out << printData[row + col * printStride] << endl;
7272 const int recvRank = 1;
7273 const int circBufInd = recvRank % 3;
7275 std::ostringstream os;
7276 os << myRank <<
": Wait on receive-size receive from Process "
7277 << recvRank << endl;
7281 wait<int>(comm, outArg(recvSizeReqs[circBufInd]));
7285 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7286 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid()) {
7288 if (!err.is_null()) {
7289 std::ostringstream os;
7290 os << myRank <<
": Result of receive-size receive from Process "
7291 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7292 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid() <<
". "
7293 "This should never happen, and suggests that the receive never "
7294 "got posted. Please report this bug to the Tpetra developers."
7309 recvDataBufs[circBufInd].resize(recvNumRows * numCols);
7311 std::ostringstream os;
7312 os << myRank <<
": Post receive-data receive from Process "
7313 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7314 << recvDataBufs[circBufInd].size() << endl;
7317 if (!recvSizeReqs[circBufInd].is_null()) {
7319 if (!err.is_null()) {
7320 std::ostringstream os;
7321 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7322 "null, before posting the receive-data receive from Process "
7323 << recvRank <<
". This should never happen. Please report "
7324 "this bug to the Tpetra developers."
7329 recvDataReqs[circBufInd] =
7330 ireceive<int, scalar_type>(recvDataBufs[circBufInd],
7331 recvRank, dataTag, comm);
7333 }
else if (myRank == 1) {
7336 std::ostringstream os;
7337 os << myRank <<
": Wait on my send-size send" << endl;
7340 wait<int>(comm, outArg(sendReqSize));
7346 for (
int p = 1; p < numProcs; ++p) {
7348 if (p + 2 < numProcs) {
7350 const int recvRank = p + 2;
7351 const int circBufInd = recvRank % 3;
7353 std::ostringstream os;
7354 os << myRank <<
": Post receive-size receive from Process "
7355 << recvRank <<
": tag = " << sizeTag << endl;
7358 if (!recvSizeReqs[circBufInd].is_null()) {
7360 if (!err.is_null()) {
7361 std::ostringstream os;
7362 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7363 <<
"null, for the receive-size receive from Process "
7364 << recvRank <<
"! This may mean that this process never "
7365 <<
"finished waiting for the receive from Process "
7366 << (recvRank - 3) <<
"." << endl;
7370 recvSizeReqs[circBufInd] =
7371 ireceive<int, size_t>(recvSizeBufs[circBufInd],
7372 recvRank, sizeTag, comm);
7375 if (p + 1 < numProcs) {
7376 const int recvRank = p + 1;
7377 const int circBufInd = recvRank % 3;
7381 std::ostringstream os;
7382 os << myRank <<
": Wait on receive-size receive from Process "
7383 << recvRank << endl;
7386 wait<int>(comm, outArg(recvSizeReqs[circBufInd]));
7390 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7391 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid()) {
7393 if (!err.is_null()) {
7394 std::ostringstream os;
7395 os << myRank <<
": Result of receive-size receive from Process "
7396 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7397 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid() <<
". "
7398 "This should never happen, and suggests that the receive never "
7399 "got posted. Please report this bug to the Tpetra developers."
7413 recvDataBufs[circBufInd].resize(recvNumRows * numCols);
7415 std::ostringstream os;
7416 os << myRank <<
": Post receive-data receive from Process "
7417 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7418 << recvDataBufs[circBufInd].size() << endl;
7421 if (!recvDataReqs[circBufInd].is_null()) {
7423 if (!err.is_null()) {
7424 std::ostringstream os;
7425 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
7426 <<
"null, for the receive-data receive from Process "
7427 << recvRank <<
"! This may mean that this process never "
7428 <<
"finished waiting for the receive from Process "
7429 << (recvRank - 3) <<
"." << endl;
7433 recvDataReqs[circBufInd] =
7434 ireceive<int, scalar_type>(recvDataBufs[circBufInd],
7435 recvRank, dataTag, comm);
7439 const int recvRank = p;
7440 const int circBufInd = recvRank % 3;
7442 std::ostringstream os;
7443 os << myRank <<
": Wait on receive-data receive from Process "
7444 << recvRank << endl;
7447 wait<int>(comm, outArg(recvDataReqs[circBufInd]));
7454 std::ostringstream os;
7455 os << myRank <<
": Write entries from Process " << recvRank
7457 *dbg << os.str() << endl;
7459 size_t printNumRows = (recvSizeBufs[circBufInd])[0];
7460 if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid()) {
7462 if (!err.is_null()) {
7463 std::ostringstream os;
7464 os << myRank <<
": Result of receive-size receive from Process "
7465 << recvRank <<
" was Teuchos::OrdinalTraits<size_t>::"
7467 << Teuchos::OrdinalTraits<size_t>::invalid()
7468 <<
". This should never happen, and suggests that its "
7469 "receive-size receive was never posted. "
7470 "Please report this bug to the Tpetra developers."
7478 if (printNumRows > 0 && recvDataBufs[circBufInd].is_null()) {
7480 if (!err.is_null()) {
7481 std::ostringstream os;
7482 os << myRank <<
": Result of receive-size receive from Proc "
7483 << recvRank <<
" was " << printNumRows <<
" > 0, but "
7485 << circBufInd <<
"] is null. This should "
7486 "never happen. Please report this bug to the Tpetra "
7498 ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd])();
7499 const size_t printStride = printNumRows;
7503 for (
size_t col = 0; col < numCols; ++col) {
7504 for (
size_t row = 0; row < printNumRows; ++row) {
7505 if (STS::isComplex) {
7506 out << STS::real(printData[row + col * printStride]) <<
" "
7507 << STS::imag(printData[row + col * printStride]) << endl;
7509 out << printData[row + col * printStride] << endl;
7513 }
else if (myRank == p) {
7516 std::ostringstream os;
7517 os << myRank <<
": Wait on my send-data send" << endl;
7520 wait<int>(comm, outArg(sendReqData));
7521 }
else if (myRank == p + 1) {
7524 std::ostringstream os;
7525 os << myRank <<
": Post send-data send: tag = " << dataTag
7529 sendReqData = isend<int, scalar_type>(sendDataBuf, 0, dataTag, comm);
7532 std::ostringstream os;
7533 os << myRank <<
": Wait on my send-size send" << endl;
7536 wait<int>(comm, outArg(sendReqSize));
7537 }
else if (myRank == p + 2) {
7540 std::ostringstream os;
7541 os << myRank <<
": Post send-size send: size = "
7542 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7545 sendReqSize = isend<int, size_t>(sendDataSize, 0, sizeTag, comm);
7550 reduceAll<int, int>(comm, REDUCE_MAX, lclErr, outArg(gblErr));
7551 TEUCHOS_TEST_FOR_EXCEPTION(
7552 gblErr == 1, std::runtime_error,
7553 "Tpetra::MatrixMarket::writeDense "
7554 "experienced some kind of error and was unable to complete.");
7558 *dbg << myRank <<
": writeDenseColumn: Done" << endl;
7571 const Teuchos::RCP<const multivector_type>&
X,
7574 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7575 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
7577 X.is_null(), std::invalid_argument,
7578 "Tpetra::MatrixMarket::"
7579 "writeDense: The input MultiVector X is null.");
7591 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7592 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
7603 const Teuchos::RCP<const multivector_type>&
X,
7604 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7605 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null) {
7607 X.is_null(), std::invalid_argument,
7608 "Tpetra::MatrixMarket::"
7609 "writeDense: The input MultiVector X is null.");
7634 Teuchos::RCP<Teuchos::FancyOStream>
err =
7635 Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
7650 const Teuchos::RCP<Teuchos::FancyOStream>&
err,
7651 const bool debug =
false) {
7653 using Teuchos::Array;
7654 using Teuchos::ArrayRCP;
7655 using Teuchos::ArrayView;
7656 using Teuchos::Comm;
7657 using Teuchos::CommRequest;
7658 using Teuchos::ireceive;
7659 using Teuchos::isend;
7661 using Teuchos::TypeNameTraits;
7662 using Teuchos::wait;
7681 sizeof(GO) >
sizeof(
int_type), std::logic_error,
7683 <<
" is too big for ptrdiff_t. sizeof(GO) = " <<
sizeof(GO)
7684 <<
" > sizeof(ptrdiff_t) = " <<
sizeof(
ptrdiff_t) <<
".");
7688 "sizeof(pid_type) = "
7689 <<
sizeof(
pid_type) <<
" > sizeof(ptrdiff_t)"
7694 const int myRank = comm.getRank();
7695 const int numProcs = comm.getSize();
7697 if (!
err.is_null()) {
7701 std::ostringstream
os;
7705 if (!
err.is_null()) {
7789 std::ostringstream
os;
7790 os <<
myRank <<
": Post receive-size receives from "
7791 "Procs 1 and 2: tag = "
7812 std::ostringstream
os;
7813 os <<
myRank <<
": Post send-size send: size = "
7824 std::ostringstream
os;
7825 os <<
myRank <<
": Not posting my send-size send yet" <<
endl;
7836 std::ostringstream
os;
7843 if (
map.isContiguous()) {
7863 std::ostringstream
os;
7864 os <<
myRank <<
": Done packing my GIDs and PIDs" <<
endl;
7884 std::ostringstream
hdr;
7888 hdr <<
"%%MatrixMarket matrix array integer general" <<
endl
7889 <<
"% Format: Version 2.0" <<
endl
7891 <<
"% This file encodes a Tpetra::Map." <<
endl
7892 <<
"% It is stored as a dense vector, with twice as many " <<
endl
7893 <<
"% entries as the global number of GIDs (global indices)." <<
endl
7894 <<
"% (GID, PID) pairs are stored contiguously, where the PID " <<
endl
7895 <<
"% is the rank of the process owning that GID." <<
endl
7896 << (2 *
map.getGlobalNumElements()) <<
" " << 1 <<
endl;
7900 std::ostringstream
os;
7922 std::ostringstream
os;
7923 os <<
myRank <<
": Wait on receive-size receive from Process "
7934 std::ostringstream
os;
7935 os <<
myRank <<
": Result of receive-size receive from Process "
7936 <<
recvRank <<
" is -1. This should never happen, and "
7937 "suggests that the receive never got posted. Please report "
7938 "this bug to the Tpetra developers."
7946 std::ostringstream
os;
7947 os <<
myRank <<
": Post receive-data receive from Process "
7953 std::ostringstream
os;
7955 "null, before posting the receive-data receive from Process "
7956 <<
recvRank <<
". This should never happen. Please report "
7957 "this bug to the Tpetra developers."
7965 }
else if (
myRank == 1) {
7968 std::ostringstream
os;
7985 std::ostringstream
os;
7986 os <<
myRank <<
": Post receive-size receive from Process "
7991 std::ostringstream
os;
7993 <<
"null, for the receive-size receive from Process "
7994 <<
recvRank <<
"! This may mean that this process never "
7995 <<
"finished waiting for the receive from Process "
8010 std::ostringstream
os;
8011 os <<
myRank <<
": Wait on receive-size receive from Process "
8021 std::ostringstream
os;
8022 os <<
myRank <<
": Result of receive-size receive from Process "
8023 <<
recvRank <<
" is -1. This should never happen, and "
8024 "suggests that the receive never got posted. Please report "
8025 "this bug to the Tpetra developers."
8033 std::ostringstream
os;
8034 os <<
myRank <<
": Post receive-data receive from Process "
8040 std::ostringstream
os;
8042 <<
"null, for the receive-data receive from Process "
8043 <<
recvRank <<
"! This may mean that this process never "
8044 <<
"finished waiting for the receive from Process "
8057 std::ostringstream
os;
8058 os <<
myRank <<
": Wait on receive-data receive from Process "
8069 std::ostringstream
os;
8070 os <<
myRank <<
": Write GIDs and PIDs from Process "
8076 std::ostringstream
os;
8077 os <<
myRank <<
": Result of receive-size receive from Process "
8078 <<
recvRank <<
" was -1. This should never happen, and "
8079 "suggests that its receive-size receive was never posted. "
8080 "Please report this bug to the Tpetra developers."
8085 std::ostringstream
os;
8086 os <<
myRank <<
": Result of receive-size receive from Proc "
8090 "never happen. Please report this bug to the Tpetra "
8105 std::ostringstream
os;
8113 std::ostringstream
os;
8121 std::ostringstream
os;
8129 std::ostringstream
os;
8130 os <<
myRank <<
": Post send-size send: size = "
8138 if (!
err.is_null()) {
8144 if (!
err.is_null()) {
8153 const int myRank =
map.getComm()->getRank();
8188 printAsComment(std::ostream&
out,
const std::string&
str) {
8194 if (!
line.empty()) {
8197 if (
line[0] ==
'%') {
8227 Teuchos::ParameterList
pl;
8253 Teuchos::ParameterList
pl;
8296 const Teuchos::ParameterList&
params) {
8299 const int myRank =
A.getDomainMap()->getComm()->getRank();
8312 "writeOperator: temporary file " <<
tmpFile <<
" already exists");
8314 if (
params.isParameter(
"precision")) {
8328 if (
params.isParameter(
"print MatrixMarket header"))
8335 std::ifstream src(
tmpFile, std::ios_base::binary);
8384 const Teuchos::ParameterList&
params) {
8385 const int myRank =
A.getDomainMap()->getComm()->getRank();
8402 std::ostringstream
tmpOut;
8404 if (
params.isParameter(
"precision") &&
params.isType<
int>(
"precision")) {
8413 if (
params.isParameter(
"print MatrixMarket header") &&
8414 params.isType<
bool>(
"print MatrixMarket header")) {
8446 writeOperatorImpl(std::ostream&
os,
8447 const operator_type&
A,
8448 const Teuchos::ParameterList&
params) {
8449 using Teuchos::Array;
8450 using Teuchos::ArrayRCP;
8457 typedef Teuchos::OrdinalTraits<LO>
TLOT;
8458 typedef Teuchos::OrdinalTraits<GO>
TGOT;
8462 const map_type& domainMap = *(
A.getDomainMap());
8465 const int myRank = comm->getRank();
8466 const size_t numProcs = comm->getSize();
8469 if (
params.isParameter(
"probing size"))
8473 GO
minColGid = domainMap.getMinAllGlobalIndex();
8474 GO
maxColGid = domainMap.getMaxAllGlobalIndex();
8485 if (
params.isParameter(
"zero-based indexing")) {
8486 if (
params.get<
bool>(
"zero-based indexing") ==
true)
8512 importGidList.reserve(numTargetMapEntries);
8514 RCP<map_type> importGidMap = rcp(
new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
8517 import_type gidImporter(allGidsMap, importGidMap);
8518 mv_type_go importedGids(importGidMap, 1);
8519 importedGids.doImport(allGids, gidImporter,
INSERT);
8522 ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
8523 RCP<const map_type> importMap = rcp(
new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm));
8526 import_type importer(rangeMap, importMap);
8528 RCP<mv_type> colsOnPid0 = rcp(
new mv_type(importMap, numMVs));
8530 RCP<mv_type> ei = rcp(
new mv_type(A.getDomainMap(), numMVs));
8531 RCP<mv_type> colsA = rcp(
new mv_type(A.getRangeMap(), numMVs));
8533 Array<GO> globalColsArray, localColsArray;
8534 globalColsArray.reserve(numMVs);
8535 localColsArray.reserve(numMVs);
8537 ArrayRCP<ArrayRCP<Scalar>> eiData(numMVs);
8538 for (
size_t i = 0; i < numMVs; ++i)
8539 eiData[i] = ei->getDataNonConst(i);
8544 for (GO k = 0; k < numChunks; ++k) {
8545 for (
size_t j = 0; j < numMVs; ++j) {
8547 GO curGlobalCol = minColGid + k * numMVs + j;
8548 globalColsArray.push_back(curGlobalCol);
8550 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8551 if (curLocalCol != TLOT::invalid()) {
8552 eiData[j][curLocalCol] = TGOT::one();
8553 localColsArray.push_back(curLocalCol);
8558 for (
size_t i = 0; i < numMVs; ++i)
8559 eiData[i] = Teuchos::null;
8561 A.apply(*ei, *colsA);
8563 colsOnPid0->doImport(*colsA, importer,
INSERT);
8566 globalNnz += writeColumns(os, *colsOnPid0, numMVs, importedGidsData(),
8567 globalColsArray, offsetToUseInPrinting);
8570 for (
size_t i = 0; i < numMVs; ++i)
8571 eiData[i] = ei->getDataNonConst(i);
8572 for (
size_t j = 0; j < numMVs; ++j) {
8573 GO curGlobalCol = minColGid + k * numMVs + j;
8575 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8576 if (curLocalCol != TLOT::invalid()) {
8577 eiData[j][curLocalCol] = TGOT::one();
8582 for (
size_t j = 0; j < numMVs; ++j) {
8583 for (
int i = 0; i < localColsArray.size(); ++i)
8584 eiData[j][localColsArray[i]] = TGOT::zero();
8586 globalColsArray.clear();
8587 localColsArray.clear();
8594 for (
int j = 0; j < rem; ++j) {
8595 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8596 globalColsArray.push_back(curGlobalCol);
8598 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8599 if (curLocalCol != TLOT::invalid()) {
8600 eiData[j][curLocalCol] = TGOT::one();
8601 localColsArray.push_back(curLocalCol);
8606 for (
size_t i = 0; i < numMVs; ++i)
8607 eiData[i] = Teuchos::null;
8609 A.apply(*ei, *colsA);
8611 colsOnPid0->doImport(*colsA, importer,
INSERT);
8613 globalNnz += writeColumns(os, *colsOnPid0, rem, importedGidsData(),
8614 globalColsArray, offsetToUseInPrinting);
8617 for (
size_t i = 0; i < numMVs; ++i)
8618 eiData[i] = ei->getDataNonConst(i);
8619 for (
int j = 0; j < rem; ++j) {
8620 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8622 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8623 if (curLocalCol != TLOT::invalid()) {
8624 eiData[j][curLocalCol] = TGOT::one();
8629 for (
int j = 0; j < rem; ++j) {
8630 for (
int i = 0; i < localColsArray.size(); ++i)
8631 eiData[j][localColsArray[i]] = TGOT::zero();
8633 globalColsArray.clear();
8634 localColsArray.clear();
8642 std::ostringstream oss;
8644 oss <<
"%%MatrixMarket matrix coordinate ";
8645 if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
8650 oss <<
" general" << std::endl;
8651 oss <<
"% Tpetra::Operator" << std::endl;
8652 std::time_t now = std::time(NULL);
8653 oss <<
"% time stamp: " << ctime(&now);
8654 oss <<
"% written from " << numProcs <<
" processes" << std::endl;
8655 size_t numRows = rangeMap->getGlobalNumElements();
8656 size_t numCols = domainMap.getGlobalNumElements();
8657 oss << numRows <<
" " << numCols <<
" " << globalNnz << std::endl;
8664 writeColumns(std::ostream& os, mv_type
const& colsA,
size_t const& numCols,
8665 Teuchos::ArrayView<const global_ordinal_type>
const& rowGids,
8666 Teuchos::Array<global_ordinal_type>
const& colsArray,
8670 typedef Teuchos::ScalarTraits<Scalar> STS;
8673 const Scalar zero = STS::zero();
8674 const size_t numRows = colsA.getGlobalLength();
8675 for (
size_t j = 0; j < numCols; ++j) {
8676 Teuchos::ArrayRCP<const Scalar>
const curCol = colsA.getData(j);
8677 const GO J = colsArray[j];
8678 for (
size_t i = 0; i < numRows; ++i) {
8679 const Scalar val = curCol[i];
8681 os << rowGids[i] + indexBase <<
" " << J + indexBase <<
" " << val << std::endl;
8705 const bool debug =
false) {
8709 using STS =
typename Teuchos::ScalarTraits<ST>;
8715 "The input matrix's communicator (Teuchos::Comm object) is null.");
8717 "The input matrix must not be GloballyIndexed and must be fillComplete.");
8720 const int myRank = comm->getRank();
8721 const int numProc = comm->getSize();
8743 out <<
"%%MatrixMarket matrix coordinate "
8744 << (STS::isComplex ?
"complex" :
"real")
8745 <<
" general" << std::endl;
8765 Teuchos::SetScientific<ST>
sci(
out);
8770 typename sparse_matrix_type::local_inds_host_view_type indices;
8771 typename sparse_matrix_type::values_host_view_type
values;
8773 for (
size_t ii = 0;
ii < indices.extent(0);
ii++) {
8774 const GO
g_col = colMap->getGlobalElement(indices(
ii));
8779 if (STS::isComplex) {
8801 template <
typename T>
8803 return obj.is_null() ? Teuchos::null :
obj->getComm();
8808 return comm.is_null() ? 0 : comm->getRank();