42 #ifndef __MatrixMarket_Tpetra_hpp 43 #define __MatrixMarket_Tpetra_hpp 56 #include "Tpetra_Details_gathervPrint.hpp" 57 #include "Tpetra_CrsMatrix.hpp" 58 #include "Tpetra_Operator.hpp" 59 #include "Tpetra_Vector.hpp" 61 #include "Teuchos_MatrixMarket_Raw_Adder.hpp" 62 #include "Teuchos_MatrixMarket_Raw_Graph_Adder.hpp" 63 #include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp" 64 #include "Teuchos_MatrixMarket_SymmetrizingGraphAdder.hpp" 65 #include "Teuchos_MatrixMarket_assignScalar.hpp" 66 #include "Teuchos_MatrixMarket_Banner.hpp" 67 #include "Teuchos_MatrixMarket_CoordDataReader.hpp" 68 #include "Teuchos_MatrixMarket_SetScientific.hpp" 163 template<
class SparseMatrixType>
168 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
183 typedef typename SparseMatrixType::global_ordinal_type
205 typedef Teuchos::Comm<int> comm_type;
209 typedef Teuchos::RCP<const comm_type> comm_ptr;
210 typedef Teuchos::RCP<const map_type> map_ptr;
211 typedef Teuchos::RCP<node_type> node_ptr;
219 typedef Teuchos::ArrayRCP<int>::size_type size_type;
232 static Teuchos::RCP<const map_type>
233 makeRangeMap (
const Teuchos::RCP<const comm_type>& pComm,
234 const Teuchos::RCP<node_type>& pNode,
238 if (pNode.is_null ()) {
239 return rcp (
new map_type (static_cast<global_size_t> (numRows),
240 static_cast<global_ordinal_type> (0),
241 pComm, GloballyDistributed));
244 return rcp (
new map_type (static_cast<global_size_t> (numRows),
245 static_cast<global_ordinal_type> (0),
246 pComm, GloballyDistributed, pNode));
278 static Teuchos::RCP<const map_type>
279 makeRowMap (
const Teuchos::RCP<const map_type>& pRowMap,
280 const Teuchos::RCP<const comm_type>& pComm,
281 const Teuchos::RCP<node_type>& pNode,
286 if (pRowMap.is_null ()) {
287 if (pNode.is_null ()) {
288 return rcp (
new map_type (static_cast<global_size_t> (numRows),
289 static_cast<global_ordinal_type> (0),
290 pComm, GloballyDistributed));
293 return rcp (
new map_type (static_cast<global_size_t> (numRows),
294 static_cast<global_ordinal_type> (0),
295 pComm, GloballyDistributed, pNode));
299 TEUCHOS_TEST_FOR_EXCEPTION
300 (! pRowMap->isDistributed () && pComm->getSize () > 1,
301 std::invalid_argument,
"The specified row map is not distributed, " 302 "but the given communicator includes more than one process (in " 303 "fact, there are " << pComm->getSize () <<
" processes).");
304 TEUCHOS_TEST_FOR_EXCEPTION
305 (pRowMap->getComm () != pComm, std::invalid_argument,
306 "The specified row Map's communicator (pRowMap->getComm()) " 307 "differs from the given separately supplied communicator pComm.");
326 static Teuchos::RCP<const map_type>
327 makeDomainMap (
const Teuchos::RCP<const map_type>& pRangeMap,
336 if (numRows == numCols) {
339 return createUniformContigMapWithNode<LO,GO,NT> (numCols,
340 pRangeMap->getComm (),
341 pRangeMap->getNode ());
418 distribute (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
419 Teuchos::ArrayRCP<size_t>& myRowPtr,
420 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
421 Teuchos::ArrayRCP<scalar_type>& myValues,
422 const Teuchos::RCP<const map_type>& pRowMap,
423 Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
424 Teuchos::ArrayRCP<size_t>& rowPtr,
425 Teuchos::ArrayRCP<global_ordinal_type>& colInd,
426 Teuchos::ArrayRCP<scalar_type>& values,
427 const bool debug=
false)
430 using Teuchos::ArrayRCP;
431 using Teuchos::ArrayView;
434 using Teuchos::CommRequest;
437 using Teuchos::receive;
442 const bool extraDebug =
false;
443 RCP<const comm_type> pComm = pRowMap->getComm ();
444 const int numProcs = pComm->getSize ();
445 const int myRank = pComm->getRank ();
446 const int rootRank = 0;
453 ArrayView<const GO> myRows = pRowMap->getNodeElementList();
454 const size_type myNumRows = myRows.size();
455 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
456 pRowMap->getNodeNumElements(),
458 "pRowMap->getNodeElementList().size() = " 460 <<
" != pRowMap->getNodeNumElements() = " 461 << pRowMap->getNodeNumElements() <<
". " 462 "Please report this bug to the Tpetra developers.");
463 TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
465 "On Proc 0: numEntriesPerRow.size() = " 466 << numEntriesPerRow.size()
467 <<
" != pRowMap->getNodeElementList().size() = " 468 << myNumRows <<
". Please report this bug to the " 469 "Tpetra developers.");
473 myNumEntriesPerRow = arcp<size_t> (myNumRows);
475 if (myRank != rootRank) {
479 send (*pComm, myNumRows, rootRank);
480 if (myNumRows != 0) {
484 send (*pComm, static_cast<int> (myNumRows),
485 myRows.getRawPtr(), rootRank);
495 receive (*pComm, rootRank,
496 static_cast<int> (myNumRows),
497 myNumEntriesPerRow.getRawPtr());
502 std::accumulate (myNumEntriesPerRow.begin(),
503 myNumEntriesPerRow.end(), 0);
509 myColInd = arcp<GO> (myNumEntries);
510 myValues = arcp<scalar_type> (myNumEntries);
511 if (myNumEntries > 0) {
514 receive (*pComm, rootRank,
515 static_cast<int> (myNumEntries),
516 myColInd.getRawPtr());
517 receive (*pComm, rootRank,
518 static_cast<int> (myNumEntries),
519 myValues.getRawPtr());
525 cerr <<
"-- Proc 0: Copying my data from global arrays" << endl;
529 for (size_type k = 0; k < myNumRows; ++k) {
530 const GO myCurRow = myRows[k];
532 myNumEntriesPerRow[k] = numEntriesInThisRow;
534 if (extraDebug && debug) {
535 cerr <<
"Proc " << pRowMap->getComm ()->getRank ()
536 <<
": myNumEntriesPerRow[0.." << (myNumRows-1) <<
"] = [";
537 for (size_type k = 0; k < myNumRows; ++k) {
538 cerr << myNumEntriesPerRow[k];
539 if (k < myNumRows-1) {
547 std::accumulate (myNumEntriesPerRow.begin(),
548 myNumEntriesPerRow.end(), 0);
550 cerr <<
"-- Proc 0: I own " << myNumRows <<
" rows and " 551 << myNumEntries <<
" entries" << endl;
553 myColInd = arcp<GO> (myNumEntries);
554 myValues = arcp<scalar_type> (myNumEntries);
562 for (size_type k = 0; k < myNumRows;
563 myCurPos += myNumEntriesPerRow[k], ++k) {
565 const GO myRow = myRows[k];
566 const size_t curPos = rowPtr[myRow];
569 if (curNumEntries > 0) {
570 ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
571 ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
572 std::copy (colIndView.begin(), colIndView.end(),
573 myColIndView.begin());
575 ArrayView<scalar_type> valuesView =
576 values (curPos, curNumEntries);
577 ArrayView<scalar_type> myValuesView =
578 myValues (myCurPos, curNumEntries);
579 std::copy (valuesView.begin(), valuesView.end(),
580 myValuesView.begin());
585 for (
int p = 1; p < numProcs; ++p) {
587 cerr <<
"-- Proc 0: Processing proc " << p << endl;
590 size_type theirNumRows = 0;
595 receive (*pComm, p, &theirNumRows);
597 cerr <<
"-- Proc 0: Proc " << p <<
" owns " 598 << theirNumRows <<
" rows" << endl;
600 if (theirNumRows != 0) {
605 ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
606 receive (*pComm, p, as<int> (theirNumRows),
607 theirRows.getRawPtr ());
616 const global_size_t numRows = pRowMap->getGlobalNumElements ();
617 const GO indexBase = pRowMap->getIndexBase ();
618 bool theirRowsValid =
true;
619 for (size_type k = 0; k < theirNumRows; ++k) {
620 if (theirRows[k] < indexBase ||
621 as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
622 theirRowsValid =
false;
625 if (! theirRowsValid) {
626 TEUCHOS_TEST_FOR_EXCEPTION(
627 ! theirRowsValid, std::logic_error,
628 "Proc " << p <<
" has at least one invalid row index. " 629 "Here are all of them: " <<
630 Teuchos::toString (theirRows ()) <<
". Valid row index " 631 "range (zero-based): [0, " << (numRows - 1) <<
"].");
646 ArrayRCP<size_t> theirNumEntriesPerRow;
647 theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
648 for (size_type k = 0; k < theirNumRows; ++k) {
649 theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
656 send (*pComm, static_cast<int> (theirNumRows),
657 theirNumEntriesPerRow.getRawPtr(), p);
661 std::accumulate (theirNumEntriesPerRow.begin(),
662 theirNumEntriesPerRow.end(), 0);
665 cerr <<
"-- Proc 0: Proc " << p <<
" owns " 666 << theirNumEntries <<
" entries" << endl;
671 if (theirNumEntries == 0) {
680 ArrayRCP<GO> theirColInd (theirNumEntries);
681 ArrayRCP<scalar_type> theirValues (theirNumEntries);
689 for (size_type k = 0; k < theirNumRows;
690 theirCurPos += theirNumEntriesPerRow[k], k++) {
692 const GO theirRow = theirRows[k];
698 if (curNumEntries > 0) {
699 ArrayView<GO> colIndView =
700 colInd (curPos, curNumEntries);
701 ArrayView<GO> theirColIndView =
702 theirColInd (theirCurPos, curNumEntries);
703 std::copy (colIndView.begin(), colIndView.end(),
704 theirColIndView.begin());
706 ArrayView<scalar_type> valuesView =
707 values (curPos, curNumEntries);
708 ArrayView<scalar_type> theirValuesView =
709 theirValues (theirCurPos, curNumEntries);
710 std::copy (valuesView.begin(), valuesView.end(),
711 theirValuesView.begin());
718 send (*pComm, static_cast<int> (theirNumEntries),
719 theirColInd.getRawPtr(), p);
720 send (*pComm, static_cast<int> (theirNumEntries),
721 theirValues.getRawPtr(), p);
724 cerr <<
"-- Proc 0: Finished with proc " << p << endl;
732 numEntriesPerRow = null;
737 if (debug && myRank == 0) {
738 cerr <<
"-- Proc 0: About to fill in myRowPtr" << endl;
746 myRowPtr = arcp<size_t> (myNumRows+1);
748 for (size_type k = 1; k < myNumRows+1; ++k) {
749 myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
751 if (extraDebug && debug) {
752 cerr <<
"Proc " << Teuchos::rank (*(pRowMap->getComm()))
753 <<
": myRowPtr[0.." << myNumRows <<
"] = [";
754 for (size_type k = 0; k < myNumRows+1; ++k) {
760 cerr <<
"]" << endl << endl;
763 if (debug && myRank == 0) {
764 cerr <<
"-- Proc 0: Done with distribute" << endl;
781 static Teuchos::RCP<sparse_matrix_type>
782 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
783 Teuchos::ArrayRCP<size_t>& myRowPtr,
784 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
785 Teuchos::ArrayRCP<scalar_type>& myValues,
786 const Teuchos::RCP<const map_type>& pRowMap,
787 const Teuchos::RCP<const map_type>& pRangeMap,
788 const Teuchos::RCP<const map_type>& pDomainMap,
789 const bool callFillComplete =
true)
791 using Teuchos::ArrayView;
805 TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
806 "makeMatrix: myRowPtr array is null. " 807 "Please report this bug to the Tpetra developers.");
808 TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
809 "makeMatrix: domain map is null. " 810 "Please report this bug to the Tpetra developers.");
811 TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
812 "makeMatrix: range map is null. " 813 "Please report this bug to the Tpetra developers.");
814 TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
815 "makeMatrix: row map is null. " 816 "Please report this bug to the Tpetra developers.");
823 RCP<sparse_matrix_type> A =
829 ArrayView<const GO> myRows = pRowMap->getNodeElementList ();
830 const size_type myNumRows = myRows.size ();
833 const GO indexBase = pRowMap->getIndexBase ();
834 for (size_type i = 0; i < myNumRows; ++i) {
835 const size_type myCurPos = myRowPtr[i];
837 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
838 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
841 for (size_type k = 0; k < curNumEntries; ++k) {
842 curColInd[k] += indexBase;
845 if (curNumEntries > 0) {
846 A->insertGlobalValues (myRows[i], curColInd, curValues);
853 myNumEntriesPerRow = null;
858 if (callFillComplete) {
859 A->fillComplete (pDomainMap, pRangeMap);
869 static Teuchos::RCP<sparse_matrix_type>
870 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
871 Teuchos::ArrayRCP<size_t>& myRowPtr,
872 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
873 Teuchos::ArrayRCP<scalar_type>& myValues,
874 const Teuchos::RCP<const map_type>& pRowMap,
875 const Teuchos::RCP<const map_type>& pRangeMap,
876 const Teuchos::RCP<const map_type>& pDomainMap,
877 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
878 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams)
880 using Teuchos::ArrayView;
894 TEUCHOS_TEST_FOR_EXCEPTION(
895 myRowPtr.is_null(), std::logic_error,
896 "makeMatrix: myRowPtr array is null. " 897 "Please report this bug to the Tpetra developers.");
898 TEUCHOS_TEST_FOR_EXCEPTION(
899 pDomainMap.is_null(), std::logic_error,
900 "makeMatrix: domain map is null. " 901 "Please report this bug to the Tpetra developers.");
902 TEUCHOS_TEST_FOR_EXCEPTION(
903 pRangeMap.is_null(), std::logic_error,
904 "makeMatrix: range map is null. " 905 "Please report this bug to the Tpetra developers.");
906 TEUCHOS_TEST_FOR_EXCEPTION(
907 pRowMap.is_null(), std::logic_error,
908 "makeMatrix: row map is null. " 909 "Please report this bug to the Tpetra developers.");
916 RCP<sparse_matrix_type> A =
922 ArrayView<const GO> myRows = pRowMap->getNodeElementList();
923 const size_type myNumRows = myRows.size();
926 const GO indexBase = pRowMap->getIndexBase ();
927 for (size_type i = 0; i < myNumRows; ++i) {
928 const size_type myCurPos = myRowPtr[i];
930 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
931 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
934 for (size_type k = 0; k < curNumEntries; ++k) {
935 curColInd[k] += indexBase;
937 if (curNumEntries > 0) {
938 A->insertGlobalValues (myRows[i], curColInd, curValues);
945 myNumEntriesPerRow = null;
950 A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
958 static Teuchos::RCP<sparse_matrix_type>
959 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
960 Teuchos::ArrayRCP<size_t>& myRowPtr,
961 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
962 Teuchos::ArrayRCP<scalar_type>& myValues,
963 const Teuchos::RCP<const map_type>& rowMap,
964 Teuchos::RCP<const map_type>& colMap,
965 const Teuchos::RCP<const map_type>& domainMap,
966 const Teuchos::RCP<const map_type>& rangeMap,
967 const bool callFillComplete =
true)
969 using Teuchos::ArrayView;
975 typedef typename ArrayView<const GO>::size_type size_type;
981 RCP<sparse_matrix_type> A;
982 if (colMap.is_null ()) {
990 ArrayView<const GO> myRows = rowMap->getNodeElementList ();
991 const size_type myNumRows = myRows.size ();
994 const GO indexBase = rowMap->getIndexBase ();
995 for (size_type i = 0; i < myNumRows; ++i) {
996 const size_type myCurPos = myRowPtr[i];
997 const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
998 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
999 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
1002 for (size_type k = 0; k < curNumEntries; ++k) {
1003 curColInd[k] += indexBase;
1005 if (curNumEntries > 0) {
1006 A->insertGlobalValues (myRows[i], curColInd, curValues);
1013 myNumEntriesPerRow = null;
1018 if (callFillComplete) {
1019 A->fillComplete (domainMap, rangeMap);
1020 if (colMap.is_null ()) {
1021 colMap = A->getColMap ();
1045 static Teuchos::RCP<const Teuchos::MatrixMarket::Banner>
1046 readBanner (std::istream& in,
1048 const bool tolerant=
false,
1049 const bool debug=
false,
1050 const bool isGraph=
false)
1052 using Teuchos::MatrixMarket::Banner;
1057 typedef Teuchos::ScalarTraits<scalar_type> STS;
1059 RCP<Banner> pBanner;
1063 const bool readFailed = ! getline(in, line);
1064 TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1065 "Failed to get Matrix Market banner line from input.");
1072 pBanner = rcp (
new Banner (line, tolerant));
1073 }
catch (std::exception& e) {
1074 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1075 "Matrix Market banner line contains syntax error(s): " 1078 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() !=
"matrix",
1079 std::invalid_argument,
"The Matrix Market file does not contain " 1080 "matrix data. Its Banner (first) line says that its object type is \"" 1081 << pBanner->matrixType() <<
"\", rather than the required \"matrix\".");
1085 TEUCHOS_TEST_FOR_EXCEPTION(
1086 ! STS::isComplex && pBanner->dataType() ==
"complex",
1087 std::invalid_argument,
1088 "The Matrix Market file contains complex-valued data, but you are " 1089 "trying to read it into a matrix containing entries of the real-" 1090 "valued Scalar type \"" 1091 << Teuchos::TypeNameTraits<scalar_type>::name() <<
"\".");
1092 TEUCHOS_TEST_FOR_EXCEPTION(
1094 pBanner->dataType() !=
"real" &&
1095 pBanner->dataType() !=
"complex" &&
1096 pBanner->dataType() !=
"integer",
1097 std::invalid_argument,
1098 "When reading Matrix Market data into a Tpetra::CrsMatrix, the " 1099 "Matrix Market file may not contain a \"pattern\" matrix. A " 1100 "pattern matrix is really just a graph with no weights. It " 1101 "should be stored in a CrsGraph, not a CrsMatrix.");
1103 TEUCHOS_TEST_FOR_EXCEPTION(
1105 pBanner->dataType() !=
"pattern",
1106 std::invalid_argument,
1107 "When reading Matrix Market data into a Tpetra::CrsGraph, the " 1108 "Matrix Market file must contain a \"pattern\" matrix.");
1135 static Teuchos::Tuple<global_ordinal_type, 3>
1136 readCoordDims (std::istream& in,
1138 const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1139 const Teuchos::RCP<const comm_type>& pComm,
1140 const bool tolerant =
false,
1141 const bool debug =
false)
1143 using Teuchos::MatrixMarket::readCoordinateDimensions;
1144 using Teuchos::Tuple;
1149 Tuple<global_ordinal_type, 3> dims;
1155 bool success =
false;
1156 if (pComm->getRank() == 0) {
1157 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() !=
"coordinate",
1158 std::invalid_argument,
"The Tpetra::CrsMatrix Matrix Market reader " 1159 "only accepts \"coordinate\" (sparse) matrix data.");
1163 success = readCoordinateDimensions (in, numRows, numCols,
1164 numNonzeros, lineNumber,
1170 dims[2] = numNonzeros;
1178 int the_success = success ? 1 : 0;
1179 Teuchos::broadcast (*pComm, 0, &the_success);
1180 success = (the_success == 1);
1185 Teuchos::broadcast (*pComm, 0, dims);
1193 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1194 "Error reading Matrix Market sparse matrix: failed to read " 1195 "coordinate matrix dimensions.");
1210 typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1212 typedef Teuchos::MatrixMarket::SymmetrizingGraphAdder<Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> > graph_adder_type;
1239 static Teuchos::RCP<adder_type>
1240 makeAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1241 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1242 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1243 const bool tolerant=
false,
1244 const bool debug=
false)
1246 if (pComm->getRank () == 0) {
1247 typedef Teuchos::MatrixMarket::Raw::Adder<
scalar_type,
1250 Teuchos::RCP<raw_adder_type> pRaw =
1251 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1253 return Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
1256 return Teuchos::null;
1285 static Teuchos::RCP<graph_adder_type>
1286 makeGraphAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1287 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1288 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1289 const bool tolerant=
false,
1290 const bool debug=
false)
1292 if (pComm->getRank () == 0) {
1293 typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
1294 Teuchos::RCP<raw_adder_type> pRaw =
1295 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1297 return Teuchos::rcp (
new graph_adder_type (pRaw, pBanner->symmType ()));
1300 return Teuchos::null;
1305 static Teuchos::RCP<sparse_graph_type>
1306 readSparseGraphHelper (std::istream& in,
1307 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1308 const Teuchos::RCP<node_type>& pNode,
1309 const Teuchos::RCP<const map_type>& rowMap,
1310 Teuchos::RCP<const map_type>& colMap,
1311 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1312 const bool tolerant=
false,
1313 const bool debug=
false)
1315 using Teuchos::MatrixMarket::Banner;
1318 using Teuchos::Tuple;
1322 const int myRank = pComm->getRank ();
1323 const int rootRank = 0;
1328 size_t lineNumber = 1;
1330 if (debug && myRank == rootRank) {
1331 cerr <<
"Matrix Market reader: readGraph:" << endl
1332 <<
"-- Reading banner line" << endl;
1341 RCP<const Banner> pBanner;
1347 int bannerIsCorrect = 1;
1348 std::ostringstream errMsg;
1350 if (myRank == rootRank) {
1353 pBanner = readBanner (in, lineNumber, tolerant, debug,
true);
1355 catch (std::exception& e) {
1356 errMsg <<
"Attempt to read the Matrix Market file's Banner line " 1357 "threw an exception: " << e.what();
1358 bannerIsCorrect = 0;
1361 if (bannerIsCorrect) {
1366 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
1367 bannerIsCorrect = 0;
1368 errMsg <<
"The Matrix Market input file must contain a " 1369 "\"coordinate\"-format sparse graph in order to create a " 1370 "Tpetra::CrsGraph object from it, but the file's matrix " 1371 "type is \"" << pBanner->matrixType() <<
"\" instead.";
1376 if (tolerant && pBanner->matrixType() ==
"array") {
1377 bannerIsCorrect = 0;
1378 errMsg <<
"Matrix Market file must contain a \"coordinate\"-" 1379 "format sparse graph in order to create a Tpetra::CrsGraph " 1380 "object from it, but the file's matrix type is \"array\" " 1381 "instead. That probably means the file contains dense matrix " 1388 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1395 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1396 std::invalid_argument, errMsg.str ());
1398 if (debug && myRank == rootRank) {
1399 cerr <<
"-- Reading dimensions line" << endl;
1407 Tuple<global_ordinal_type, 3> dims =
1408 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1410 if (debug && myRank == rootRank) {
1411 cerr <<
"-- Making Adder for collecting graph data" << endl;
1418 RCP<graph_adder_type> pAdder =
1419 makeGraphAdder (pComm, pBanner, dims, tolerant, debug);
1421 if (debug && myRank == rootRank) {
1422 cerr <<
"-- Reading graph data" << endl;
1432 int readSuccess = 1;
1433 std::ostringstream errMsg;
1434 if (myRank == rootRank) {
1437 typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
1439 reader_type reader (pAdder);
1442 std::pair<bool, std::vector<size_t> > results =
1443 reader.read (in, lineNumber, tolerant, debug);
1444 readSuccess = results.first ? 1 : 0;
1446 catch (std::exception& e) {
1451 broadcast (*pComm, rootRank, ptr (&readSuccess));
1460 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1461 "Failed to read the Matrix Market sparse graph file: " 1465 if (debug && myRank == rootRank) {
1466 cerr <<
"-- Successfully read the Matrix Market data" << endl;
1479 if (debug && myRank == rootRank) {
1480 cerr <<
"-- Tolerant mode: rebroadcasting graph dimensions" 1482 <<
"----- Dimensions before: " 1483 << dims[0] <<
" x " << dims[1]
1487 Tuple<global_ordinal_type, 2> updatedDims;
1488 if (myRank == rootRank) {
1495 std::max (dims[0], pAdder->getAdder()->numRows());
1496 updatedDims[1] = pAdder->getAdder()->numCols();
1498 broadcast (*pComm, rootRank, updatedDims);
1499 dims[0] = updatedDims[0];
1500 dims[1] = updatedDims[1];
1501 if (debug && myRank == rootRank) {
1502 cerr <<
"----- Dimensions after: " << dims[0] <<
" x " 1516 if (myRank == rootRank) {
1523 if (dims[0] < pAdder->getAdder ()->numRows ()) {
1527 broadcast (*pComm, 0, ptr (&dimsMatch));
1528 if (dimsMatch == 0) {
1535 Tuple<global_ordinal_type, 2> addersDims;
1536 if (myRank == rootRank) {
1537 addersDims[0] = pAdder->getAdder()->numRows();
1538 addersDims[1] = pAdder->getAdder()->numCols();
1540 broadcast (*pComm, 0, addersDims);
1541 TEUCHOS_TEST_FOR_EXCEPTION(
1542 dimsMatch == 0, std::runtime_error,
1543 "The graph metadata says that the graph is " << dims[0] <<
" x " 1544 << dims[1] <<
", but the actual data says that the graph is " 1545 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the " 1546 "data includes more rows than reported in the metadata. This " 1547 "is not allowed when parsing in strict mode. Parse the graph in " 1548 "tolerant mode to ignore the metadata when it disagrees with the " 1554 RCP<map_type> proc0Map;
1556 if(Teuchos::is_null(rowMap)) {
1560 indexBase = rowMap->getIndexBase();
1562 if(myRank == rootRank) {
1563 proc0Map = rcp(
new map_type(dims[0],dims[0],indexBase,pComm,pNode));
1566 proc0Map = rcp(
new map_type(dims[0],0,indexBase,pComm,pNode));
1570 RCP<sparse_graph_type> proc0Graph =
1572 if(myRank == rootRank) {
1573 typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
1576 const std::vector<element_type>& entries =
1577 pAdder->getAdder()->getEntries();
1580 for(
size_t curPos=0; curPos<entries.size(); curPos++) {
1581 const element_type& curEntry = entries[curPos];
1584 Teuchos::ArrayView<const global_ordinal_type> colView(&curCol,1);
1585 proc0Graph->insertGlobalIndices(curRow,colView);
1588 proc0Graph->fillComplete();
1590 RCP<sparse_graph_type> distGraph;
1591 if(Teuchos::is_null(rowMap))
1594 RCP<map_type> distMap =
1595 rcp(
new map_type(dims[0],0,pComm,GloballyDistributed,pNode));
1601 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1602 import_type importer (proc0Map, distMap);
1605 distGraph->doImport(*proc0Graph,importer,
INSERT);
1611 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1612 import_type importer (proc0Map, rowMap);
1615 distGraph->doImport(*proc0Graph,importer,
INSERT);
1645 static Teuchos::RCP<sparse_graph_type>
1647 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1648 const bool callFillComplete=
true,
1649 const bool tolerant=
false,
1650 const bool debug=
false)
1652 return readSparseGraph (filename, pComm, Teuchos::null, callFillComplete, tolerant, debug);
1656 static Teuchos::RCP<sparse_graph_type>
1658 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1659 const Teuchos::RCP<node_type>& pNode,
1660 const bool callFillComplete=
true,
1661 const bool tolerant=
false,
1662 const bool debug=
false)
1664 const int myRank = pComm->getRank ();
1669 in.open (filename.c_str ());
1677 return readSparseGraph (in, pComm, pNode, callFillComplete, tolerant, debug);
1711 static Teuchos::RCP<sparse_graph_type>
1713 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1714 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1715 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1716 const bool tolerant=
false,
1717 const bool debug=
false)
1720 constructorParams, fillCompleteParams,
1725 static Teuchos::RCP<sparse_graph_type>
1727 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1728 const Teuchos::RCP<node_type>& pNode,
1729 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1730 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1731 const bool tolerant=
false,
1732 const bool debug=
false)
1735 if (pComm->getRank () == 0) {
1736 in.open (filename.c_str ());
1739 fillCompleteParams, tolerant, debug);
1779 static Teuchos::RCP<sparse_graph_type>
1781 const Teuchos::RCP<const map_type>& rowMap,
1782 Teuchos::RCP<const map_type>& colMap,
1783 const Teuchos::RCP<const map_type>& domainMap,
1784 const Teuchos::RCP<const map_type>& rangeMap,
1785 const bool callFillComplete=
true,
1786 const bool tolerant=
false,
1787 const bool debug=
false)
1789 using Teuchos::broadcast;
1790 using Teuchos::Comm;
1791 using Teuchos::outArg;
1794 TEUCHOS_TEST_FOR_EXCEPTION(
1795 rowMap.is_null (), std::invalid_argument,
1796 "Row Map must be nonnull.");
1798 RCP<const Comm<int> > comm = rowMap->getComm ();
1799 const int myRank = comm->getRank ();
1809 in.open (filename.c_str ());
1816 broadcast<int, int> (*comm, 0, outArg (opened));
1817 TEUCHOS_TEST_FOR_EXCEPTION(
1818 opened == 0, std::runtime_error,
1819 "readSparseGraph: Failed to open file \"" << filename <<
"\" on " 1822 callFillComplete, tolerant, debug);
1850 static Teuchos::RCP<sparse_graph_type>
1852 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1853 const bool callFillComplete=
true,
1854 const bool tolerant=
false,
1855 const bool debug=
false)
1857 return readSparseGraph (in, pComm, Teuchos::null, callFillComplete, tolerant, debug);
1861 static Teuchos::RCP<sparse_graph_type>
1863 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1864 const Teuchos::RCP<node_type>& pNode,
1865 const bool callFillComplete=
true,
1866 const bool tolerant=
false,
1867 const bool debug=
false)
1869 Teuchos::RCP<sparse_graph_type> graph = readSparseGraphHelper(in, pComm, pNode, Teuchos::null, Teuchos::null, Teuchos::null,tolerant,debug);
1870 if(callFillComplete)
1871 graph->FillComplete();
1903 static Teuchos::RCP<sparse_graph_type>
1905 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1906 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1907 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1908 const bool tolerant=
false,
1909 const bool debug=
false)
1912 fillCompleteParams, tolerant, debug);
1916 static Teuchos::RCP<sparse_graph_type>
1918 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1919 const Teuchos::RCP<node_type>& pNode,
1920 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1921 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1922 const bool tolerant=
false,
1923 const bool debug=
false)
1925 Teuchos::RCP<sparse_graph_type> graph = readSparseGraphHelper(in, pComm, pNode,
1926 Teuchos::null, Teuchos::null, constructorParams, tolerant, debug);
1927 graph->FillComplete(fillCompleteParams);
1971 static Teuchos::RCP<sparse_graph_type>
1973 const Teuchos::RCP<const map_type>& rowMap,
1974 Teuchos::RCP<const map_type>& colMap,
1975 const Teuchos::RCP<const map_type>& domainMap,
1976 const Teuchos::RCP<const map_type>& rangeMap,
1977 const bool callFillComplete=
true,
1978 const bool tolerant=
false,
1979 const bool debug=
false)
1981 Teuchos::RCP<sparse_graph_type> graph = readSparseGraphHelper(in, rowMap->getComm(),
1982 rowMap->getNode(), rowMap, colMap, Teuchos::null, tolerant, debug);
1983 if(callFillComplete)
1984 graph->fillComplete();
2011 static Teuchos::RCP<sparse_matrix_type>
2013 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2014 const bool callFillComplete=
true,
2015 const bool tolerant=
false,
2016 const bool debug=
false)
2018 return readSparseFile (filename, pComm, Teuchos::null, callFillComplete, tolerant, debug);
2022 static Teuchos::RCP<sparse_matrix_type>
2024 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2025 const Teuchos::RCP<node_type>& pNode,
2026 const bool callFillComplete=
true,
2027 const bool tolerant=
false,
2028 const bool debug=
false)
2030 const int myRank = pComm->getRank ();
2035 in.open (filename.c_str ());
2043 return readSparse (in, pComm, pNode, callFillComplete, tolerant, debug);
2077 static Teuchos::RCP<sparse_matrix_type>
2079 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2080 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2081 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2082 const bool tolerant=
false,
2083 const bool debug=
false)
2086 constructorParams, fillCompleteParams,
2091 static Teuchos::RCP<sparse_matrix_type>
2093 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2094 const Teuchos::RCP<node_type>& pNode,
2095 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2096 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2097 const bool tolerant=
false,
2098 const bool debug=
false)
2101 if (pComm->getRank () == 0) {
2102 in.open (filename.c_str ());
2104 return readSparse (in, pComm, pNode, constructorParams,
2105 fillCompleteParams, tolerant, debug);
2145 static Teuchos::RCP<sparse_matrix_type>
2147 const Teuchos::RCP<const map_type>& rowMap,
2148 Teuchos::RCP<const map_type>& colMap,
2149 const Teuchos::RCP<const map_type>& domainMap,
2150 const Teuchos::RCP<const map_type>& rangeMap,
2151 const bool callFillComplete=
true,
2152 const bool tolerant=
false,
2153 const bool debug=
false)
2155 using Teuchos::broadcast;
2156 using Teuchos::Comm;
2157 using Teuchos::outArg;
2160 TEUCHOS_TEST_FOR_EXCEPTION(
2161 rowMap.is_null (), std::invalid_argument,
2162 "Row Map must be nonnull.");
2164 RCP<const Comm<int> > comm = rowMap->getComm ();
2165 const int myRank = comm->getRank ();
2175 in.open (filename.c_str ());
2182 broadcast<int, int> (*comm, 0, outArg (opened));
2183 TEUCHOS_TEST_FOR_EXCEPTION(
2184 opened == 0, std::runtime_error,
2185 "readSparseFile: Failed to open file \"" << filename <<
"\" on " 2187 return readSparse (in, rowMap, colMap, domainMap, rangeMap,
2188 callFillComplete, tolerant, debug);
2216 static Teuchos::RCP<sparse_matrix_type>
2218 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2219 const bool callFillComplete=
true,
2220 const bool tolerant=
false,
2221 const bool debug=
false)
2223 return readSparse (in, pComm, Teuchos::null, callFillComplete, tolerant, debug);
2227 static Teuchos::RCP<sparse_matrix_type>
2229 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2230 const Teuchos::RCP<node_type>& pNode,
2231 const bool callFillComplete=
true,
2232 const bool tolerant=
false,
2233 const bool debug=
false)
2235 using Teuchos::MatrixMarket::Banner;
2236 using Teuchos::arcp;
2237 using Teuchos::ArrayRCP;
2238 using Teuchos::broadcast;
2239 using Teuchos::null;
2242 using Teuchos::REDUCE_MAX;
2243 using Teuchos::reduceAll;
2244 using Teuchos::Tuple;
2247 typedef Teuchos::ScalarTraits<scalar_type> STS;
2249 const bool extraDebug =
false;
2250 const int myRank = pComm->getRank ();
2251 const int rootRank = 0;
2256 size_t lineNumber = 1;
2258 if (debug && myRank == rootRank) {
2259 cerr <<
"Matrix Market reader: readSparse:" << endl
2260 <<
"-- Reading banner line" << endl;
2269 RCP<const Banner> pBanner;
2275 int bannerIsCorrect = 1;
2276 std::ostringstream errMsg;
2278 if (myRank == rootRank) {
2281 pBanner = readBanner (in, lineNumber, tolerant, debug);
2283 catch (std::exception& e) {
2284 errMsg <<
"Attempt to read the Matrix Market file's Banner line " 2285 "threw an exception: " << e.what();
2286 bannerIsCorrect = 0;
2289 if (bannerIsCorrect) {
2294 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2295 bannerIsCorrect = 0;
2296 errMsg <<
"The Matrix Market input file must contain a " 2297 "\"coordinate\"-format sparse matrix in order to create a " 2298 "Tpetra::CrsMatrix object from it, but the file's matrix " 2299 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2304 if (tolerant && pBanner->matrixType() ==
"array") {
2305 bannerIsCorrect = 0;
2306 errMsg <<
"Matrix Market file must contain a \"coordinate\"-" 2307 "format sparse matrix in order to create a Tpetra::CrsMatrix " 2308 "object from it, but the file's matrix type is \"array\" " 2309 "instead. That probably means the file contains dense matrix " 2316 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2323 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2324 std::invalid_argument, errMsg.str ());
2326 if (debug && myRank == rootRank) {
2327 cerr <<
"-- Reading dimensions line" << endl;
2335 Tuple<global_ordinal_type, 3> dims =
2336 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2338 if (debug && myRank == rootRank) {
2339 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2344 RCP<adder_type> pAdder =
2345 makeAdder (pComm, pBanner, dims, tolerant, debug);
2347 if (debug && myRank == rootRank) {
2348 cerr <<
"-- Reading matrix data" << endl;
2358 int readSuccess = 1;
2359 std::ostringstream errMsg;
2360 if (myRank == rootRank) {
2363 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2365 reader_type reader (pAdder);
2368 std::pair<bool, std::vector<size_t> > results =
2369 reader.read (in, lineNumber, tolerant, debug);
2370 readSuccess = results.first ? 1 : 0;
2372 catch (std::exception& e) {
2377 broadcast (*pComm, rootRank, ptr (&readSuccess));
2386 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2387 "Failed to read the Matrix Market sparse matrix file: " 2391 if (debug && myRank == rootRank) {
2392 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2405 if (debug && myRank == rootRank) {
2406 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions" 2408 <<
"----- Dimensions before: " 2409 << dims[0] <<
" x " << dims[1]
2413 Tuple<global_ordinal_type, 2> updatedDims;
2414 if (myRank == rootRank) {
2421 std::max (dims[0], pAdder->getAdder()->numRows());
2422 updatedDims[1] = pAdder->getAdder()->numCols();
2424 broadcast (*pComm, rootRank, updatedDims);
2425 dims[0] = updatedDims[0];
2426 dims[1] = updatedDims[1];
2427 if (debug && myRank == rootRank) {
2428 cerr <<
"----- Dimensions after: " << dims[0] <<
" x " 2442 if (myRank == rootRank) {
2449 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2453 broadcast (*pComm, 0, ptr (&dimsMatch));
2454 if (dimsMatch == 0) {
2461 Tuple<global_ordinal_type, 2> addersDims;
2462 if (myRank == rootRank) {
2463 addersDims[0] = pAdder->getAdder()->numRows();
2464 addersDims[1] = pAdder->getAdder()->numCols();
2466 broadcast (*pComm, 0, addersDims);
2467 TEUCHOS_TEST_FOR_EXCEPTION(
2468 dimsMatch == 0, std::runtime_error,
2469 "The matrix metadata says that the matrix is " << dims[0] <<
" x " 2470 << dims[1] <<
", but the actual data says that the matrix is " 2471 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the " 2472 "data includes more rows than reported in the metadata. This " 2473 "is not allowed when parsing in strict mode. Parse the matrix in " 2474 "tolerant mode to ignore the metadata when it disagrees with the " 2479 if (debug && myRank == rootRank) {
2480 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
2492 ArrayRCP<size_t> numEntriesPerRow;
2493 ArrayRCP<size_t> rowPtr;
2494 ArrayRCP<global_ordinal_type> colInd;
2495 ArrayRCP<scalar_type> values;
2500 int mergeAndConvertSucceeded = 1;
2501 std::ostringstream errMsg;
2503 if (myRank == rootRank) {
2505 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2515 const size_type numRows = dims[0];
2518 pAdder->getAdder()->merge ();
2521 const std::vector<element_type>& entries =
2522 pAdder->getAdder()->getEntries();
2525 const size_t numEntries = (size_t)entries.size();
2528 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
2529 <<
" rows and numEntries=" << numEntries
2530 <<
" entries." << endl;
2536 numEntriesPerRow = arcp<size_t> (numRows);
2537 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2538 rowPtr = arcp<size_t> (numRows+1);
2539 std::fill (rowPtr.begin(), rowPtr.end(), 0);
2540 colInd = arcp<global_ordinal_type> (numEntries);
2541 values = arcp<scalar_type> (numEntries);
2548 for (curPos = 0; curPos < numEntries; ++curPos) {
2549 const element_type& curEntry = entries[curPos];
2551 TEUCHOS_TEST_FOR_EXCEPTION(
2552 curRow < prvRow, std::logic_error,
2553 "Row indices are out of order, even though they are supposed " 2554 "to be sorted. curRow = " << curRow <<
", prvRow = " 2555 << prvRow <<
", at curPos = " << curPos <<
". Please report " 2556 "this bug to the Tpetra developers.");
2557 if (curRow > prvRow) {
2563 numEntriesPerRow[curRow]++;
2564 colInd[curPos] = curEntry.colIndex();
2565 values[curPos] = curEntry.value();
2570 rowPtr[numRows] = numEntries;
2572 catch (std::exception& e) {
2573 mergeAndConvertSucceeded = 0;
2574 errMsg <<
"Failed to merge sparse matrix entries and convert to " 2575 "CSR format: " << e.what();
2578 if (debug && mergeAndConvertSucceeded) {
2580 const size_type numRows = dims[0];
2581 const size_type maxToDisplay = 100;
2583 cerr <<
"----- Proc 0: numEntriesPerRow[0.." 2584 << (numEntriesPerRow.size()-1) <<
"] ";
2585 if (numRows > maxToDisplay) {
2586 cerr <<
"(only showing first and last few entries) ";
2590 if (numRows > maxToDisplay) {
2591 for (size_type k = 0; k < 2; ++k) {
2592 cerr << numEntriesPerRow[k] <<
" ";
2595 for (size_type k = numRows-2; k < numRows-1; ++k) {
2596 cerr << numEntriesPerRow[k] <<
" ";
2600 for (size_type k = 0; k < numRows-1; ++k) {
2601 cerr << numEntriesPerRow[k] <<
" ";
2604 cerr << numEntriesPerRow[numRows-1];
2606 cerr <<
"]" << endl;
2608 cerr <<
"----- Proc 0: rowPtr ";
2609 if (numRows > maxToDisplay) {
2610 cerr <<
"(only showing first and last few entries) ";
2613 if (numRows > maxToDisplay) {
2614 for (size_type k = 0; k < 2; ++k) {
2615 cerr << rowPtr[k] <<
" ";
2618 for (size_type k = numRows-2; k < numRows; ++k) {
2619 cerr << rowPtr[k] <<
" ";
2623 for (size_type k = 0; k < numRows; ++k) {
2624 cerr << rowPtr[k] <<
" ";
2627 cerr << rowPtr[numRows] <<
"]" << endl;
2638 if (debug && myRank == rootRank) {
2639 cerr <<
"-- Making range, domain, and row maps" << endl;
2646 RCP<const map_type> pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
2647 RCP<const map_type> pDomainMap =
2648 makeDomainMap (pRangeMap, dims[0], dims[1]);
2649 RCP<const map_type> pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
2651 if (debug && myRank == rootRank) {
2652 cerr <<
"-- Distributing the matrix data" << endl;
2666 ArrayRCP<size_t> myNumEntriesPerRow;
2667 ArrayRCP<size_t> myRowPtr;
2668 ArrayRCP<global_ordinal_type> myColInd;
2669 ArrayRCP<scalar_type> myValues;
2671 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
2672 numEntriesPerRow, rowPtr, colInd, values, debug);
2674 if (debug && myRank == rootRank) {
2675 cerr <<
"-- Inserting matrix entries on each processor";
2676 if (callFillComplete) {
2677 cerr <<
" and calling fillComplete()";
2688 RCP<sparse_matrix_type> pMatrix =
2689 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
2690 pRowMap, pRangeMap, pDomainMap, callFillComplete);
2695 int localIsNull = pMatrix.is_null () ? 1 : 0;
2696 int globalIsNull = 0;
2697 reduceAll (*pComm, REDUCE_MAX, localIsNull, ptr (&globalIsNull));
2698 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
2699 "Reader::makeMatrix() returned a null pointer on at least one " 2700 "process. Please report this bug to the Tpetra developers.");
2703 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
2704 "Reader::makeMatrix() returned a null pointer. " 2705 "Please report this bug to the Tpetra developers.");
2717 if (callFillComplete) {
2718 const int numProcs = pComm->getSize ();
2720 if (extraDebug && debug) {
2722 pRangeMap->getGlobalNumElements ();
2724 pDomainMap->getGlobalNumElements ();
2725 if (myRank == rootRank) {
2726 cerr <<
"-- Matrix is " 2727 << globalNumRows <<
" x " << globalNumCols
2728 <<
" with " << pMatrix->getGlobalNumEntries()
2729 <<
" entries, and index base " 2730 << pMatrix->getIndexBase() <<
"." << endl;
2733 for (
int p = 0; p < numProcs; ++p) {
2735 cerr <<
"-- Proc " << p <<
" owns " 2736 << pMatrix->getNodeNumCols() <<
" columns, and " 2737 << pMatrix->getNodeNumEntries() <<
" entries." << endl;
2744 if (debug && myRank == rootRank) {
2745 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data" 2780 static Teuchos::RCP<sparse_matrix_type>
2782 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2783 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2784 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2785 const bool tolerant=
false,
2786 const bool debug=
false)
2788 return readSparse (in, pComm, Teuchos::null, constructorParams,
2789 fillCompleteParams, tolerant, debug);
2793 static Teuchos::RCP<sparse_matrix_type>
2795 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2796 const Teuchos::RCP<node_type>& pNode,
2797 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2798 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2799 const bool tolerant=
false,
2800 const bool debug=
false)
2802 using Teuchos::MatrixMarket::Banner;
2803 using Teuchos::arcp;
2804 using Teuchos::ArrayRCP;
2805 using Teuchos::broadcast;
2806 using Teuchos::null;
2809 using Teuchos::reduceAll;
2810 using Teuchos::Tuple;
2813 typedef Teuchos::ScalarTraits<scalar_type> STS;
2815 const bool extraDebug =
false;
2816 const int myRank = pComm->getRank ();
2817 const int rootRank = 0;
2822 size_t lineNumber = 1;
2824 if (debug && myRank == rootRank) {
2825 cerr <<
"Matrix Market reader: readSparse:" << endl
2826 <<
"-- Reading banner line" << endl;
2835 RCP<const Banner> pBanner;
2841 int bannerIsCorrect = 1;
2842 std::ostringstream errMsg;
2844 if (myRank == rootRank) {
2847 pBanner = readBanner (in, lineNumber, tolerant, debug);
2849 catch (std::exception& e) {
2850 errMsg <<
"Attempt to read the Matrix Market file's Banner line " 2851 "threw an exception: " << e.what();
2852 bannerIsCorrect = 0;
2855 if (bannerIsCorrect) {
2860 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2861 bannerIsCorrect = 0;
2862 errMsg <<
"The Matrix Market input file must contain a " 2863 "\"coordinate\"-format sparse matrix in order to create a " 2864 "Tpetra::CrsMatrix object from it, but the file's matrix " 2865 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2870 if (tolerant && pBanner->matrixType() ==
"array") {
2871 bannerIsCorrect = 0;
2872 errMsg <<
"Matrix Market file must contain a \"coordinate\"-" 2873 "format sparse matrix in order to create a Tpetra::CrsMatrix " 2874 "object from it, but the file's matrix type is \"array\" " 2875 "instead. That probably means the file contains dense matrix " 2882 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2889 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2890 std::invalid_argument, errMsg.str ());
2892 if (debug && myRank == rootRank) {
2893 cerr <<
"-- Reading dimensions line" << endl;
2901 Tuple<global_ordinal_type, 3> dims =
2902 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2904 if (debug && myRank == rootRank) {
2905 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2910 RCP<adder_type> pAdder =
2911 makeAdder (pComm, pBanner, dims, tolerant, debug);
2913 if (debug && myRank == rootRank) {
2914 cerr <<
"-- Reading matrix data" << endl;
2924 int readSuccess = 1;
2925 std::ostringstream errMsg;
2926 if (myRank == rootRank) {
2929 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2931 reader_type reader (pAdder);
2934 std::pair<bool, std::vector<size_t> > results =
2935 reader.read (in, lineNumber, tolerant, debug);
2936 readSuccess = results.first ? 1 : 0;
2938 catch (std::exception& e) {
2943 broadcast (*pComm, rootRank, ptr (&readSuccess));
2952 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2953 "Failed to read the Matrix Market sparse matrix file: " 2957 if (debug && myRank == rootRank) {
2958 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2971 if (debug && myRank == rootRank) {
2972 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions" 2974 <<
"----- Dimensions before: " 2975 << dims[0] <<
" x " << dims[1]
2979 Tuple<global_ordinal_type, 2> updatedDims;
2980 if (myRank == rootRank) {
2987 std::max (dims[0], pAdder->getAdder()->numRows());
2988 updatedDims[1] = pAdder->getAdder()->numCols();
2990 broadcast (*pComm, rootRank, updatedDims);
2991 dims[0] = updatedDims[0];
2992 dims[1] = updatedDims[1];
2993 if (debug && myRank == rootRank) {
2994 cerr <<
"----- Dimensions after: " << dims[0] <<
" x " 3008 if (myRank == rootRank) {
3015 if (dims[0] < pAdder->getAdder ()->numRows ()) {
3019 broadcast (*pComm, 0, ptr (&dimsMatch));
3020 if (dimsMatch == 0) {
3027 Tuple<global_ordinal_type, 2> addersDims;
3028 if (myRank == rootRank) {
3029 addersDims[0] = pAdder->getAdder()->numRows();
3030 addersDims[1] = pAdder->getAdder()->numCols();
3032 broadcast (*pComm, 0, addersDims);
3033 TEUCHOS_TEST_FOR_EXCEPTION(
3034 dimsMatch == 0, std::runtime_error,
3035 "The matrix metadata says that the matrix is " << dims[0] <<
" x " 3036 << dims[1] <<
", but the actual data says that the matrix is " 3037 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the " 3038 "data includes more rows than reported in the metadata. This " 3039 "is not allowed when parsing in strict mode. Parse the matrix in " 3040 "tolerant mode to ignore the metadata when it disagrees with the " 3045 if (debug && myRank == rootRank) {
3046 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
3058 ArrayRCP<size_t> numEntriesPerRow;
3059 ArrayRCP<size_t> rowPtr;
3060 ArrayRCP<global_ordinal_type> colInd;
3061 ArrayRCP<scalar_type> values;
3066 int mergeAndConvertSucceeded = 1;
3067 std::ostringstream errMsg;
3069 if (myRank == rootRank) {
3071 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3081 const size_type numRows = dims[0];
3084 pAdder->getAdder()->merge ();
3087 const std::vector<element_type>& entries =
3088 pAdder->getAdder()->getEntries();
3091 const size_t numEntries = (size_t)entries.size();
3094 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
3095 <<
" rows and numEntries=" << numEntries
3096 <<
" entries." << endl;
3102 numEntriesPerRow = arcp<size_t> (numRows);
3103 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3104 rowPtr = arcp<size_t> (numRows+1);
3105 std::fill (rowPtr.begin(), rowPtr.end(), 0);
3106 colInd = arcp<global_ordinal_type> (numEntries);
3107 values = arcp<scalar_type> (numEntries);
3114 for (curPos = 0; curPos < numEntries; ++curPos) {
3115 const element_type& curEntry = entries[curPos];
3117 TEUCHOS_TEST_FOR_EXCEPTION(
3118 curRow < prvRow, std::logic_error,
3119 "Row indices are out of order, even though they are supposed " 3120 "to be sorted. curRow = " << curRow <<
", prvRow = " 3121 << prvRow <<
", at curPos = " << curPos <<
". Please report " 3122 "this bug to the Tpetra developers.");
3123 if (curRow > prvRow) {
3129 numEntriesPerRow[curRow]++;
3130 colInd[curPos] = curEntry.colIndex();
3131 values[curPos] = curEntry.value();
3136 rowPtr[numRows] = numEntries;
3138 catch (std::exception& e) {
3139 mergeAndConvertSucceeded = 0;
3140 errMsg <<
"Failed to merge sparse matrix entries and convert to " 3141 "CSR format: " << e.what();
3144 if (debug && mergeAndConvertSucceeded) {
3146 const size_type numRows = dims[0];
3147 const size_type maxToDisplay = 100;
3149 cerr <<
"----- Proc 0: numEntriesPerRow[0.." 3150 << (numEntriesPerRow.size()-1) <<
"] ";
3151 if (numRows > maxToDisplay) {
3152 cerr <<
"(only showing first and last few entries) ";
3156 if (numRows > maxToDisplay) {
3157 for (size_type k = 0; k < 2; ++k) {
3158 cerr << numEntriesPerRow[k] <<
" ";
3161 for (size_type k = numRows-2; k < numRows-1; ++k) {
3162 cerr << numEntriesPerRow[k] <<
" ";
3166 for (size_type k = 0; k < numRows-1; ++k) {
3167 cerr << numEntriesPerRow[k] <<
" ";
3170 cerr << numEntriesPerRow[numRows-1];
3172 cerr <<
"]" << endl;
3174 cerr <<
"----- Proc 0: rowPtr ";
3175 if (numRows > maxToDisplay) {
3176 cerr <<
"(only showing first and last few entries) ";
3179 if (numRows > maxToDisplay) {
3180 for (size_type k = 0; k < 2; ++k) {
3181 cerr << rowPtr[k] <<
" ";
3184 for (size_type k = numRows-2; k < numRows; ++k) {
3185 cerr << rowPtr[k] <<
" ";
3189 for (size_type k = 0; k < numRows; ++k) {
3190 cerr << rowPtr[k] <<
" ";
3193 cerr << rowPtr[numRows] <<
"]" << endl;
3204 if (debug && myRank == rootRank) {
3205 cerr <<
"-- Making range, domain, and row maps" << endl;
3212 RCP<const map_type> pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
3213 RCP<const map_type> pDomainMap =
3214 makeDomainMap (pRangeMap, dims[0], dims[1]);
3215 RCP<const map_type> pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
3217 if (debug && myRank == rootRank) {
3218 cerr <<
"-- Distributing the matrix data" << endl;
3232 ArrayRCP<size_t> myNumEntriesPerRow;
3233 ArrayRCP<size_t> myRowPtr;
3234 ArrayRCP<global_ordinal_type> myColInd;
3235 ArrayRCP<scalar_type> myValues;
3237 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
3238 numEntriesPerRow, rowPtr, colInd, values, debug);
3240 if (debug && myRank == rootRank) {
3241 cerr <<
"-- Inserting matrix entries on each process " 3242 "and calling fillComplete()" << endl;
3251 Teuchos::RCP<sparse_matrix_type> pMatrix =
3252 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
3253 pRowMap, pRangeMap, pDomainMap, constructorParams,
3254 fillCompleteParams);
3259 int localIsNull = pMatrix.is_null () ? 1 : 0;
3260 int globalIsNull = 0;
3261 reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
3262 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
3263 "Reader::makeMatrix() returned a null pointer on at least one " 3264 "process. Please report this bug to the Tpetra developers.");
3267 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
3268 "Reader::makeMatrix() returned a null pointer. " 3269 "Please report this bug to the Tpetra developers.");
3278 if (extraDebug && debug) {
3279 const int numProcs = pComm->getSize ();
3281 pRangeMap->getGlobalNumElements();
3283 pDomainMap->getGlobalNumElements();
3284 if (myRank == rootRank) {
3285 cerr <<
"-- Matrix is " 3286 << globalNumRows <<
" x " << globalNumCols
3287 <<
" with " << pMatrix->getGlobalNumEntries()
3288 <<
" entries, and index base " 3289 << pMatrix->getIndexBase() <<
"." << endl;
3292 for (
int p = 0; p < numProcs; ++p) {
3294 cerr <<
"-- Proc " << p <<
" owns " 3295 << pMatrix->getNodeNumCols() <<
" columns, and " 3296 << pMatrix->getNodeNumEntries() <<
" entries." << endl;
3302 if (debug && myRank == rootRank) {
3303 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data" 3349 static Teuchos::RCP<sparse_matrix_type>
3351 const Teuchos::RCP<const map_type>& rowMap,
3352 Teuchos::RCP<const map_type>& colMap,
3353 const Teuchos::RCP<const map_type>& domainMap,
3354 const Teuchos::RCP<const map_type>& rangeMap,
3355 const bool callFillComplete=
true,
3356 const bool tolerant=
false,
3357 const bool debug=
false)
3359 using Teuchos::MatrixMarket::Banner;
3360 using Teuchos::arcp;
3361 using Teuchos::ArrayRCP;
3362 using Teuchos::ArrayView;
3364 using Teuchos::broadcast;
3365 using Teuchos::Comm;
3366 using Teuchos::null;
3369 using Teuchos::reduceAll;
3370 using Teuchos::Tuple;
3373 typedef Teuchos::ScalarTraits<scalar_type> STS;
3375 RCP<const Comm<int> > pComm = rowMap->getComm ();
3376 const int myRank = pComm->getRank ();
3377 const int rootRank = 0;
3378 const bool extraDebug =
false;
3383 TEUCHOS_TEST_FOR_EXCEPTION(
3384 rowMap.is_null (), std::invalid_argument,
3385 "Row Map must be nonnull.");
3386 TEUCHOS_TEST_FOR_EXCEPTION(
3387 rangeMap.is_null (), std::invalid_argument,
3388 "Range Map must be nonnull.");
3389 TEUCHOS_TEST_FOR_EXCEPTION(
3390 domainMap.is_null (), std::invalid_argument,
3391 "Domain Map must be nonnull.");
3392 TEUCHOS_TEST_FOR_EXCEPTION(
3393 rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
3394 std::invalid_argument,
3395 "The specified row Map's communicator (rowMap->getComm())" 3396 "differs from the given separately supplied communicator pComm.");
3397 TEUCHOS_TEST_FOR_EXCEPTION(
3398 domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
3399 std::invalid_argument,
3400 "The specified domain Map's communicator (domainMap->getComm())" 3401 "differs from the given separately supplied communicator pComm.");
3402 TEUCHOS_TEST_FOR_EXCEPTION(
3403 rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
3404 std::invalid_argument,
3405 "The specified range Map's communicator (rangeMap->getComm())" 3406 "differs from the given separately supplied communicator pComm.");
3411 size_t lineNumber = 1;
3413 if (debug && myRank == rootRank) {
3414 cerr <<
"Matrix Market reader: readSparse:" << endl
3415 <<
"-- Reading banner line" << endl;
3424 RCP<const Banner> pBanner;
3430 int bannerIsCorrect = 1;
3431 std::ostringstream errMsg;
3433 if (myRank == rootRank) {
3436 pBanner = readBanner (in, lineNumber, tolerant, debug);
3438 catch (std::exception& e) {
3439 errMsg <<
"Attempt to read the Matrix Market file's Banner line " 3440 "threw an exception: " << e.what();
3441 bannerIsCorrect = 0;
3444 if (bannerIsCorrect) {
3449 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
3450 bannerIsCorrect = 0;
3451 errMsg <<
"The Matrix Market input file must contain a " 3452 "\"coordinate\"-format sparse matrix in order to create a " 3453 "Tpetra::CrsMatrix object from it, but the file's matrix " 3454 "type is \"" << pBanner->matrixType() <<
"\" instead.";
3459 if (tolerant && pBanner->matrixType() ==
"array") {
3460 bannerIsCorrect = 0;
3461 errMsg <<
"Matrix Market file must contain a \"coordinate\"-" 3462 "format sparse matrix in order to create a Tpetra::CrsMatrix " 3463 "object from it, but the file's matrix type is \"array\" " 3464 "instead. That probably means the file contains dense matrix " 3471 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
3478 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
3479 std::invalid_argument, errMsg.str ());
3481 if (debug && myRank == rootRank) {
3482 cerr <<
"-- Reading dimensions line" << endl;
3490 Tuple<global_ordinal_type, 3> dims =
3491 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
3493 if (debug && myRank == rootRank) {
3494 cerr <<
"-- Making Adder for collecting matrix data" << endl;
3501 RCP<adder_type> pAdder =
3502 makeAdder (pComm, pBanner, dims, tolerant, debug);
3504 if (debug && myRank == rootRank) {
3505 cerr <<
"-- Reading matrix data" << endl;
3515 int readSuccess = 1;
3516 std::ostringstream errMsg;
3517 if (myRank == rootRank) {
3520 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3522 reader_type reader (pAdder);
3525 std::pair<bool, std::vector<size_t> > results =
3526 reader.read (in, lineNumber, tolerant, debug);
3527 readSuccess = results.first ? 1 : 0;
3529 catch (std::exception& e) {
3534 broadcast (*pComm, rootRank, ptr (&readSuccess));
3543 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
3544 "Failed to read the Matrix Market sparse matrix file: " 3548 if (debug && myRank == rootRank) {
3549 cerr <<
"-- Successfully read the Matrix Market data" << endl;
3562 if (debug && myRank == rootRank) {
3563 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions" 3565 <<
"----- Dimensions before: " 3566 << dims[0] <<
" x " << dims[1]
3570 Tuple<global_ordinal_type, 2> updatedDims;
3571 if (myRank == rootRank) {
3578 std::max (dims[0], pAdder->getAdder()->numRows());
3579 updatedDims[1] = pAdder->getAdder()->numCols();
3581 broadcast (*pComm, rootRank, updatedDims);
3582 dims[0] = updatedDims[0];
3583 dims[1] = updatedDims[1];
3584 if (debug && myRank == rootRank) {
3585 cerr <<
"----- Dimensions after: " << dims[0] <<
" x " 3599 if (myRank == rootRank) {
3606 if (dims[0] < pAdder->getAdder ()->numRows ()) {
3610 broadcast (*pComm, 0, ptr (&dimsMatch));
3611 if (dimsMatch == 0) {
3618 Tuple<global_ordinal_type, 2> addersDims;
3619 if (myRank == rootRank) {
3620 addersDims[0] = pAdder->getAdder()->numRows();
3621 addersDims[1] = pAdder->getAdder()->numCols();
3623 broadcast (*pComm, 0, addersDims);
3624 TEUCHOS_TEST_FOR_EXCEPTION(
3625 dimsMatch == 0, std::runtime_error,
3626 "The matrix metadata says that the matrix is " << dims[0] <<
" x " 3627 << dims[1] <<
", but the actual data says that the matrix is " 3628 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the " 3629 "data includes more rows than reported in the metadata. This " 3630 "is not allowed when parsing in strict mode. Parse the matrix in " 3631 "tolerant mode to ignore the metadata when it disagrees with the " 3636 if (debug && myRank == rootRank) {
3637 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
3649 ArrayRCP<size_t> numEntriesPerRow;
3650 ArrayRCP<size_t> rowPtr;
3651 ArrayRCP<global_ordinal_type> colInd;
3652 ArrayRCP<scalar_type> values;
3657 int mergeAndConvertSucceeded = 1;
3658 std::ostringstream errMsg;
3660 if (myRank == rootRank) {
3662 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3672 const size_type numRows = dims[0];
3675 pAdder->getAdder()->merge ();
3678 const std::vector<element_type>& entries =
3679 pAdder->getAdder()->getEntries();
3682 const size_t numEntries = (size_t)entries.size();
3685 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
3686 <<
" rows and numEntries=" << numEntries
3687 <<
" entries." << endl;
3693 numEntriesPerRow = arcp<size_t> (numRows);
3694 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3695 rowPtr = arcp<size_t> (numRows+1);
3696 std::fill (rowPtr.begin(), rowPtr.end(), 0);
3697 colInd = arcp<global_ordinal_type> (numEntries);
3698 values = arcp<scalar_type> (numEntries);
3705 for (curPos = 0; curPos < numEntries; ++curPos) {
3706 const element_type& curEntry = entries[curPos];
3708 TEUCHOS_TEST_FOR_EXCEPTION(
3709 curRow < prvRow, std::logic_error,
3710 "Row indices are out of order, even though they are supposed " 3711 "to be sorted. curRow = " << curRow <<
", prvRow = " 3712 << prvRow <<
", at curPos = " << curPos <<
". Please report " 3713 "this bug to the Tpetra developers.");
3714 if (curRow > prvRow) {
3720 numEntriesPerRow[curRow]++;
3721 colInd[curPos] = curEntry.colIndex();
3722 values[curPos] = curEntry.value();
3727 rowPtr[numRows] = numEntries;
3729 catch (std::exception& e) {
3730 mergeAndConvertSucceeded = 0;
3731 errMsg <<
"Failed to merge sparse matrix entries and convert to " 3732 "CSR format: " << e.what();
3735 if (debug && mergeAndConvertSucceeded) {
3737 const size_type numRows = dims[0];
3738 const size_type maxToDisplay = 100;
3740 cerr <<
"----- Proc 0: numEntriesPerRow[0.." 3741 << (numEntriesPerRow.size()-1) <<
"] ";
3742 if (numRows > maxToDisplay) {
3743 cerr <<
"(only showing first and last few entries) ";
3747 if (numRows > maxToDisplay) {
3748 for (size_type k = 0; k < 2; ++k) {
3749 cerr << numEntriesPerRow[k] <<
" ";
3752 for (size_type k = numRows-2; k < numRows-1; ++k) {
3753 cerr << numEntriesPerRow[k] <<
" ";
3757 for (size_type k = 0; k < numRows-1; ++k) {
3758 cerr << numEntriesPerRow[k] <<
" ";
3761 cerr << numEntriesPerRow[numRows-1];
3763 cerr <<
"]" << endl;
3765 cerr <<
"----- Proc 0: rowPtr ";
3766 if (numRows > maxToDisplay) {
3767 cerr <<
"(only showing first and last few entries) ";
3770 if (numRows > maxToDisplay) {
3771 for (size_type k = 0; k < 2; ++k) {
3772 cerr << rowPtr[k] <<
" ";
3775 for (size_type k = numRows-2; k < numRows; ++k) {
3776 cerr << rowPtr[k] <<
" ";
3780 for (size_type k = 0; k < numRows; ++k) {
3781 cerr << rowPtr[k] <<
" ";
3784 cerr << rowPtr[numRows] <<
"]" << endl;
3786 cerr <<
"----- Proc 0: colInd = [";
3787 for (
size_t k = 0; k < rowPtr[numRows]; ++k) {
3788 cerr << colInd[k] <<
" ";
3790 cerr <<
"]" << endl;
3804 if (debug && myRank == rootRank) {
3805 cerr <<
"-- Verifying Maps" << endl;
3807 TEUCHOS_TEST_FOR_EXCEPTION(
3808 as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
3809 std::invalid_argument,
3810 "The range Map has " << rangeMap->getGlobalNumElements ()
3811 <<
" entries, but the matrix has a global number of rows " << dims[0]
3813 TEUCHOS_TEST_FOR_EXCEPTION(
3814 as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
3815 std::invalid_argument,
3816 "The domain Map has " << domainMap->getGlobalNumElements ()
3817 <<
" entries, but the matrix has a global number of columns " 3821 RCP<Teuchos::FancyOStream> err = debug ?
3822 Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
3824 RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
3825 ArrayView<const global_ordinal_type> myRows =
3826 gatherRowMap->getNodeElementList ();
3827 const size_type myNumRows = myRows.size ();
3830 ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
3831 for (size_type i_ = 0; i_ < myNumRows; i_++) {
3832 gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_]-indexBase];
3838 RCP<sparse_matrix_type> A_proc0 =
3841 if (myRank == rootRank) {
3843 cerr <<
"-- Proc 0: Filling gather matrix" << endl;
3846 cerr <<
"---- Rows: " << Teuchos::toString (myRows) << endl;
3850 for (size_type i_ = 0; i_ < myNumRows; ++i_) {
3851 size_type i = myRows[i_] - indexBase;
3853 const size_type curPos = as<size_type> (rowPtr[i]);
3855 ArrayView<global_ordinal_type> curColInd =
3856 colInd.view (curPos, curNumEntries);
3857 ArrayView<scalar_type> curValues =
3858 values.view (curPos, curNumEntries);
3861 for (size_type k = 0; k < curNumEntries; ++k) {
3862 curColInd[k] += indexBase;
3865 cerr <<
"------ Columns: " << Teuchos::toString (curColInd) << endl;
3866 cerr <<
"------ Values: " << Teuchos::toString (curValues) << endl;
3869 if (curNumEntries > 0) {
3870 A_proc0->insertGlobalValues (myRows[i_], curColInd, curValues);
3876 numEntriesPerRow = null;
3882 RCP<sparse_matrix_type> A;
3883 if (colMap.is_null ()) {
3889 export_type exp (gatherRowMap, rowMap);
3890 A->doExport (*A_proc0, exp,
INSERT);
3892 if (callFillComplete) {
3893 A->fillComplete (domainMap, rangeMap);
3905 if (callFillComplete) {
3906 const int numProcs = pComm->getSize ();
3908 if (extraDebug && debug) {
3909 const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
3910 const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
3911 if (myRank == rootRank) {
3912 cerr <<
"-- Matrix is " 3913 << globalNumRows <<
" x " << globalNumCols
3914 <<
" with " << A->getGlobalNumEntries()
3915 <<
" entries, and index base " 3916 << A->getIndexBase() <<
"." << endl;
3919 for (
int p = 0; p < numProcs; ++p) {
3921 cerr <<
"-- Proc " << p <<
" owns " 3922 << A->getNodeNumCols() <<
" columns, and " 3923 << A->getNodeNumEntries() <<
" entries." << endl;
3930 if (debug && myRank == rootRank) {
3931 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data" 3966 static Teuchos::RCP<multivector_type>
3968 const Teuchos::RCP<const comm_type>& comm,
3969 Teuchos::RCP<const map_type>& map,
3970 const bool tolerant=
false,
3971 const bool debug=
false)
3974 if (comm->getRank () == 0) {
3975 in.open (filename.c_str ());
3977 return readDense (in, comm, map, tolerant, debug);
3981 static Teuchos::RCP<multivector_type>
3983 const Teuchos::RCP<const comm_type>& comm,
3984 const Teuchos::RCP<node_type>& node,
3985 Teuchos::RCP<const map_type>& map,
3986 const bool tolerant=
false,
3987 const bool debug=
false)
3990 if (comm->getRank () == 0) {
3991 in.open (filename.c_str ());
3993 return readDense (in, comm, node, map, tolerant, debug);
4025 static Teuchos::RCP<vector_type>
4027 const Teuchos::RCP<const comm_type>& comm,
4028 Teuchos::RCP<const map_type>& map,
4029 const bool tolerant=
false,
4030 const bool debug=
false)
4033 if (comm->getRank () == 0) {
4034 in.open (filename.c_str ());
4036 return readVector (in, comm, map, tolerant, debug);
4041 static Teuchos::RCP<vector_type>
4043 const Teuchos::RCP<const comm_type>& comm,
4044 const Teuchos::RCP<node_type>& node,
4045 Teuchos::RCP<const map_type>& map,
4046 const bool tolerant=
false,
4047 const bool debug=
false)
4050 if (comm->getRank () == 0) {
4051 in.open (filename.c_str ());
4053 return readVector (in, comm, node, map, tolerant, debug);
4123 static Teuchos::RCP<multivector_type>
4125 const Teuchos::RCP<const comm_type>& comm,
4126 Teuchos::RCP<const map_type>& map,
4127 const bool tolerant=
false,
4128 const bool debug=
false)
4130 return readDense (in, comm, Teuchos::null, map, tolerant, debug);
4134 static Teuchos::RCP<multivector_type>
4136 const Teuchos::RCP<const comm_type>& comm,
4137 const Teuchos::RCP<node_type>& node,
4138 Teuchos::RCP<const map_type>& map,
4139 const bool tolerant=
false,
4140 const bool debug=
false)
4142 Teuchos::RCP<Teuchos::FancyOStream> err =
4143 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4144 return readDenseImpl<scalar_type> (in, comm, node, map,
4145 err, tolerant, debug);
4149 static Teuchos::RCP<vector_type>
4151 const Teuchos::RCP<const comm_type>& comm,
4152 Teuchos::RCP<const map_type>& map,
4153 const bool tolerant=
false,
4154 const bool debug=
false)
4156 Teuchos::RCP<Teuchos::FancyOStream> err =
4157 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4158 return readVectorImpl<scalar_type> (in, comm, Teuchos::null, map,
4159 err, tolerant, debug);
4163 static Teuchos::RCP<vector_type>
4165 const Teuchos::RCP<const comm_type>& comm,
4166 const Teuchos::RCP<node_type>& node,
4167 Teuchos::RCP<const map_type>& map,
4168 const bool tolerant=
false,
4169 const bool debug=
false)
4171 Teuchos::RCP<Teuchos::FancyOStream> err =
4172 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4173 return readVectorImpl<scalar_type> (in, comm, node, map,
4174 err, tolerant, debug);
4197 static Teuchos::RCP<const map_type>
4199 const Teuchos::RCP<const comm_type>& comm,
4200 const bool tolerant=
false,
4201 const bool debug=
false)
4203 return readMapFile (filename, comm, Teuchos::null, tolerant, debug);
4208 static Teuchos::RCP<const map_type>
4210 const Teuchos::RCP<const comm_type>& comm,
4211 const Teuchos::RCP<node_type>& node,
4212 const bool tolerant=
false,
4213 const bool debug=
false)
4215 using Teuchos::inOutArg;
4216 using Teuchos::broadcast;
4220 if (comm->getRank () == 0) {
4221 in.open (filename.c_str ());
4226 broadcast<int, int> (*comm, 0, inOutArg (success));
4227 TEUCHOS_TEST_FOR_EXCEPTION(
4228 success == 0, std::runtime_error,
4229 "Tpetra::MatrixMarket::Reader::readMapFile: " 4230 "Failed to read file \"" << filename <<
"\" on Process 0.");
4231 return readMap (in, comm, node, tolerant, debug);
4235 template<
class MultiVectorScalarType>
4240 readDenseImpl (std::istream& in,
4241 const Teuchos::RCP<const comm_type>& comm,
4242 const Teuchos::RCP<node_type>& node,
4243 Teuchos::RCP<const map_type>& map,
4244 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4245 const bool tolerant=
false,
4246 const bool debug=
false)
4248 using Teuchos::MatrixMarket::Banner;
4249 using Teuchos::MatrixMarket::checkCommentLine;
4250 using Teuchos::ArrayRCP;
4252 using Teuchos::broadcast;
4253 using Teuchos::outArg;
4255 using Teuchos::Tuple;
4257 typedef MultiVectorScalarType ST;
4261 typedef Teuchos::ScalarTraits<ST> STS;
4262 typedef typename STS::magnitudeType MT;
4263 typedef Teuchos::ScalarTraits<MT> STM;
4268 const int myRank = comm->getRank ();
4270 if (! err.is_null ()) {
4274 *err << myRank <<
": readDenseImpl" << endl;
4276 if (! err.is_null ()) {
4310 size_t lineNumber = 1;
4313 std::ostringstream exMsg;
4314 int localBannerReadSuccess = 1;
4315 int localDimsReadSuccess = 1;
4320 *err << myRank <<
": readDenseImpl: Reading banner line (dense)" << endl;
4326 RCP<const Banner> pBanner;
4328 pBanner = readBanner (in, lineNumber, tolerant, debug);
4329 }
catch (std::exception& e) {
4331 localBannerReadSuccess = 0;
4334 if (localBannerReadSuccess) {
4335 if (pBanner->matrixType () !=
"array") {
4336 exMsg <<
"The Matrix Market file does not contain dense matrix " 4337 "data. Its banner (first) line says that its matrix type is \"" 4338 << pBanner->matrixType () <<
"\", rather that the required " 4340 localBannerReadSuccess = 0;
4341 }
else if (pBanner->dataType() ==
"pattern") {
4342 exMsg <<
"The Matrix Market file's banner (first) " 4343 "line claims that the matrix's data type is \"pattern\". This does " 4344 "not make sense for a dense matrix, yet the file reports the matrix " 4345 "as dense. The only valid data types for a dense matrix are " 4346 "\"real\", \"complex\", and \"integer\".";
4347 localBannerReadSuccess = 0;
4351 dims[2] = encodeDataType (pBanner->dataType ());
4357 if (localBannerReadSuccess) {
4359 *err << myRank <<
": readDenseImpl: Reading dimensions line (dense)" << endl;
4368 bool commentLine =
true;
4370 while (commentLine) {
4373 if (in.eof () || in.fail ()) {
4374 exMsg <<
"Unable to get array dimensions line (at all) from line " 4375 << lineNumber <<
" of input stream. The input stream " 4376 <<
"claims that it is " 4377 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4378 localDimsReadSuccess = 0;
4381 if (getline (in, line)) {
4388 size_t start = 0, size = 0;
4389 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4396 std::istringstream istr (line);
4400 if (istr.eof () || istr.fail ()) {
4401 exMsg <<
"Unable to read any data from line " << lineNumber
4402 <<
" of input; the line should contain the matrix dimensions " 4403 <<
"\"<numRows> <numCols>\".";
4404 localDimsReadSuccess = 0;
4409 exMsg <<
"Failed to get number of rows from line " 4410 << lineNumber <<
" of input; the line should contains the " 4411 <<
"matrix dimensions \"<numRows> <numCols>\".";
4412 localDimsReadSuccess = 0;
4414 dims[0] = theNumRows;
4416 exMsg <<
"No more data after number of rows on line " 4417 << lineNumber <<
" of input; the line should contain the " 4418 <<
"matrix dimensions \"<numRows> <numCols>\".";
4419 localDimsReadSuccess = 0;
4424 exMsg <<
"Failed to get number of columns from line " 4425 << lineNumber <<
" of input; the line should contain " 4426 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4427 localDimsReadSuccess = 0;
4429 dims[1] = theNumCols;
4440 Tuple<GO, 5> bannerDimsReadResult;
4442 bannerDimsReadResult[0] = dims[0];
4443 bannerDimsReadResult[1] = dims[1];
4444 bannerDimsReadResult[2] = dims[2];
4445 bannerDimsReadResult[3] = localBannerReadSuccess;
4446 bannerDimsReadResult[4] = localDimsReadSuccess;
4450 broadcast (*comm, 0, bannerDimsReadResult);
4452 TEUCHOS_TEST_FOR_EXCEPTION(
4453 bannerDimsReadResult[3] == 0, std::runtime_error,
4454 "Failed to read banner line: " << exMsg.str ());
4455 TEUCHOS_TEST_FOR_EXCEPTION(
4456 bannerDimsReadResult[4] == 0, std::runtime_error,
4457 "Failed to read matrix dimensions line: " << exMsg.str ());
4459 dims[0] = bannerDimsReadResult[0];
4460 dims[1] = bannerDimsReadResult[1];
4461 dims[2] = bannerDimsReadResult[2];
4466 const size_t numCols =
static_cast<size_t> (dims[1]);
4471 RCP<const map_type> proc0Map;
4472 if (map.is_null ()) {
4476 if (node.is_null ()) {
4477 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4479 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
4481 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4483 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4484 comm, map->getNode ());
4487 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4491 RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
4497 int localReadDataSuccess = 1;
4501 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)" 4506 TEUCHOS_TEST_FOR_EXCEPTION(
4507 ! X->isConstantStride (), std::logic_error,
4508 "Can't get a 1-D view of the entries of the MultiVector X on " 4509 "Process 0, because the stride between the columns of X is not " 4510 "constant. This shouldn't happen because we just created X and " 4511 "haven't filled it in yet. Please report this bug to the Tpetra " 4518 ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4519 TEUCHOS_TEST_FOR_EXCEPTION(
4520 as<global_size_t> (X_view.size ()) < numRows * numCols,
4522 "The view of X has size " << X_view <<
" which is not enough to " 4523 "accommodate the expected number of entries numRows*numCols = " 4524 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". " 4525 "Please report this bug to the Tpetra developers.");
4526 const size_t stride = X->getStride ();
4533 const bool isComplex = (dims[2] == 1);
4534 size_type count = 0, curRow = 0, curCol = 0;
4537 while (getline (in, line)) {
4541 size_t start = 0, size = 0;
4542 const bool commentLine =
4543 checkCommentLine (line, start, size, lineNumber, tolerant);
4544 if (! commentLine) {
4550 if (count >= X_view.size()) {
4555 TEUCHOS_TEST_FOR_EXCEPTION(
4556 count >= X_view.size(),
4558 "The Matrix Market input stream has more data in it than " 4559 "its metadata reported. Current line number is " 4560 << lineNumber <<
".");
4569 const size_t pos = line.substr (start, size).find (
':');
4570 if (pos != std::string::npos) {
4574 std::istringstream istr (line.substr (start, size));
4578 if (istr.eof() || istr.fail()) {
4585 TEUCHOS_TEST_FOR_EXCEPTION(
4586 ! tolerant && (istr.eof() || istr.fail()),
4588 "Line " << lineNumber <<
" of the Matrix Market file is " 4589 "empty, or we cannot read from it for some other reason.");
4592 ST val = STS::zero();
4595 MT real = STM::zero(), imag = STM::zero();
4612 TEUCHOS_TEST_FOR_EXCEPTION(
4613 ! tolerant && istr.eof(), std::runtime_error,
4614 "Failed to get the real part of a complex-valued matrix " 4615 "entry from line " << lineNumber <<
" of the Matrix Market " 4621 }
else if (istr.eof()) {
4622 TEUCHOS_TEST_FOR_EXCEPTION(
4623 ! tolerant && istr.eof(), std::runtime_error,
4624 "Missing imaginary part of a complex-valued matrix entry " 4625 "on line " << lineNumber <<
" of the Matrix Market file.");
4631 TEUCHOS_TEST_FOR_EXCEPTION(
4632 ! tolerant && istr.fail(), std::runtime_error,
4633 "Failed to get the imaginary part of a complex-valued " 4634 "matrix entry from line " << lineNumber <<
" of the " 4635 "Matrix Market file.");
4642 TEUCHOS_TEST_FOR_EXCEPTION(
4643 ! tolerant && istr.fail(), std::runtime_error,
4644 "Failed to get a real-valued matrix entry from line " 4645 << lineNumber <<
" of the Matrix Market file.");
4648 if (istr.fail() && tolerant) {
4654 TEUCHOS_TEST_FOR_EXCEPTION(
4655 ! tolerant && istr.fail(), std::runtime_error,
4656 "Failed to read matrix data from line " << lineNumber
4657 <<
" of the Matrix Market file.");
4660 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4662 curRow = count % numRows;
4663 curCol = count / numRows;
4664 X_view[curRow + curCol*stride] = val;
4669 TEUCHOS_TEST_FOR_EXCEPTION(
4670 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
4672 "The Matrix Market metadata reports that the dense matrix is " 4673 << numRows <<
" x " << numCols <<
", and thus has " 4674 << numRows*numCols <<
" total entries, but we only found " << count
4675 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
4676 }
catch (std::exception& e) {
4678 localReadDataSuccess = 0;
4683 *err << myRank <<
": readDenseImpl: done reading data" << endl;
4687 int globalReadDataSuccess = localReadDataSuccess;
4688 broadcast (*comm, 0, outArg (globalReadDataSuccess));
4689 TEUCHOS_TEST_FOR_EXCEPTION(
4690 globalReadDataSuccess == 0, std::runtime_error,
4691 "Failed to read the multivector's data: " << exMsg.str ());
4696 if (comm->getSize () == 1 && map.is_null ()) {
4698 if (! err.is_null ()) {
4702 *err << myRank <<
": readDenseImpl: done" << endl;
4704 if (! err.is_null ()) {
4711 *err << myRank <<
": readDenseImpl: Creating target MV" << endl;
4715 RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
4718 *err << myRank <<
": readDenseImpl: Creating Export" << endl;
4724 Export<LO, GO, NT> exporter (proc0Map, map, err);
4727 *err << myRank <<
": readDenseImpl: Exporting" << endl;
4730 Y->doExport (*X, exporter,
INSERT);
4732 if (! err.is_null ()) {
4736 *err << myRank <<
": readDenseImpl: done" << endl;
4738 if (! err.is_null ()) {
4747 template<
class VectorScalarType>
4752 readVectorImpl (std::istream& in,
4753 const Teuchos::RCP<const comm_type>& comm,
4754 const Teuchos::RCP<node_type>& node,
4755 Teuchos::RCP<const map_type>& map,
4756 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4757 const bool tolerant=
false,
4758 const bool debug=
false)
4760 using Teuchos::MatrixMarket::Banner;
4761 using Teuchos::MatrixMarket::checkCommentLine;
4763 using Teuchos::broadcast;
4764 using Teuchos::outArg;
4766 using Teuchos::Tuple;
4768 typedef VectorScalarType ST;
4772 typedef Teuchos::ScalarTraits<ST> STS;
4773 typedef typename STS::magnitudeType MT;
4774 typedef Teuchos::ScalarTraits<MT> STM;
4779 const int myRank = comm->getRank ();
4781 if (! err.is_null ()) {
4785 *err << myRank <<
": readVectorImpl" << endl;
4787 if (! err.is_null ()) {
4821 size_t lineNumber = 1;
4824 std::ostringstream exMsg;
4825 int localBannerReadSuccess = 1;
4826 int localDimsReadSuccess = 1;
4831 *err << myRank <<
": readVectorImpl: Reading banner line (dense)" << endl;
4837 RCP<const Banner> pBanner;
4839 pBanner = readBanner (in, lineNumber, tolerant, debug);
4840 }
catch (std::exception& e) {
4842 localBannerReadSuccess = 0;
4845 if (localBannerReadSuccess) {
4846 if (pBanner->matrixType () !=
"array") {
4847 exMsg <<
"The Matrix Market file does not contain dense matrix " 4848 "data. Its banner (first) line says that its matrix type is \"" 4849 << pBanner->matrixType () <<
"\", rather that the required " 4851 localBannerReadSuccess = 0;
4852 }
else if (pBanner->dataType() ==
"pattern") {
4853 exMsg <<
"The Matrix Market file's banner (first) " 4854 "line claims that the matrix's data type is \"pattern\". This does " 4855 "not make sense for a dense matrix, yet the file reports the matrix " 4856 "as dense. The only valid data types for a dense matrix are " 4857 "\"real\", \"complex\", and \"integer\".";
4858 localBannerReadSuccess = 0;
4862 dims[2] = encodeDataType (pBanner->dataType ());
4868 if (localBannerReadSuccess) {
4870 *err << myRank <<
": readVectorImpl: Reading dimensions line (dense)" << endl;
4879 bool commentLine =
true;
4881 while (commentLine) {
4884 if (in.eof () || in.fail ()) {
4885 exMsg <<
"Unable to get array dimensions line (at all) from line " 4886 << lineNumber <<
" of input stream. The input stream " 4887 <<
"claims that it is " 4888 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4889 localDimsReadSuccess = 0;
4892 if (getline (in, line)) {
4899 size_t start = 0, size = 0;
4900 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4907 std::istringstream istr (line);
4911 if (istr.eof () || istr.fail ()) {
4912 exMsg <<
"Unable to read any data from line " << lineNumber
4913 <<
" of input; the line should contain the matrix dimensions " 4914 <<
"\"<numRows> <numCols>\".";
4915 localDimsReadSuccess = 0;
4920 exMsg <<
"Failed to get number of rows from line " 4921 << lineNumber <<
" of input; the line should contains the " 4922 <<
"matrix dimensions \"<numRows> <numCols>\".";
4923 localDimsReadSuccess = 0;
4925 dims[0] = theNumRows;
4927 exMsg <<
"No more data after number of rows on line " 4928 << lineNumber <<
" of input; the line should contain the " 4929 <<
"matrix dimensions \"<numRows> <numCols>\".";
4930 localDimsReadSuccess = 0;
4935 exMsg <<
"Failed to get number of columns from line " 4936 << lineNumber <<
" of input; the line should contain " 4937 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4938 localDimsReadSuccess = 0;
4940 dims[1] = theNumCols;
4950 exMsg <<
"File does not contain a 1D Vector";
4951 localDimsReadSuccess = 0;
4957 Tuple<GO, 5> bannerDimsReadResult;
4959 bannerDimsReadResult[0] = dims[0];
4960 bannerDimsReadResult[1] = dims[1];
4961 bannerDimsReadResult[2] = dims[2];
4962 bannerDimsReadResult[3] = localBannerReadSuccess;
4963 bannerDimsReadResult[4] = localDimsReadSuccess;
4968 broadcast (*comm, 0, bannerDimsReadResult);
4970 TEUCHOS_TEST_FOR_EXCEPTION(
4971 bannerDimsReadResult[3] == 0, std::runtime_error,
4972 "Failed to read banner line: " << exMsg.str ());
4973 TEUCHOS_TEST_FOR_EXCEPTION(
4974 bannerDimsReadResult[4] == 0, std::runtime_error,
4975 "Failed to read matrix dimensions line: " << exMsg.str ());
4977 dims[0] = bannerDimsReadResult[0];
4978 dims[1] = bannerDimsReadResult[1];
4979 dims[2] = bannerDimsReadResult[2];
4984 const size_t numCols =
static_cast<size_t> (dims[1]);
4989 RCP<const map_type> proc0Map;
4990 if (map.is_null ()) {
4994 if (node.is_null ()) {
4995 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4997 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
4999 const size_t localNumRows = (myRank == 0) ? numRows : 0;
5001 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
5002 comm, map->getNode ());
5005 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
5009 RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
5015 int localReadDataSuccess = 1;
5019 *err << myRank <<
": readVectorImpl: Reading matrix data (dense)" 5024 TEUCHOS_TEST_FOR_EXCEPTION(
5025 ! X->isConstantStride (), std::logic_error,
5026 "Can't get a 1-D view of the entries of the MultiVector X on " 5027 "Process 0, because the stride between the columns of X is not " 5028 "constant. This shouldn't happen because we just created X and " 5029 "haven't filled it in yet. Please report this bug to the Tpetra " 5036 Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst ();
5037 TEUCHOS_TEST_FOR_EXCEPTION(
5038 as<global_size_t> (X_view.size ()) < numRows * numCols,
5040 "The view of X has size " << X_view <<
" which is not enough to " 5041 "accommodate the expected number of entries numRows*numCols = " 5042 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". " 5043 "Please report this bug to the Tpetra developers.");
5044 const size_t stride = X->getStride ();
5051 const bool isComplex = (dims[2] == 1);
5052 size_type count = 0, curRow = 0, curCol = 0;
5055 while (getline (in, line)) {
5059 size_t start = 0, size = 0;
5060 const bool commentLine =
5061 checkCommentLine (line, start, size, lineNumber, tolerant);
5062 if (! commentLine) {
5068 if (count >= X_view.size()) {
5073 TEUCHOS_TEST_FOR_EXCEPTION(
5074 count >= X_view.size(),
5076 "The Matrix Market input stream has more data in it than " 5077 "its metadata reported. Current line number is " 5078 << lineNumber <<
".");
5087 const size_t pos = line.substr (start, size).find (
':');
5088 if (pos != std::string::npos) {
5092 std::istringstream istr (line.substr (start, size));
5096 if (istr.eof() || istr.fail()) {
5103 TEUCHOS_TEST_FOR_EXCEPTION(
5104 ! tolerant && (istr.eof() || istr.fail()),
5106 "Line " << lineNumber <<
" of the Matrix Market file is " 5107 "empty, or we cannot read from it for some other reason.");
5110 ST val = STS::zero();
5113 MT real = STM::zero(), imag = STM::zero();
5130 TEUCHOS_TEST_FOR_EXCEPTION(
5131 ! tolerant && istr.eof(), std::runtime_error,
5132 "Failed to get the real part of a complex-valued matrix " 5133 "entry from line " << lineNumber <<
" of the Matrix Market " 5139 }
else if (istr.eof()) {
5140 TEUCHOS_TEST_FOR_EXCEPTION(
5141 ! tolerant && istr.eof(), std::runtime_error,
5142 "Missing imaginary part of a complex-valued matrix entry " 5143 "on line " << lineNumber <<
" of the Matrix Market file.");
5149 TEUCHOS_TEST_FOR_EXCEPTION(
5150 ! tolerant && istr.fail(), std::runtime_error,
5151 "Failed to get the imaginary part of a complex-valued " 5152 "matrix entry from line " << lineNumber <<
" of the " 5153 "Matrix Market file.");
5160 TEUCHOS_TEST_FOR_EXCEPTION(
5161 ! tolerant && istr.fail(), std::runtime_error,
5162 "Failed to get a real-valued matrix entry from line " 5163 << lineNumber <<
" of the Matrix Market file.");
5166 if (istr.fail() && tolerant) {
5172 TEUCHOS_TEST_FOR_EXCEPTION(
5173 ! tolerant && istr.fail(), std::runtime_error,
5174 "Failed to read matrix data from line " << lineNumber
5175 <<
" of the Matrix Market file.");
5178 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
5180 curRow = count % numRows;
5181 curCol = count / numRows;
5182 X_view[curRow + curCol*stride] = val;
5187 TEUCHOS_TEST_FOR_EXCEPTION(
5188 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
5190 "The Matrix Market metadata reports that the dense matrix is " 5191 << numRows <<
" x " << numCols <<
", and thus has " 5192 << numRows*numCols <<
" total entries, but we only found " << count
5193 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
5194 }
catch (std::exception& e) {
5196 localReadDataSuccess = 0;
5201 *err << myRank <<
": readVectorImpl: done reading data" << endl;
5205 int globalReadDataSuccess = localReadDataSuccess;
5206 broadcast (*comm, 0, outArg (globalReadDataSuccess));
5207 TEUCHOS_TEST_FOR_EXCEPTION(
5208 globalReadDataSuccess == 0, std::runtime_error,
5209 "Failed to read the multivector's data: " << exMsg.str ());
5214 if (comm->getSize () == 1 && map.is_null ()) {
5216 if (! err.is_null ()) {
5220 *err << myRank <<
": readVectorImpl: done" << endl;
5222 if (! err.is_null ()) {
5229 *err << myRank <<
": readVectorImpl: Creating target MV" << endl;
5233 RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
5236 *err << myRank <<
": readVectorImpl: Creating Export" << endl;
5242 Export<LO, GO, NT> exporter (proc0Map, map, err);
5245 *err << myRank <<
": readVectorImpl: Exporting" << endl;
5248 Y->doExport (*X, exporter,
INSERT);
5250 if (! err.is_null ()) {
5254 *err << myRank <<
": readVectorImpl: done" << endl;
5256 if (! err.is_null ()) {
5285 static Teuchos::RCP<const map_type>
5287 const Teuchos::RCP<const comm_type>& comm,
5288 const bool tolerant=
false,
5289 const bool debug=
false)
5291 Teuchos::RCP<Teuchos::FancyOStream> err =
5292 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
5293 return readMap (in, comm, err, tolerant, debug);
5298 static Teuchos::RCP<const map_type>
5300 const Teuchos::RCP<const comm_type>& comm,
5301 const Teuchos::RCP<node_type>& node,
5302 const bool tolerant=
false,
5303 const bool debug=
false)
5305 Teuchos::RCP<Teuchos::FancyOStream> err =
5306 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
5307 return readMap (in, comm, node, err, tolerant, debug);
5335 static Teuchos::RCP<const map_type>
5337 const Teuchos::RCP<const comm_type>& comm,
5338 const Teuchos::RCP<Teuchos::FancyOStream>& err,
5339 const bool tolerant=
false,
5340 const bool debug=
false)
5342 return readMap (in, comm, Teuchos::null, err, tolerant, debug);
5347 static Teuchos::RCP<const map_type>
5349 const Teuchos::RCP<const comm_type>& comm,
5350 const Teuchos::RCP<node_type>& node,
5351 const Teuchos::RCP<Teuchos::FancyOStream>& err,
5352 const bool tolerant=
false,
5353 const bool debug=
false)
5355 using Teuchos::arcp;
5356 using Teuchos::Array;
5357 using Teuchos::ArrayRCP;
5359 using Teuchos::broadcast;
5360 using Teuchos::Comm;
5361 using Teuchos::CommRequest;
5362 using Teuchos::inOutArg;
5363 using Teuchos::ireceive;
5364 using Teuchos::outArg;
5366 using Teuchos::receive;
5367 using Teuchos::reduceAll;
5368 using Teuchos::REDUCE_MIN;
5369 using Teuchos::isend;
5370 using Teuchos::SerialComm;
5371 using Teuchos::toString;
5372 using Teuchos::wait;
5375 typedef ptrdiff_t int_type;
5381 const int numProcs = comm->getSize ();
5382 const int myRank = comm->getRank ();
5384 if (err.is_null ()) {
5388 std::ostringstream os;
5389 os << myRank <<
": readMap: " << endl;
5392 if (err.is_null ()) {
5398 const int sizeTag = 1339;
5400 const int dataTag = 1340;
5406 RCP<CommRequest<int> > sizeReq;
5407 RCP<CommRequest<int> > dataReq;
5412 ArrayRCP<int_type> numGidsToRecv (1);
5413 numGidsToRecv[0] = 0;
5415 sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
5418 int readSuccess = 1;
5419 std::ostringstream exMsg;
5428 RCP<const Comm<int> > proc0Comm (
new SerialComm<int> ());
5430 RCP<const map_type> dataMap;
5434 data = readDenseImpl<GO> (in, proc0Comm, node, dataMap,
5435 err, tolerant, debug);
5437 if (data.is_null ()) {
5439 exMsg <<
"readDenseImpl() returned null." << endl;
5441 }
catch (std::exception& e) {
5443 exMsg << e.what () << endl;
5449 std::map<int, Array<GO> > pid2gids;
5454 int_type globalNumGIDs = 0;
5464 if (myRank == 0 && readSuccess == 1) {
5465 if (data->getNumVectors () == 2) {
5466 ArrayRCP<const GO> GIDs = data->getData (0);
5467 ArrayRCP<const GO> PIDs = data->getData (1);
5468 globalNumGIDs = GIDs.size ();
5472 if (globalNumGIDs > 0) {
5473 const int pid =
static_cast<int> (PIDs[0]);
5475 if (pid < 0 || pid >= numProcs) {
5477 exMsg <<
"Tpetra::MatrixMarket::readMap: " 5478 <<
"Encountered invalid PID " << pid <<
"." << endl;
5481 const GO gid = GIDs[0];
5482 pid2gids[pid].push_back (gid);
5486 if (readSuccess == 1) {
5489 for (size_type k = 1; k < globalNumGIDs; ++k) {
5490 const int pid =
static_cast<int> (PIDs[k]);
5491 if (pid < 0 || pid >= numProcs) {
5493 exMsg <<
"Tpetra::MatrixMarket::readMap: " 5494 <<
"Encountered invalid PID " << pid <<
"." << endl;
5497 const int_type gid = GIDs[k];
5498 pid2gids[pid].push_back (gid);
5499 if (gid < indexBase) {
5506 else if (data->getNumVectors () == 1) {
5507 if (data->getGlobalLength () % 2 !=
static_cast<GST
> (0)) {
5509 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data has the " 5510 "wrong format (for Map format 2.0). The global number of rows " 5511 "in the MultiVector must be even (divisible by 2)." << endl;
5514 ArrayRCP<const GO> theData = data->getData (0);
5515 globalNumGIDs =
static_cast<GO
> (data->getGlobalLength ()) /
5516 static_cast<GO> (2);
5520 if (globalNumGIDs > 0) {
5521 const int pid =
static_cast<int> (theData[1]);
5522 if (pid < 0 || pid >= numProcs) {
5524 exMsg <<
"Tpetra::MatrixMarket::readMap: " 5525 <<
"Encountered invalid PID " << pid <<
"." << endl;
5528 const GO gid = theData[0];
5529 pid2gids[pid].push_back (gid);
5535 for (int_type k = 1; k < globalNumGIDs; ++k) {
5536 const int pid =
static_cast<int> (theData[2*k + 1]);
5537 if (pid < 0 || pid >= numProcs) {
5539 exMsg <<
"Tpetra::MatrixMarket::readMap: " 5540 <<
"Encountered invalid PID " << pid <<
"." << endl;
5543 const GO gid = theData[2*k];
5544 pid2gids[pid].push_back (gid);
5545 if (gid < indexBase) {
5554 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data must have " 5555 "either 1 column (for the new Map format 2.0) or 2 columns (for " 5556 "the old Map format 1.0).";
5563 int_type readResults[3];
5564 readResults[0] =
static_cast<int_type
> (indexBase);
5565 readResults[1] =
static_cast<int_type
> (globalNumGIDs);
5566 readResults[2] =
static_cast<int_type
> (readSuccess);
5567 broadcast<int, int_type> (*comm, 0, 3, readResults);
5569 indexBase =
static_cast<GO
> (readResults[0]);
5570 globalNumGIDs =
static_cast<int_type
> (readResults[1]);
5571 readSuccess =
static_cast<int> (readResults[2]);
5577 TEUCHOS_TEST_FOR_EXCEPTION(
5578 readSuccess != 1, std::runtime_error,
5579 "Tpetra::MatrixMarket::readMap: Reading the Map failed with the " 5580 "following exception message: " << exMsg.str ());
5584 for (
int p = 1; p < numProcs; ++p) {
5585 ArrayRCP<int_type> numGidsToSend (1);
5587 typename std::map<int, Array<GO> >::const_iterator it = pid2gids.find (p);
5588 if (it == pid2gids.end ()) {
5589 numGidsToSend[0] = 0;
5591 numGidsToSend[0] = it->second.size ();
5593 sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
5594 wait<int> (*comm, outArg (sizeReq));
5599 wait<int> (*comm, outArg (sizeReq));
5604 ArrayRCP<GO> myGids;
5605 int_type myNumGids = 0;
5607 GO* myGidsRaw = NULL;
5609 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
5610 if (it != pid2gids.end ()) {
5611 myGidsRaw = it->second.getRawPtr ();
5612 myNumGids = it->second.size ();
5614 myGids = arcp<GO> (myGidsRaw, 0, myNumGids,
false);
5618 myNumGids = numGidsToRecv[0];
5619 myGids = arcp<GO> (myNumGids);
5626 if (myNumGids > 0) {
5627 dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
5631 for (
int p = 1; p < numProcs; ++p) {
5633 ArrayRCP<GO> sendGids;
5634 GO* sendGidsRaw = NULL;
5635 int_type numSendGids = 0;
5637 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
5638 if (it != pid2gids.end ()) {
5639 numSendGids = it->second.size ();
5640 sendGidsRaw = it->second.getRawPtr ();
5641 sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids,
false);
5644 if (numSendGids > 0) {
5645 dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
5647 wait<int> (*comm, outArg (dataReq));
5649 else if (myRank == p) {
5651 wait<int> (*comm, outArg (dataReq));
5656 std::ostringstream os;
5657 os << myRank <<
": readMap: creating Map" << endl;
5660 const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
5661 RCP<const map_type> newMap;
5668 std::ostringstream errStrm;
5670 if (node.is_null ()) {
5671 newMap = rcp (
new map_type (INVALID, myGids (), indexBase, comm));
5674 newMap = rcp (
new map_type (INVALID, myGids (), indexBase, comm, node));
5677 catch (std::exception& e) {
5679 errStrm <<
"Process " << comm->getRank () <<
" threw an exception: " 5680 << e.what () << std::endl;
5684 errStrm <<
"Process " << comm->getRank () <<
" threw an exception " 5685 "not a subclass of std::exception" << std::endl;
5687 Teuchos::reduceAll<int, int> (*comm, Teuchos::REDUCE_MIN,
5688 lclSuccess, Teuchos::outArg (gblSuccess));
5689 if (gblSuccess != 1) {
5692 TEUCHOS_TEST_FOR_EXCEPTION(gblSuccess != 1, std::runtime_error,
"Map constructor failed!");
5694 if (err.is_null ()) {
5698 std::ostringstream os;
5699 os << myRank <<
": readMap: done" << endl;
5702 if (err.is_null ()) {
5721 encodeDataType (
const std::string& dataType)
5723 if (dataType ==
"real" || dataType ==
"integer") {
5725 }
else if (dataType ==
"complex") {
5727 }
else if (dataType ==
"pattern") {
5733 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
5734 "Unrecognized Matrix Market data type \"" << dataType
5735 <<
"\". We should never get here. " 5736 "Please report this bug to the Tpetra developers.");
5769 template<
class SparseMatrixType>
5774 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
5835 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5836 const std::string& matrixName,
5837 const std::string& matrixDescription,
5838 const bool debug=
false)
5840 TEUCHOS_TEST_FOR_EXCEPTION(
5841 pMatrix.is_null (), std::invalid_argument,
5842 "The input matrix is null.");
5843 Teuchos::RCP<const Teuchos::Comm<int> > comm = pMatrix->getComm ();
5844 TEUCHOS_TEST_FOR_EXCEPTION(
5845 comm.is_null (), std::invalid_argument,
5846 "The input matrix's communicator (Teuchos::Comm object) is null.");
5847 const int myRank = comm->getRank ();
5852 out.open (filename.c_str ());
5854 writeSparse (out, pMatrix, matrixName, matrixDescription, debug);
5882 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5883 const bool debug=
false)
5920 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5921 const std::string& matrixName,
5922 const std::string& matrixDescription,
5923 const bool debug=
false)
5925 using Teuchos::ArrayView;
5926 using Teuchos::Comm;
5927 using Teuchos::FancyOStream;
5928 using Teuchos::getFancyOStream;
5929 using Teuchos::null;
5931 using Teuchos::rcpFromRef;
5937 typedef typename Teuchos::ScalarTraits<ST> STS;
5938 typedef typename ArrayView<const LO>::const_iterator lo_iter;
5939 typedef typename ArrayView<const GO>::const_iterator go_iter;
5940 typedef typename ArrayView<const ST>::const_iterator st_iter;
5942 TEUCHOS_TEST_FOR_EXCEPTION(
5943 pMatrix.is_null (), std::invalid_argument,
5944 "The input matrix is null.");
5950 Teuchos::MatrixMarket::details::SetScientific<ST> sci (out);
5953 RCP<const Comm<int> > comm = pMatrix->getComm ();
5954 TEUCHOS_TEST_FOR_EXCEPTION(
5955 comm.is_null (), std::invalid_argument,
5956 "The input matrix's communicator (Teuchos::Comm object) is null.");
5957 const int myRank = comm->getRank ();
5960 RCP<FancyOStream> err =
5961 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
5963 std::ostringstream os;
5964 os << myRank <<
": writeSparse" << endl;
5967 os <<
"-- " << myRank <<
": past barrier" << endl;
5972 const bool debugPrint = debug && myRank == 0;
5974 RCP<const map_type> rowMap = pMatrix->getRowMap ();
5975 RCP<const map_type> colMap = pMatrix->getColMap ();
5976 RCP<const map_type> domainMap = pMatrix->getDomainMap ();
5977 RCP<const map_type> rangeMap = pMatrix->getRangeMap ();
5979 const global_size_t numRows = rangeMap->getGlobalNumElements ();
5980 const global_size_t numCols = domainMap->getGlobalNumElements ();
5982 if (debug && myRank == 0) {
5983 std::ostringstream os;
5984 os <<
"-- Input sparse matrix is:" 5985 <<
"---- " << numRows <<
" x " << numCols <<
" with " 5986 << pMatrix->getGlobalNumEntries() <<
" entries;" << endl
5988 << (pMatrix->isGloballyIndexed() ?
"Globally" :
"Locally")
5989 <<
" indexed." << endl
5990 <<
"---- Its row map has " << rowMap->getGlobalNumElements ()
5991 <<
" elements." << endl
5992 <<
"---- Its col map has " << colMap->getGlobalNumElements ()
5993 <<
" elements." << endl;
5998 const size_t localNumRows = (myRank == 0) ? numRows : 0;
5999 RCP<node_type> node = rowMap->getNode();
6001 std::ostringstream os;
6002 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6005 RCP<const map_type> gatherRowMap =
6006 Details::computeGatherMap (rowMap, err, debug);
6014 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6015 RCP<const map_type> gatherColMap =
6016 Details::computeGatherMap (colMap, err, debug);
6020 import_type importer (rowMap, gatherRowMap);
6025 RCP<sparse_matrix_type> newMatrix =
6027 static_cast<size_t> (0)));
6029 newMatrix->doImport (*pMatrix, importer,
INSERT);
6034 RCP<const map_type> gatherDomainMap =
6035 rcp (
new map_type (numCols, localNumCols,
6036 domainMap->getIndexBase (),
6038 RCP<const map_type> gatherRangeMap =
6039 rcp (
new map_type (numRows, localNumRows,
6040 rangeMap->getIndexBase (),
6042 newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
6046 cerr <<
"-- Output sparse matrix is:" 6047 <<
"---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
6049 << newMatrix->getDomainMap ()->getGlobalNumElements ()
6051 << newMatrix->getGlobalNumEntries () <<
" entries;" << endl
6053 << (newMatrix->isGloballyIndexed () ?
"Globally" :
"Locally")
6054 <<
" indexed." << endl
6055 <<
"---- Its row map has " 6056 << newMatrix->getRowMap ()->getGlobalNumElements ()
6057 <<
" elements, with index base " 6058 << newMatrix->getRowMap ()->getIndexBase () <<
"." << endl
6059 <<
"---- Its col map has " 6060 << newMatrix->getColMap ()->getGlobalNumElements ()
6061 <<
" elements, with index base " 6062 << newMatrix->getColMap ()->getIndexBase () <<
"." << endl
6063 <<
"---- Element count of output matrix's column Map may differ " 6064 <<
"from that of the input matrix's column Map, if some columns " 6065 <<
"of the matrix contain no entries." << endl;
6078 out <<
"%%MatrixMarket matrix coordinate " 6079 << (STS::isComplex ?
"complex" :
"real")
6080 <<
" general" << endl;
6083 if (matrixName !=
"") {
6084 printAsComment (out, matrixName);
6086 if (matrixDescription !=
"") {
6087 printAsComment (out, matrixDescription);
6097 out << newMatrix->getRangeMap ()->getGlobalNumElements () <<
" " 6098 << newMatrix->getDomainMap ()->getGlobalNumElements () <<
" " 6099 << newMatrix->getGlobalNumEntries () << endl;
6104 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6105 const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
6113 if (newMatrix->isGloballyIndexed()) {
6116 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6117 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6118 for (GO globalRowIndex = minAllGlobalIndex;
6119 globalRowIndex <= maxAllGlobalIndex;
6121 ArrayView<const GO> ind;
6122 ArrayView<const ST> val;
6123 newMatrix->getGlobalRowView (globalRowIndex, ind, val);
6124 go_iter indIter = ind.begin ();
6125 st_iter valIter = val.begin ();
6126 for (; indIter != ind.end() && valIter != val.end();
6127 ++indIter, ++valIter) {
6128 const GO globalColIndex = *indIter;
6131 out << (globalRowIndex + 1 - rowIndexBase) <<
" " 6132 << (globalColIndex + 1 - colIndexBase) <<
" ";
6133 if (STS::isComplex) {
6134 out << STS::real (*valIter) <<
" " << STS::imag (*valIter);
6142 typedef Teuchos::OrdinalTraits<GO> OTG;
6143 for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
6144 localRowIndex <= gatherRowMap->getMaxLocalIndex();
6147 const GO globalRowIndex =
6148 gatherRowMap->getGlobalElement (localRowIndex);
6149 TEUCHOS_TEST_FOR_EXCEPTION(
6150 globalRowIndex == OTG::invalid(), std::logic_error,
6151 "Failed to convert the supposed local row index " 6152 << localRowIndex <<
" into a global row index. " 6153 "Please report this bug to the Tpetra developers.");
6154 ArrayView<const LO> ind;
6155 ArrayView<const ST> val;
6156 newMatrix->getLocalRowView (localRowIndex, ind, val);
6157 lo_iter indIter = ind.begin ();
6158 st_iter valIter = val.begin ();
6159 for (; indIter != ind.end() && valIter != val.end();
6160 ++indIter, ++valIter) {
6162 const GO globalColIndex =
6163 newMatrix->getColMap()->getGlobalElement (*indIter);
6164 TEUCHOS_TEST_FOR_EXCEPTION(
6165 globalColIndex == OTG::invalid(), std::logic_error,
6166 "On local row " << localRowIndex <<
" of the sparse matrix: " 6167 "Failed to convert the supposed local column index " 6168 << *indIter <<
" into a global column index. Please report " 6169 "this bug to the Tpetra developers.");
6172 out << (globalRowIndex + 1 - rowIndexBase) <<
" " 6173 << (globalColIndex + 1 - colIndexBase) <<
" ";
6174 if (STS::isComplex) {
6175 out << STS::real (*valIter) <<
" " << STS::imag (*valIter);
6220 const std::string& graphName,
6221 const std::string& graphDescription,
6222 const bool debug=
false)
6224 using Teuchos::ArrayView;
6225 using Teuchos::Comm;
6226 using Teuchos::FancyOStream;
6227 using Teuchos::getFancyOStream;
6228 using Teuchos::null;
6230 using Teuchos::rcpFromRef;
6241 if (rowMap.is_null ()) {
6244 auto comm = rowMap->getComm ();
6245 if (comm.is_null ()) {
6248 const int myRank = comm->getRank ();
6251 RCP<FancyOStream> err =
6252 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6254 std::ostringstream os;
6255 os << myRank <<
": writeSparseGraph" << endl;
6258 os <<
"-- " << myRank <<
": past barrier" << endl;
6263 const bool debugPrint = debug && myRank == 0;
6270 const global_size_t numRows = rangeMap->getGlobalNumElements ();
6271 const global_size_t numCols = domainMap->getGlobalNumElements ();
6273 if (debug && myRank == 0) {
6274 std::ostringstream os;
6275 os <<
"-- Input sparse graph is:" 6276 <<
"---- " << numRows <<
" x " << numCols <<
" with " 6280 <<
" indexed." << endl
6281 <<
"---- Its row Map has " << rowMap->getGlobalNumElements ()
6282 <<
" elements." << endl
6283 <<
"---- Its col Map has " << colMap->getGlobalNumElements ()
6284 <<
" elements." << endl;
6289 const size_t localNumRows = (myRank == 0) ? numRows : 0;
6291 std::ostringstream os;
6292 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6295 auto gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
6303 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6304 auto gatherColMap = Details::computeGatherMap (colMap, err, debug);
6313 static_cast<size_t> (0));
6320 RCP<const map_type> gatherDomainMap =
6321 rcp (
new map_type (numCols, localNumCols,
6322 domainMap->getIndexBase (),
6324 RCP<const map_type> gatherRangeMap =
6325 rcp (
new map_type (numRows, localNumRows,
6326 rangeMap->getIndexBase (),
6328 newGraph.
fillComplete (gatherDomainMap, gatherRangeMap);
6332 cerr <<
"-- Output sparse graph is:" 6333 <<
"---- " << newGraph.
getRangeMap ()->getGlobalNumElements ()
6340 <<
" indexed." << endl
6341 <<
"---- Its row map has " 6342 << newGraph.
getRowMap ()->getGlobalNumElements ()
6343 <<
" elements, with index base " 6344 << newGraph.
getRowMap ()->getIndexBase () <<
"." << endl
6345 <<
"---- Its col map has " 6346 << newGraph.
getColMap ()->getGlobalNumElements ()
6347 <<
" elements, with index base " 6348 << newGraph.
getColMap ()->getIndexBase () <<
"." << endl
6349 <<
"---- Element count of output graph's column Map may differ " 6350 <<
"from that of the input matrix's column Map, if some columns " 6351 <<
"of the matrix contain no entries." << endl;
6365 out <<
"%%MatrixMarket matrix coordinate pattern general" << endl;
6368 if (graphName !=
"") {
6369 printAsComment (out, graphName);
6371 if (graphDescription !=
"") {
6372 printAsComment (out, graphDescription);
6383 out << newGraph.
getRangeMap ()->getGlobalNumElements () <<
" " 6384 << newGraph.
getDomainMap ()->getGlobalNumElements () <<
" " 6390 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6391 const GO colIndexBase = newGraph.
getColMap()->getIndexBase ();
6402 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6403 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6404 for (GO globalRowIndex = minAllGlobalIndex;
6405 globalRowIndex <= maxAllGlobalIndex;
6407 ArrayView<const GO> ind;
6409 for (
auto indIter = ind.begin (); indIter != ind.end (); ++indIter) {
6410 const GO globalColIndex = *indIter;
6413 out << (globalRowIndex + 1 - rowIndexBase) <<
" " 6414 << (globalColIndex + 1 - colIndexBase) <<
" ";
6420 typedef Teuchos::OrdinalTraits<GO> OTG;
6421 for (LO localRowIndex = gatherRowMap->getMinLocalIndex ();
6422 localRowIndex <= gatherRowMap->getMaxLocalIndex ();
6425 const GO globalRowIndex =
6426 gatherRowMap->getGlobalElement (localRowIndex);
6427 TEUCHOS_TEST_FOR_EXCEPTION
6428 (globalRowIndex == OTG::invalid (), std::logic_error,
"Failed " 6429 "to convert the supposed local row index " << localRowIndex <<
6430 " into a global row index. Please report this bug to the " 6431 "Tpetra developers.");
6432 ArrayView<const LO> ind;
6434 for (
auto indIter = ind.begin (); indIter != ind.end (); ++indIter) {
6436 const GO globalColIndex =
6437 newGraph.
getColMap ()->getGlobalElement (*indIter);
6438 TEUCHOS_TEST_FOR_EXCEPTION(
6439 globalColIndex == OTG::invalid(), std::logic_error,
6440 "On local row " << localRowIndex <<
" of the sparse graph: " 6441 "Failed to convert the supposed local column index " 6442 << *indIter <<
" into a global column index. Please report " 6443 "this bug to the Tpetra developers.");
6446 out << (globalRowIndex + 1 - rowIndexBase) <<
" " 6447 << (globalColIndex + 1 - colIndexBase) <<
" ";
6463 const bool debug=
false)
6505 const std::string& graphName,
6506 const std::string& graphDescription,
6507 const bool debug=
false)
6510 if (comm.is_null ()) {
6518 const int myRank = comm->getRank ();
6523 out.open (filename.c_str ());
6538 const bool debug=
false)
6553 const Teuchos::RCP<const crs_graph_type>& pGraph,
6554 const std::string& graphName,
6555 const std::string& graphDescription,
6556 const bool debug=
false)
6572 const Teuchos::RCP<const crs_graph_type>& pGraph,
6573 const bool debug=
false)
6602 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6603 const bool debug=
false)
6639 const std::string& matrixName,
6640 const std::string& matrixDescription,
6641 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6642 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6644 const int myRank = X.
getMap ().is_null () ? 0 :
6645 (X.
getMap ()->getComm ().is_null () ? 0 :
6646 X.
getMap ()->getComm ()->getRank ());
6650 out.open (filename.c_str());
6653 writeDense (out, X, matrixName, matrixDescription, err, dbg);
6666 const Teuchos::RCP<const multivector_type>& X,
6667 const std::string& matrixName,
6668 const std::string& matrixDescription,
6669 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6670 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6672 TEUCHOS_TEST_FOR_EXCEPTION(
6673 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 6674 "writeDenseFile: The input MultiVector X is null.");
6675 writeDenseFile (filename, *X, matrixName, matrixDescription, err, dbg);
6686 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6687 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6699 const Teuchos::RCP<const multivector_type>& X,
6700 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6701 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6703 TEUCHOS_TEST_FOR_EXCEPTION(
6704 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 6705 "writeDenseFile: The input MultiVector X is null.");
6743 const std::string& matrixName,
6744 const std::string& matrixDescription,
6745 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6746 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6748 using Teuchos::Comm;
6749 using Teuchos::outArg;
6750 using Teuchos::REDUCE_MAX;
6751 using Teuchos::reduceAll;
6755 RCP<const Comm<int> > comm = X.
getMap ().is_null () ?
6756 Teuchos::null : X.
getMap ()->getComm ();
6757 const int myRank = comm.is_null () ? 0 : comm->getRank ();
6762 const bool debug = ! dbg.is_null ();
6765 std::ostringstream os;
6766 os << myRank <<
": writeDense" << endl;
6771 writeDenseHeader (out, X, matrixName, matrixDescription, err, dbg);
6776 for (
size_t j = 0; j < numVecs; ++j) {
6777 writeDenseColumn (out, * (X.
getVector (j)), err, dbg);
6782 std::ostringstream os;
6783 os << myRank <<
": writeDense: Done" << endl;
6817 writeDenseHeader (std::ostream& out,
6819 const std::string& matrixName,
6820 const std::string& matrixDescription,
6821 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6822 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6824 using Teuchos::Comm;
6825 using Teuchos::outArg;
6827 using Teuchos::REDUCE_MAX;
6828 using Teuchos::reduceAll;
6831 typedef Teuchos::ScalarTraits<scalar_type> STS;
6832 const char prefix[] =
"Tpetra::MatrixMarket::writeDenseHeader: ";
6834 RCP<const Comm<int> > comm = X.getMap ().is_null () ?
6835 Teuchos::null : X.getMap ()->getComm ();
6836 const int myRank = comm.is_null () ? 0 : comm->getRank ();
6843 const bool debug = ! dbg.is_null ();
6847 std::ostringstream os;
6848 os << myRank <<
": writeDenseHeader" << endl;
6867 std::ostringstream hdr;
6869 std::string dataType;
6870 if (STS::isComplex) {
6871 dataType =
"complex";
6872 }
else if (STS::isOrdinal) {
6873 dataType =
"integer";
6877 hdr <<
"%%MatrixMarket matrix array " << dataType <<
" general" 6882 if (matrixName !=
"") {
6883 printAsComment (hdr, matrixName);
6885 if (matrixDescription !=
"") {
6886 printAsComment (hdr, matrixDescription);
6889 hdr << X.getGlobalLength () <<
" " << X.getNumVectors () << endl;
6893 }
catch (std::exception& e) {
6894 if (! err.is_null ()) {
6895 *err << prefix <<
"While writing the Matrix Market header, " 6896 "Process 0 threw an exception: " << e.what () << endl;
6905 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
6906 TEUCHOS_TEST_FOR_EXCEPTION(
6907 gblErr == 1, std::runtime_error, prefix <<
"Some error occurred " 6908 "which prevented this method from completing.");
6912 *dbg << myRank <<
": writeDenseHeader: Done" << endl;
6935 writeDenseColumn (std::ostream& out,
6937 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6938 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6940 using Teuchos::arcp;
6941 using Teuchos::Array;
6942 using Teuchos::ArrayRCP;
6943 using Teuchos::ArrayView;
6944 using Teuchos::Comm;
6945 using Teuchos::CommRequest;
6946 using Teuchos::ireceive;
6947 using Teuchos::isend;
6948 using Teuchos::outArg;
6949 using Teuchos::REDUCE_MAX;
6950 using Teuchos::reduceAll;
6952 using Teuchos::TypeNameTraits;
6953 using Teuchos::wait;
6956 typedef Teuchos::ScalarTraits<scalar_type> STS;
6958 const Comm<int>& comm = * (X.getMap ()->getComm ());
6959 const int myRank = comm.getRank ();
6960 const int numProcs = comm.getSize ();
6967 const bool debug = ! dbg.is_null ();
6971 std::ostringstream os;
6972 os << myRank <<
": writeDenseColumn" << endl;
6981 Teuchos::MatrixMarket::details::SetScientific<scalar_type> sci (out);
6983 const size_t myNumRows = X.getLocalLength ();
6984 const size_t numCols = X.getNumVectors ();
6987 const int sizeTag = 1337;
6988 const int dataTag = 1338;
7049 RCP<CommRequest<int> > sendReqSize, sendReqData;
7055 Array<ArrayRCP<size_t> > recvSizeBufs (3);
7056 Array<ArrayRCP<scalar_type> > recvDataBufs (3);
7057 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7058 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7061 ArrayRCP<size_t> sendDataSize (1);
7062 sendDataSize[0] = myNumRows;
7066 std::ostringstream os;
7067 os << myRank <<
": Post receive-size receives from " 7068 "Procs 1 and 2: tag = " << sizeTag << endl;
7072 recvSizeBufs[0].resize (1);
7076 (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7077 recvSizeBufs[1].resize (1);
7078 (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7079 recvSizeBufs[2].resize (1);
7080 (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7083 ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
7087 ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
7090 else if (myRank == 1 || myRank == 2) {
7092 std::ostringstream os;
7093 os << myRank <<
": Post send-size send: size = " 7094 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7101 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7105 std::ostringstream os;
7106 os << myRank <<
": Not posting my send-size send yet" << endl;
7115 std::ostringstream os;
7116 os << myRank <<
": Pack my entries" << endl;
7119 ArrayRCP<scalar_type> sendDataBuf;
7121 sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
7122 X.get1dCopy (sendDataBuf (), myNumRows);
7124 catch (std::exception& e) {
7126 if (! err.is_null ()) {
7127 std::ostringstream os;
7128 os <<
"Process " << myRank <<
": Attempt to pack my MultiVector " 7129 "entries threw an exception: " << e.what () << endl;
7134 std::ostringstream os;
7135 os << myRank <<
": Done packing my entries" << endl;
7144 *dbg << myRank <<
": Post send-data send: tag = " << dataTag
7147 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7155 std::ostringstream os;
7156 os << myRank <<
": Write my entries" << endl;
7162 const size_t printNumRows = myNumRows;
7163 ArrayView<const scalar_type> printData = sendDataBuf ();
7164 const size_t printStride = printNumRows;
7165 if (static_cast<size_t> (printData.size ()) < printStride * numCols) {
7167 if (! err.is_null ()) {
7168 std::ostringstream os;
7169 os <<
"Process " << myRank <<
": My MultiVector data's size " 7170 << printData.size () <<
" does not match my local dimensions " 7171 << printStride <<
" x " << numCols <<
"." << endl;
7179 for (
size_t col = 0; col < numCols; ++col) {
7180 for (
size_t row = 0; row < printNumRows; ++row) {
7181 if (STS::isComplex) {
7182 out << STS::real (printData[row + col * printStride]) <<
" " 7183 << STS::imag (printData[row + col * printStride]) << endl;
7185 out << printData[row + col * printStride] << endl;
7194 const int recvRank = 1;
7195 const int circBufInd = recvRank % 3;
7197 std::ostringstream os;
7198 os << myRank <<
": Wait on receive-size receive from Process " 7199 << recvRank << endl;
7203 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7207 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7208 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7210 if (! err.is_null ()) {
7211 std::ostringstream os;
7212 os << myRank <<
": Result of receive-size receive from Process " 7213 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() " 7214 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". " 7215 "This should never happen, and suggests that the receive never " 7216 "got posted. Please report this bug to the Tpetra developers." 7231 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7233 std::ostringstream os;
7234 os << myRank <<
": Post receive-data receive from Process " 7235 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 7236 << recvDataBufs[circBufInd].size () << endl;
7239 if (! recvSizeReqs[circBufInd].is_null ()) {
7241 if (! err.is_null ()) {
7242 std::ostringstream os;
7243 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 7244 "null, before posting the receive-data receive from Process " 7245 << recvRank <<
". This should never happen. Please report " 7246 "this bug to the Tpetra developers." << endl;
7250 recvDataReqs[circBufInd] =
7251 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7252 recvRank, dataTag, comm);
7255 else if (myRank == 1) {
7258 std::ostringstream os;
7259 os << myRank <<
": Wait on my send-size send" << endl;
7262 wait<int> (comm, outArg (sendReqSize));
7268 for (
int p = 1; p < numProcs; ++p) {
7270 if (p + 2 < numProcs) {
7272 const int recvRank = p + 2;
7273 const int circBufInd = recvRank % 3;
7275 std::ostringstream os;
7276 os << myRank <<
": Post receive-size receive from Process " 7277 << recvRank <<
": tag = " << sizeTag << endl;
7280 if (! recvSizeReqs[circBufInd].is_null ()) {
7282 if (! err.is_null ()) {
7283 std::ostringstream os;
7284 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 7285 <<
"null, for the receive-size receive from Process " 7286 << recvRank <<
"! This may mean that this process never " 7287 <<
"finished waiting for the receive from Process " 7288 << (recvRank - 3) <<
"." << endl;
7292 recvSizeReqs[circBufInd] =
7293 ireceive<int, size_t> (recvSizeBufs[circBufInd],
7294 recvRank, sizeTag, comm);
7297 if (p + 1 < numProcs) {
7298 const int recvRank = p + 1;
7299 const int circBufInd = recvRank % 3;
7303 std::ostringstream os;
7304 os << myRank <<
": Wait on receive-size receive from Process " 7305 << recvRank << endl;
7308 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7312 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7313 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7315 if (! err.is_null ()) {
7316 std::ostringstream os;
7317 os << myRank <<
": Result of receive-size receive from Process " 7318 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() " 7319 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". " 7320 "This should never happen, and suggests that the receive never " 7321 "got posted. Please report this bug to the Tpetra developers." 7335 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7337 std::ostringstream os;
7338 os << myRank <<
": Post receive-data receive from Process " 7339 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 7340 << recvDataBufs[circBufInd].size () << endl;
7343 if (! recvDataReqs[circBufInd].is_null ()) {
7345 if (! err.is_null ()) {
7346 std::ostringstream os;
7347 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not " 7348 <<
"null, for the receive-data receive from Process " 7349 << recvRank <<
"! This may mean that this process never " 7350 <<
"finished waiting for the receive from Process " 7351 << (recvRank - 3) <<
"." << endl;
7355 recvDataReqs[circBufInd] =
7356 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7357 recvRank, dataTag, comm);
7361 const int recvRank = p;
7362 const int circBufInd = recvRank % 3;
7364 std::ostringstream os;
7365 os << myRank <<
": Wait on receive-data receive from Process " 7366 << recvRank << endl;
7369 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
7376 std::ostringstream os;
7377 os << myRank <<
": Write entries from Process " << recvRank
7379 *dbg << os.str () << endl;
7381 size_t printNumRows = (recvSizeBufs[circBufInd])[0];
7382 if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7384 if (! err.is_null ()) {
7385 std::ostringstream os;
7386 os << myRank <<
": Result of receive-size receive from Process " 7387 << recvRank <<
" was Teuchos::OrdinalTraits<size_t>::" 7388 "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
7389 <<
". This should never happen, and suggests that its " 7390 "receive-size receive was never posted. " 7391 "Please report this bug to the Tpetra developers." << endl;
7398 if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
7400 if (! err.is_null ()) {
7401 std::ostringstream os;
7402 os << myRank <<
": Result of receive-size receive from Proc " 7403 << recvRank <<
" was " << printNumRows <<
" > 0, but " 7404 "recvDataBufs[" << circBufInd <<
"] is null. This should " 7405 "never happen. Please report this bug to the Tpetra " 7406 "developers." << endl;
7416 ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
7417 const size_t printStride = printNumRows;
7421 for (
size_t col = 0; col < numCols; ++col) {
7422 for (
size_t row = 0; row < printNumRows; ++row) {
7423 if (STS::isComplex) {
7424 out << STS::real (printData[row + col * printStride]) <<
" " 7425 << STS::imag (printData[row + col * printStride]) << endl;
7427 out << printData[row + col * printStride] << endl;
7432 else if (myRank == p) {
7435 std::ostringstream os;
7436 os << myRank <<
": Wait on my send-data send" << endl;
7439 wait<int> (comm, outArg (sendReqData));
7441 else if (myRank == p + 1) {
7444 std::ostringstream os;
7445 os << myRank <<
": Post send-data send: tag = " << dataTag
7449 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7452 std::ostringstream os;
7453 os << myRank <<
": Wait on my send-size send" << endl;
7456 wait<int> (comm, outArg (sendReqSize));
7458 else if (myRank == p + 2) {
7461 std::ostringstream os;
7462 os << myRank <<
": Post send-size send: size = " 7463 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7466 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7471 reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
7472 TEUCHOS_TEST_FOR_EXCEPTION(
7473 gblErr == 1, std::runtime_error,
"Tpetra::MatrixMarket::writeDense " 7474 "experienced some kind of error and was unable to complete.");
7478 *dbg << myRank <<
": writeDenseColumn: Done" << endl;
7492 const Teuchos::RCP<const multivector_type>& X,
7493 const std::string& matrixName,
7494 const std::string& matrixDescription,
7495 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7496 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7498 TEUCHOS_TEST_FOR_EXCEPTION(
7499 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 7500 "writeDense: The input MultiVector X is null.");
7501 writeDense (out, *X, matrixName, matrixDescription, err, dbg);
7512 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7513 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7525 const Teuchos::RCP<const multivector_type>& X,
7526 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7527 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7529 TEUCHOS_TEST_FOR_EXCEPTION(
7530 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 7531 "writeDense: The input MultiVector X is null.");
7557 Teuchos::RCP<Teuchos::FancyOStream> err =
7558 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
7573 const Teuchos::RCP<Teuchos::FancyOStream>& err,
7574 const bool debug=
false)
7576 using Teuchos::Array;
7577 using Teuchos::ArrayRCP;
7578 using Teuchos::ArrayView;
7579 using Teuchos::Comm;
7580 using Teuchos::CommRequest;
7581 using Teuchos::ireceive;
7582 using Teuchos::isend;
7584 using Teuchos::TypeNameTraits;
7585 using Teuchos::wait;
7588 typedef int pid_type;
7603 typedef ptrdiff_t int_type;
7604 TEUCHOS_TEST_FOR_EXCEPTION(
7605 sizeof (GO) >
sizeof (int_type), std::logic_error,
7606 "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
7607 <<
" is too big for ptrdiff_t. sizeof(GO) = " <<
sizeof (GO)
7608 <<
" > sizeof(ptrdiff_t) = " <<
sizeof (ptrdiff_t) <<
".");
7609 TEUCHOS_TEST_FOR_EXCEPTION(
7610 sizeof (pid_type) >
sizeof (int_type), std::logic_error,
7611 "The (MPI) process rank type pid_type=" <<
7612 TypeNameTraits<pid_type>::name () <<
" is too big for ptrdiff_t. " 7613 "sizeof(pid_type) = " <<
sizeof (pid_type) <<
" > sizeof(ptrdiff_t)" 7614 " = " <<
sizeof (ptrdiff_t) <<
".");
7616 const Comm<int>& comm = * (map.
getComm ());
7617 const int myRank = comm.getRank ();
7618 const int numProcs = comm.getSize ();
7620 if (! err.is_null ()) {
7624 std::ostringstream os;
7625 os << myRank <<
": writeMap" << endl;
7628 if (! err.is_null ()) {
7635 const int sizeTag = 1337;
7636 const int dataTag = 1338;
7695 RCP<CommRequest<int> > sendReqSize, sendReqData;
7701 Array<ArrayRCP<int_type> > recvSizeBufs (3);
7702 Array<ArrayRCP<int_type> > recvDataBufs (3);
7703 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7704 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7707 ArrayRCP<int_type> sendDataSize (1);
7708 sendDataSize[0] = myNumRows;
7712 std::ostringstream os;
7713 os << myRank <<
": Post receive-size receives from " 7714 "Procs 1 and 2: tag = " << sizeTag << endl;
7718 recvSizeBufs[0].resize (1);
7719 (recvSizeBufs[0])[0] = -1;
7720 recvSizeBufs[1].resize (1);
7721 (recvSizeBufs[1])[0] = -1;
7722 recvSizeBufs[2].resize (1);
7723 (recvSizeBufs[2])[0] = -1;
7726 ireceive<int, int_type> (recvSizeBufs[1], 1, sizeTag, comm);
7730 ireceive<int, int_type> (recvSizeBufs[2], 2, sizeTag, comm);
7733 else if (myRank == 1 || myRank == 2) {
7735 std::ostringstream os;
7736 os << myRank <<
": Post send-size send: size = " 7737 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7744 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
7748 std::ostringstream os;
7749 os << myRank <<
": Not posting my send-size send yet" << endl;
7760 std::ostringstream os;
7761 os << myRank <<
": Pack my GIDs and PIDs" << endl;
7765 ArrayRCP<int_type> sendDataBuf (myNumRows * 2);
7768 const int_type myMinGblIdx =
7770 for (
size_t k = 0; k < myNumRows; ++k) {
7771 const int_type gid = myMinGblIdx +
static_cast<int_type
> (k);
7772 const int_type pid =
static_cast<int_type
> (myRank);
7773 sendDataBuf[2*k] = gid;
7774 sendDataBuf[2*k+1] = pid;
7779 for (
size_t k = 0; k < myNumRows; ++k) {
7780 const int_type gid =
static_cast<int_type
> (myGblInds[k]);
7781 const int_type pid =
static_cast<int_type
> (myRank);
7782 sendDataBuf[2*k] = gid;
7783 sendDataBuf[2*k+1] = pid;
7788 std::ostringstream os;
7789 os << myRank <<
": Done packing my GIDs and PIDs" << endl;
7796 *err << myRank <<
": Post send-data send: tag = " << dataTag
7799 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
7804 *err << myRank <<
": Write MatrixMarket header" << endl;
7809 std::ostringstream hdr;
7813 hdr <<
"%%MatrixMarket matrix array integer general" << endl
7814 <<
"% Format: Version 2.0" << endl
7816 <<
"% This file encodes a Tpetra::Map." << endl
7817 <<
"% It is stored as a dense vector, with twice as many " << endl
7818 <<
"% entries as the global number of GIDs (global indices)." << endl
7819 <<
"% (GID, PID) pairs are stored contiguously, where the PID " << endl
7820 <<
"% is the rank of the process owning that GID." << endl
7825 std::ostringstream os;
7826 os << myRank <<
": Write my GIDs and PIDs" << endl;
7832 const int_type printNumRows = myNumRows;
7833 ArrayView<const int_type> printData = sendDataBuf ();
7834 for (int_type k = 0; k < printNumRows; ++k) {
7835 const int_type gid = printData[2*k];
7836 const int_type pid = printData[2*k+1];
7837 out << gid << endl << pid << endl;
7843 const int recvRank = 1;
7844 const int circBufInd = recvRank % 3;
7846 std::ostringstream os;
7847 os << myRank <<
": Wait on receive-size receive from Process " 7848 << recvRank << endl;
7852 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7856 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
7857 if (debug && recvNumRows == -1) {
7858 std::ostringstream os;
7859 os << myRank <<
": Result of receive-size receive from Process " 7860 << recvRank <<
" is -1. This should never happen, and " 7861 "suggests that the receive never got posted. Please report " 7862 "this bug to the Tpetra developers." << endl;
7867 recvDataBufs[circBufInd].resize (recvNumRows * 2);
7869 std::ostringstream os;
7870 os << myRank <<
": Post receive-data receive from Process " 7871 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 7872 << recvDataBufs[circBufInd].size () << endl;
7875 if (! recvSizeReqs[circBufInd].is_null ()) {
7876 std::ostringstream os;
7877 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 7878 "null, before posting the receive-data receive from Process " 7879 << recvRank <<
". This should never happen. Please report " 7880 "this bug to the Tpetra developers." << endl;
7883 recvDataReqs[circBufInd] =
7884 ireceive<int, int_type> (recvDataBufs[circBufInd],
7885 recvRank, dataTag, comm);
7888 else if (myRank == 1) {
7891 std::ostringstream os;
7892 os << myRank <<
": Wait on my send-size send" << endl;
7895 wait<int> (comm, outArg (sendReqSize));
7901 for (
int p = 1; p < numProcs; ++p) {
7903 if (p + 2 < numProcs) {
7905 const int recvRank = p + 2;
7906 const int circBufInd = recvRank % 3;
7908 std::ostringstream os;
7909 os << myRank <<
": Post receive-size receive from Process " 7910 << recvRank <<
": tag = " << sizeTag << endl;
7913 if (! recvSizeReqs[circBufInd].is_null ()) {
7914 std::ostringstream os;
7915 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 7916 <<
"null, for the receive-size receive from Process " 7917 << recvRank <<
"! This may mean that this process never " 7918 <<
"finished waiting for the receive from Process " 7919 << (recvRank - 3) <<
"." << endl;
7922 recvSizeReqs[circBufInd] =
7923 ireceive<int, int_type> (recvSizeBufs[circBufInd],
7924 recvRank, sizeTag, comm);
7927 if (p + 1 < numProcs) {
7928 const int recvRank = p + 1;
7929 const int circBufInd = recvRank % 3;
7933 std::ostringstream os;
7934 os << myRank <<
": Wait on receive-size receive from Process " 7935 << recvRank << endl;
7938 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7942 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
7943 if (debug && recvNumRows == -1) {
7944 std::ostringstream os;
7945 os << myRank <<
": Result of receive-size receive from Process " 7946 << recvRank <<
" is -1. This should never happen, and " 7947 "suggests that the receive never got posted. Please report " 7948 "this bug to the Tpetra developers." << endl;
7953 recvDataBufs[circBufInd].resize (recvNumRows * 2);
7955 std::ostringstream os;
7956 os << myRank <<
": Post receive-data receive from Process " 7957 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 7958 << recvDataBufs[circBufInd].size () << endl;
7961 if (! recvDataReqs[circBufInd].is_null ()) {
7962 std::ostringstream os;
7963 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not " 7964 <<
"null, for the receive-data receive from Process " 7965 << recvRank <<
"! This may mean that this process never " 7966 <<
"finished waiting for the receive from Process " 7967 << (recvRank - 3) <<
"." << endl;
7970 recvDataReqs[circBufInd] =
7971 ireceive<int, int_type> (recvDataBufs[circBufInd],
7972 recvRank, dataTag, comm);
7976 const int recvRank = p;
7977 const int circBufInd = recvRank % 3;
7979 std::ostringstream os;
7980 os << myRank <<
": Wait on receive-data receive from Process " 7981 << recvRank << endl;
7984 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
7991 std::ostringstream os;
7992 os << myRank <<
": Write GIDs and PIDs from Process " 7993 << recvRank << endl;
7994 *err << os.str () << endl;
7996 const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
7997 if (debug && printNumRows == -1) {
7998 std::ostringstream os;
7999 os << myRank <<
": Result of receive-size receive from Process " 8000 << recvRank <<
" was -1. This should never happen, and " 8001 "suggests that its receive-size receive was never posted. " 8002 "Please report this bug to the Tpetra developers." << endl;
8005 if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
8006 std::ostringstream os;
8007 os << myRank <<
": Result of receive-size receive from Proc " 8008 << recvRank <<
" was " << printNumRows <<
" > 0, but " 8009 "recvDataBufs[" << circBufInd <<
"] is null. This should " 8010 "never happen. Please report this bug to the Tpetra " 8011 "developers." << endl;
8014 ArrayView<const int_type> printData = (recvDataBufs[circBufInd]) ();
8015 for (int_type k = 0; k < printNumRows; ++k) {
8016 const int_type gid = printData[2*k];
8017 const int_type pid = printData[2*k+1];
8018 out << gid << endl << pid << endl;
8021 else if (myRank == p) {
8024 std::ostringstream os;
8025 os << myRank <<
": Wait on my send-data send" << endl;
8028 wait<int> (comm, outArg (sendReqData));
8030 else if (myRank == p + 1) {
8033 std::ostringstream os;
8034 os << myRank <<
": Post send-data send: tag = " << dataTag
8038 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
8041 std::ostringstream os;
8042 os << myRank <<
": Wait on my send-size send" << endl;
8045 wait<int> (comm, outArg (sendReqSize));
8047 else if (myRank == p + 2) {
8050 std::ostringstream os;
8051 os << myRank <<
": Post send-size send: size = " 8052 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
8055 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
8059 if (! err.is_null ()) {
8063 *err << myRank <<
": writeMap: Done" << endl;
8065 if (! err.is_null ()) {
8075 const int myRank = map.
getComm ()->getRank ();
8078 out.open (filename.c_str());
8111 printAsComment (std::ostream& out,
const std::string& str)
8114 std::istringstream inpstream (str);
8117 while (getline (inpstream, line)) {
8118 if (! line.empty()) {
8121 if (line[0] ==
'%') {
8122 out << line << endl;
8125 out <<
"%% " << line << endl;
8153 Teuchos::ParameterList pl;
8179 Teuchos::ParameterList pl;
8222 const Teuchos::ParameterList& params)
8225 std::string tmpFile =
"__TMP__" + fileName;
8226 const int myRank = A.
getDomainMap()->getComm()->getRank();
8227 bool precisionChanged=
false;
8237 if (std::ifstream(tmpFile))
8238 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
8239 "writeOperator: temporary file " << tmpFile <<
" already exists");
8240 out.open(tmpFile.c_str());
8241 if (params.isParameter(
"precision")) {
8242 oldPrecision = out.precision(params.get<
int>(
"precision"));
8243 precisionChanged=
true;
8247 const std::string header = writeOperatorImpl(out, A, params);
8250 if (precisionChanged)
8251 out.precision(oldPrecision);
8253 out.open(fileName.c_str(), std::ios::binary);
8254 bool printMatrixMarketHeader =
true;
8255 if (params.isParameter(
"print MatrixMarket header"))
8256 printMatrixMarketHeader = params.get<
bool>(
"print MatrixMarket header");
8257 if (printMatrixMarketHeader && myRank == 0) {
8262 std::ifstream src(tmpFile, std::ios_base::binary);
8266 remove(tmpFile.c_str());
8311 const Teuchos::ParameterList& params)
8313 const int myRank = A.
getDomainMap ()->getComm ()->getRank ();
8330 std::ostringstream tmpOut;
8332 if (params.isParameter (
"precision") && params.isType<
int> (
"precision")) {
8333 (void) tmpOut.precision (params.get<
int> (
"precision"));
8337 const std::string header = writeOperatorImpl (tmpOut, A, params);
8340 bool printMatrixMarketHeader =
true;
8341 if (params.isParameter (
"print MatrixMarket header") &&
8342 params.isType<
bool> (
"print MatrixMarket header")) {
8343 printMatrixMarketHeader = params.get<
bool> (
"print MatrixMarket header");
8345 if (printMatrixMarketHeader && myRank == 0) {
8361 out << tmpOut.str ();
8375 writeOperatorImpl (std::ostream& os,
8376 const operator_type& A,
8377 const Teuchos::ParameterList& params)
8381 using Teuchos::ArrayRCP;
8382 using Teuchos::Array;
8387 typedef Teuchos::OrdinalTraits<LO> TLOT;
8388 typedef Teuchos::OrdinalTraits<GO> TGOT;
8392 const map_type& domainMap = *(A.getDomainMap());
8393 RCP<const map_type> rangeMap = A.getRangeMap();
8394 RCP<const Teuchos::Comm<int> > comm = rangeMap->getComm();
8395 const int myRank = comm->getRank();
8396 const size_t numProcs = comm->getSize();
8399 if (params.isParameter(
"probing size"))
8400 numMVs = params.get<
int>(
"probing size");
8403 GO minColGid = domainMap.getMinAllGlobalIndex();
8404 GO maxColGid = domainMap.getMaxAllGlobalIndex();
8409 GO numGlobElts = maxColGid - minColGid + TGOT::one();
8410 GO numChunks = numGlobElts / numMVs;
8411 GO rem = numGlobElts % numMVs;
8412 GO indexBase = rangeMap->getIndexBase();
8414 int offsetToUseInPrinting = 1 - indexBase;
8415 if (params.isParameter(
"zero-based indexing")) {
8416 if (params.get<
bool>(
"zero-based indexing") ==
true)
8417 offsetToUseInPrinting = -indexBase;
8421 size_t numLocalRangeEntries = rangeMap->getNodeNumElements();
8424 RCP<const map_type> allGidsMap = rcp(
new map_type(TGOT::invalid(), numLocalRangeEntries,
8427 mv_type_go allGids(allGidsMap,1);
8428 Teuchos::ArrayRCP<GO> allGidsData = allGids.getDataNonConst(0);
8430 for (
size_t i=0; i<numLocalRangeEntries; i++)
8431 allGidsData[i] = rangeMap->getGlobalElement(i);
8432 allGidsData = Teuchos::null;
8435 GO numTargetMapEntries=TGOT::zero();
8436 Teuchos::Array<GO> importGidList;
8438 numTargetMapEntries = rangeMap->getGlobalNumElements();
8439 importGidList.reserve(numTargetMapEntries);
8440 for (GO j=0; j<numTargetMapEntries; ++j) importGidList.push_back(j + indexBase);
8442 importGidList.reserve(numTargetMapEntries);
8444 RCP<map_type> importGidMap = rcp(
new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
8447 import_type gidImporter(allGidsMap, importGidMap);
8448 mv_type_go importedGids(importGidMap, 1);
8449 importedGids.doImport(allGids, gidImporter,
INSERT);
8452 ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
8453 RCP<const map_type> importMap = rcp(
new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm) );
8456 import_type importer(rangeMap, importMap);
8458 RCP<mv_type> colsOnPid0 = rcp(
new mv_type(importMap,numMVs));
8460 RCP<mv_type> ei = rcp(
new mv_type(A.getDomainMap(),numMVs));
8461 RCP<mv_type> colsA = rcp(
new mv_type(A.getRangeMap(),numMVs));
8463 Array<GO> globalColsArray, localColsArray;
8464 globalColsArray.reserve(numMVs);
8465 localColsArray.reserve(numMVs);
8467 ArrayRCP<ArrayRCP<Scalar> > eiData(numMVs);
8468 for (
size_t i=0; i<numMVs; ++i)
8469 eiData[i] = ei->getDataNonConst(i);
8474 for (GO k=0; k<numChunks; ++k) {
8475 for (
size_t j=0; j<numMVs; ++j ) {
8477 GO curGlobalCol = minColGid + k*numMVs + j;
8478 globalColsArray.push_back(curGlobalCol);
8480 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8481 if (curLocalCol != TLOT::invalid()) {
8482 eiData[j][curLocalCol] = TGOT::one();
8483 localColsArray.push_back(curLocalCol);
8489 A.apply(*ei,*colsA);
8491 colsOnPid0->doImport(*colsA,importer,
INSERT);
8494 globalNnz += writeColumns(os,*colsOnPid0, numMVs, importedGidsData(),
8495 globalColsArray, offsetToUseInPrinting);
8498 for (
size_t j=0; j<numMVs; ++j ) {
8499 for (
int i=0; i<localColsArray.size(); ++i)
8500 eiData[j][localColsArray[i]] = TGOT::zero();
8502 globalColsArray.clear();
8503 localColsArray.clear();
8511 for (
int j=0; j<rem; ++j ) {
8512 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8513 globalColsArray.push_back(curGlobalCol);
8515 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8516 if (curLocalCol != TLOT::invalid()) {
8517 eiData[j][curLocalCol] = TGOT::one();
8518 localColsArray.push_back(curLocalCol);
8524 A.apply(*ei,*colsA);
8526 colsOnPid0->doImport(*colsA,importer,
INSERT);
8528 globalNnz += writeColumns(os,*colsOnPid0, rem, importedGidsData(),
8529 globalColsArray, offsetToUseInPrinting);
8532 for (
int j=0; j<rem; ++j ) {
8533 for (
int i=0; i<localColsArray.size(); ++i)
8534 eiData[j][localColsArray[i]] = TGOT::zero();
8536 globalColsArray.clear();
8537 localColsArray.clear();
8546 std::ostringstream oss;
8548 oss <<
"%%MatrixMarket matrix coordinate ";
8549 if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
8554 oss <<
" general" << std::endl;
8555 oss <<
"% Tpetra::Operator" << std::endl;
8556 std::time_t now = std::time(NULL);
8557 oss <<
"% time stamp: " << ctime(&now);
8558 oss <<
"% written from " << numProcs <<
" processes" << std::endl;
8559 size_t numRows = rangeMap->getGlobalNumElements();
8560 size_t numCols = domainMap.getGlobalNumElements();
8561 oss << numRows <<
" " << numCols <<
" " << globalNnz << std::endl;
8568 writeColumns(std::ostream& os, mv_type
const &colsA,
size_t const &numCols,
8569 Teuchos::ArrayView<const global_ordinal_type>
const &rowGids,
8570 Teuchos::Array<global_ordinal_type>
const &colsArray,
8575 typedef Teuchos::ScalarTraits<Scalar> STS;
8578 const Scalar zero = STS::zero();
8579 const size_t numRows = colsA.getGlobalLength();
8580 for (
size_t j=0; j<numCols; ++j) {
8581 Teuchos::ArrayRCP<const Scalar>
const curCol = colsA.getData(j);
8582 const GO J = colsArray[j];
8583 for (
size_t i=0; i<numRows; ++i) {
8584 const Scalar val = curCol[i];
8586 os << rowGids[i]+indexBase <<
" " << J+indexBase <<
" " << val << std::endl;
8603 #endif // __MatrixMarket_Tpetra_hpp static Teuchos::RCP< vector_type > readVectorFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read a Vector from the given Matrix Market file.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market input stream.
static void writeSparseGraphFile(const std::string &filename, const Teuchos::RCP< const crs_graph_type > &pGraph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), with no comments...
static void writeOperator(std::ostream &out, const operator_type &A)
Write a Tpetra::Operator to an output stream.
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market file.
static Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< node_type > &node, const bool tolerant=false, const bool debug=false)
Variant of readMap (above) that takes an explicit Node instance.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
SparseMatrixType::node_type node_type
The fourth template parameter of CrsMatrix and MultiVector.
static void writeDense(std::ostream &out, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
global_size_t getGlobalNumEntries() const
Returns the global number of entries in the graph.
Teuchos::RCP< const map_type > getRowMap() const
Returns the Map that describes the row distribution in this graph.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of readSparseFile above that takes a Node object.
SparseMatrixType::local_ordinal_type local_ordinal_type
SparseMatrixType::global_ordinal_type global_ordinal_type
Teuchos::RCP< const map_type > getRangeMap() const
Returns the Map associated with the domain of this graph.
One or more distributed dense vectors.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream.
SparseMatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the sparse matrix.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
bool isContiguous() const
True if this Map is distributed contiguously, else false.
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of readSparseFile above that takes a Node object.
Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > getVector(const size_t j) const
Return a Vector which is a const view of column j.
static void writeDenseFile(const std::string &filename, const Teuchos::RCP< const multivector_type > &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
SparseMatrixType::scalar_type scalar_type
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of the above readSparse() method that takes a Kokkos Node.
void getLocalRowView(LocalOrdinal LocalRow, Teuchos::ArrayView< const LocalOrdinal > &indices) const
Get a const, non-persisting view of the given local row's local column indices, as a Teuchos::ArrayVi...
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market file, with provided Maps.
void fillComplete(const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Tell the graph that you are done changing its structure.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the communicator.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
Specialization of Tpetra::MultiVector that matches SparseMatrixType.
global_size_t getGlobalNumElements() const
The number of elements in this Map.
static Teuchos::RCP< const map_type > readMapFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given Matrix Market file.
static void writeMap(std::ostream &out, const map_type &map, const bool debug=false)
Print the Map to the given output stream.
static void writeDenseFile(const std::string &filename, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static void writeDenseFile(const std::string &filename, const Teuchos::RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static Teuchos::RCP< const map_type > readMapFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< node_type > &node, const bool tolerant=false, const bool debug=false)
Variant of readMapFile (above) that takes an explicit Node instance.
size_t getNodeNumElements() const
The number of elements belonging to the calling process.
SparseMatrixType::global_ordinal_type global_ordinal_type
Type of indices as read from the Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream, with provided Maps. ...
static Teuchos::RCP< multivector_type > readDense(std::istream &in, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market input stream.
static Teuchos::RCP< multivector_type > readDenseFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< node_type > &node, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Variant of readDenseMatrix (see above) that takes a Node.
void gathervPrint(std::ostream &out, const std::string &s, const Teuchos::Comm< int > &comm)
On Process 0 in the given communicator, print strings from each process in that communicator, in rank order.
Teuchos::RCP< const map_type > getDomainMap() const
Returns the Map associated with the domain of this graph.
static void writeSparseGraph(std::ostream &out, const crs_graph_type &graph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given output stream, with no comments...
size_t global_size_t
Global size_t object.
static Teuchos::RCP< vector_type > readVector(std::istream &in, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream.
SparseMatrixType::node_type node_type
The Kokkos Node type; fourth template parameter of Tpetra::CrsMatrix.
static void writeOperator(const std::string &fileName, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to a file, with options.
static Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Specialization of Tpetra::Map that matches SparseMatrixType.
Insert new values that don't currently exist.
SparseMatrixType::scalar_type scalar_type
Type of the entries of the sparse matrix.
Scalar scalar_type
This class' first template parameter; the type of each entry in the MultiVector.
static void writeMapFile(const std::string &filename, const map_type &map)
Write the Map to the given file.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of the above readSparseGraph() method that takes a Kokkos Node.
static Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< node_type > &node, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool tolerant=false, const bool debug=false)
Variant of readMap (above) that takes an explicit Node instance.
Abstract interface for operators (e.g., matrices and preconditioners).
static void writeSparseGraphFile(const std::string &filename, const Teuchos::RCP< const crs_graph_type > &pGraph, const std::string &graphName, const std::string &graphDescription, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), taking the graph by T...
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparseGraph() above that takes a Node object.
CrsGraph< local_ordinal_type, global_ordinal_type, node_type > crs_graph_type
Specialization of Tpetra::CrsGraph that matches SparseMatrixType.
static void writeDense(std::ostream &out, const Teuchos::RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file, with provided Maps.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
From a distributed map build a map with all GIDs on the root node.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market file.
static void writeMap(std::ostream &out, const map_type &map, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool debug=false)
Print the Map to the given output stream out.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
void getGlobalRowView(GlobalOrdinal GlobalRow, Teuchos::ArrayView< const GlobalOrdinal > &Indices) const
Get a const, non-persisting view of the given global row's global column indices, as a Teuchos::Array...
static Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream, with optional debugging output stream...
static Teuchos::RCP< multivector_type > readDenseFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market file.
static void writeSparseGraphFile(const std::string &filename, const crs_graph_type &graph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), with no comments...
static void writeDense(std::ostream &out, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static void writeSparse(std::ostream &out, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a NONOWNING view of the global indices owned by this process.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
The MultiVector specialization associated with SparseMatrixType.
GlobalOrdinal getMinGlobalIndex() const
The minimum global index owned by the calling process.
static void writeSparseGraphFile(const std::string &filename, const crs_graph_type &graph, const std::string &graphName, const std::string &graphDescription, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename).
static Teuchos::RCP< vector_type > readVectorFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< node_type > &node, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Like readVectorFile() (see above), but with a supplied Node object.
bool isGloballyIndexed() const
If graph indices are in the global range, this function returns true. Otherwise, this function return...
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > vector_type
The Vector specialization associated with SparseMatrixType.
Describes a parallel distribution of objects over processes.
Teuchos::RCP< const map_type > getColMap() const
Returns the Map that describes the column distribution in this graph.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market input stream.
size_t getNumVectors() const
Number of columns in the multivector.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparse() above that takes a Node object.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream.
A distributed dense vector.
static void writeSparseGraph(std::ostream &out, const crs_graph_type &graph, const std::string &graphName, const std::string &graphDescription, const bool debug=false)
Print the sparse graph in Matrix Market format to the given output stream.
CrsGraph< local_ordinal_type, global_ordinal_type, node_type > sparse_graph_type
The CrsGraph specialization associated with SparseMatrixType.
SparseMatrixType sparse_matrix_type
This class' template parameter; a specialization of CrsMatrix.
SparseMatrixType sparse_matrix_type
Template parameter of this class; specialization of CrsMatrix.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream, with provided Maps. ...
static void writeDenseFile(const std::string &filename, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparseGraph that takes a Node object.
static Teuchos::RCP< multivector_type > readDense(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< node_type > &node, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Variant of readDense (see above) that takes a Node.
static void writeOperator(std::ostream &out, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to an output stream, with options.
Matrix Market file reader for CrsMatrix and MultiVector.
static Teuchos::RCP< vector_type > readVector(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< node_type > &node, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream, with a supplied Node.
Matrix Market file writer for CrsMatrix and MultiVector.
static void writeSparse(std::ostream &out, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
static void writeDense(std::ostream &out, const Teuchos::RCP< const multivector_type > &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and or description.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparseFile that takes a Node object.
Matrix Market file readers and writers for sparse and dense matrices (as CrsMatrix resp...
static void writeOperator(const std::string &fileName, operator_type const &A)
Write a Tpetra::Operator to a file.
void doImport(const SrcDistObject &source, const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, CombineMode CM)
Import data into this object using an Import object ("forward mode").