42 #ifndef __MatrixMarket_Tpetra_hpp
43 #define __MatrixMarket_Tpetra_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_SetScientific.hpp"
69 #include "Teuchos_TimeMonitor.hpp"
72 #include "mmio_Tpetra.h"
74 #include "Tpetra_Distribution.hpp"
115 namespace MatrixMarket {
171 template<
class SparseMatrixType>
176 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
191 typedef typename SparseMatrixType::global_ordinal_type
224 typedef Teuchos::ArrayRCP<int>::size_type size_type;
236 static Teuchos::RCP<const map_type>
241 return rcp (
new map_type (static_cast<global_size_t> (numRows),
242 static_cast<global_ordinal_type> (0),
243 pComm, GloballyDistributed));
273 static Teuchos::RCP<const map_type>
274 makeRowMap (
const Teuchos::RCP<const map_type>& pRowMap,
280 if (pRowMap.is_null ()) {
281 return rcp (
new map_type (static_cast<global_size_t> (numRows),
282 static_cast<global_ordinal_type> (0),
283 pComm, GloballyDistributed));
286 TEUCHOS_TEST_FOR_EXCEPTION
287 (! pRowMap->isDistributed () && pComm->getSize () > 1,
288 std::invalid_argument,
"The specified row map is not distributed, "
289 "but the given communicator includes more than one process (in "
290 "fact, there are " << pComm->getSize () <<
" processes).");
291 TEUCHOS_TEST_FOR_EXCEPTION
292 (pRowMap->getComm () != pComm, std::invalid_argument,
293 "The specified row Map's communicator (pRowMap->getComm()) "
294 "differs from the given separately supplied communicator pComm.");
313 static Teuchos::RCP<const map_type>
314 makeDomainMap (
const Teuchos::RCP<const map_type>& pRangeMap,
323 if (numRows == numCols) {
326 return createUniformContigMapWithNode<LO,GO,NT> (numCols,
327 pRangeMap->getComm ());
404 distribute (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
405 Teuchos::ArrayRCP<size_t>& myRowPtr,
406 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
407 Teuchos::ArrayRCP<scalar_type>& myValues,
408 const Teuchos::RCP<const map_type>& pRowMap,
409 Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
410 Teuchos::ArrayRCP<size_t>& rowPtr,
411 Teuchos::ArrayRCP<global_ordinal_type>& colInd,
412 Teuchos::ArrayRCP<scalar_type>& values,
413 const bool debug=
false)
416 using Teuchos::ArrayRCP;
417 using Teuchos::ArrayView;
420 using Teuchos::CommRequest;
423 using Teuchos::receive;
428 const bool extraDebug =
false;
430 const int numProcs = pComm->getSize ();
431 const int myRank = pComm->getRank ();
432 const int rootRank = 0;
439 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
440 const size_type myNumRows = myRows.size();
441 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
442 pRowMap->getLocalNumElements(),
444 "pRowMap->getLocalElementList().size() = "
446 <<
" != pRowMap->getLocalNumElements() = "
447 << pRowMap->getLocalNumElements() <<
". "
448 "Please report this bug to the Tpetra developers.");
449 TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
451 "On Proc 0: numEntriesPerRow.size() = "
452 << numEntriesPerRow.size()
453 <<
" != pRowMap->getLocalElementList().size() = "
454 << myNumRows <<
". Please report this bug to the "
455 "Tpetra developers.");
459 myNumEntriesPerRow = arcp<size_t> (myNumRows);
461 if (myRank != rootRank) {
465 send (*pComm, myNumRows, rootRank);
466 if (myNumRows != 0) {
470 send (*pComm, static_cast<int> (myNumRows),
471 myRows.getRawPtr(), rootRank);
481 receive (*pComm, rootRank,
482 static_cast<int> (myNumRows),
483 myNumEntriesPerRow.getRawPtr());
488 std::accumulate (myNumEntriesPerRow.begin(),
489 myNumEntriesPerRow.end(), 0);
495 myColInd = arcp<GO> (myNumEntries);
496 myValues = arcp<scalar_type> (myNumEntries);
497 if (myNumEntries > 0) {
500 receive (*pComm, rootRank,
501 static_cast<int> (myNumEntries),
502 myColInd.getRawPtr());
503 receive (*pComm, rootRank,
504 static_cast<int> (myNumEntries),
505 myValues.getRawPtr());
511 cerr <<
"-- Proc 0: Copying my data from global arrays" << endl;
515 for (size_type k = 0; k < myNumRows; ++k) {
516 const GO myCurRow = myRows[k];
518 myNumEntriesPerRow[k] = numEntriesInThisRow;
520 if (extraDebug && debug) {
521 cerr <<
"Proc " << pRowMap->getComm ()->getRank ()
522 <<
": myNumEntriesPerRow[0.." << (myNumRows-1) <<
"] = [";
523 for (size_type k = 0; k < myNumRows; ++k) {
524 cerr << myNumEntriesPerRow[k];
525 if (k < myNumRows-1) {
533 std::accumulate (myNumEntriesPerRow.begin(),
534 myNumEntriesPerRow.end(), 0);
536 cerr <<
"-- Proc 0: I own " << myNumRows <<
" rows and "
537 << myNumEntries <<
" entries" << endl;
539 myColInd = arcp<GO> (myNumEntries);
540 myValues = arcp<scalar_type> (myNumEntries);
548 for (size_type k = 0; k < myNumRows;
549 myCurPos += myNumEntriesPerRow[k], ++k) {
551 const GO myRow = myRows[k];
552 const size_t curPos = rowPtr[myRow];
555 if (curNumEntries > 0) {
556 ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
557 ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
558 std::copy (colIndView.begin(), colIndView.end(),
559 myColIndView.begin());
561 ArrayView<scalar_type> valuesView =
562 values (curPos, curNumEntries);
563 ArrayView<scalar_type> myValuesView =
564 myValues (myCurPos, curNumEntries);
565 std::copy (valuesView.begin(), valuesView.end(),
566 myValuesView.begin());
571 for (
int p = 1; p < numProcs; ++p) {
573 cerr <<
"-- Proc 0: Processing proc " << p << endl;
576 size_type theirNumRows = 0;
581 receive (*pComm, p, &theirNumRows);
583 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
584 << theirNumRows <<
" rows" << endl;
586 if (theirNumRows != 0) {
591 ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
592 receive (*pComm, p, as<int> (theirNumRows),
593 theirRows.getRawPtr ());
602 const global_size_t numRows = pRowMap->getGlobalNumElements ();
603 const GO indexBase = pRowMap->getIndexBase ();
604 bool theirRowsValid =
true;
605 for (size_type k = 0; k < theirNumRows; ++k) {
606 if (theirRows[k] < indexBase ||
607 as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
608 theirRowsValid =
false;
611 if (! theirRowsValid) {
612 TEUCHOS_TEST_FOR_EXCEPTION(
613 ! theirRowsValid, std::logic_error,
614 "Proc " << p <<
" has at least one invalid row index. "
615 "Here are all of them: " <<
616 Teuchos::toString (theirRows ()) <<
". Valid row index "
617 "range (zero-based): [0, " << (numRows - 1) <<
"].");
632 ArrayRCP<size_t> theirNumEntriesPerRow;
633 theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
634 for (size_type k = 0; k < theirNumRows; ++k) {
635 theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
642 send (*pComm, static_cast<int> (theirNumRows),
643 theirNumEntriesPerRow.getRawPtr(), p);
647 std::accumulate (theirNumEntriesPerRow.begin(),
648 theirNumEntriesPerRow.end(), 0);
651 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
652 << theirNumEntries <<
" entries" << endl;
657 if (theirNumEntries == 0) {
666 ArrayRCP<GO> theirColInd (theirNumEntries);
667 ArrayRCP<scalar_type> theirValues (theirNumEntries);
675 for (size_type k = 0; k < theirNumRows;
676 theirCurPos += theirNumEntriesPerRow[k], k++) {
678 const GO theirRow = theirRows[k];
684 if (curNumEntries > 0) {
685 ArrayView<GO> colIndView =
686 colInd (curPos, curNumEntries);
687 ArrayView<GO> theirColIndView =
688 theirColInd (theirCurPos, curNumEntries);
689 std::copy (colIndView.begin(), colIndView.end(),
690 theirColIndView.begin());
692 ArrayView<scalar_type> valuesView =
693 values (curPos, curNumEntries);
694 ArrayView<scalar_type> theirValuesView =
695 theirValues (theirCurPos, curNumEntries);
696 std::copy (valuesView.begin(), valuesView.end(),
697 theirValuesView.begin());
704 send (*pComm, static_cast<int> (theirNumEntries),
705 theirColInd.getRawPtr(), p);
706 send (*pComm, static_cast<int> (theirNumEntries),
707 theirValues.getRawPtr(), p);
710 cerr <<
"-- Proc 0: Finished with proc " << p << endl;
718 numEntriesPerRow = null;
723 if (debug && myRank == 0) {
724 cerr <<
"-- Proc 0: About to fill in myRowPtr" << endl;
732 myRowPtr = arcp<size_t> (myNumRows+1);
734 for (size_type k = 1; k < myNumRows+1; ++k) {
735 myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
737 if (extraDebug && debug) {
738 cerr <<
"Proc " << Teuchos::rank (*(pRowMap->getComm()))
739 <<
": myRowPtr[0.." << myNumRows <<
"] = [";
740 for (size_type k = 0; k < myNumRows+1; ++k) {
746 cerr <<
"]" << endl << endl;
749 if (debug && myRank == 0) {
750 cerr <<
"-- Proc 0: Done with distribute" << endl;
767 static Teuchos::RCP<sparse_matrix_type>
768 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
769 Teuchos::ArrayRCP<size_t>& myRowPtr,
770 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
771 Teuchos::ArrayRCP<scalar_type>& myValues,
772 const Teuchos::RCP<const map_type>& pRowMap,
773 const Teuchos::RCP<const map_type>& pRangeMap,
774 const Teuchos::RCP<const map_type>& pDomainMap,
775 const bool callFillComplete =
true)
777 using Teuchos::ArrayView;
791 TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
792 "makeMatrix: myRowPtr array is null. "
793 "Please report this bug to the Tpetra developers.");
794 TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
795 "makeMatrix: domain map is null. "
796 "Please report this bug to the Tpetra developers.");
797 TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
798 "makeMatrix: range map is null. "
799 "Please report this bug to the Tpetra developers.");
800 TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
801 "makeMatrix: row map is null. "
802 "Please report this bug to the Tpetra developers.");
806 RCP<sparse_matrix_type> A =
811 ArrayView<const GO> myRows = pRowMap->getLocalElementList ();
812 const size_type myNumRows = myRows.size ();
815 const GO indexBase = pRowMap->getIndexBase ();
816 for (size_type i = 0; i < myNumRows; ++i) {
817 const size_type myCurPos = myRowPtr[i];
819 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
820 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
823 for (size_type k = 0; k < curNumEntries; ++k) {
824 curColInd[k] += indexBase;
827 if (curNumEntries > 0) {
828 A->insertGlobalValues (myRows[i], curColInd, curValues);
835 myNumEntriesPerRow = null;
840 if (callFillComplete) {
841 A->fillComplete (pDomainMap, pRangeMap);
851 static Teuchos::RCP<sparse_matrix_type>
852 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
853 Teuchos::ArrayRCP<size_t>& myRowPtr,
854 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
855 Teuchos::ArrayRCP<scalar_type>& myValues,
856 const Teuchos::RCP<const map_type>& pRowMap,
857 const Teuchos::RCP<const map_type>& pRangeMap,
858 const Teuchos::RCP<const map_type>& pDomainMap,
859 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
860 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams)
862 using Teuchos::ArrayView;
876 TEUCHOS_TEST_FOR_EXCEPTION(
877 myRowPtr.is_null(), std::logic_error,
878 "makeMatrix: myRowPtr array is null. "
879 "Please report this bug to the Tpetra developers.");
880 TEUCHOS_TEST_FOR_EXCEPTION(
881 pDomainMap.is_null(), std::logic_error,
882 "makeMatrix: domain map is null. "
883 "Please report this bug to the Tpetra developers.");
884 TEUCHOS_TEST_FOR_EXCEPTION(
885 pRangeMap.is_null(), std::logic_error,
886 "makeMatrix: range map is null. "
887 "Please report this bug to the Tpetra developers.");
888 TEUCHOS_TEST_FOR_EXCEPTION(
889 pRowMap.is_null(), std::logic_error,
890 "makeMatrix: row map is null. "
891 "Please report this bug to the Tpetra developers.");
895 RCP<sparse_matrix_type> A =
901 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
902 const size_type myNumRows = myRows.size();
905 const GO indexBase = pRowMap->getIndexBase ();
906 for (size_type i = 0; i < myNumRows; ++i) {
907 const size_type myCurPos = myRowPtr[i];
909 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
910 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
913 for (size_type k = 0; k < curNumEntries; ++k) {
914 curColInd[k] += indexBase;
916 if (curNumEntries > 0) {
917 A->insertGlobalValues (myRows[i], curColInd, curValues);
924 myNumEntriesPerRow = null;
929 A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
937 static Teuchos::RCP<sparse_matrix_type>
938 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
939 Teuchos::ArrayRCP<size_t>& myRowPtr,
940 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
941 Teuchos::ArrayRCP<scalar_type>& myValues,
942 const Teuchos::RCP<const map_type>& rowMap,
943 Teuchos::RCP<const map_type>& colMap,
944 const Teuchos::RCP<const map_type>& domainMap,
945 const Teuchos::RCP<const map_type>& rangeMap,
946 const bool callFillComplete =
true)
948 using Teuchos::ArrayView;
957 RCP<sparse_matrix_type> A;
958 if (colMap.is_null ()) {
966 ArrayView<const GO> myRows = rowMap->getLocalElementList ();
967 const size_type myNumRows = myRows.size ();
970 const GO indexBase = rowMap->getIndexBase ();
971 for (size_type i = 0; i < myNumRows; ++i) {
972 const size_type myCurPos = myRowPtr[i];
973 const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
974 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
975 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
978 for (size_type k = 0; k < curNumEntries; ++k) {
979 curColInd[k] += indexBase;
981 if (curNumEntries > 0) {
982 A->insertGlobalValues (myRows[i], curColInd, curValues);
989 myNumEntriesPerRow = null;
994 if (callFillComplete) {
995 A->fillComplete (domainMap, rangeMap);
996 if (colMap.is_null ()) {
997 colMap = A->getColMap ();
1021 static Teuchos::RCP<const Teuchos::MatrixMarket::Banner>
1022 readBanner (std::istream& in,
1024 const bool tolerant=
false,
1026 const bool isGraph=
false)
1028 using Teuchos::MatrixMarket::Banner;
1033 typedef Teuchos::ScalarTraits<scalar_type> STS;
1035 RCP<Banner> pBanner;
1039 const bool readFailed = ! getline(in, line);
1040 TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1041 "Failed to get Matrix Market banner line from input.");
1048 pBanner = rcp (
new Banner (line, tolerant));
1049 }
catch (std::exception& e) {
1051 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1052 "Matrix Market banner line contains syntax error(s): "
1056 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() !=
"matrix",
1057 std::invalid_argument,
"The Matrix Market file does not contain "
1058 "matrix data. Its Banner (first) line says that its object type is \""
1059 << pBanner->matrixType() <<
"\", rather than the required \"matrix\".");
1063 TEUCHOS_TEST_FOR_EXCEPTION(
1064 ! STS::isComplex && pBanner->dataType() ==
"complex",
1065 std::invalid_argument,
1066 "The Matrix Market file contains complex-valued data, but you are "
1067 "trying to read it into a matrix containing entries of the real-"
1068 "valued Scalar type \""
1069 << Teuchos::TypeNameTraits<scalar_type>::name() <<
"\".");
1070 TEUCHOS_TEST_FOR_EXCEPTION(
1072 pBanner->dataType() !=
"real" &&
1073 pBanner->dataType() !=
"complex" &&
1074 pBanner->dataType() !=
"integer",
1075 std::invalid_argument,
1076 "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
1077 "Matrix Market file may not contain a \"pattern\" matrix. A "
1078 "pattern matrix is really just a graph with no weights. It "
1079 "should be stored in a CrsGraph, not a CrsMatrix.");
1081 TEUCHOS_TEST_FOR_EXCEPTION(
1083 pBanner->dataType() !=
"pattern",
1084 std::invalid_argument,
1085 "When reading Matrix Market data into a Tpetra::CrsGraph, the "
1086 "Matrix Market file must contain a \"pattern\" matrix.");
1113 static Teuchos::Tuple<global_ordinal_type, 3>
1114 readCoordDims (std::istream& in,
1116 const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1118 const bool tolerant =
false,
1121 using Teuchos::MatrixMarket::readCoordinateDimensions;
1122 using Teuchos::Tuple;
1127 Tuple<global_ordinal_type, 3> dims;
1133 bool success =
false;
1134 if (pComm->getRank() == 0) {
1135 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() !=
"coordinate",
1136 std::invalid_argument,
"The Tpetra::CrsMatrix Matrix Market reader "
1137 "only accepts \"coordinate\" (sparse) matrix data.");
1141 success = readCoordinateDimensions (in, numRows, numCols,
1142 numNonzeros, lineNumber,
1148 dims[2] = numNonzeros;
1156 int the_success = success ? 1 : 0;
1157 Teuchos::broadcast (*pComm, 0, &the_success);
1158 success = (the_success == 1);
1163 Teuchos::broadcast (*pComm, 0, dims);
1171 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1172 "Error reading Matrix Market sparse matrix: failed to read "
1173 "coordinate matrix dimensions.");
1188 typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1190 typedef Teuchos::MatrixMarket::SymmetrizingGraphAdder<Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> > graph_adder_type;
1217 static Teuchos::RCP<adder_type>
1218 makeAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1219 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1220 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1221 const bool tolerant=
false,
1222 const bool debug=
false)
1224 if (pComm->getRank () == 0) {
1225 typedef Teuchos::MatrixMarket::Raw::Adder<
scalar_type,
1228 Teuchos::RCP<raw_adder_type> pRaw =
1229 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1231 return Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
1234 return Teuchos::null;
1263 static Teuchos::RCP<graph_adder_type>
1264 makeGraphAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1265 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1266 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1267 const bool tolerant=
false,
1268 const bool debug=
false)
1270 if (pComm->getRank () == 0) {
1271 typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
1272 Teuchos::RCP<raw_adder_type> pRaw =
1273 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1275 return Teuchos::rcp (
new graph_adder_type (pRaw, pBanner->symmType ()));
1278 return Teuchos::null;
1283 static Teuchos::RCP<sparse_graph_type>
1284 readSparseGraphHelper (std::istream& in,
1285 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1286 const Teuchos::RCP<const map_type>& rowMap,
1287 Teuchos::RCP<const map_type>& colMap,
1288 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1289 const bool tolerant,
1292 using Teuchos::MatrixMarket::Banner;
1295 using Teuchos::Tuple;
1299 const int myRank = pComm->getRank ();
1300 const int rootRank = 0;
1305 size_t lineNumber = 1;
1307 if (debug && myRank == rootRank) {
1308 cerr <<
"Matrix Market reader: readGraph:" << endl
1309 <<
"-- Reading banner line" << endl;
1318 RCP<const Banner> pBanner;
1324 int bannerIsCorrect = 1;
1325 std::ostringstream errMsg;
1327 if (myRank == rootRank) {
1330 pBanner = readBanner (in, lineNumber, tolerant, debug,
true);
1332 catch (std::exception& e) {
1333 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
1334 "threw an exception: " << e.what();
1335 bannerIsCorrect = 0;
1338 if (bannerIsCorrect) {
1343 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
1344 bannerIsCorrect = 0;
1345 errMsg <<
"The Matrix Market input file must contain a "
1346 "\"coordinate\"-format sparse graph in order to create a "
1347 "Tpetra::CrsGraph object from it, but the file's matrix "
1348 "type is \"" << pBanner->matrixType() <<
"\" instead.";
1353 if (tolerant && pBanner->matrixType() ==
"array") {
1354 bannerIsCorrect = 0;
1355 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
1356 "format sparse graph in order to create a Tpetra::CrsGraph "
1357 "object from it, but the file's matrix type is \"array\" "
1358 "instead. That probably means the file contains dense matrix "
1365 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1372 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1373 std::invalid_argument, errMsg.str ());
1375 if (debug && myRank == rootRank) {
1376 cerr <<
"-- Reading dimensions line" << endl;
1384 Tuple<global_ordinal_type, 3> dims =
1385 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1387 if (debug && myRank == rootRank) {
1388 cerr <<
"-- Making Adder for collecting graph data" << endl;
1395 RCP<graph_adder_type> pAdder =
1396 makeGraphAdder (pComm, pBanner, dims, tolerant, debug);
1398 if (debug && myRank == rootRank) {
1399 cerr <<
"-- Reading graph data" << endl;
1409 int readSuccess = 1;
1410 std::ostringstream errMsg;
1411 if (myRank == rootRank) {
1414 typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
1416 reader_type reader (pAdder);
1419 std::pair<bool, std::vector<size_t> > results =
1420 reader.read (in, lineNumber, tolerant, debug);
1421 readSuccess = results.first ? 1 : 0;
1423 catch (std::exception& e) {
1428 broadcast (*pComm, rootRank, ptr (&readSuccess));
1437 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1438 "Failed to read the Matrix Market sparse graph file: "
1442 if (debug && myRank == rootRank) {
1443 cerr <<
"-- Successfully read the Matrix Market data" << endl;
1456 if (debug && myRank == rootRank) {
1457 cerr <<
"-- Tolerant mode: rebroadcasting graph dimensions"
1459 <<
"----- Dimensions before: "
1460 << dims[0] <<
" x " << dims[1]
1464 Tuple<global_ordinal_type, 2> updatedDims;
1465 if (myRank == rootRank) {
1472 std::max (dims[0], pAdder->getAdder()->numRows());
1473 updatedDims[1] = pAdder->getAdder()->numCols();
1475 broadcast (*pComm, rootRank, updatedDims);
1476 dims[0] = updatedDims[0];
1477 dims[1] = updatedDims[1];
1478 if (debug && myRank == rootRank) {
1479 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
1493 if (myRank == rootRank) {
1500 if (dims[0] < pAdder->getAdder ()->numRows ()) {
1504 broadcast (*pComm, 0, ptr (&dimsMatch));
1505 if (dimsMatch == 0) {
1512 Tuple<global_ordinal_type, 2> addersDims;
1513 if (myRank == rootRank) {
1514 addersDims[0] = pAdder->getAdder()->numRows();
1515 addersDims[1] = pAdder->getAdder()->numCols();
1517 broadcast (*pComm, 0, addersDims);
1518 TEUCHOS_TEST_FOR_EXCEPTION(
1519 dimsMatch == 0, std::runtime_error,
1520 "The graph metadata says that the graph is " << dims[0] <<
" x "
1521 << dims[1] <<
", but the actual data says that the graph is "
1522 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
1523 "data includes more rows than reported in the metadata. This "
1524 "is not allowed when parsing in strict mode. Parse the graph in "
1525 "tolerant mode to ignore the metadata when it disagrees with the "
1531 RCP<map_type> proc0Map;
1533 if(Teuchos::is_null(rowMap)) {
1537 indexBase = rowMap->getIndexBase();
1539 if(myRank == rootRank) {
1540 proc0Map = rcp(
new map_type(dims[0],dims[0],indexBase,pComm));
1543 proc0Map = rcp(
new map_type(dims[0],0,indexBase,pComm));
1547 std::map<global_ordinal_type, size_t> numEntriesPerRow_map;
1548 if (myRank == rootRank) {
1549 const auto& entries = pAdder()->getAdder()->getEntries();
1554 for (
const auto& entry : entries) {
1556 ++numEntriesPerRow_map[gblRow];
1560 Teuchos::Array<size_t> numEntriesPerRow (proc0Map->getLocalNumElements ());
1561 for (
const auto& ent : numEntriesPerRow_map) {
1563 numEntriesPerRow[lclRow] = ent.second;
1570 std::map<global_ordinal_type, size_t> empty_map;
1571 std::swap (numEntriesPerRow_map, empty_map);
1574 RCP<sparse_graph_type> proc0Graph =
1576 constructorParams));
1577 if(myRank == rootRank) {
1578 typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
1581 const std::vector<element_type>& entries =
1582 pAdder->getAdder()->getEntries();
1585 for(
size_t curPos=0; curPos<entries.size(); curPos++) {
1586 const element_type& curEntry = entries[curPos];
1589 Teuchos::ArrayView<const global_ordinal_type> colView(&curCol,1);
1590 proc0Graph->insertGlobalIndices(curRow,colView);
1593 proc0Graph->fillComplete();
1595 RCP<sparse_graph_type> distGraph;
1596 if(Teuchos::is_null(rowMap))
1599 RCP<map_type> distMap =
1600 rcp(
new map_type(dims[0],0,pComm,GloballyDistributed));
1606 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1607 import_type importer (proc0Map, distMap);
1610 distGraph->doImport(*proc0Graph,importer,
INSERT);
1616 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1617 import_type importer (proc0Map, rowMap);
1620 distGraph->doImport(*proc0Graph,importer,
INSERT);
1650 static Teuchos::RCP<sparse_graph_type>
1653 const bool callFillComplete=
true,
1654 const bool tolerant=
false,
1655 const bool debug=
false)
1696 static Teuchos::RCP<sparse_graph_type>
1698 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1699 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1700 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1701 const bool tolerant=
false,
1702 const bool debug=
false)
1708 fillCompleteParams, tolerant, debug);
1752 static Teuchos::RCP<sparse_graph_type>
1754 const Teuchos::RCP<const map_type>& rowMap,
1755 Teuchos::RCP<const map_type>& colMap,
1756 const Teuchos::RCP<const map_type>& domainMap,
1757 const Teuchos::RCP<const map_type>& rangeMap,
1758 const bool callFillComplete=
true,
1759 const bool tolerant=
false,
1760 const bool debug=
false)
1762 TEUCHOS_TEST_FOR_EXCEPTION
1763 (rowMap.is_null (), std::invalid_argument,
1764 "Input rowMap must be nonnull.");
1766 if (comm.is_null ()) {
1769 return Teuchos::null;
1775 callFillComplete, tolerant, debug);
1803 static Teuchos::RCP<sparse_graph_type>
1805 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1806 const bool callFillComplete=
true,
1807 const bool tolerant=
false,
1808 const bool debug=
false)
1810 Teuchos::RCP<const map_type> fakeRowMap;
1811 Teuchos::RCP<const map_type> fakeColMap;
1812 Teuchos::RCP<Teuchos::ParameterList> fakeCtorParams;
1814 Teuchos::RCP<sparse_graph_type> graph =
1815 readSparseGraphHelper (in, pComm,
1816 fakeRowMap, fakeColMap,
1817 fakeCtorParams, tolerant, debug);
1818 if (callFillComplete) {
1819 graph->fillComplete ();
1854 static Teuchos::RCP<sparse_graph_type>
1856 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1857 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1858 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1859 const bool tolerant=
false,
1860 const bool debug=
false)
1862 Teuchos::RCP<const map_type> fakeRowMap;
1863 Teuchos::RCP<const map_type> fakeColMap;
1864 Teuchos::RCP<sparse_graph_type> graph =
1865 readSparseGraphHelper (in, pComm,
1866 fakeRowMap, fakeColMap,
1867 constructorParams, tolerant, debug);
1868 graph->fillComplete (fillCompleteParams);
1913 static Teuchos::RCP<sparse_graph_type>
1915 const Teuchos::RCP<const map_type>& rowMap,
1916 Teuchos::RCP<const map_type>& colMap,
1917 const Teuchos::RCP<const map_type>& domainMap,
1918 const Teuchos::RCP<const map_type>& rangeMap,
1919 const bool callFillComplete=
true,
1920 const bool tolerant=
false,
1921 const bool debug=
false)
1923 Teuchos::RCP<sparse_graph_type> graph =
1924 readSparseGraphHelper (in, rowMap->getComm (),
1925 rowMap, colMap, Teuchos::null, tolerant,
1927 if (callFillComplete) {
1928 graph->fillComplete (domainMap, rangeMap);
1933 #include "MatrixMarket_TpetraNew.hpp"
1958 static Teuchos::RCP<sparse_matrix_type>
1961 const bool callFillComplete=
true,
1962 const bool tolerant=
false,
1963 const bool debug=
false)
1967 return readSparse (in, comm, callFillComplete, tolerant, debug);
2002 static Teuchos::RCP<sparse_matrix_type>
2005 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2006 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2007 const bool tolerant=
false,
2008 const bool debug=
false)
2012 return readSparse (in, comm, constructorParams,
2013 fillCompleteParams, tolerant, debug);
2054 static Teuchos::RCP<sparse_matrix_type>
2056 const Teuchos::RCP<const map_type>& rowMap,
2057 Teuchos::RCP<const map_type>& colMap,
2058 const Teuchos::RCP<const map_type>& domainMap,
2059 const Teuchos::RCP<const map_type>& rangeMap,
2060 const bool callFillComplete=
true,
2061 const bool tolerant=
false,
2062 const bool debug=
false)
2064 TEUCHOS_TEST_FOR_EXCEPTION(
2065 rowMap.is_null (), std::invalid_argument,
2066 "Row Map must be nonnull.");
2072 return readSparse (in, rowMap, colMap, domainMap, rangeMap,
2073 callFillComplete, tolerant, debug);
2101 static Teuchos::RCP<sparse_matrix_type>
2103 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2104 const bool callFillComplete=
true,
2105 const bool tolerant=
false,
2106 const bool debug=
false)
2108 using Teuchos::MatrixMarket::Banner;
2109 using Teuchos::arcp;
2110 using Teuchos::ArrayRCP;
2111 using Teuchos::broadcast;
2112 using Teuchos::null;
2115 using Teuchos::REDUCE_MAX;
2116 using Teuchos::reduceAll;
2117 using Teuchos::Tuple;
2120 typedef Teuchos::ScalarTraits<scalar_type> STS;
2122 const bool extraDebug =
false;
2123 const int myRank = pComm->getRank ();
2124 const int rootRank = 0;
2129 size_t lineNumber = 1;
2131 if (debug && myRank == rootRank) {
2132 cerr <<
"Matrix Market reader: readSparse:" << endl
2133 <<
"-- Reading banner line" << endl;
2142 RCP<const Banner> pBanner;
2148 int bannerIsCorrect = 1;
2149 std::ostringstream errMsg;
2151 if (myRank == rootRank) {
2154 pBanner = readBanner (in, lineNumber, tolerant, debug);
2156 catch (std::exception& e) {
2157 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2158 "threw an exception: " << e.what();
2159 bannerIsCorrect = 0;
2162 if (bannerIsCorrect) {
2167 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2168 bannerIsCorrect = 0;
2169 errMsg <<
"The Matrix Market input file must contain a "
2170 "\"coordinate\"-format sparse matrix in order to create a "
2171 "Tpetra::CrsMatrix object from it, but the file's matrix "
2172 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2177 if (tolerant && pBanner->matrixType() ==
"array") {
2178 bannerIsCorrect = 0;
2179 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2180 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2181 "object from it, but the file's matrix type is \"array\" "
2182 "instead. That probably means the file contains dense matrix "
2189 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2196 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2197 std::invalid_argument, errMsg.str ());
2199 if (debug && myRank == rootRank) {
2200 cerr <<
"-- Reading dimensions line" << endl;
2208 Tuple<global_ordinal_type, 3> dims =
2209 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2211 if (debug && myRank == rootRank) {
2212 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2217 RCP<adder_type> pAdder =
2218 makeAdder (pComm, pBanner, dims, tolerant, debug);
2220 if (debug && myRank == rootRank) {
2221 cerr <<
"-- Reading matrix data" << endl;
2231 int readSuccess = 1;
2232 std::ostringstream errMsg;
2233 if (myRank == rootRank) {
2236 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2238 reader_type reader (pAdder);
2241 std::pair<bool, std::vector<size_t> > results =
2242 reader.read (in, lineNumber, tolerant, debug);
2243 readSuccess = results.first ? 1 : 0;
2245 catch (std::exception& e) {
2250 broadcast (*pComm, rootRank, ptr (&readSuccess));
2259 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2260 "Failed to read the Matrix Market sparse matrix file: "
2264 if (debug && myRank == rootRank) {
2265 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2278 if (debug && myRank == rootRank) {
2279 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2281 <<
"----- Dimensions before: "
2282 << dims[0] <<
" x " << dims[1]
2286 Tuple<global_ordinal_type, 2> updatedDims;
2287 if (myRank == rootRank) {
2294 std::max (dims[0], pAdder->getAdder()->numRows());
2295 updatedDims[1] = pAdder->getAdder()->numCols();
2297 broadcast (*pComm, rootRank, updatedDims);
2298 dims[0] = updatedDims[0];
2299 dims[1] = updatedDims[1];
2300 if (debug && myRank == rootRank) {
2301 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
2315 if (myRank == rootRank) {
2322 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2326 broadcast (*pComm, 0, ptr (&dimsMatch));
2327 if (dimsMatch == 0) {
2334 Tuple<global_ordinal_type, 2> addersDims;
2335 if (myRank == rootRank) {
2336 addersDims[0] = pAdder->getAdder()->numRows();
2337 addersDims[1] = pAdder->getAdder()->numCols();
2339 broadcast (*pComm, 0, addersDims);
2340 TEUCHOS_TEST_FOR_EXCEPTION(
2341 dimsMatch == 0, std::runtime_error,
2342 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
2343 << dims[1] <<
", but the actual data says that the matrix is "
2344 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
2345 "data includes more rows than reported in the metadata. This "
2346 "is not allowed when parsing in strict mode. Parse the matrix in "
2347 "tolerant mode to ignore the metadata when it disagrees with the "
2352 if (debug && myRank == rootRank) {
2353 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
2365 ArrayRCP<size_t> numEntriesPerRow;
2366 ArrayRCP<size_t> rowPtr;
2367 ArrayRCP<global_ordinal_type> colInd;
2368 ArrayRCP<scalar_type> values;
2373 int mergeAndConvertSucceeded = 1;
2374 std::ostringstream errMsg;
2376 if (myRank == rootRank) {
2378 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2388 const size_type numRows = dims[0];
2391 pAdder->getAdder()->merge ();
2394 const std::vector<element_type>& entries =
2395 pAdder->getAdder()->getEntries();
2398 const size_t numEntries = (size_t)entries.size();
2401 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
2402 <<
" rows and numEntries=" << numEntries
2403 <<
" entries." << endl;
2409 numEntriesPerRow = arcp<size_t> (numRows);
2410 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2411 rowPtr = arcp<size_t> (numRows+1);
2412 std::fill (rowPtr.begin(), rowPtr.end(), 0);
2413 colInd = arcp<global_ordinal_type> (numEntries);
2414 values = arcp<scalar_type> (numEntries);
2421 for (curPos = 0; curPos < numEntries; ++curPos) {
2422 const element_type& curEntry = entries[curPos];
2424 TEUCHOS_TEST_FOR_EXCEPTION(
2425 curRow < prvRow, std::logic_error,
2426 "Row indices are out of order, even though they are supposed "
2427 "to be sorted. curRow = " << curRow <<
", prvRow = "
2428 << prvRow <<
", at curPos = " << curPos <<
". Please report "
2429 "this bug to the Tpetra developers.");
2430 if (curRow > prvRow) {
2436 numEntriesPerRow[curRow]++;
2437 colInd[curPos] = curEntry.colIndex();
2438 values[curPos] = curEntry.value();
2443 rowPtr[numRows] = numEntries;
2445 catch (std::exception& e) {
2446 mergeAndConvertSucceeded = 0;
2447 errMsg <<
"Failed to merge sparse matrix entries and convert to "
2448 "CSR format: " << e.what();
2451 if (debug && mergeAndConvertSucceeded) {
2453 const size_type numRows = dims[0];
2454 const size_type maxToDisplay = 100;
2456 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
2457 << (numEntriesPerRow.size()-1) <<
"] ";
2458 if (numRows > maxToDisplay) {
2459 cerr <<
"(only showing first and last few entries) ";
2463 if (numRows > maxToDisplay) {
2464 for (size_type k = 0; k < 2; ++k) {
2465 cerr << numEntriesPerRow[k] <<
" ";
2468 for (size_type k = numRows-2; k < numRows-1; ++k) {
2469 cerr << numEntriesPerRow[k] <<
" ";
2473 for (size_type k = 0; k < numRows-1; ++k) {
2474 cerr << numEntriesPerRow[k] <<
" ";
2477 cerr << numEntriesPerRow[numRows-1];
2479 cerr <<
"]" << endl;
2481 cerr <<
"----- Proc 0: rowPtr ";
2482 if (numRows > maxToDisplay) {
2483 cerr <<
"(only showing first and last few entries) ";
2486 if (numRows > maxToDisplay) {
2487 for (size_type k = 0; k < 2; ++k) {
2488 cerr << rowPtr[k] <<
" ";
2491 for (size_type k = numRows-2; k < numRows; ++k) {
2492 cerr << rowPtr[k] <<
" ";
2496 for (size_type k = 0; k < numRows; ++k) {
2497 cerr << rowPtr[k] <<
" ";
2500 cerr << rowPtr[numRows] <<
"]" << endl;
2511 if (debug && myRank == rootRank) {
2512 cerr <<
"-- Making range, domain, and row maps" << endl;
2519 RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
2520 RCP<const map_type> pDomainMap =
2521 makeDomainMap (pRangeMap, dims[0], dims[1]);
2522 RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
2524 if (debug && myRank == rootRank) {
2525 cerr <<
"-- Distributing the matrix data" << endl;
2539 ArrayRCP<size_t> myNumEntriesPerRow;
2540 ArrayRCP<size_t> myRowPtr;
2541 ArrayRCP<global_ordinal_type> myColInd;
2542 ArrayRCP<scalar_type> myValues;
2544 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
2545 numEntriesPerRow, rowPtr, colInd, values, debug);
2547 if (debug && myRank == rootRank) {
2548 cerr <<
"-- Inserting matrix entries on each processor";
2549 if (callFillComplete) {
2550 cerr <<
" and calling fillComplete()";
2561 RCP<sparse_matrix_type> pMatrix =
2562 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
2563 pRowMap, pRangeMap, pDomainMap, callFillComplete);
2568 int localIsNull = pMatrix.is_null () ? 1 : 0;
2569 int globalIsNull = 0;
2570 reduceAll (*pComm, REDUCE_MAX, localIsNull, ptr (&globalIsNull));
2571 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
2572 "Reader::makeMatrix() returned a null pointer on at least one "
2573 "process. Please report this bug to the Tpetra developers.");
2576 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
2577 "Reader::makeMatrix() returned a null pointer. "
2578 "Please report this bug to the Tpetra developers.");
2590 if (callFillComplete) {
2591 const int numProcs = pComm->getSize ();
2593 if (extraDebug && debug) {
2595 pRangeMap->getGlobalNumElements ();
2597 pDomainMap->getGlobalNumElements ();
2598 if (myRank == rootRank) {
2599 cerr <<
"-- Matrix is "
2600 << globalNumRows <<
" x " << globalNumCols
2601 <<
" with " << pMatrix->getGlobalNumEntries()
2602 <<
" entries, and index base "
2603 << pMatrix->getIndexBase() <<
"." << endl;
2606 for (
int p = 0; p < numProcs; ++p) {
2608 cerr <<
"-- Proc " << p <<
" owns "
2609 << pMatrix->getLocalNumCols() <<
" columns, and "
2610 << pMatrix->getLocalNumEntries() <<
" entries." << endl;
2617 if (debug && myRank == rootRank) {
2618 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
2653 static Teuchos::RCP<sparse_matrix_type>
2655 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2656 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2657 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2658 const bool tolerant=
false,
2659 const bool debug=
false)
2661 using Teuchos::MatrixMarket::Banner;
2662 using Teuchos::arcp;
2663 using Teuchos::ArrayRCP;
2664 using Teuchos::broadcast;
2665 using Teuchos::null;
2668 using Teuchos::reduceAll;
2669 using Teuchos::Tuple;
2672 typedef Teuchos::ScalarTraits<scalar_type> STS;
2674 const bool extraDebug =
false;
2675 const int myRank = pComm->getRank ();
2676 const int rootRank = 0;
2681 size_t lineNumber = 1;
2683 if (debug && myRank == rootRank) {
2684 cerr <<
"Matrix Market reader: readSparse:" << endl
2685 <<
"-- Reading banner line" << endl;
2694 RCP<const Banner> pBanner;
2700 int bannerIsCorrect = 1;
2701 std::ostringstream errMsg;
2703 if (myRank == rootRank) {
2706 pBanner = readBanner (in, lineNumber, tolerant, debug);
2708 catch (std::exception& e) {
2709 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2710 "threw an exception: " << e.what();
2711 bannerIsCorrect = 0;
2714 if (bannerIsCorrect) {
2719 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2720 bannerIsCorrect = 0;
2721 errMsg <<
"The Matrix Market input file must contain a "
2722 "\"coordinate\"-format sparse matrix in order to create a "
2723 "Tpetra::CrsMatrix object from it, but the file's matrix "
2724 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2729 if (tolerant && pBanner->matrixType() ==
"array") {
2730 bannerIsCorrect = 0;
2731 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2732 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2733 "object from it, but the file's matrix type is \"array\" "
2734 "instead. That probably means the file contains dense matrix "
2741 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2748 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2749 std::invalid_argument, errMsg.str ());
2751 if (debug && myRank == rootRank) {
2752 cerr <<
"-- Reading dimensions line" << endl;
2760 Tuple<global_ordinal_type, 3> dims =
2761 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2763 if (debug && myRank == rootRank) {
2764 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2769 RCP<adder_type> pAdder =
2770 makeAdder (pComm, pBanner, dims, tolerant, debug);
2772 if (debug && myRank == rootRank) {
2773 cerr <<
"-- Reading matrix data" << endl;
2783 int readSuccess = 1;
2784 std::ostringstream errMsg;
2785 if (myRank == rootRank) {
2788 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2790 reader_type reader (pAdder);
2793 std::pair<bool, std::vector<size_t> > results =
2794 reader.read (in, lineNumber, tolerant, debug);
2795 readSuccess = results.first ? 1 : 0;
2797 catch (std::exception& e) {
2802 broadcast (*pComm, rootRank, ptr (&readSuccess));
2811 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2812 "Failed to read the Matrix Market sparse matrix file: "
2816 if (debug && myRank == rootRank) {
2817 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2830 if (debug && myRank == rootRank) {
2831 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2833 <<
"----- Dimensions before: "
2834 << dims[0] <<
" x " << dims[1]
2838 Tuple<global_ordinal_type, 2> updatedDims;
2839 if (myRank == rootRank) {
2846 std::max (dims[0], pAdder->getAdder()->numRows());
2847 updatedDims[1] = pAdder->getAdder()->numCols();
2849 broadcast (*pComm, rootRank, updatedDims);
2850 dims[0] = updatedDims[0];
2851 dims[1] = updatedDims[1];
2852 if (debug && myRank == rootRank) {
2853 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
2867 if (myRank == rootRank) {
2874 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2878 broadcast (*pComm, 0, ptr (&dimsMatch));
2879 if (dimsMatch == 0) {
2886 Tuple<global_ordinal_type, 2> addersDims;
2887 if (myRank == rootRank) {
2888 addersDims[0] = pAdder->getAdder()->numRows();
2889 addersDims[1] = pAdder->getAdder()->numCols();
2891 broadcast (*pComm, 0, addersDims);
2892 TEUCHOS_TEST_FOR_EXCEPTION(
2893 dimsMatch == 0, std::runtime_error,
2894 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
2895 << dims[1] <<
", but the actual data says that the matrix is "
2896 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
2897 "data includes more rows than reported in the metadata. This "
2898 "is not allowed when parsing in strict mode. Parse the matrix in "
2899 "tolerant mode to ignore the metadata when it disagrees with the "
2904 if (debug && myRank == rootRank) {
2905 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
2917 ArrayRCP<size_t> numEntriesPerRow;
2918 ArrayRCP<size_t> rowPtr;
2919 ArrayRCP<global_ordinal_type> colInd;
2920 ArrayRCP<scalar_type> values;
2925 int mergeAndConvertSucceeded = 1;
2926 std::ostringstream errMsg;
2928 if (myRank == rootRank) {
2930 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2940 const size_type numRows = dims[0];
2943 pAdder->getAdder()->merge ();
2946 const std::vector<element_type>& entries =
2947 pAdder->getAdder()->getEntries();
2950 const size_t numEntries = (size_t)entries.size();
2953 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
2954 <<
" rows and numEntries=" << numEntries
2955 <<
" entries." << endl;
2961 numEntriesPerRow = arcp<size_t> (numRows);
2962 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2963 rowPtr = arcp<size_t> (numRows+1);
2964 std::fill (rowPtr.begin(), rowPtr.end(), 0);
2965 colInd = arcp<global_ordinal_type> (numEntries);
2966 values = arcp<scalar_type> (numEntries);
2973 for (curPos = 0; curPos < numEntries; ++curPos) {
2974 const element_type& curEntry = entries[curPos];
2976 TEUCHOS_TEST_FOR_EXCEPTION(
2977 curRow < prvRow, std::logic_error,
2978 "Row indices are out of order, even though they are supposed "
2979 "to be sorted. curRow = " << curRow <<
", prvRow = "
2980 << prvRow <<
", at curPos = " << curPos <<
". Please report "
2981 "this bug to the Tpetra developers.");
2982 if (curRow > prvRow) {
2988 numEntriesPerRow[curRow]++;
2989 colInd[curPos] = curEntry.colIndex();
2990 values[curPos] = curEntry.value();
2995 rowPtr[numRows] = numEntries;
2997 catch (std::exception& e) {
2998 mergeAndConvertSucceeded = 0;
2999 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3000 "CSR format: " << e.what();
3003 if (debug && mergeAndConvertSucceeded) {
3005 const size_type numRows = dims[0];
3006 const size_type maxToDisplay = 100;
3008 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3009 << (numEntriesPerRow.size()-1) <<
"] ";
3010 if (numRows > maxToDisplay) {
3011 cerr <<
"(only showing first and last few entries) ";
3015 if (numRows > maxToDisplay) {
3016 for (size_type k = 0; k < 2; ++k) {
3017 cerr << numEntriesPerRow[k] <<
" ";
3020 for (size_type k = numRows-2; k < numRows-1; ++k) {
3021 cerr << numEntriesPerRow[k] <<
" ";
3025 for (size_type k = 0; k < numRows-1; ++k) {
3026 cerr << numEntriesPerRow[k] <<
" ";
3029 cerr << numEntriesPerRow[numRows-1];
3031 cerr <<
"]" << endl;
3033 cerr <<
"----- Proc 0: rowPtr ";
3034 if (numRows > maxToDisplay) {
3035 cerr <<
"(only showing first and last few entries) ";
3038 if (numRows > maxToDisplay) {
3039 for (size_type k = 0; k < 2; ++k) {
3040 cerr << rowPtr[k] <<
" ";
3043 for (size_type k = numRows-2; k < numRows; ++k) {
3044 cerr << rowPtr[k] <<
" ";
3048 for (size_type k = 0; k < numRows; ++k) {
3049 cerr << rowPtr[k] <<
" ";
3052 cerr << rowPtr[numRows] <<
"]" << endl;
3063 if (debug && myRank == rootRank) {
3064 cerr <<
"-- Making range, domain, and row maps" << endl;
3071 RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
3072 RCP<const map_type> pDomainMap =
3073 makeDomainMap (pRangeMap, dims[0], dims[1]);
3074 RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
3076 if (debug && myRank == rootRank) {
3077 cerr <<
"-- Distributing the matrix data" << endl;
3091 ArrayRCP<size_t> myNumEntriesPerRow;
3092 ArrayRCP<size_t> myRowPtr;
3093 ArrayRCP<global_ordinal_type> myColInd;
3094 ArrayRCP<scalar_type> myValues;
3096 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
3097 numEntriesPerRow, rowPtr, colInd, values, debug);
3099 if (debug && myRank == rootRank) {
3100 cerr <<
"-- Inserting matrix entries on each process "
3101 "and calling fillComplete()" << endl;
3110 Teuchos::RCP<sparse_matrix_type> pMatrix =
3111 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
3112 pRowMap, pRangeMap, pDomainMap, constructorParams,
3113 fillCompleteParams);
3118 int localIsNull = pMatrix.is_null () ? 1 : 0;
3119 int globalIsNull = 0;
3120 reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
3121 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
3122 "Reader::makeMatrix() returned a null pointer on at least one "
3123 "process. Please report this bug to the Tpetra developers.");
3126 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
3127 "Reader::makeMatrix() returned a null pointer. "
3128 "Please report this bug to the Tpetra developers.");
3137 if (extraDebug && debug) {
3138 const int numProcs = pComm->getSize ();
3140 pRangeMap->getGlobalNumElements();
3142 pDomainMap->getGlobalNumElements();
3143 if (myRank == rootRank) {
3144 cerr <<
"-- Matrix is "
3145 << globalNumRows <<
" x " << globalNumCols
3146 <<
" with " << pMatrix->getGlobalNumEntries()
3147 <<
" entries, and index base "
3148 << pMatrix->getIndexBase() <<
"." << endl;
3151 for (
int p = 0; p < numProcs; ++p) {
3153 cerr <<
"-- Proc " << p <<
" owns "
3154 << pMatrix->getLocalNumCols() <<
" columns, and "
3155 << pMatrix->getLocalNumEntries() <<
" entries." << endl;
3161 if (debug && myRank == rootRank) {
3162 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3209 static Teuchos::RCP<sparse_matrix_type>
3211 const Teuchos::RCP<const map_type>& rowMap,
3212 Teuchos::RCP<const map_type>& colMap,
3213 const Teuchos::RCP<const map_type>& domainMap,
3214 const Teuchos::RCP<const map_type>& rangeMap,
3215 const bool callFillComplete=
true,
3216 const bool tolerant=
false,
3217 const bool debug=
false)
3219 using Teuchos::MatrixMarket::Banner;
3220 using Teuchos::arcp;
3221 using Teuchos::ArrayRCP;
3222 using Teuchos::ArrayView;
3224 using Teuchos::broadcast;
3225 using Teuchos::Comm;
3226 using Teuchos::null;
3229 using Teuchos::reduceAll;
3230 using Teuchos::Tuple;
3233 typedef Teuchos::ScalarTraits<scalar_type> STS;
3235 RCP<const Comm<int> > pComm = rowMap->getComm ();
3236 const int myRank = pComm->getRank ();
3237 const int rootRank = 0;
3238 const bool extraDebug =
false;
3243 TEUCHOS_TEST_FOR_EXCEPTION(
3244 rowMap.is_null (), std::invalid_argument,
3245 "Row Map must be nonnull.");
3246 TEUCHOS_TEST_FOR_EXCEPTION(
3247 rangeMap.is_null (), std::invalid_argument,
3248 "Range Map must be nonnull.");
3249 TEUCHOS_TEST_FOR_EXCEPTION(
3250 domainMap.is_null (), std::invalid_argument,
3251 "Domain Map must be nonnull.");
3252 TEUCHOS_TEST_FOR_EXCEPTION(
3253 rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
3254 std::invalid_argument,
3255 "The specified row Map's communicator (rowMap->getComm())"
3256 "differs from the given separately supplied communicator pComm.");
3257 TEUCHOS_TEST_FOR_EXCEPTION(
3258 domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
3259 std::invalid_argument,
3260 "The specified domain Map's communicator (domainMap->getComm())"
3261 "differs from the given separately supplied communicator pComm.");
3262 TEUCHOS_TEST_FOR_EXCEPTION(
3263 rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
3264 std::invalid_argument,
3265 "The specified range Map's communicator (rangeMap->getComm())"
3266 "differs from the given separately supplied communicator pComm.");
3271 size_t lineNumber = 1;
3273 if (debug && myRank == rootRank) {
3274 cerr <<
"Matrix Market reader: readSparse:" << endl
3275 <<
"-- Reading banner line" << endl;
3284 RCP<const Banner> pBanner;
3290 int bannerIsCorrect = 1;
3291 std::ostringstream errMsg;
3293 if (myRank == rootRank) {
3296 pBanner = readBanner (in, lineNumber, tolerant, debug);
3298 catch (std::exception& e) {
3299 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
3300 "threw an exception: " << e.what();
3301 bannerIsCorrect = 0;
3304 if (bannerIsCorrect) {
3309 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
3310 bannerIsCorrect = 0;
3311 errMsg <<
"The Matrix Market input file must contain a "
3312 "\"coordinate\"-format sparse matrix in order to create a "
3313 "Tpetra::CrsMatrix object from it, but the file's matrix "
3314 "type is \"" << pBanner->matrixType() <<
"\" instead.";
3319 if (tolerant && pBanner->matrixType() ==
"array") {
3320 bannerIsCorrect = 0;
3321 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
3322 "format sparse matrix in order to create a Tpetra::CrsMatrix "
3323 "object from it, but the file's matrix type is \"array\" "
3324 "instead. That probably means the file contains dense matrix "
3331 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
3338 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
3339 std::invalid_argument, errMsg.str ());
3341 if (debug && myRank == rootRank) {
3342 cerr <<
"-- Reading dimensions line" << endl;
3350 Tuple<global_ordinal_type, 3> dims =
3351 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
3353 if (debug && myRank == rootRank) {
3354 cerr <<
"-- Making Adder for collecting matrix data" << endl;
3361 RCP<adder_type> pAdder =
3362 makeAdder (pComm, pBanner, dims, tolerant, debug);
3364 if (debug && myRank == rootRank) {
3365 cerr <<
"-- Reading matrix data" << endl;
3375 int readSuccess = 1;
3376 std::ostringstream errMsg;
3377 if (myRank == rootRank) {
3380 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3382 reader_type reader (pAdder);
3385 std::pair<bool, std::vector<size_t> > results =
3386 reader.read (in, lineNumber, tolerant, debug);
3387 readSuccess = results.first ? 1 : 0;
3389 catch (std::exception& e) {
3394 broadcast (*pComm, rootRank, ptr (&readSuccess));
3403 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
3404 "Failed to read the Matrix Market sparse matrix file: "
3408 if (debug && myRank == rootRank) {
3409 cerr <<
"-- Successfully read the Matrix Market data" << endl;
3422 if (debug && myRank == rootRank) {
3423 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
3425 <<
"----- Dimensions before: "
3426 << dims[0] <<
" x " << dims[1]
3430 Tuple<global_ordinal_type, 2> updatedDims;
3431 if (myRank == rootRank) {
3438 std::max (dims[0], pAdder->getAdder()->numRows());
3439 updatedDims[1] = pAdder->getAdder()->numCols();
3441 broadcast (*pComm, rootRank, updatedDims);
3442 dims[0] = updatedDims[0];
3443 dims[1] = updatedDims[1];
3444 if (debug && myRank == rootRank) {
3445 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
3459 if (myRank == rootRank) {
3466 if (dims[0] < pAdder->getAdder ()->numRows ()) {
3470 broadcast (*pComm, 0, ptr (&dimsMatch));
3471 if (dimsMatch == 0) {
3478 Tuple<global_ordinal_type, 2> addersDims;
3479 if (myRank == rootRank) {
3480 addersDims[0] = pAdder->getAdder()->numRows();
3481 addersDims[1] = pAdder->getAdder()->numCols();
3483 broadcast (*pComm, 0, addersDims);
3484 TEUCHOS_TEST_FOR_EXCEPTION(
3485 dimsMatch == 0, std::runtime_error,
3486 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
3487 << dims[1] <<
", but the actual data says that the matrix is "
3488 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
3489 "data includes more rows than reported in the metadata. This "
3490 "is not allowed when parsing in strict mode. Parse the matrix in "
3491 "tolerant mode to ignore the metadata when it disagrees with the "
3496 if (debug && myRank == rootRank) {
3497 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
3509 ArrayRCP<size_t> numEntriesPerRow;
3510 ArrayRCP<size_t> rowPtr;
3511 ArrayRCP<global_ordinal_type> colInd;
3512 ArrayRCP<scalar_type> values;
3517 int mergeAndConvertSucceeded = 1;
3518 std::ostringstream errMsg;
3520 if (myRank == rootRank) {
3522 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3532 const size_type numRows = dims[0];
3535 pAdder->getAdder()->merge ();
3538 const std::vector<element_type>& entries =
3539 pAdder->getAdder()->getEntries();
3542 const size_t numEntries = (size_t)entries.size();
3545 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
3546 <<
" rows and numEntries=" << numEntries
3547 <<
" entries." << endl;
3553 numEntriesPerRow = arcp<size_t> (numRows);
3554 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3555 rowPtr = arcp<size_t> (numRows+1);
3556 std::fill (rowPtr.begin(), rowPtr.end(), 0);
3557 colInd = arcp<global_ordinal_type> (numEntries);
3558 values = arcp<scalar_type> (numEntries);
3565 for (curPos = 0; curPos < numEntries; ++curPos) {
3566 const element_type& curEntry = entries[curPos];
3568 TEUCHOS_TEST_FOR_EXCEPTION(
3569 curRow < prvRow, std::logic_error,
3570 "Row indices are out of order, even though they are supposed "
3571 "to be sorted. curRow = " << curRow <<
", prvRow = "
3572 << prvRow <<
", at curPos = " << curPos <<
". Please report "
3573 "this bug to the Tpetra developers.");
3574 if (curRow > prvRow) {
3577 numEntriesPerRow[curRow]++;
3578 colInd[curPos] = curEntry.colIndex();
3579 values[curPos] = curEntry.value();
3584 rowPtr[row] = numEntriesPerRow[row-1] + rowPtr[row-1];
3587 catch (std::exception& e) {
3588 mergeAndConvertSucceeded = 0;
3589 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3590 "CSR format: " << e.what();
3593 if (debug && mergeAndConvertSucceeded) {
3595 const size_type numRows = dims[0];
3596 const size_type maxToDisplay = 100;
3598 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3599 << (numEntriesPerRow.size()-1) <<
"] ";
3600 if (numRows > maxToDisplay) {
3601 cerr <<
"(only showing first and last few entries) ";
3605 if (numRows > maxToDisplay) {
3606 for (size_type k = 0; k < 2; ++k) {
3607 cerr << numEntriesPerRow[k] <<
" ";
3610 for (size_type k = numRows-2; k < numRows-1; ++k) {
3611 cerr << numEntriesPerRow[k] <<
" ";
3615 for (size_type k = 0; k < numRows-1; ++k) {
3616 cerr << numEntriesPerRow[k] <<
" ";
3619 cerr << numEntriesPerRow[numRows-1];
3621 cerr <<
"]" << endl;
3623 cerr <<
"----- Proc 0: rowPtr ";
3624 if (numRows > maxToDisplay) {
3625 cerr <<
"(only showing first and last few entries) ";
3628 if (numRows > maxToDisplay) {
3629 for (size_type k = 0; k < 2; ++k) {
3630 cerr << rowPtr[k] <<
" ";
3633 for (size_type k = numRows-2; k < numRows; ++k) {
3634 cerr << rowPtr[k] <<
" ";
3638 for (size_type k = 0; k < numRows; ++k) {
3639 cerr << rowPtr[k] <<
" ";
3642 cerr << rowPtr[numRows] <<
"]" << endl;
3644 cerr <<
"----- Proc 0: colInd = [";
3645 for (
size_t k = 0; k < rowPtr[numRows]; ++k) {
3646 cerr << colInd[k] <<
" ";
3648 cerr <<
"]" << endl;
3662 if (debug && myRank == rootRank) {
3663 cerr <<
"-- Verifying Maps" << endl;
3665 TEUCHOS_TEST_FOR_EXCEPTION(
3666 as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
3667 std::invalid_argument,
3668 "The range Map has " << rangeMap->getGlobalNumElements ()
3669 <<
" entries, but the matrix has a global number of rows " << dims[0]
3671 TEUCHOS_TEST_FOR_EXCEPTION(
3672 as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
3673 std::invalid_argument,
3674 "The domain Map has " << domainMap->getGlobalNumElements ()
3675 <<
" entries, but the matrix has a global number of columns "
3679 RCP<Teuchos::FancyOStream> err = debug ?
3680 Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
3682 RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
3683 ArrayView<const global_ordinal_type> myRows =
3684 gatherRowMap->getLocalElementList ();
3685 const size_type myNumRows = myRows.size ();
3688 ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
3689 for (size_type i_ = 0; i_ < myNumRows; i_++) {
3690 gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_]-indexBase];
3696 RCP<sparse_matrix_type> A_proc0 =
3698 if (myRank == rootRank) {
3700 cerr <<
"-- Proc 0: Filling gather matrix" << endl;
3703 cerr <<
"---- Rows: " << Teuchos::toString (myRows) << endl;
3707 for (size_type i_ = 0; i_ < myNumRows; ++i_) {
3708 size_type i = myRows[i_] - indexBase;
3710 const size_type curPos = as<size_type> (rowPtr[i]);
3712 ArrayView<global_ordinal_type> curColInd =
3713 colInd.view (curPos, curNumEntries);
3714 ArrayView<scalar_type> curValues =
3715 values.view (curPos, curNumEntries);
3718 for (size_type k = 0; k < curNumEntries; ++k) {
3719 curColInd[k] += indexBase;
3722 cerr <<
"------ Columns: " << Teuchos::toString (curColInd) << endl;
3723 cerr <<
"------ Values: " << Teuchos::toString (curValues) << endl;
3726 if (curNumEntries > 0) {
3727 A_proc0->insertGlobalValues (myRows[i_], curColInd, curValues);
3733 numEntriesPerRow = null;
3739 RCP<sparse_matrix_type> A;
3741 export_type exp (gatherRowMap, rowMap);
3748 mv_type_go target_nnzPerRow(rowMap,1);
3749 mv_type_go source_nnzPerRow(gatherRowMap,1);
3750 Teuchos::ArrayRCP<GO> srcData = source_nnzPerRow.getDataNonConst(0);
3751 for (
int i=0; i<myNumRows; i++)
3752 srcData[i] = gatherNumEntriesPerRow[i];
3753 srcData = Teuchos::null;
3755 Teuchos::ArrayRCP<GO> targetData = target_nnzPerRow.getDataNonConst(0);
3756 ArrayRCP<size_t> targetData_size_t = arcp<size_t>(targetData.size());
3757 for (
int i=0; i<targetData.size(); i++)
3758 targetData_size_t[i] = targetData[i];
3760 if (colMap.is_null ()) {
3765 A->doExport (*A_proc0, exp,
INSERT);
3766 if (callFillComplete) {
3767 A->fillComplete (domainMap, rangeMap);
3779 if (callFillComplete) {
3780 const int numProcs = pComm->getSize ();
3782 if (extraDebug && debug) {
3783 const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
3784 const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
3785 if (myRank == rootRank) {
3786 cerr <<
"-- Matrix is "
3787 << globalNumRows <<
" x " << globalNumCols
3788 <<
" with " << A->getGlobalNumEntries()
3789 <<
" entries, and index base "
3790 << A->getIndexBase() <<
"." << endl;
3793 for (
int p = 0; p < numProcs; ++p) {
3795 cerr <<
"-- Proc " << p <<
" owns "
3796 << A->getLocalNumCols() <<
" columns, and "
3797 << A->getLocalNumEntries() <<
" entries." << endl;
3804 if (debug && myRank == rootRank) {
3805 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3840 static Teuchos::RCP<multivector_type>
3843 Teuchos::RCP<const map_type>& map,
3844 const bool tolerant=
false,
3845 const bool debug=
false,
3846 const bool binary=
false)
3850 return readDense (in, comm, map, tolerant, debug, binary);
3865 const std::string& filename,
const bool safe =
true,
3866 std::ios_base::openmode mode = std::ios_base::in
3872 int all_should_stop =
false;
3875 if (comm->getRank() == 0) {
3876 in.open(filename, mode);
3877 all_should_stop = !in && safe;
3881 if(comm) Teuchos::broadcast(*comm, 0, &all_should_stop);
3883 TEUCHOS_TEST_FOR_EXCEPTION(
3886 "Could not open input file '" + filename +
"' on root rank 0."
3922 static Teuchos::RCP<vector_type>
3925 Teuchos::RCP<const map_type>& map,
3926 const bool tolerant=
false,
3927 const bool debug=
false)
3931 return readVector (in, comm, map, tolerant, debug);
4003 static Teuchos::RCP<multivector_type>
4006 Teuchos::RCP<const map_type>& map,
4007 const bool tolerant=
false,
4008 const bool debug=
false,
4009 const bool binary=
false)
4011 Teuchos::RCP<Teuchos::FancyOStream> err =
4012 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4013 return readDenseImpl<scalar_type> (in, comm, map, err, tolerant, debug, binary);
4018 static Teuchos::RCP<vector_type>
4021 Teuchos::RCP<const map_type>& map,
4022 const bool tolerant=
false,
4023 const bool debug=
false)
4025 Teuchos::RCP<Teuchos::FancyOStream> err =
4026 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4027 return readVectorImpl<scalar_type> (in, comm, map, err, tolerant, debug);
4052 static Teuchos::RCP<const map_type>
4055 const bool tolerant=
false,
4056 const bool debug=
false,
4057 const bool binary=
false)
4061 return readMap (in, comm, tolerant, debug, binary);
4066 template<
class MultiVectorScalarType>
4071 readDenseImpl (std::istream& in,
4073 Teuchos::RCP<const map_type>& map,
4074 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4075 const bool tolerant=
false,
4076 const bool debug=
false,
4077 const bool binary=
false)
4079 using Teuchos::MatrixMarket::Banner;
4080 using Teuchos::MatrixMarket::checkCommentLine;
4081 using Teuchos::ArrayRCP;
4083 using Teuchos::broadcast;
4084 using Teuchos::outArg;
4086 using Teuchos::Tuple;
4088 typedef MultiVectorScalarType ST;
4092 typedef Teuchos::ScalarTraits<ST> STS;
4093 typedef typename STS::magnitudeType MT;
4094 typedef Teuchos::ScalarTraits<MT> STM;
4099 const int myRank = comm->getRank ();
4101 if (! err.is_null ()) {
4105 *err << myRank <<
": readDenseImpl" << endl;
4107 if (! err.is_null ()) {
4142 size_t lineNumber = 1;
4145 std::ostringstream exMsg;
4146 int localBannerReadSuccess = 1;
4147 int localDimsReadSuccess = 1;
4153 *err << myRank <<
": readDenseImpl: Reading banner line (dense)" << endl;
4159 RCP<const Banner> pBanner;
4161 pBanner = readBanner (in, lineNumber, tolerant, debug);
4162 }
catch (std::exception& e) {
4164 localBannerReadSuccess = 0;
4167 if (localBannerReadSuccess) {
4168 if (pBanner->matrixType () !=
"array") {
4169 exMsg <<
"The Matrix Market file does not contain dense matrix "
4170 "data. Its banner (first) line says that its matrix type is \""
4171 << pBanner->matrixType () <<
"\", rather that the required "
4173 localBannerReadSuccess = 0;
4174 }
else if (pBanner->dataType() ==
"pattern") {
4175 exMsg <<
"The Matrix Market file's banner (first) "
4176 "line claims that the matrix's data type is \"pattern\". This does "
4177 "not make sense for a dense matrix, yet the file reports the matrix "
4178 "as dense. The only valid data types for a dense matrix are "
4179 "\"real\", \"complex\", and \"integer\".";
4180 localBannerReadSuccess = 0;
4184 dims[2] = encodeDataType (pBanner->dataType ());
4190 if (localBannerReadSuccess) {
4192 *err << myRank <<
": readDenseImpl: Reading dimensions line (dense)" << endl;
4201 bool commentLine =
true;
4203 while (commentLine) {
4206 if (in.eof () || in.fail ()) {
4207 exMsg <<
"Unable to get array dimensions line (at all) from line "
4208 << lineNumber <<
" of input stream. The input stream "
4209 <<
"claims that it is "
4210 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4211 localDimsReadSuccess = 0;
4214 if (getline (in, line)) {
4221 size_t start = 0, size = 0;
4222 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4229 std::istringstream istr (line);
4233 if (istr.eof () || istr.fail ()) {
4234 exMsg <<
"Unable to read any data from line " << lineNumber
4235 <<
" of input; the line should contain the matrix dimensions "
4236 <<
"\"<numRows> <numCols>\".";
4237 localDimsReadSuccess = 0;
4242 exMsg <<
"Failed to get number of rows from line "
4243 << lineNumber <<
" of input; the line should contains the "
4244 <<
"matrix dimensions \"<numRows> <numCols>\".";
4245 localDimsReadSuccess = 0;
4247 dims[0] = theNumRows;
4249 exMsg <<
"No more data after number of rows on line "
4250 << lineNumber <<
" of input; the line should contain the "
4251 <<
"matrix dimensions \"<numRows> <numCols>\".";
4252 localDimsReadSuccess = 0;
4257 exMsg <<
"Failed to get number of columns from line "
4258 << lineNumber <<
" of input; the line should contain "
4259 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4260 localDimsReadSuccess = 0;
4262 dims[1] = theNumCols;
4271 in.read(reinterpret_cast<char*>(&numRows),
sizeof(numRows));
4272 in.read(reinterpret_cast<char*>(&numCols),
sizeof(numCols));
4273 dims[0] = Teuchos::as<GO>(numRows);
4274 dims[1] = Teuchos::as<GO>(numCols);
4275 if ((
typeid(ST) ==
typeid(double)) || Teuchos::ScalarTraits<ST>::isOrdinal) {
4277 }
else if (Teuchos::ScalarTraits<ST>::isComplex) {
4280 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
4281 "Unrecognized Matrix Market data type. "
4282 "We should never get here. "
4283 "Please report this bug to the Tpetra developers.");
4285 localDimsReadSuccess =
true;
4286 localDimsReadSuccess =
true;
4293 Tuple<GO, 5> bannerDimsReadResult;
4295 bannerDimsReadResult[0] = dims[0];
4296 bannerDimsReadResult[1] = dims[1];
4297 bannerDimsReadResult[2] = dims[2];
4298 bannerDimsReadResult[3] = localBannerReadSuccess;
4299 bannerDimsReadResult[4] = localDimsReadSuccess;
4303 broadcast (*comm, 0, bannerDimsReadResult);
4305 TEUCHOS_TEST_FOR_EXCEPTION(
4306 bannerDimsReadResult[3] == 0, std::runtime_error,
4307 "Failed to read banner line: " << exMsg.str ());
4308 TEUCHOS_TEST_FOR_EXCEPTION(
4309 bannerDimsReadResult[4] == 0, std::runtime_error,
4310 "Failed to read matrix dimensions line: " << exMsg.str ());
4312 dims[0] = bannerDimsReadResult[0];
4313 dims[1] = bannerDimsReadResult[1];
4314 dims[2] = bannerDimsReadResult[2];
4319 const size_t numCols =
static_cast<size_t> (dims[1]);
4324 RCP<const map_type> proc0Map;
4325 if (map.is_null ()) {
4329 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4330 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4332 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4336 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4340 RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
4346 int localReadDataSuccess = 1;
4351 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4356 TEUCHOS_TEST_FOR_EXCEPTION(
4357 ! X->isConstantStride (), std::logic_error,
4358 "Can't get a 1-D view of the entries of the MultiVector X on "
4359 "Process 0, because the stride between the columns of X is not "
4360 "constant. This shouldn't happen because we just created X and "
4361 "haven't filled it in yet. Please report this bug to the Tpetra "
4368 ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4369 TEUCHOS_TEST_FOR_EXCEPTION(
4370 as<global_size_t> (X_view.size ()) < numRows * numCols,
4372 "The view of X has size " << X_view.size() <<
" which is not enough to "
4373 "accommodate the expected number of entries numRows*numCols = "
4374 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
4375 "Please report this bug to the Tpetra developers.");
4376 const size_t stride = X->getStride ();
4383 const bool isComplex = (dims[2] == 1);
4384 size_type count = 0, curRow = 0, curCol = 0;
4387 while (getline (in, line)) {
4391 size_t start = 0, size = 0;
4392 const bool commentLine =
4393 checkCommentLine (line, start, size, lineNumber, tolerant);
4394 if (! commentLine) {
4400 if (count >= X_view.size()) {
4405 TEUCHOS_TEST_FOR_EXCEPTION(
4406 count >= X_view.size(),
4408 "The Matrix Market input stream has more data in it than "
4409 "its metadata reported. Current line number is "
4410 << lineNumber <<
".");
4419 const size_t pos = line.substr (start, size).find (
':');
4420 if (pos != std::string::npos) {
4424 std::istringstream istr (line.substr (start, size));
4428 if (istr.eof() || istr.fail()) {
4435 TEUCHOS_TEST_FOR_EXCEPTION(
4436 ! tolerant && (istr.eof() || istr.fail()),
4438 "Line " << lineNumber <<
" of the Matrix Market file is "
4439 "empty, or we cannot read from it for some other reason.");
4442 ST val = STS::zero();
4445 MT real = STM::zero(), imag = STM::zero();
4462 TEUCHOS_TEST_FOR_EXCEPTION(
4463 ! tolerant && istr.eof(), std::runtime_error,
4464 "Failed to get the real part of a complex-valued matrix "
4465 "entry from line " << lineNumber <<
" of the Matrix Market "
4471 }
else if (istr.eof()) {
4472 TEUCHOS_TEST_FOR_EXCEPTION(
4473 ! tolerant && istr.eof(), std::runtime_error,
4474 "Missing imaginary part of a complex-valued matrix entry "
4475 "on line " << lineNumber <<
" of the Matrix Market file.");
4481 TEUCHOS_TEST_FOR_EXCEPTION(
4482 ! tolerant && istr.fail(), std::runtime_error,
4483 "Failed to get the imaginary part of a complex-valued "
4484 "matrix entry from line " << lineNumber <<
" of the "
4485 "Matrix Market file.");
4492 TEUCHOS_TEST_FOR_EXCEPTION(
4493 ! tolerant && istr.fail(), std::runtime_error,
4494 "Failed to get a real-valued matrix entry from line "
4495 << lineNumber <<
" of the Matrix Market file.");
4498 if (istr.fail() && tolerant) {
4504 TEUCHOS_TEST_FOR_EXCEPTION(
4505 ! tolerant && istr.fail(), std::runtime_error,
4506 "Failed to read matrix data from line " << lineNumber
4507 <<
" of the Matrix Market file.");
4510 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4512 curRow = count % numRows;
4513 curCol = count / numRows;
4514 X_view[curRow + curCol*stride] = val;
4519 TEUCHOS_TEST_FOR_EXCEPTION(
4520 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
4522 "The Matrix Market metadata reports that the dense matrix is "
4523 << numRows <<
" x " << numCols <<
", and thus has "
4524 << numRows*numCols <<
" total entries, but we only found " << count
4525 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
4526 }
catch (std::exception& e) {
4528 localReadDataSuccess = 0;
4533 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4538 TEUCHOS_TEST_FOR_EXCEPTION(
4539 ! X->isConstantStride (), std::logic_error,
4540 "Can't get a 1-D view of the entries of the MultiVector X on "
4541 "Process 0, because the stride between the columns of X is not "
4542 "constant. This shouldn't happen because we just created X and "
4543 "haven't filled it in yet. Please report this bug to the Tpetra "
4550 auto X_view = X->getLocalViewHost (Access::OverwriteAll);
4552 TEUCHOS_TEST_FOR_EXCEPTION(
4553 as<global_size_t> (X_view.extent(0)) < numRows,
4555 "The view of X has " << X_view.extent(0) <<
" rows which is not enough to "
4556 "accommodate the expected number of entries numRows = "
4558 "Please report this bug to the Tpetra developers.");
4559 TEUCHOS_TEST_FOR_EXCEPTION(
4560 as<global_size_t> (X_view.extent(1)) < numCols,
4562 "The view of X has " << X_view.extent(1) <<
" colums which is not enough to "
4563 "accommodate the expected number of entries numRows = "
4565 "Please report this bug to the Tpetra developers.");
4567 for (
size_t curRow = 0; curRow < numRows; ++curRow) {
4568 for (
size_t curCol = 0; curCol < numCols; ++curCol) {
4569 if (Teuchos::ScalarTraits<ST>::isOrdinal){
4571 in.read(reinterpret_cast<char*>(&val),
sizeof(val));
4572 X_view(curRow, curCol) = val;
4575 in.read(reinterpret_cast<char*>(&val),
sizeof(val));
4576 X_view(curRow, curCol) = val;
4584 *err << myRank <<
": readDenseImpl: done reading data" << endl;
4588 int globalReadDataSuccess = localReadDataSuccess;
4589 broadcast (*comm, 0, outArg (globalReadDataSuccess));
4590 TEUCHOS_TEST_FOR_EXCEPTION(
4591 globalReadDataSuccess == 0, std::runtime_error,
4592 "Failed to read the multivector's data: " << exMsg.str ());
4597 if (comm->getSize () == 1 && map.is_null ()) {
4599 if (! err.is_null ()) {
4603 *err << myRank <<
": readDenseImpl: done" << endl;
4605 if (! err.is_null ()) {
4612 *err << myRank <<
": readDenseImpl: Creating target MV" << endl;
4616 RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
4619 *err << myRank <<
": readDenseImpl: Creating Export" << endl;
4625 Export<LO, GO, NT> exporter (proc0Map, map, err);
4628 *err << myRank <<
": readDenseImpl: Exporting" << endl;
4631 Y->doExport (*X, exporter,
INSERT);
4633 if (! err.is_null ()) {
4637 *err << myRank <<
": readDenseImpl: done" << endl;
4639 if (! err.is_null ()) {
4648 template<
class VectorScalarType>
4653 readVectorImpl (std::istream& in,
4655 Teuchos::RCP<const map_type>& map,
4656 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4657 const bool tolerant=
false,
4658 const bool debug=
false)
4660 using Teuchos::MatrixMarket::Banner;
4661 using Teuchos::MatrixMarket::checkCommentLine;
4663 using Teuchos::broadcast;
4664 using Teuchos::outArg;
4666 using Teuchos::Tuple;
4668 typedef VectorScalarType ST;
4672 typedef Teuchos::ScalarTraits<ST> STS;
4673 typedef typename STS::magnitudeType MT;
4674 typedef Teuchos::ScalarTraits<MT> STM;
4679 const int myRank = comm->getRank ();
4681 if (! err.is_null ()) {
4685 *err << myRank <<
": readVectorImpl" << endl;
4687 if (! err.is_null ()) {
4721 size_t lineNumber = 1;
4724 std::ostringstream exMsg;
4725 int localBannerReadSuccess = 1;
4726 int localDimsReadSuccess = 1;
4731 *err << myRank <<
": readVectorImpl: Reading banner line (dense)" << endl;
4737 RCP<const Banner> pBanner;
4739 pBanner = readBanner (in, lineNumber, tolerant, debug);
4740 }
catch (std::exception& e) {
4742 localBannerReadSuccess = 0;
4745 if (localBannerReadSuccess) {
4746 if (pBanner->matrixType () !=
"array") {
4747 exMsg <<
"The Matrix Market file does not contain dense matrix "
4748 "data. Its banner (first) line says that its matrix type is \""
4749 << pBanner->matrixType () <<
"\", rather that the required "
4751 localBannerReadSuccess = 0;
4752 }
else if (pBanner->dataType() ==
"pattern") {
4753 exMsg <<
"The Matrix Market file's banner (first) "
4754 "line claims that the matrix's data type is \"pattern\". This does "
4755 "not make sense for a dense matrix, yet the file reports the matrix "
4756 "as dense. The only valid data types for a dense matrix are "
4757 "\"real\", \"complex\", and \"integer\".";
4758 localBannerReadSuccess = 0;
4762 dims[2] = encodeDataType (pBanner->dataType ());
4768 if (localBannerReadSuccess) {
4770 *err << myRank <<
": readVectorImpl: Reading dimensions line (dense)" << endl;
4779 bool commentLine =
true;
4781 while (commentLine) {
4784 if (in.eof () || in.fail ()) {
4785 exMsg <<
"Unable to get array dimensions line (at all) from line "
4786 << lineNumber <<
" of input stream. The input stream "
4787 <<
"claims that it is "
4788 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4789 localDimsReadSuccess = 0;
4792 if (getline (in, line)) {
4799 size_t start = 0, size = 0;
4800 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4807 std::istringstream istr (line);
4811 if (istr.eof () || istr.fail ()) {
4812 exMsg <<
"Unable to read any data from line " << lineNumber
4813 <<
" of input; the line should contain the matrix dimensions "
4814 <<
"\"<numRows> <numCols>\".";
4815 localDimsReadSuccess = 0;
4820 exMsg <<
"Failed to get number of rows from line "
4821 << lineNumber <<
" of input; the line should contains the "
4822 <<
"matrix dimensions \"<numRows> <numCols>\".";
4823 localDimsReadSuccess = 0;
4825 dims[0] = theNumRows;
4827 exMsg <<
"No more data after number of rows on line "
4828 << lineNumber <<
" of input; the line should contain the "
4829 <<
"matrix dimensions \"<numRows> <numCols>\".";
4830 localDimsReadSuccess = 0;
4835 exMsg <<
"Failed to get number of columns from line "
4836 << lineNumber <<
" of input; the line should contain "
4837 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4838 localDimsReadSuccess = 0;
4840 dims[1] = theNumCols;
4850 exMsg <<
"File does not contain a 1D Vector";
4851 localDimsReadSuccess = 0;
4857 Tuple<GO, 5> bannerDimsReadResult;
4859 bannerDimsReadResult[0] = dims[0];
4860 bannerDimsReadResult[1] = dims[1];
4861 bannerDimsReadResult[2] = dims[2];
4862 bannerDimsReadResult[3] = localBannerReadSuccess;
4863 bannerDimsReadResult[4] = localDimsReadSuccess;
4868 broadcast (*comm, 0, bannerDimsReadResult);
4870 TEUCHOS_TEST_FOR_EXCEPTION(
4871 bannerDimsReadResult[3] == 0, std::runtime_error,
4872 "Failed to read banner line: " << exMsg.str ());
4873 TEUCHOS_TEST_FOR_EXCEPTION(
4874 bannerDimsReadResult[4] == 0, std::runtime_error,
4875 "Failed to read matrix dimensions line: " << exMsg.str ());
4877 dims[0] = bannerDimsReadResult[0];
4878 dims[1] = bannerDimsReadResult[1];
4879 dims[2] = bannerDimsReadResult[2];
4884 const size_t numCols =
static_cast<size_t> (dims[1]);
4889 RCP<const map_type> proc0Map;
4890 if (map.is_null ()) {
4894 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4895 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4897 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4901 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4905 RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
4911 int localReadDataSuccess = 1;
4915 *err << myRank <<
": readVectorImpl: Reading matrix data (dense)"
4920 TEUCHOS_TEST_FOR_EXCEPTION(
4921 ! X->isConstantStride (), std::logic_error,
4922 "Can't get a 1-D view of the entries of the MultiVector X on "
4923 "Process 0, because the stride between the columns of X is not "
4924 "constant. This shouldn't happen because we just created X and "
4925 "haven't filled it in yet. Please report this bug to the Tpetra "
4932 Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4933 TEUCHOS_TEST_FOR_EXCEPTION(
4934 as<global_size_t> (X_view.size ()) < numRows * numCols,
4936 "The view of X has size " << X_view <<
" which is not enough to "
4937 "accommodate the expected number of entries numRows*numCols = "
4938 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
4939 "Please report this bug to the Tpetra developers.");
4940 const size_t stride = X->getStride ();
4947 const bool isComplex = (dims[2] == 1);
4948 size_type count = 0, curRow = 0, curCol = 0;
4951 while (getline (in, line)) {
4955 size_t start = 0, size = 0;
4956 const bool commentLine =
4957 checkCommentLine (line, start, size, lineNumber, tolerant);
4958 if (! commentLine) {
4964 if (count >= X_view.size()) {
4969 TEUCHOS_TEST_FOR_EXCEPTION(
4970 count >= X_view.size(),
4972 "The Matrix Market input stream has more data in it than "
4973 "its metadata reported. Current line number is "
4974 << lineNumber <<
".");
4983 const size_t pos = line.substr (start, size).find (
':');
4984 if (pos != std::string::npos) {
4988 std::istringstream istr (line.substr (start, size));
4992 if (istr.eof() || istr.fail()) {
4999 TEUCHOS_TEST_FOR_EXCEPTION(
5000 ! tolerant && (istr.eof() || istr.fail()),
5002 "Line " << lineNumber <<
" of the Matrix Market file is "
5003 "empty, or we cannot read from it for some other reason.");
5006 ST val = STS::zero();
5009 MT real = STM::zero(), imag = STM::zero();
5026 TEUCHOS_TEST_FOR_EXCEPTION(
5027 ! tolerant && istr.eof(), std::runtime_error,
5028 "Failed to get the real part of a complex-valued matrix "
5029 "entry from line " << lineNumber <<
" of the Matrix Market "
5035 }
else if (istr.eof()) {
5036 TEUCHOS_TEST_FOR_EXCEPTION(
5037 ! tolerant && istr.eof(), std::runtime_error,
5038 "Missing imaginary part of a complex-valued matrix entry "
5039 "on line " << lineNumber <<
" of the Matrix Market file.");
5045 TEUCHOS_TEST_FOR_EXCEPTION(
5046 ! tolerant && istr.fail(), std::runtime_error,
5047 "Failed to get the imaginary part of a complex-valued "
5048 "matrix entry from line " << lineNumber <<
" of the "
5049 "Matrix Market file.");
5056 TEUCHOS_TEST_FOR_EXCEPTION(
5057 ! tolerant && istr.fail(), std::runtime_error,
5058 "Failed to get a real-valued matrix entry from line "
5059 << lineNumber <<
" of the Matrix Market file.");
5062 if (istr.fail() && tolerant) {
5068 TEUCHOS_TEST_FOR_EXCEPTION(
5069 ! tolerant && istr.fail(), std::runtime_error,
5070 "Failed to read matrix data from line " << lineNumber
5071 <<
" of the Matrix Market file.");
5074 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
5076 curRow = count % numRows;
5077 curCol = count / numRows;
5078 X_view[curRow + curCol*stride] = val;
5083 TEUCHOS_TEST_FOR_EXCEPTION(
5084 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
5086 "The Matrix Market metadata reports that the dense matrix is "
5087 << numRows <<
" x " << numCols <<
", and thus has "
5088 << numRows*numCols <<
" total entries, but we only found " << count
5089 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
5090 }
catch (std::exception& e) {
5092 localReadDataSuccess = 0;
5097 *err << myRank <<
": readVectorImpl: done reading data" << endl;
5101 int globalReadDataSuccess = localReadDataSuccess;
5102 broadcast (*comm, 0, outArg (globalReadDataSuccess));
5103 TEUCHOS_TEST_FOR_EXCEPTION(
5104 globalReadDataSuccess == 0, std::runtime_error,
5105 "Failed to read the multivector's data: " << exMsg.str ());
5110 if (comm->getSize () == 1 && map.is_null ()) {
5112 if (! err.is_null ()) {
5116 *err << myRank <<
": readVectorImpl: done" << endl;
5118 if (! err.is_null ()) {
5125 *err << myRank <<
": readVectorImpl: Creating target MV" << endl;
5129 RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
5132 *err << myRank <<
": readVectorImpl: Creating Export" << endl;
5138 Export<LO, GO, NT> exporter (proc0Map, map, err);
5141 *err << myRank <<
": readVectorImpl: Exporting" << endl;
5144 Y->doExport (*X, exporter,
INSERT);
5146 if (! err.is_null ()) {
5150 *err << myRank <<
": readVectorImpl: done" << endl;
5152 if (! err.is_null ()) {
5181 static Teuchos::RCP<const map_type>
5184 const bool tolerant=
false,
5185 const bool debug=
false,
5186 const bool binary=
false)
5188 Teuchos::RCP<Teuchos::FancyOStream> err =
5189 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
5190 return readMap (in, comm, err, tolerant, debug, binary);
5220 static Teuchos::RCP<const map_type>
5223 const Teuchos::RCP<Teuchos::FancyOStream>& err,
5224 const bool tolerant=
false,
5225 const bool debug=
false,
5226 const bool binary=
false)
5228 using Teuchos::arcp;
5229 using Teuchos::Array;
5230 using Teuchos::ArrayRCP;
5232 using Teuchos::broadcast;
5233 using Teuchos::Comm;
5234 using Teuchos::CommRequest;
5235 using Teuchos::inOutArg;
5236 using Teuchos::ireceive;
5237 using Teuchos::outArg;
5239 using Teuchos::receive;
5240 using Teuchos::reduceAll;
5241 using Teuchos::REDUCE_MIN;
5242 using Teuchos::isend;
5243 using Teuchos::SerialComm;
5244 using Teuchos::toString;
5245 using Teuchos::wait;
5248 typedef ptrdiff_t int_type;
5254 const int numProcs = comm->getSize ();
5255 const int myRank = comm->getRank ();
5257 if (err.is_null ()) {
5261 std::ostringstream os;
5262 os << myRank <<
": readMap: " << endl;
5265 if (err.is_null ()) {
5271 const int sizeTag = 1339;
5273 const int dataTag = 1340;
5279 RCP<CommRequest<int> > sizeReq;
5280 RCP<CommRequest<int> > dataReq;
5285 ArrayRCP<int_type> numGidsToRecv (1);
5286 numGidsToRecv[0] = 0;
5288 sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
5291 int readSuccess = 1;
5292 std::ostringstream exMsg;
5301 RCP<const Comm<int> > proc0Comm (
new SerialComm<int> ());
5303 RCP<const map_type> dataMap;
5307 data = readDenseImpl<GO> (in, proc0Comm, dataMap, err, tolerant, debug, binary);
5309 if (data.is_null ()) {
5311 exMsg <<
"readDenseImpl() returned null." << endl;
5313 }
catch (std::exception& e) {
5315 exMsg << e.what () << endl;
5321 std::map<int, Array<GO> > pid2gids;
5326 int_type globalNumGIDs = 0;
5336 if (myRank == 0 && readSuccess == 1) {
5337 if (data->getNumVectors () == 2) {
5338 ArrayRCP<const GO> GIDs = data->getData (0);
5339 ArrayRCP<const GO> PIDs = data->getData (1);
5340 globalNumGIDs = GIDs.size ();
5344 if (globalNumGIDs > 0) {
5345 const int pid =
static_cast<int> (PIDs[0]);
5347 if (pid < 0 || pid >= numProcs) {
5349 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5350 <<
"Encountered invalid PID " << pid <<
"." << endl;
5353 const GO gid = GIDs[0];
5354 pid2gids[pid].push_back (gid);
5358 if (readSuccess == 1) {
5361 for (size_type k = 1; k < globalNumGIDs; ++k) {
5362 const int pid =
static_cast<int> (PIDs[k]);
5363 if (pid < 0 || pid >= numProcs) {
5365 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5366 <<
"Encountered invalid PID " << pid <<
"." << endl;
5369 const int_type gid = GIDs[k];
5370 pid2gids[pid].push_back (gid);
5371 if (gid < indexBase) {
5378 else if (data->getNumVectors () == 1) {
5379 if (data->getGlobalLength () % 2 !=
static_cast<GST
> (0)) {
5381 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data has the "
5382 "wrong format (for Map format 2.0). The global number of rows "
5383 "in the MultiVector must be even (divisible by 2)." << endl;
5386 ArrayRCP<const GO> theData = data->getData (0);
5387 globalNumGIDs =
static_cast<GO
> (data->getGlobalLength ()) /
5388 static_cast<GO> (2);
5392 if (globalNumGIDs > 0) {
5393 const int pid =
static_cast<int> (theData[1]);
5394 if (pid < 0 || pid >= numProcs) {
5396 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5397 <<
"Encountered invalid PID " << pid <<
"." << endl;
5400 const GO gid = theData[0];
5401 pid2gids[pid].push_back (gid);
5407 for (int_type k = 1; k < globalNumGIDs; ++k) {
5408 const int pid =
static_cast<int> (theData[2*k + 1]);
5409 if (pid < 0 || pid >= numProcs) {
5411 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5412 <<
"Encountered invalid PID " << pid <<
"." << endl;
5415 const GO gid = theData[2*k];
5416 pid2gids[pid].push_back (gid);
5417 if (gid < indexBase) {
5426 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data must have "
5427 "either 1 column (for the new Map format 2.0) or 2 columns (for "
5428 "the old Map format 1.0).";
5435 int_type readResults[3];
5436 readResults[0] =
static_cast<int_type
> (indexBase);
5437 readResults[1] =
static_cast<int_type
> (globalNumGIDs);
5438 readResults[2] =
static_cast<int_type
> (readSuccess);
5439 broadcast<int, int_type> (*comm, 0, 3, readResults);
5441 indexBase =
static_cast<GO
> (readResults[0]);
5442 globalNumGIDs =
static_cast<int_type
> (readResults[1]);
5443 readSuccess =
static_cast<int> (readResults[2]);
5449 TEUCHOS_TEST_FOR_EXCEPTION(
5450 readSuccess != 1, std::runtime_error,
5451 "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
5452 "following exception message: " << exMsg.str ());
5456 for (
int p = 1; p < numProcs; ++p) {
5457 ArrayRCP<int_type> numGidsToSend (1);
5459 auto it = pid2gids.find (p);
5460 if (it == pid2gids.end ()) {
5461 numGidsToSend[0] = 0;
5463 numGidsToSend[0] = it->second.size ();
5465 sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
5466 wait<int> (*comm, outArg (sizeReq));
5471 wait<int> (*comm, outArg (sizeReq));
5476 ArrayRCP<GO> myGids;
5477 int_type myNumGids = 0;
5479 GO* myGidsRaw = NULL;
5481 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
5482 if (it != pid2gids.end ()) {
5483 myGidsRaw = it->second.getRawPtr ();
5484 myNumGids = it->second.size ();
5486 myGids = arcp<GO> (myGidsRaw, 0, myNumGids,
false);
5490 myNumGids = numGidsToRecv[0];
5491 myGids = arcp<GO> (myNumGids);
5498 if (myNumGids > 0) {
5499 dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
5503 for (
int p = 1; p < numProcs; ++p) {
5505 ArrayRCP<GO> sendGids;
5506 GO* sendGidsRaw = NULL;
5507 int_type numSendGids = 0;
5509 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
5510 if (it != pid2gids.end ()) {
5511 numSendGids = it->second.size ();
5512 sendGidsRaw = it->second.getRawPtr ();
5513 sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids,
false);
5516 if (numSendGids > 0) {
5517 dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
5519 wait<int> (*comm, outArg (dataReq));
5521 else if (myRank == p) {
5523 wait<int> (*comm, outArg (dataReq));
5528 std::ostringstream os;
5529 os << myRank <<
": readMap: creating Map" << endl;
5532 const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
5533 RCP<const map_type> newMap;
5540 std::ostringstream errStrm;
5542 newMap = rcp (
new map_type (INVALID, myGids (), indexBase, comm));
5544 catch (std::exception& e) {
5546 errStrm <<
"Process " << comm->getRank () <<
" threw an exception: "
5547 << e.what () << std::endl;
5551 errStrm <<
"Process " << comm->getRank () <<
" threw an exception "
5552 "not a subclass of std::exception" << std::endl;
5554 Teuchos::reduceAll<int, int> (*comm, Teuchos::REDUCE_MIN,
5555 lclSuccess, Teuchos::outArg (gblSuccess));
5556 if (gblSuccess != 1) {
5559 TEUCHOS_TEST_FOR_EXCEPTION(gblSuccess != 1, std::runtime_error,
"Map constructor failed!");
5561 if (err.is_null ()) {
5565 std::ostringstream os;
5566 os << myRank <<
": readMap: done" << endl;
5569 if (err.is_null ()) {
5589 encodeDataType (
const std::string& dataType)
5591 if (dataType ==
"real" || dataType ==
"integer") {
5593 }
else if (dataType ==
"complex") {
5595 }
else if (dataType ==
"pattern") {
5601 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
5602 "Unrecognized Matrix Market data type \"" << dataType
5603 <<
"\". We should never get here. "
5604 "Please report this bug to the Tpetra developers.");
5639 static Teuchos::RCP<sparse_matrix_type>
5641 const std::string& filename_suffix,
5642 const Teuchos::RCP<const map_type>& rowMap,
5643 Teuchos::RCP<const map_type>& colMap,
5644 const Teuchos::RCP<const map_type>& domainMap,
5645 const Teuchos::RCP<const map_type>& rangeMap,
5646 const bool callFillComplete=
true,
5647 const bool tolerant=
false,
5648 const int ranksToReadAtOnce=8,
5649 const bool debug=
false)
5654 using STS =
typename Teuchos::ScalarTraits<ST>;
5656 using Teuchos::ArrayRCP;
5657 using Teuchos::arcp;
5664 TEUCHOS_TEST_FOR_EXCEPTION(
5665 rowMap.is_null (), std::invalid_argument,
5666 "Row Map must be nonnull.");
5667 Teuchos::RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
5668 TEUCHOS_TEST_FOR_EXCEPTION
5669 (comm.is_null (), std::invalid_argument,
5670 "The input row map's communicator (Teuchos::Comm object) is null.");
5671 TEUCHOS_TEST_FOR_EXCEPTION(
5672 rangeMap.is_null (), std::invalid_argument,
5673 "Range Map must be nonnull.");
5674 TEUCHOS_TEST_FOR_EXCEPTION(
5675 domainMap.is_null (), std::invalid_argument,
5676 "Domain Map must be nonnull.");
5677 TEUCHOS_TEST_FOR_EXCEPTION(
5678 domainMap->getComm().getRawPtr() != comm.getRawPtr(),
5679 std::invalid_argument,
5680 "The specified domain Map's communicator (domainMap->getComm())"
5681 "differs from row Map's communicator");
5682 TEUCHOS_TEST_FOR_EXCEPTION(
5683 rangeMap->getComm().getRawPtr() != comm.getRawPtr(),
5684 std::invalid_argument,
5685 "The specified range Map's communicator (rangeMap->getComm())"
5686 "differs from row Map's communicator");
5689 const int myRank = comm->getRank();
5690 const int numProc = comm->getSize();
5691 std::string filename = filename_prefix + std::to_string(myRank) + filename_suffix;
5694 int rank_limit = std::min(std::max(ranksToReadAtOnce,1),numProc);
5697 ArrayRCP<size_t> numEntriesPerRow;
5698 ArrayRCP<size_t> rowPtr;
5699 ArrayRCP<global_ordinal_type> colInd;
5700 ArrayRCP<scalar_type> values;
5701 std::ostringstream errMsg;
5708 bool success =
true;
5709 int bannerIsCorrect = 1, readSuccess = 1;
5710 LO numRows, numCols, numNonzeros;
5711 for(
int base_rank = 0; base_rank < numProc; base_rank += rank_limit) {
5712 int stop = std::min(base_rank+rank_limit,numProc);
5715 if(base_rank <= myRank && myRank < stop) {
5717 std::ifstream in(filename);
5718 using Teuchos::MatrixMarket::Banner;
5719 size_t lineNumber = 1;
5720 RCP<const Banner> pBanner;
5722 pBanner = readBanner (in, lineNumber, tolerant, debug);
5724 catch (std::exception& e) {
5725 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
5726 "threw an exception: " << e.what();
5727 bannerIsCorrect = 0;
5729 if (bannerIsCorrect) {
5734 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
5735 bannerIsCorrect = 0;
5736 errMsg <<
"The Matrix Market input file must contain a "
5737 "\"coordinate\"-format sparse matrix in order to create a "
5738 "Tpetra::CrsMatrix object from it, but the file's matrix "
5739 "type is \"" << pBanner->matrixType() <<
"\" instead.";
5744 if (tolerant && pBanner->matrixType() ==
"array") {
5745 bannerIsCorrect = 0;
5746 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
5747 "format sparse matrix in order to create a Tpetra::CrsMatrix "
5748 "object from it, but the file's matrix type is \"array\" "
5749 "instead. That probably means the file contains dense matrix "
5755 using Teuchos::MatrixMarket::readCoordinateDimensions;
5756 success = readCoordinateDimensions (in, numRows, numCols,
5757 numNonzeros, lineNumber,
5761 TEUCHOS_TEST_FOR_EXCEPTION(numRows != (LO)rowMap->getLocalNumElements(), std::invalid_argument,
5762 "# rows in file does not match rowmap.");
5763 TEUCHOS_TEST_FOR_EXCEPTION(!colMap.is_null() && numCols != (LO)colMap->getLocalNumElements(), std::invalid_argument,
5764 "# rows in file does not match colmap.");
5768 typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type,global_ordinal_type> raw_adder_type;
5769 bool tolerant_required =
true;
5770 Teuchos::RCP<raw_adder_type> pRaw =
5771 Teuchos::rcp (
new raw_adder_type (numRows,numCols,numNonzeros,tolerant_required,debug));
5772 RCP<adder_type> pAdder = Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
5775 std::cerr <<
"-- Reading matrix data" << std::endl;
5780 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
5782 reader_type reader (pAdder);
5785 std::pair<bool, std::vector<size_t> > results = reader.read (in, lineNumber, tolerant_required, debug);
5788 readSuccess = results.first ? 1 : 0;
5790 catch (std::exception& e) {
5797 typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,global_ordinal_type> element_type;
5800 pAdder->getAdder()->merge ();
5803 const std::vector<element_type>& entries = pAdder->getAdder()->getEntries();
5806 const size_t numEntries = (size_t)entries.size();
5809 std::cerr <<
"----- Proc "<<myRank<<
": Matrix has numRows=" << numRows
5810 <<
" rows and numEntries=" << numEntries
5811 <<
" entries." << std::endl;
5818 numEntriesPerRow = arcp<size_t> (numRows);
5819 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
5820 rowPtr = arcp<size_t> (numRows+1);
5821 std::fill (rowPtr.begin(), rowPtr.end(), 0);
5822 colInd = arcp<global_ordinal_type> (numEntries);
5823 values = arcp<scalar_type> (numEntries);
5829 LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
5831 LO indexBase = rowMap->getIndexBase();
5832 for (curPos = 0; curPos < numEntries; ++curPos) {
5833 const element_type& curEntry = entries[curPos];
5835 LO l_curRow = rowMap->getLocalElement(curRow);
5838 TEUCHOS_TEST_FOR_EXCEPTION(l_curRow == INVALID,std::logic_error,
5839 "Current global row "<< curRow <<
" is invalid.");
5841 TEUCHOS_TEST_FOR_EXCEPTION(l_curRow < l_prvRow, std::logic_error,
5842 "Row indices are out of order, even though they are supposed "
5843 "to be sorted. curRow = " << l_curRow <<
", prvRow = "
5844 << l_prvRow <<
", at curPos = " << curPos <<
". Please report "
5845 "this bug to the Tpetra developers.");
5846 if (l_curRow > l_prvRow) {
5847 for (LO r = l_prvRow+1; r <= l_curRow; ++r) {
5850 l_prvRow = l_curRow;
5852 numEntriesPerRow[l_curRow]++;
5853 colInd[curPos] = curEntry.colIndex() + indexBase;
5854 values[curPos] = curEntry.value();
5860 rowPtr[numRows] = numEntries;
5871 RCP<sparse_matrix_type> A;
5872 if(colMap.is_null()) {
5874 for(
size_t i=0; i<rowMap->getLocalNumElements(); i++) {
5875 GO g_row = rowMap->getGlobalElement(i);
5876 size_t start = rowPtr[i];
5877 size_t size = rowPtr[i+1] - rowPtr[i];
5879 A->insertGlobalValues(g_row,size,&values[start],&colInd[start]);
5884 throw std::runtime_error(
"Reading with a column map is not yet implemented");
5886 RCP<const map_type> myDomainMap = domainMap.is_null() ? rowMap : domainMap;
5887 RCP<const map_type> myRangeMap = rangeMap.is_null() ? rowMap : rangeMap;
5889 A->fillComplete(myDomainMap,myRangeMap);
5893 TEUCHOS_TEST_FOR_EXCEPTION(success ==
false, std::runtime_error,
5930 template<
class SparseMatrixType>
5935 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
6000 const std::string& matrixName,
6001 const std::string& matrixDescription,
6002 const bool debug=
false)
6005 TEUCHOS_TEST_FOR_EXCEPTION
6006 (comm.is_null (), std::invalid_argument,
6007 "The input matrix's communicator (Teuchos::Comm object) is null.");
6008 const int myRank = comm->getRank ();
6010 auto out = Writer::openOutFileOnRankZero(comm, filename, myRank,
true);
6012 writeSparse (out, matrix, matrixName, matrixDescription, debug);
6021 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6022 const std::string& matrixName,
6023 const std::string& matrixDescription,
6024 const bool debug=
false)
6026 TEUCHOS_TEST_FOR_EXCEPTION
6027 (pMatrix.is_null (), std::invalid_argument,
6028 "The input matrix is null.");
6030 matrixDescription, debug);
6055 const bool debug=
false)
6063 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6064 const bool debug=
false)
6102 const std::string& matrixName,
6103 const std::string& matrixDescription,
6104 const bool debug=
false)
6106 using Teuchos::ArrayView;
6107 using Teuchos::Comm;
6108 using Teuchos::FancyOStream;
6109 using Teuchos::getFancyOStream;
6110 using Teuchos::null;
6112 using Teuchos::rcpFromRef;
6118 using STS =
typename Teuchos::ScalarTraits<ST>;
6124 Teuchos::SetScientific<ST> sci (out);
6128 TEUCHOS_TEST_FOR_EXCEPTION(
6129 comm.is_null (), std::invalid_argument,
6130 "The input matrix's communicator (Teuchos::Comm object) is null.");
6131 const int myRank = comm->getRank ();
6134 RCP<FancyOStream> err =
6135 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6137 std::ostringstream os;
6138 os << myRank <<
": writeSparse" << endl;
6141 os <<
"-- " << myRank <<
": past barrier" << endl;
6146 const bool debugPrint = debug && myRank == 0;
6148 RCP<const map_type> rowMap = matrix.getRowMap ();
6149 RCP<const map_type> colMap = matrix.getColMap ();
6150 RCP<const map_type> domainMap = matrix.getDomainMap ();
6151 RCP<const map_type> rangeMap = matrix.getRangeMap ();
6153 const global_size_t numRows = rangeMap->getGlobalNumElements ();
6154 const global_size_t numCols = domainMap->getGlobalNumElements ();
6156 if (debug && myRank == 0) {
6157 std::ostringstream os;
6158 os <<
"-- Input sparse matrix is:"
6159 <<
"---- " << numRows <<
" x " << numCols << endl
6161 << (matrix.isGloballyIndexed() ?
"Globally" :
"Locally")
6162 <<
" indexed." << endl
6163 <<
"---- Its row map has " << rowMap->getGlobalNumElements ()
6164 <<
" elements." << endl
6165 <<
"---- Its col map has " << colMap->getGlobalNumElements ()
6166 <<
" elements." << endl;
6171 const size_t localNumRows = (myRank == 0) ? numRows : 0;
6173 std::ostringstream os;
6174 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6177 RCP<const map_type> gatherRowMap =
6178 Details::computeGatherMap (rowMap, err, debug);
6186 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6187 RCP<const map_type> gatherColMap =
6188 Details::computeGatherMap (colMap, err, debug);
6192 import_type importer (rowMap, gatherRowMap);
6197 RCP<sparse_matrix_type> newMatrix =
6199 static_cast<size_t> (0)));
6201 newMatrix->doImport (matrix, importer,
INSERT);
6206 RCP<const map_type> gatherDomainMap =
6207 rcp (
new map_type (numCols, localNumCols,
6208 domainMap->getIndexBase (),
6210 RCP<const map_type> gatherRangeMap =
6211 rcp (
new map_type (numRows, localNumRows,
6212 rangeMap->getIndexBase (),
6214 newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
6218 cerr <<
"-- Output sparse matrix is:"
6219 <<
"---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
6221 << newMatrix->getDomainMap ()->getGlobalNumElements ()
6223 << newMatrix->getGlobalNumEntries () <<
" entries;" << endl
6225 << (newMatrix->isGloballyIndexed () ?
"Globally" :
"Locally")
6226 <<
" indexed." << endl
6227 <<
"---- Its row map has "
6228 << newMatrix->getRowMap ()->getGlobalNumElements ()
6229 <<
" elements, with index base "
6230 << newMatrix->getRowMap ()->getIndexBase () <<
"." << endl
6231 <<
"---- Its col map has "
6232 << newMatrix->getColMap ()->getGlobalNumElements ()
6233 <<
" elements, with index base "
6234 << newMatrix->getColMap ()->getIndexBase () <<
"." << endl
6235 <<
"---- Element count of output matrix's column Map may differ "
6236 <<
"from that of the input matrix's column Map, if some columns "
6237 <<
"of the matrix contain no entries." << endl;
6250 out <<
"%%MatrixMarket matrix coordinate "
6251 << (STS::isComplex ?
"complex" :
"real")
6252 <<
" general" << endl;
6255 if (matrixName !=
"") {
6256 printAsComment (out, matrixName);
6258 if (matrixDescription !=
"") {
6259 printAsComment (out, matrixDescription);
6269 out << newMatrix->getRangeMap ()->getGlobalNumElements () <<
" "
6270 << newMatrix->getDomainMap ()->getGlobalNumElements () <<
" "
6271 << newMatrix->getGlobalNumEntries () << endl;
6276 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6277 const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
6285 if (newMatrix->isGloballyIndexed()) {
6288 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6289 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6290 for (GO globalRowIndex = minAllGlobalIndex;
6291 globalRowIndex <= maxAllGlobalIndex;
6293 typename sparse_matrix_type::global_inds_host_view_type ind;
6294 typename sparse_matrix_type::values_host_view_type val;
6295 newMatrix->getGlobalRowView (globalRowIndex, ind, val);
6296 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6297 const GO globalColIndex = ind(ii);
6300 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6301 << (globalColIndex + 1 - colIndexBase) <<
" ";
6302 if (STS::isComplex) {
6303 out << STS::real (val(ii)) <<
" " << STS::imag (val(ii));
6312 using OTG = Teuchos::OrdinalTraits<GO>;
6313 for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
6314 localRowIndex <= gatherRowMap->getMaxLocalIndex();
6317 const GO globalRowIndex =
6318 gatherRowMap->getGlobalElement (localRowIndex);
6319 TEUCHOS_TEST_FOR_EXCEPTION(
6320 globalRowIndex == OTG::invalid(), std::logic_error,
6321 "Failed to convert the supposed local row index "
6322 << localRowIndex <<
" into a global row index. "
6323 "Please report this bug to the Tpetra developers.");
6324 typename sparse_matrix_type::local_inds_host_view_type ind;
6325 typename sparse_matrix_type::values_host_view_type val;
6326 newMatrix->getLocalRowView (localRowIndex, ind, val);
6327 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6329 const GO globalColIndex =
6330 newMatrix->getColMap()->getGlobalElement (ind(ii));
6331 TEUCHOS_TEST_FOR_EXCEPTION(
6332 globalColIndex == OTG::invalid(), std::logic_error,
6333 "On local row " << localRowIndex <<
" of the sparse matrix: "
6334 "Failed to convert the supposed local column index "
6335 << ind(ii) <<
" into a global column index. Please report "
6336 "this bug to the Tpetra developers.");
6339 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6340 << (globalColIndex + 1 - colIndexBase) <<
" ";
6341 if (STS::isComplex) {
6342 out << STS::real (val(ii)) <<
" " << STS::imag (val(ii));
6356 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6357 const std::string& matrixName,
6358 const std::string& matrixDescription,
6359 const bool debug=
false)
6361 TEUCHOS_TEST_FOR_EXCEPTION
6362 (pMatrix.is_null (), std::invalid_argument,
6363 "The input matrix is null.");
6364 writeSparse (out, *pMatrix, matrixName, matrixDescription, debug);
6400 const std::string& graphName,
6401 const std::string& graphDescription,
6402 const bool debug=
false)
6404 using Teuchos::ArrayView;
6405 using Teuchos::Comm;
6406 using Teuchos::FancyOStream;
6407 using Teuchos::getFancyOStream;
6408 using Teuchos::null;
6410 using Teuchos::rcpFromRef;
6421 if (rowMap.is_null ()) {
6424 auto comm = rowMap->getComm ();
6425 if (comm.is_null ()) {
6428 const int myRank = comm->getRank ();
6431 RCP<FancyOStream> err =
6432 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6434 std::ostringstream os;
6435 os << myRank <<
": writeSparseGraph" << endl;
6438 os <<
"-- " << myRank <<
": past barrier" << endl;
6443 const bool debugPrint = debug && myRank == 0;
6450 const global_size_t numRows = rangeMap->getGlobalNumElements ();
6451 const global_size_t numCols = domainMap->getGlobalNumElements ();
6453 if (debug && myRank == 0) {
6454 std::ostringstream os;
6455 os <<
"-- Input sparse graph is:"
6456 <<
"---- " << numRows <<
" x " << numCols <<
" with "
6460 <<
" indexed." << endl
6461 <<
"---- Its row Map has " << rowMap->getGlobalNumElements ()
6462 <<
" elements." << endl
6463 <<
"---- Its col Map has " << colMap->getGlobalNumElements ()
6464 <<
" elements." << endl;
6469 const size_t localNumRows = (myRank == 0) ? numRows : 0;
6471 std::ostringstream os;
6472 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6475 auto gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
6483 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6484 auto gatherColMap = Details::computeGatherMap (colMap, err, debug);
6493 static_cast<size_t> (0));
6500 RCP<const map_type> gatherDomainMap =
6501 rcp (
new map_type (numCols, localNumCols,
6502 domainMap->getIndexBase (),
6504 RCP<const map_type> gatherRangeMap =
6505 rcp (
new map_type (numRows, localNumRows,
6506 rangeMap->getIndexBase (),
6508 newGraph.
fillComplete (gatherDomainMap, gatherRangeMap);
6512 cerr <<
"-- Output sparse graph is:"
6513 <<
"---- " << newGraph.
getRangeMap ()->getGlobalNumElements ()
6520 <<
" indexed." << endl
6521 <<
"---- Its row map has "
6522 << newGraph.
getRowMap ()->getGlobalNumElements ()
6523 <<
" elements, with index base "
6524 << newGraph.
getRowMap ()->getIndexBase () <<
"." << endl
6525 <<
"---- Its col map has "
6526 << newGraph.
getColMap ()->getGlobalNumElements ()
6527 <<
" elements, with index base "
6528 << newGraph.
getColMap ()->getIndexBase () <<
"." << endl
6529 <<
"---- Element count of output graph's column Map may differ "
6530 <<
"from that of the input matrix's column Map, if some columns "
6531 <<
"of the matrix contain no entries." << endl;
6545 out <<
"%%MatrixMarket matrix coordinate pattern general" << endl;
6548 if (graphName !=
"") {
6549 printAsComment (out, graphName);
6551 if (graphDescription !=
"") {
6552 printAsComment (out, graphDescription);
6563 out << newGraph.
getRangeMap ()->getGlobalNumElements () <<
" "
6564 << newGraph.
getDomainMap ()->getGlobalNumElements () <<
" "
6570 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6571 const GO colIndexBase = newGraph.
getColMap()->getIndexBase ();
6582 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6583 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6584 for (GO globalRowIndex = minAllGlobalIndex;
6585 globalRowIndex <= maxAllGlobalIndex;
6587 typename crs_graph_type::global_inds_host_view_type ind;
6589 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6590 const GO globalColIndex = ind(ii);
6593 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6594 << (globalColIndex + 1 - colIndexBase) <<
" ";
6600 typedef Teuchos::OrdinalTraits<GO> OTG;
6601 for (LO localRowIndex = gatherRowMap->getMinLocalIndex ();
6602 localRowIndex <= gatherRowMap->getMaxLocalIndex ();
6605 const GO globalRowIndex =
6606 gatherRowMap->getGlobalElement (localRowIndex);
6607 TEUCHOS_TEST_FOR_EXCEPTION
6608 (globalRowIndex == OTG::invalid (), std::logic_error,
"Failed "
6609 "to convert the supposed local row index " << localRowIndex <<
6610 " into a global row index. Please report this bug to the "
6611 "Tpetra developers.");
6612 typename crs_graph_type::local_inds_host_view_type ind;
6614 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6616 const GO globalColIndex =
6617 newGraph.
getColMap ()->getGlobalElement (ind(ii));
6618 TEUCHOS_TEST_FOR_EXCEPTION(
6619 globalColIndex == OTG::invalid(), std::logic_error,
6620 "On local row " << localRowIndex <<
" of the sparse graph: "
6621 "Failed to convert the supposed local column index "
6622 << ind(ii) <<
" into a global column index. Please report "
6623 "this bug to the Tpetra developers.");
6626 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6627 << (globalColIndex + 1 - colIndexBase) <<
" ";
6643 const bool debug=
false)
6685 const std::string& graphName,
6686 const std::string& graphDescription,
6687 const bool debug=
false)
6690 if (comm.is_null ()) {
6698 const int myRank = comm->getRank ();
6700 auto out = Writer::openOutFileOnRankZero(comm, filename, myRank,
true);
6715 const bool debug=
false)
6730 const Teuchos::RCP<const crs_graph_type>& pGraph,
6731 const std::string& graphName,
6732 const std::string& graphDescription,
6733 const bool debug=
false)
6749 const Teuchos::RCP<const crs_graph_type>& pGraph,
6750 const bool debug=
false)
6780 const bool debug=
false)
6788 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6789 const bool debug=
false)
6825 const std::string& matrixName,
6826 const std::string& matrixDescription,
6827 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6828 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6833 auto out = Writer::openOutFileOnRankZero(comm, filename, myRank,
true);
6835 writeDense (out, X, matrixName, matrixDescription, err, dbg);
6848 const Teuchos::RCP<const multivector_type>& X,
6849 const std::string& matrixName,
6850 const std::string& matrixDescription,
6851 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6852 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6854 TEUCHOS_TEST_FOR_EXCEPTION(
6855 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
6856 "writeDenseFile: The input MultiVector X is null.");
6857 writeDenseFile (filename, *X, matrixName, matrixDescription, err, dbg);
6868 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6869 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6881 const Teuchos::RCP<const multivector_type>& X,
6882 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6883 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6885 TEUCHOS_TEST_FOR_EXCEPTION(
6886 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
6887 "writeDenseFile: The input MultiVector X is null.");
6925 const std::string& matrixName,
6926 const std::string& matrixDescription,
6927 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6928 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6930 using Teuchos::Comm;
6931 using Teuchos::outArg;
6932 using Teuchos::REDUCE_MAX;
6933 using Teuchos::reduceAll;
6942 const bool debug = ! dbg.is_null ();
6945 std::ostringstream os;
6946 os << myRank <<
": writeDense" << endl;
6951 writeDenseHeader (out, X, matrixName, matrixDescription, err, dbg);
6956 for (
size_t j = 0; j < numVecs; ++j) {
6957 writeDenseColumn (out, * (X.
getVector (j)), err, dbg);
6962 std::ostringstream os;
6963 os << myRank <<
": writeDense: Done" << endl;
6975 static std::ofstream openOutFileOnRankZero(
6977 const std::string& filename,
const int rank,
const bool safe =
true,
6978 const std::ios_base::openmode mode = std::ios_base::out
6984 int all_should_stop = 0;
6988 out.open(filename, mode);
6989 all_should_stop = !out && safe;
6993 if(comm) Teuchos::broadcast(*comm, 0, &all_should_stop);
6995 TEUCHOS_TEST_FOR_EXCEPTION(
6998 "Could not open output file '" + filename +
"' on root rank 0."
7030 writeDenseHeader (std::ostream& out,
7032 const std::string& matrixName,
7033 const std::string& matrixDescription,
7034 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7035 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7037 using Teuchos::Comm;
7038 using Teuchos::outArg;
7040 using Teuchos::REDUCE_MAX;
7041 using Teuchos::reduceAll;
7043 typedef Teuchos::ScalarTraits<scalar_type> STS;
7044 const char prefix[] =
"Tpetra::MatrixMarket::writeDenseHeader: ";
7054 const bool debug = ! dbg.is_null ();
7058 std::ostringstream os;
7059 os << myRank <<
": writeDenseHeader" << endl;
7078 std::ostringstream hdr;
7080 std::string dataType;
7081 if (STS::isComplex) {
7082 dataType =
"complex";
7083 }
else if (STS::isOrdinal) {
7084 dataType =
"integer";
7088 hdr <<
"%%MatrixMarket matrix array " << dataType <<
" general"
7093 if (matrixName !=
"") {
7094 printAsComment (hdr, matrixName);
7096 if (matrixDescription !=
"") {
7097 printAsComment (hdr, matrixDescription);
7100 hdr << X.getGlobalLength () <<
" " << X.getNumVectors () << endl;
7104 }
catch (std::exception& e) {
7105 if (! err.is_null ()) {
7106 *err << prefix <<
"While writing the Matrix Market header, "
7107 "Process 0 threw an exception: " << e.what () << endl;
7116 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
7117 TEUCHOS_TEST_FOR_EXCEPTION(
7118 gblErr == 1, std::runtime_error, prefix <<
"Some error occurred "
7119 "which prevented this method from completing.");
7123 *dbg << myRank <<
": writeDenseHeader: Done" << endl;
7146 writeDenseColumn (std::ostream& out,
7148 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7149 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7151 using Teuchos::arcp;
7152 using Teuchos::Array;
7153 using Teuchos::ArrayRCP;
7154 using Teuchos::ArrayView;
7155 using Teuchos::Comm;
7156 using Teuchos::CommRequest;
7157 using Teuchos::ireceive;
7158 using Teuchos::isend;
7159 using Teuchos::outArg;
7160 using Teuchos::REDUCE_MAX;
7161 using Teuchos::reduceAll;
7163 using Teuchos::TypeNameTraits;
7164 using Teuchos::wait;
7166 typedef Teuchos::ScalarTraits<scalar_type> STS;
7168 const Comm<int>& comm = * (X.getMap ()->getComm ());
7169 const int myRank = comm.getRank ();
7170 const int numProcs = comm.getSize ();
7177 const bool debug = ! dbg.is_null ();
7181 std::ostringstream os;
7182 os << myRank <<
": writeDenseColumn" << endl;
7191 Teuchos::SetScientific<scalar_type> sci (out);
7193 const size_t myNumRows = X.getLocalLength ();
7194 const size_t numCols = X.getNumVectors ();
7197 const int sizeTag = 1337;
7198 const int dataTag = 1338;
7259 RCP<CommRequest<int> > sendReqSize, sendReqData;
7265 Array<ArrayRCP<size_t> > recvSizeBufs (3);
7266 Array<ArrayRCP<scalar_type> > recvDataBufs (3);
7267 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7268 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7271 ArrayRCP<size_t> sendDataSize (1);
7272 sendDataSize[0] = myNumRows;
7276 std::ostringstream os;
7277 os << myRank <<
": Post receive-size receives from "
7278 "Procs 1 and 2: tag = " << sizeTag << endl;
7282 recvSizeBufs[0].resize (1);
7286 (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7287 recvSizeBufs[1].resize (1);
7288 (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7289 recvSizeBufs[2].resize (1);
7290 (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7293 ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
7297 ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
7300 else if (myRank == 1 || myRank == 2) {
7302 std::ostringstream os;
7303 os << myRank <<
": Post send-size send: size = "
7304 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7311 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7315 std::ostringstream os;
7316 os << myRank <<
": Not posting my send-size send yet" << endl;
7325 std::ostringstream os;
7326 os << myRank <<
": Pack my entries" << endl;
7329 ArrayRCP<scalar_type> sendDataBuf;
7331 sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
7332 X.get1dCopy (sendDataBuf (), myNumRows);
7334 catch (std::exception& e) {
7336 if (! err.is_null ()) {
7337 std::ostringstream os;
7338 os <<
"Process " << myRank <<
": Attempt to pack my MultiVector "
7339 "entries threw an exception: " << e.what () << endl;
7344 std::ostringstream os;
7345 os << myRank <<
": Done packing my entries" << endl;
7354 *dbg << myRank <<
": Post send-data send: tag = " << dataTag
7357 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7365 std::ostringstream os;
7366 os << myRank <<
": Write my entries" << endl;
7372 const size_t printNumRows = myNumRows;
7373 ArrayView<const scalar_type> printData = sendDataBuf ();
7374 const size_t printStride = printNumRows;
7375 if (static_cast<size_t> (printData.size ()) < printStride * numCols) {
7377 if (! err.is_null ()) {
7378 std::ostringstream os;
7379 os <<
"Process " << myRank <<
": My MultiVector data's size "
7380 << printData.size () <<
" does not match my local dimensions "
7381 << printStride <<
" x " << numCols <<
"." << endl;
7389 for (
size_t col = 0; col < numCols; ++col) {
7390 for (
size_t row = 0; row < printNumRows; ++row) {
7391 if (STS::isComplex) {
7392 out << STS::real (printData[row + col * printStride]) <<
" "
7393 << STS::imag (printData[row + col * printStride]) << endl;
7395 out << printData[row + col * printStride] << endl;
7404 const int recvRank = 1;
7405 const int circBufInd = recvRank % 3;
7407 std::ostringstream os;
7408 os << myRank <<
": Wait on receive-size receive from Process "
7409 << recvRank << endl;
7413 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7417 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7418 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7420 if (! err.is_null ()) {
7421 std::ostringstream os;
7422 os << myRank <<
": Result of receive-size receive from Process "
7423 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7424 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7425 "This should never happen, and suggests that the receive never "
7426 "got posted. Please report this bug to the Tpetra developers."
7441 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7443 std::ostringstream os;
7444 os << myRank <<
": Post receive-data receive from Process "
7445 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7446 << recvDataBufs[circBufInd].size () << endl;
7449 if (! recvSizeReqs[circBufInd].is_null ()) {
7451 if (! err.is_null ()) {
7452 std::ostringstream os;
7453 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7454 "null, before posting the receive-data receive from Process "
7455 << recvRank <<
". This should never happen. Please report "
7456 "this bug to the Tpetra developers." << endl;
7460 recvDataReqs[circBufInd] =
7461 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7462 recvRank, dataTag, comm);
7465 else if (myRank == 1) {
7468 std::ostringstream os;
7469 os << myRank <<
": Wait on my send-size send" << endl;
7472 wait<int> (comm, outArg (sendReqSize));
7478 for (
int p = 1; p < numProcs; ++p) {
7480 if (p + 2 < numProcs) {
7482 const int recvRank = p + 2;
7483 const int circBufInd = recvRank % 3;
7485 std::ostringstream os;
7486 os << myRank <<
": Post receive-size receive from Process "
7487 << recvRank <<
": tag = " << sizeTag << endl;
7490 if (! recvSizeReqs[circBufInd].is_null ()) {
7492 if (! err.is_null ()) {
7493 std::ostringstream os;
7494 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7495 <<
"null, for the receive-size receive from Process "
7496 << recvRank <<
"! This may mean that this process never "
7497 <<
"finished waiting for the receive from Process "
7498 << (recvRank - 3) <<
"." << endl;
7502 recvSizeReqs[circBufInd] =
7503 ireceive<int, size_t> (recvSizeBufs[circBufInd],
7504 recvRank, sizeTag, comm);
7507 if (p + 1 < numProcs) {
7508 const int recvRank = p + 1;
7509 const int circBufInd = recvRank % 3;
7513 std::ostringstream os;
7514 os << myRank <<
": Wait on receive-size receive from Process "
7515 << recvRank << endl;
7518 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7522 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7523 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7525 if (! err.is_null ()) {
7526 std::ostringstream os;
7527 os << myRank <<
": Result of receive-size receive from Process "
7528 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7529 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7530 "This should never happen, and suggests that the receive never "
7531 "got posted. Please report this bug to the Tpetra developers."
7545 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7547 std::ostringstream os;
7548 os << myRank <<
": Post receive-data receive from Process "
7549 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7550 << recvDataBufs[circBufInd].size () << endl;
7553 if (! recvDataReqs[circBufInd].is_null ()) {
7555 if (! err.is_null ()) {
7556 std::ostringstream os;
7557 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
7558 <<
"null, for the receive-data receive from Process "
7559 << recvRank <<
"! This may mean that this process never "
7560 <<
"finished waiting for the receive from Process "
7561 << (recvRank - 3) <<
"." << endl;
7565 recvDataReqs[circBufInd] =
7566 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7567 recvRank, dataTag, comm);
7571 const int recvRank = p;
7572 const int circBufInd = recvRank % 3;
7574 std::ostringstream os;
7575 os << myRank <<
": Wait on receive-data receive from Process "
7576 << recvRank << endl;
7579 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
7586 std::ostringstream os;
7587 os << myRank <<
": Write entries from Process " << recvRank
7589 *dbg << os.str () << endl;
7591 size_t printNumRows = (recvSizeBufs[circBufInd])[0];
7592 if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7594 if (! err.is_null ()) {
7595 std::ostringstream os;
7596 os << myRank <<
": Result of receive-size receive from Process "
7597 << recvRank <<
" was Teuchos::OrdinalTraits<size_t>::"
7598 "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
7599 <<
". This should never happen, and suggests that its "
7600 "receive-size receive was never posted. "
7601 "Please report this bug to the Tpetra developers." << endl;
7608 if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
7610 if (! err.is_null ()) {
7611 std::ostringstream os;
7612 os << myRank <<
": Result of receive-size receive from Proc "
7613 << recvRank <<
" was " << printNumRows <<
" > 0, but "
7614 "recvDataBufs[" << circBufInd <<
"] is null. This should "
7615 "never happen. Please report this bug to the Tpetra "
7616 "developers." << endl;
7626 ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
7627 const size_t printStride = printNumRows;
7631 for (
size_t col = 0; col < numCols; ++col) {
7632 for (
size_t row = 0; row < printNumRows; ++row) {
7633 if (STS::isComplex) {
7634 out << STS::real (printData[row + col * printStride]) <<
" "
7635 << STS::imag (printData[row + col * printStride]) << endl;
7637 out << printData[row + col * printStride] << endl;
7642 else if (myRank == p) {
7645 std::ostringstream os;
7646 os << myRank <<
": Wait on my send-data send" << endl;
7649 wait<int> (comm, outArg (sendReqData));
7651 else if (myRank == p + 1) {
7654 std::ostringstream os;
7655 os << myRank <<
": Post send-data send: tag = " << dataTag
7659 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7662 std::ostringstream os;
7663 os << myRank <<
": Wait on my send-size send" << endl;
7666 wait<int> (comm, outArg (sendReqSize));
7668 else if (myRank == p + 2) {
7671 std::ostringstream os;
7672 os << myRank <<
": Post send-size send: size = "
7673 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7676 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7681 reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
7682 TEUCHOS_TEST_FOR_EXCEPTION(
7683 gblErr == 1, std::runtime_error,
"Tpetra::MatrixMarket::writeDense "
7684 "experienced some kind of error and was unable to complete.");
7688 *dbg << myRank <<
": writeDenseColumn: Done" << endl;
7702 const Teuchos::RCP<const multivector_type>& X,
7703 const std::string& matrixName,
7704 const std::string& matrixDescription,
7705 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7706 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7708 TEUCHOS_TEST_FOR_EXCEPTION(
7709 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7710 "writeDense: The input MultiVector X is null.");
7711 writeDense (out, *X, matrixName, matrixDescription, err, dbg);
7722 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7723 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7735 const Teuchos::RCP<const multivector_type>& X,
7736 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7737 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7739 TEUCHOS_TEST_FOR_EXCEPTION(
7740 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7741 "writeDense: The input MultiVector X is null.");
7767 Teuchos::RCP<Teuchos::FancyOStream> err =
7768 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
7783 const Teuchos::RCP<Teuchos::FancyOStream>& err,
7784 const bool debug=
false)
7786 using Teuchos::Array;
7787 using Teuchos::ArrayRCP;
7788 using Teuchos::ArrayView;
7789 using Teuchos::Comm;
7790 using Teuchos::CommRequest;
7791 using Teuchos::ireceive;
7792 using Teuchos::isend;
7794 using Teuchos::TypeNameTraits;
7795 using Teuchos::wait;
7798 typedef int pid_type;
7813 typedef ptrdiff_t int_type;
7814 TEUCHOS_TEST_FOR_EXCEPTION(
7815 sizeof (GO) >
sizeof (int_type), std::logic_error,
7816 "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
7817 <<
" is too big for ptrdiff_t. sizeof(GO) = " <<
sizeof (GO)
7818 <<
" > sizeof(ptrdiff_t) = " <<
sizeof (ptrdiff_t) <<
".");
7819 TEUCHOS_TEST_FOR_EXCEPTION(
7820 sizeof (pid_type) >
sizeof (int_type), std::logic_error,
7821 "The (MPI) process rank type pid_type=" <<
7822 TypeNameTraits<pid_type>::name () <<
" is too big for ptrdiff_t. "
7823 "sizeof(pid_type) = " <<
sizeof (pid_type) <<
" > sizeof(ptrdiff_t)"
7824 " = " <<
sizeof (ptrdiff_t) <<
".");
7826 const Comm<int>& comm = * (map.
getComm ());
7827 const int myRank = comm.getRank ();
7828 const int numProcs = comm.getSize ();
7830 if (! err.is_null ()) {
7834 std::ostringstream os;
7835 os << myRank <<
": writeMap" << endl;
7838 if (! err.is_null ()) {
7845 const int sizeTag = 1337;
7846 const int dataTag = 1338;
7905 RCP<CommRequest<int> > sendReqSize, sendReqData;
7911 Array<ArrayRCP<int_type> > recvSizeBufs (3);
7912 Array<ArrayRCP<int_type> > recvDataBufs (3);
7913 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7914 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7917 ArrayRCP<int_type> sendDataSize (1);
7918 sendDataSize[0] = myNumRows;
7922 std::ostringstream os;
7923 os << myRank <<
": Post receive-size receives from "
7924 "Procs 1 and 2: tag = " << sizeTag << endl;
7928 recvSizeBufs[0].resize (1);
7929 (recvSizeBufs[0])[0] = -1;
7930 recvSizeBufs[1].resize (1);
7931 (recvSizeBufs[1])[0] = -1;
7932 recvSizeBufs[2].resize (1);
7933 (recvSizeBufs[2])[0] = -1;
7936 ireceive<int, int_type> (recvSizeBufs[1], 1, sizeTag, comm);
7940 ireceive<int, int_type> (recvSizeBufs[2], 2, sizeTag, comm);
7943 else if (myRank == 1 || myRank == 2) {
7945 std::ostringstream os;
7946 os << myRank <<
": Post send-size send: size = "
7947 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7954 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
7958 std::ostringstream os;
7959 os << myRank <<
": Not posting my send-size send yet" << endl;
7970 std::ostringstream os;
7971 os << myRank <<
": Pack my GIDs and PIDs" << endl;
7975 ArrayRCP<int_type> sendDataBuf (myNumRows * 2);
7978 const int_type myMinGblIdx =
7980 for (
size_t k = 0; k < myNumRows; ++k) {
7981 const int_type gid = myMinGblIdx +
static_cast<int_type
> (k);
7982 const int_type pid =
static_cast<int_type
> (myRank);
7983 sendDataBuf[2*k] = gid;
7984 sendDataBuf[2*k+1] = pid;
7989 for (
size_t k = 0; k < myNumRows; ++k) {
7990 const int_type gid =
static_cast<int_type
> (myGblInds[k]);
7991 const int_type pid =
static_cast<int_type
> (myRank);
7992 sendDataBuf[2*k] = gid;
7993 sendDataBuf[2*k+1] = pid;
7998 std::ostringstream os;
7999 os << myRank <<
": Done packing my GIDs and PIDs" << endl;
8006 *err << myRank <<
": Post send-data send: tag = " << dataTag
8009 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
8014 *err << myRank <<
": Write MatrixMarket header" << endl;
8019 std::ostringstream hdr;
8023 hdr <<
"%%MatrixMarket matrix array integer general" << endl
8024 <<
"% Format: Version 2.0" << endl
8026 <<
"% This file encodes a Tpetra::Map." << endl
8027 <<
"% It is stored as a dense vector, with twice as many " << endl
8028 <<
"% entries as the global number of GIDs (global indices)." << endl
8029 <<
"% (GID, PID) pairs are stored contiguously, where the PID " << endl
8030 <<
"% is the rank of the process owning that GID." << endl
8035 std::ostringstream os;
8036 os << myRank <<
": Write my GIDs and PIDs" << endl;
8042 const int_type printNumRows = myNumRows;
8043 ArrayView<const int_type> printData = sendDataBuf ();
8044 for (int_type k = 0; k < printNumRows; ++k) {
8045 const int_type gid = printData[2*k];
8046 const int_type pid = printData[2*k+1];
8047 out << gid << endl << pid << endl;
8053 const int recvRank = 1;
8054 const int circBufInd = recvRank % 3;
8056 std::ostringstream os;
8057 os << myRank <<
": Wait on receive-size receive from Process "
8058 << recvRank << endl;
8062 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
8066 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
8067 if (debug && recvNumRows == -1) {
8068 std::ostringstream os;
8069 os << myRank <<
": Result of receive-size receive from Process "
8070 << recvRank <<
" is -1. This should never happen, and "
8071 "suggests that the receive never got posted. Please report "
8072 "this bug to the Tpetra developers." << endl;
8077 recvDataBufs[circBufInd].resize (recvNumRows * 2);
8079 std::ostringstream os;
8080 os << myRank <<
": Post receive-data receive from Process "
8081 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
8082 << recvDataBufs[circBufInd].size () << endl;
8085 if (! recvSizeReqs[circBufInd].is_null ()) {
8086 std::ostringstream os;
8087 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
8088 "null, before posting the receive-data receive from Process "
8089 << recvRank <<
". This should never happen. Please report "
8090 "this bug to the Tpetra developers." << endl;
8093 recvDataReqs[circBufInd] =
8094 ireceive<int, int_type> (recvDataBufs[circBufInd],
8095 recvRank, dataTag, comm);
8098 else if (myRank == 1) {
8101 std::ostringstream os;
8102 os << myRank <<
": Wait on my send-size send" << endl;
8105 wait<int> (comm, outArg (sendReqSize));
8111 for (
int p = 1; p < numProcs; ++p) {
8113 if (p + 2 < numProcs) {
8115 const int recvRank = p + 2;
8116 const int circBufInd = recvRank % 3;
8118 std::ostringstream os;
8119 os << myRank <<
": Post receive-size receive from Process "
8120 << recvRank <<
": tag = " << sizeTag << endl;
8123 if (! recvSizeReqs[circBufInd].is_null ()) {
8124 std::ostringstream os;
8125 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
8126 <<
"null, for the receive-size receive from Process "
8127 << recvRank <<
"! This may mean that this process never "
8128 <<
"finished waiting for the receive from Process "
8129 << (recvRank - 3) <<
"." << endl;
8132 recvSizeReqs[circBufInd] =
8133 ireceive<int, int_type> (recvSizeBufs[circBufInd],
8134 recvRank, sizeTag, comm);
8137 if (p + 1 < numProcs) {
8138 const int recvRank = p + 1;
8139 const int circBufInd = recvRank % 3;
8143 std::ostringstream os;
8144 os << myRank <<
": Wait on receive-size receive from Process "
8145 << recvRank << endl;
8148 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
8152 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
8153 if (debug && recvNumRows == -1) {
8154 std::ostringstream os;
8155 os << myRank <<
": Result of receive-size receive from Process "
8156 << recvRank <<
" is -1. This should never happen, and "
8157 "suggests that the receive never got posted. Please report "
8158 "this bug to the Tpetra developers." << endl;
8163 recvDataBufs[circBufInd].resize (recvNumRows * 2);
8165 std::ostringstream os;
8166 os << myRank <<
": Post receive-data receive from Process "
8167 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
8168 << recvDataBufs[circBufInd].size () << endl;
8171 if (! recvDataReqs[circBufInd].is_null ()) {
8172 std::ostringstream os;
8173 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
8174 <<
"null, for the receive-data receive from Process "
8175 << recvRank <<
"! This may mean that this process never "
8176 <<
"finished waiting for the receive from Process "
8177 << (recvRank - 3) <<
"." << endl;
8180 recvDataReqs[circBufInd] =
8181 ireceive<int, int_type> (recvDataBufs[circBufInd],
8182 recvRank, dataTag, comm);
8186 const int recvRank = p;
8187 const int circBufInd = recvRank % 3;
8189 std::ostringstream os;
8190 os << myRank <<
": Wait on receive-data receive from Process "
8191 << recvRank << endl;
8194 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
8201 std::ostringstream os;
8202 os << myRank <<
": Write GIDs and PIDs from Process "
8203 << recvRank << endl;
8204 *err << os.str () << endl;
8206 const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
8207 if (debug && printNumRows == -1) {
8208 std::ostringstream os;
8209 os << myRank <<
": Result of receive-size receive from Process "
8210 << recvRank <<
" was -1. This should never happen, and "
8211 "suggests that its receive-size receive was never posted. "
8212 "Please report this bug to the Tpetra developers." << endl;
8215 if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
8216 std::ostringstream os;
8217 os << myRank <<
": Result of receive-size receive from Proc "
8218 << recvRank <<
" was " << printNumRows <<
" > 0, but "
8219 "recvDataBufs[" << circBufInd <<
"] is null. This should "
8220 "never happen. Please report this bug to the Tpetra "
8221 "developers." << endl;
8224 ArrayView<const int_type> printData = (recvDataBufs[circBufInd]) ();
8225 for (int_type k = 0; k < printNumRows; ++k) {
8226 const int_type gid = printData[2*k];
8227 const int_type pid = printData[2*k+1];
8228 out << gid << endl << pid << endl;
8231 else if (myRank == p) {
8234 std::ostringstream os;
8235 os << myRank <<
": Wait on my send-data send" << endl;
8238 wait<int> (comm, outArg (sendReqData));
8240 else if (myRank == p + 1) {
8243 std::ostringstream os;
8244 os << myRank <<
": Post send-data send: tag = " << dataTag
8248 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
8251 std::ostringstream os;
8252 os << myRank <<
": Wait on my send-size send" << endl;
8255 wait<int> (comm, outArg (sendReqSize));
8257 else if (myRank == p + 2) {
8260 std::ostringstream os;
8261 os << myRank <<
": Post send-size send: size = "
8262 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
8265 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
8269 if (! err.is_null ()) {
8273 *err << myRank <<
": writeMap: Done" << endl;
8275 if (! err.is_null ()) {
8285 const int myRank = map.
getComm ()->getRank ();
8287 auto out = Writer::openOutFileOnRankZero(map.
getComm(), filename, myRank,
true);
8320 printAsComment (std::ostream& out,
const std::string& str)
8323 std::istringstream inpstream (str);
8326 while (getline (inpstream, line)) {
8327 if (! line.empty()) {
8330 if (line[0] ==
'%') {
8331 out << line << endl;
8334 out <<
"%% " << line << endl;
8362 Teuchos::ParameterList pl;
8388 Teuchos::ParameterList pl;
8431 const Teuchos::ParameterList& params)
8434 std::string tmpFile =
"__TMP__" + fileName;
8435 const int myRank = A.
getDomainMap()->getComm()->getRank();
8436 bool precisionChanged=
false;
8446 if (std::ifstream(tmpFile))
8447 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
8448 "writeOperator: temporary file " << tmpFile <<
" already exists");
8449 out.open(tmpFile.c_str());
8450 if (params.isParameter(
"precision")) {
8451 oldPrecision = out.precision(params.get<
int>(
"precision"));
8452 precisionChanged=
true;
8456 const std::string header = writeOperatorImpl(out, A, params);
8459 if (precisionChanged)
8460 out.precision(oldPrecision);
8462 out.open(fileName.c_str(), std::ios::binary);
8463 bool printMatrixMarketHeader =
true;
8464 if (params.isParameter(
"print MatrixMarket header"))
8465 printMatrixMarketHeader = params.get<
bool>(
"print MatrixMarket header");
8466 if (printMatrixMarketHeader && myRank == 0) {
8471 std::ifstream src(tmpFile, std::ios_base::binary);
8475 remove(tmpFile.c_str());
8520 const Teuchos::ParameterList& params)
8522 const int myRank = A.
getDomainMap ()->getComm ()->getRank ();
8539 std::ostringstream tmpOut;
8541 if (params.isParameter (
"precision") && params.isType<
int> (
"precision")) {
8542 (void) tmpOut.precision (params.get<
int> (
"precision"));
8546 const std::string header = writeOperatorImpl (tmpOut, A, params);
8549 bool printMatrixMarketHeader =
true;
8550 if (params.isParameter (
"print MatrixMarket header") &&
8551 params.isType<
bool> (
"print MatrixMarket header")) {
8552 printMatrixMarketHeader = params.get<
bool> (
"print MatrixMarket header");
8554 if (printMatrixMarketHeader && myRank == 0) {
8570 out << tmpOut.str ();
8584 writeOperatorImpl (std::ostream& os,
8585 const operator_type& A,
8586 const Teuchos::ParameterList& params)
8590 using Teuchos::ArrayRCP;
8591 using Teuchos::Array;
8595 typedef scalar_type Scalar;
8596 typedef Teuchos::OrdinalTraits<LO> TLOT;
8597 typedef Teuchos::OrdinalTraits<GO> TGOT;
8601 const map_type& domainMap = *(A.getDomainMap());
8602 RCP<const map_type> rangeMap = A.getRangeMap();
8604 const int myRank = comm->getRank();
8605 const size_t numProcs = comm->getSize();
8608 if (params.isParameter(
"probing size"))
8609 numMVs = params.get<
int>(
"probing size");
8612 GO minColGid = domainMap.getMinAllGlobalIndex();
8613 GO maxColGid = domainMap.getMaxAllGlobalIndex();
8618 GO numGlobElts = maxColGid - minColGid + TGOT::one();
8619 GO numChunks = numGlobElts / numMVs;
8620 GO rem = numGlobElts % numMVs;
8621 GO indexBase = rangeMap->getIndexBase();
8623 int offsetToUseInPrinting = 1 - indexBase;
8624 if (params.isParameter(
"zero-based indexing")) {
8625 if (params.get<
bool>(
"zero-based indexing") ==
true)
8626 offsetToUseInPrinting = -indexBase;
8630 size_t numLocalRangeEntries = rangeMap->getLocalNumElements();
8633 RCP<const map_type> allGidsMap = rcp(
new map_type(TGOT::invalid(), numLocalRangeEntries,
8636 mv_type_go allGids(allGidsMap,1);
8637 Teuchos::ArrayRCP<GO> allGidsData = allGids.getDataNonConst(0);
8639 for (
size_t i=0; i<numLocalRangeEntries; i++)
8640 allGidsData[i] = rangeMap->getGlobalElement(i);
8641 allGidsData = Teuchos::null;
8644 GO numTargetMapEntries=TGOT::zero();
8645 Teuchos::Array<GO> importGidList;
8647 numTargetMapEntries = rangeMap->getGlobalNumElements();
8648 importGidList.reserve(numTargetMapEntries);
8649 for (GO j=0; j<numTargetMapEntries; ++j) importGidList.push_back(j + indexBase);
8651 importGidList.reserve(numTargetMapEntries);
8653 RCP<map_type> importGidMap = rcp(
new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
8656 import_type gidImporter(allGidsMap, importGidMap);
8657 mv_type_go importedGids(importGidMap, 1);
8658 importedGids.doImport(allGids, gidImporter,
INSERT);
8661 ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
8662 RCP<const map_type> importMap = rcp(
new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm) );
8665 import_type importer(rangeMap, importMap);
8667 RCP<mv_type> colsOnPid0 = rcp(
new mv_type(importMap,numMVs));
8669 RCP<mv_type> ei = rcp(
new mv_type(A.getDomainMap(),numMVs));
8670 RCP<mv_type> colsA = rcp(
new mv_type(A.getRangeMap(),numMVs));
8672 Array<GO> globalColsArray, localColsArray;
8673 globalColsArray.reserve(numMVs);
8674 localColsArray.reserve(numMVs);
8676 ArrayRCP<ArrayRCP<Scalar> > eiData(numMVs);
8677 for (
size_t i=0; i<numMVs; ++i)
8678 eiData[i] = ei->getDataNonConst(i);
8683 for (GO k=0; k<numChunks; ++k) {
8684 for (
size_t j=0; j<numMVs; ++j ) {
8686 GO curGlobalCol = minColGid + k*numMVs + j;
8687 globalColsArray.push_back(curGlobalCol);
8689 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8690 if (curLocalCol != TLOT::invalid()) {
8691 eiData[j][curLocalCol] = TGOT::one();
8692 localColsArray.push_back(curLocalCol);
8697 for (
size_t i=0; i<numMVs; ++i)
8698 eiData[i] = Teuchos::null;
8700 A.apply(*ei,*colsA);
8702 colsOnPid0->doImport(*colsA,importer,
INSERT);
8705 globalNnz += writeColumns(os,*colsOnPid0, numMVs, importedGidsData(),
8706 globalColsArray, offsetToUseInPrinting);
8709 for (
size_t i=0; i<numMVs; ++i)
8710 eiData[i] = ei->getDataNonConst(i);
8711 for (
size_t j=0; j<numMVs; ++j ) {
8712 GO curGlobalCol = minColGid + k*numMVs + j;
8714 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8715 if (curLocalCol != TLOT::invalid()) {
8716 eiData[j][curLocalCol] = TGOT::one();
8721 for (
size_t j=0; j<numMVs; ++j ) {
8722 for (
int i=0; i<localColsArray.size(); ++i)
8723 eiData[j][localColsArray[i]] = TGOT::zero();
8725 globalColsArray.clear();
8726 localColsArray.clear();
8734 for (
int j=0; j<rem; ++j ) {
8735 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8736 globalColsArray.push_back(curGlobalCol);
8738 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8739 if (curLocalCol != TLOT::invalid()) {
8740 eiData[j][curLocalCol] = TGOT::one();
8741 localColsArray.push_back(curLocalCol);
8746 for (
size_t i=0; i<numMVs; ++i)
8747 eiData[i] = Teuchos::null;
8749 A.apply(*ei,*colsA);
8751 colsOnPid0->doImport(*colsA,importer,
INSERT);
8753 globalNnz += writeColumns(os,*colsOnPid0, rem, importedGidsData(),
8754 globalColsArray, offsetToUseInPrinting);
8757 for (
size_t i=0; i<numMVs; ++i)
8758 eiData[i] = ei->getDataNonConst(i);
8759 for (
int j=0; j<rem; ++j ) {
8760 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8762 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8763 if (curLocalCol != TLOT::invalid()) {
8764 eiData[j][curLocalCol] = TGOT::one();
8769 for (
int j=0; j<rem; ++j ) {
8770 for (
int i=0; i<localColsArray.size(); ++i)
8771 eiData[j][localColsArray[i]] = TGOT::zero();
8773 globalColsArray.clear();
8774 localColsArray.clear();
8783 std::ostringstream oss;
8785 oss <<
"%%MatrixMarket matrix coordinate ";
8786 if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
8791 oss <<
" general" << std::endl;
8792 oss <<
"% Tpetra::Operator" << std::endl;
8793 std::time_t now = std::time(NULL);
8794 oss <<
"% time stamp: " << ctime(&now);
8795 oss <<
"% written from " << numProcs <<
" processes" << std::endl;
8796 size_t numRows = rangeMap->getGlobalNumElements();
8797 size_t numCols = domainMap.getGlobalNumElements();
8798 oss << numRows <<
" " << numCols <<
" " << globalNnz << std::endl;
8805 writeColumns(std::ostream& os, mv_type
const &colsA,
size_t const &numCols,
8806 Teuchos::ArrayView<const global_ordinal_type>
const &rowGids,
8807 Teuchos::Array<global_ordinal_type>
const &colsArray,
8811 typedef scalar_type Scalar;
8812 typedef Teuchos::ScalarTraits<Scalar> STS;
8815 const Scalar zero = STS::zero();
8816 const size_t numRows = colsA.getGlobalLength();
8817 for (
size_t j=0; j<numCols; ++j) {
8818 Teuchos::ArrayRCP<const Scalar>
const curCol = colsA.getData(j);
8819 const GO J = colsArray[j];
8820 for (
size_t i=0; i<numRows; ++i) {
8821 const Scalar val = curCol[i];
8823 os << rowGids[i]+indexBase <<
" " << J+indexBase <<
" " << val << std::endl;
8846 const std::string& filename_suffix,
8848 const std::string& matrixName,
8849 const std::string& matrixDescription,
8850 const int ranksToWriteAtOnce=8,
8851 const bool debug=
false) {
8856 using STS =
typename Teuchos::ScalarTraits<ST>;
8861 TEUCHOS_TEST_FOR_EXCEPTION
8862 (comm.is_null (), std::invalid_argument,
8863 "The input matrix's communicator (Teuchos::Comm object) is null.");
8864 TEUCHOS_TEST_FOR_EXCEPTION
8865 (matrix.isGloballyIndexed() || !matrix.isFillComplete(), std::invalid_argument,
8866 "The input matrix must not be GloballyIndexed and must be fillComplete.");
8869 const int myRank = comm->getRank();
8870 const int numProc = comm->getSize();
8871 std::string filename = filename_prefix + std::to_string(myRank) + filename_suffix;
8872 RCP<const map_type> rowMap = matrix.getRowMap();
8873 RCP<const map_type> colMap = matrix.getColMap();
8874 size_t local_nnz = matrix.getLocalNumEntries();
8875 size_t local_num_rows = rowMap->getLocalNumElements();
8876 size_t local_num_cols = colMap->getLocalNumElements();
8877 const GO rowIndexBase = rowMap->getIndexBase();
8878 const GO colIndexBase = colMap->getIndexBase();
8881 int rank_limit = std::min(std::max(ranksToWriteAtOnce,1),numProc);
8884 for(
int base_rank = 0; base_rank < numProc; base_rank += rank_limit) {
8885 int stop = std::min(base_rank+rank_limit,numProc);
8887 if(base_rank <= myRank && myRank < stop) {
8889 std::ofstream out(filename);
8892 out <<
"%%MatrixMarket matrix coordinate "
8893 << (STS::isComplex ?
"complex" :
"real")
8894 <<
" general" << std::endl;
8897 if (matrixName !=
"") {
8898 printAsComment (out, matrixName);
8900 if (matrixDescription !=
"") {
8901 printAsComment (out, matrixDescription);
8907 out << local_num_rows <<
" " << local_num_cols <<
" " << local_nnz <<std::endl;
8914 Teuchos::SetScientific<ST> sci (out);
8916 for(
size_t l_row = 0; l_row < local_num_rows; l_row++) {
8917 GO g_row = rowMap->getGlobalElement(l_row);
8919 typename sparse_matrix_type::local_inds_host_view_type indices;
8920 typename sparse_matrix_type::values_host_view_type values;
8921 matrix.getLocalRowView(l_row, indices, values);
8922 for (
size_t ii = 0; ii < indices.extent(0); ii++) {
8923 const GO g_col = colMap->getGlobalElement(indices(ii));
8926 out << (g_row + 1 - rowIndexBase) <<
" "
8927 << (g_col + 1 - colIndexBase) <<
" ";
8928 if (STS::isComplex) {
8929 out << STS::real(values(ii)) <<
" " << STS::imag(values(ii));
8950 template <
typename T>
8953 return obj.is_null() ? Teuchos::null : obj->getComm();
8959 return comm.is_null() ? 0 : comm->getRank();
8967 #endif // __MatrixMarket_Tpetra_hpp
static std::ifstream openInFileOnRankZero(const trcp_tcomm_t &comm, const std::string &filename, const bool safe=true, std::ios_base::openmode mode=std::ios_base::in)
Open an input file stream safely on rank zero.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Teuchos::RCP< const map_type > getRowMap() const override
Returns the Map that describes the row distribution in this graph.
static Teuchos::RCP< multivector_type > readDense(std::istream &in, const trcp_tcomm_t &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read dense matrix (as a MultiVector) from the given Matrix Market input stream.
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 Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const trcp_tcomm_t &comm, 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.
static Teuchos::RCP< vector_type > readVectorFile(const std::string &filename, const trcp_tcomm_t &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read a Vector from the given Matrix Market file.
static void writeOperator(std::ostream &out, const operator_type &A)
Write a Tpetra::Operator to an output stream.
Teuchos::RCP< const map_type > getColMap() const override
Returns the Map that describes the column distribution in this graph.
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)
Only for backwards compatibility; prefer the overload above.
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.
bool isContiguous() const
True if this Map is distributed contiguously, else false.
void doImport(const SrcDistObject &source, const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, const CombineMode CM, const bool restrictedMode=false)
Import data into this object using an Import object ("forward mode").
size_t getNumVectors() const
Number of columns in the multivector.
Declaration of a function that prints strings from each process.
size_t getLocalNumElements() const
The number of elements belonging to the calling process.
SparseMatrixType::local_ordinal_type local_ordinal_type
SparseMatrixType::global_ordinal_type global_ordinal_type
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.
Teuchos::RCP< const map_type > getDomainMap() const override
Returns the Map associated with the domain of this graph.
SparseMatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the sparse matrix.
static Teuchos::RCP< const map_type > readMapFile(const std::string &filename, const trcp_tcomm_t &comm, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read Map (as a MultiVector) from the given Matrix Market file.
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
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
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.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
Specialization of Tpetra::MultiVector that matches SparseMatrixType.
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const trcp_tcomm_t &comm, 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 trcp_tcomm_t &comm, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read Map (as a MultiVector) from the given input stream.
bool isGloballyIndexed() const override
Whether the graph's column indices are stored as global indices.
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 writeSparse(std::ostream &out, const sparse_matrix_type &matrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
static trcp_tcomm_t getComm(const Teuchos::RCP< T > &obj)
Return obj MPI communicator or Teuchos::null.
static Teuchos::RCP< const map_type > readMap(std::istream &in, const trcp_tcomm_t &comm, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read Map (as a MultiVector) from the given input stream, with optional debugging output stream...
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.
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. ...
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.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const trcp_tcomm_t &comm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
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...
Teuchos::ArrayView< const global_ordinal_type > getLocalElementList() const
Return a NONOWNING view of the global indices owned by this process.
static Teuchos::RCP< sparse_matrix_type > readSparsePerRank(const std::string &filename_prefix, const std::string &filename_suffix, 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 int ranksToReadAtOnce=8, const bool debug=false)
Read a Tpetra::CrsMatrix from a file per rank setup.
size_t global_size_t
Global size_t object.
SparseMatrixType::node_type node_type
The Kokkos Node type; fourth template parameter of Tpetra::CrsMatrix.
static void writeSparsePerRank(const std::string &filename_prefix, const std::string &filename_suffix, const sparse_matrix_type &matrix, const std::string &matrixName, const std::string &matrixDescription, const int ranksToWriteAtOnce=8, const bool debug=false)
Write a Tpetra::CrsMatrix to a file per rank.
static void writeOperator(const std::string &fileName, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to a file, with options.
Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Specialization of Tpetra::Map that matches SparseMatrixType.
static Teuchos::RCP< multivector_type > readDenseFile(const std::string &filename, const trcp_tcomm_t &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read dense matrix (as a MultiVector) from the given Matrix Market file.
Insert new values that don't currently exist.
SparseMatrixType::scalar_type scalar_type
Type of the entries of the sparse matrix.
static void writeMapFile(const std::string &filename, const map_type &map)
Write the Map to the given file.
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...
Teuchos::RCP< const Teuchos::Comm< int >> trcp_tcomm_t
Type of the MPI communicator.
static void writeSparse(std::ostream &out, const sparse_matrix_type &matrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
CrsGraph< local_ordinal_type, global_ordinal_type, node_type > crs_graph_type
Specialization of Tpetra::CrsGraph that matches SparseMatrixType.
static Teuchos::RCP< vector_type > readVector(std::istream &in, const trcp_tcomm_t &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream.
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...
static void writeSparseFile(const std::string &filename, const sparse_matrix_type &matrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
From a distributed map build a map with all GIDs on the root node.
static int getRank(const trcp_tcomm_t &comm)
Return MPI rank or 0.
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 void writeSparseFile(const std::string &filename, const sparse_matrix_type &matrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
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.
Teuchos::RCP< const Teuchos::Comm< int >> trcp_tcomm_t
Type of the MPI communicator.
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.
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...
Teuchos::RCP< const map_type > getRangeMap() const override
Returns the Map associated with the domain of this graph.
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)
Only for backwards compatibility; prefer the overload above.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Only for backwards compatibility; prefer the overload above.
Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getVector(const size_t j) const
Return a Vector which is a const view of column j.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
The MultiVector specialization associated with SparseMatrixType.
void start()
Start the deep_copy counter.
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).
global_ordinal_type getMinGlobalIndex() const
The minimum global index owned by the calling process.
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.
A parallel distribution of indices over processes.
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.
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.
void getLocalRowView(const LocalOrdinal lclRow, local_inds_host_view_type &lclColInds) const override
Get a const view of the given local row's local column indices.
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.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const override
Returns the communicator.
SparseMatrixType sparse_matrix_type
This class' template parameter; a specialization of CrsMatrix.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
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 void writeOperator(std::ostream &out, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to an output stream, with options.
global_size_t getGlobalNumEntries() const override
Returns the global number of entries in the graph.
Matrix Market file reader for CrsMatrix and MultiVector.
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)
Only for backwards compatibility; prefer the overload above.
void getGlobalRowView(const global_ordinal_type gblRow, global_inds_host_view_type &gblColInds) const override
Get a const view of the given global row's global column indices.
void stop()
Stop the deep_copy counter.
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.
global_size_t getGlobalNumElements() const
The number of elements in this Map.
static void writeOperator(const std::string &fileName, operator_type const &A)
Write a Tpetra::Operator to a file.