42 #ifndef __MatrixMarket_Tpetra_hpp
43 #define __MatrixMarket_Tpetra_hpp
56 #include "Tpetra_Details_gathervPrint.hpp"
57 #include "Tpetra_CrsMatrix.hpp"
58 #include "Tpetra_Operator.hpp"
59 #include "Tpetra_Vector.hpp"
61 #include "Teuchos_MatrixMarket_Raw_Adder.hpp"
62 #include "Teuchos_MatrixMarket_Raw_Graph_Adder.hpp"
63 #include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp"
64 #include "Teuchos_MatrixMarket_SymmetrizingGraphAdder.hpp"
65 #include "Teuchos_MatrixMarket_assignScalar.hpp"
66 #include "Teuchos_MatrixMarket_Banner.hpp"
67 #include "Teuchos_MatrixMarket_CoordDataReader.hpp"
68 #include "Teuchos_SetScientific.hpp"
108 namespace MatrixMarket {
164 template<
class SparseMatrixType>
169 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
184 typedef typename SparseMatrixType::global_ordinal_type
206 typedef Teuchos::Comm<int> comm_type;
209 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
211 typedef Teuchos::RCP<const comm_type> comm_ptr TPETRA_DEPRECATED;
212 typedef Teuchos::RCP<const map_type> map_ptr TPETRA_DEPRECATED;
213 typedef Teuchos::RCP<node_type> node_ptr TPETRA_DEPRECATED;
214 #endif // TPETRA_ENABLE_DEPRECATED_CODE
222 typedef Teuchos::ArrayRCP<int>::size_type size_type;
234 static Teuchos::RCP<const map_type>
235 makeRangeMap (
const Teuchos::RCP<const comm_type>& pComm,
239 return rcp (
new map_type (static_cast<global_size_t> (numRows),
240 static_cast<global_ordinal_type> (0),
241 pComm, GloballyDistributed));
271 static Teuchos::RCP<const map_type>
272 makeRowMap (
const Teuchos::RCP<const map_type>& pRowMap,
273 const Teuchos::RCP<const comm_type>& pComm,
278 if (pRowMap.is_null ()) {
279 return rcp (
new map_type (static_cast<global_size_t> (numRows),
280 static_cast<global_ordinal_type> (0),
281 pComm, GloballyDistributed));
284 TEUCHOS_TEST_FOR_EXCEPTION
285 (! pRowMap->isDistributed () && pComm->getSize () > 1,
286 std::invalid_argument,
"The specified row map is not distributed, "
287 "but the given communicator includes more than one process (in "
288 "fact, there are " << pComm->getSize () <<
" processes).");
289 TEUCHOS_TEST_FOR_EXCEPTION
290 (pRowMap->getComm () != pComm, std::invalid_argument,
291 "The specified row Map's communicator (pRowMap->getComm()) "
292 "differs from the given separately supplied communicator pComm.");
311 static Teuchos::RCP<const map_type>
312 makeDomainMap (
const Teuchos::RCP<const map_type>& pRangeMap,
321 if (numRows == numCols) {
324 return createUniformContigMapWithNode<LO,GO,NT> (numCols,
325 pRangeMap->getComm ());
402 distribute (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
403 Teuchos::ArrayRCP<size_t>& myRowPtr,
404 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
405 Teuchos::ArrayRCP<scalar_type>& myValues,
406 const Teuchos::RCP<const map_type>& pRowMap,
407 Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
408 Teuchos::ArrayRCP<size_t>& rowPtr,
409 Teuchos::ArrayRCP<global_ordinal_type>& colInd,
410 Teuchos::ArrayRCP<scalar_type>& values,
411 const bool debug=
false)
414 using Teuchos::ArrayRCP;
415 using Teuchos::ArrayView;
418 using Teuchos::CommRequest;
421 using Teuchos::receive;
426 const bool extraDebug =
false;
427 RCP<const comm_type> pComm = pRowMap->getComm ();
428 const int numProcs = pComm->getSize ();
429 const int myRank = pComm->getRank ();
430 const int rootRank = 0;
437 ArrayView<const GO> myRows = pRowMap->getNodeElementList();
438 const size_type myNumRows = myRows.size();
439 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
440 pRowMap->getNodeNumElements(),
442 "pRowMap->getNodeElementList().size() = "
444 <<
" != pRowMap->getNodeNumElements() = "
445 << pRowMap->getNodeNumElements() <<
". "
446 "Please report this bug to the Tpetra developers.");
447 TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
449 "On Proc 0: numEntriesPerRow.size() = "
450 << numEntriesPerRow.size()
451 <<
" != pRowMap->getNodeElementList().size() = "
452 << myNumRows <<
". Please report this bug to the "
453 "Tpetra developers.");
457 myNumEntriesPerRow = arcp<size_t> (myNumRows);
459 if (myRank != rootRank) {
463 send (*pComm, myNumRows, rootRank);
464 if (myNumRows != 0) {
468 send (*pComm, static_cast<int> (myNumRows),
469 myRows.getRawPtr(), rootRank);
479 receive (*pComm, rootRank,
480 static_cast<int> (myNumRows),
481 myNumEntriesPerRow.getRawPtr());
486 std::accumulate (myNumEntriesPerRow.begin(),
487 myNumEntriesPerRow.end(), 0);
493 myColInd = arcp<GO> (myNumEntries);
494 myValues = arcp<scalar_type> (myNumEntries);
495 if (myNumEntries > 0) {
498 receive (*pComm, rootRank,
499 static_cast<int> (myNumEntries),
500 myColInd.getRawPtr());
501 receive (*pComm, rootRank,
502 static_cast<int> (myNumEntries),
503 myValues.getRawPtr());
509 cerr <<
"-- Proc 0: Copying my data from global arrays" << endl;
513 for (size_type k = 0; k < myNumRows; ++k) {
514 const GO myCurRow = myRows[k];
516 myNumEntriesPerRow[k] = numEntriesInThisRow;
518 if (extraDebug && debug) {
519 cerr <<
"Proc " << pRowMap->getComm ()->getRank ()
520 <<
": myNumEntriesPerRow[0.." << (myNumRows-1) <<
"] = [";
521 for (size_type k = 0; k < myNumRows; ++k) {
522 cerr << myNumEntriesPerRow[k];
523 if (k < myNumRows-1) {
531 std::accumulate (myNumEntriesPerRow.begin(),
532 myNumEntriesPerRow.end(), 0);
534 cerr <<
"-- Proc 0: I own " << myNumRows <<
" rows and "
535 << myNumEntries <<
" entries" << endl;
537 myColInd = arcp<GO> (myNumEntries);
538 myValues = arcp<scalar_type> (myNumEntries);
546 for (size_type k = 0; k < myNumRows;
547 myCurPos += myNumEntriesPerRow[k], ++k) {
549 const GO myRow = myRows[k];
550 const size_t curPos = rowPtr[myRow];
553 if (curNumEntries > 0) {
554 ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
555 ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
556 std::copy (colIndView.begin(), colIndView.end(),
557 myColIndView.begin());
559 ArrayView<scalar_type> valuesView =
560 values (curPos, curNumEntries);
561 ArrayView<scalar_type> myValuesView =
562 myValues (myCurPos, curNumEntries);
563 std::copy (valuesView.begin(), valuesView.end(),
564 myValuesView.begin());
569 for (
int p = 1; p < numProcs; ++p) {
571 cerr <<
"-- Proc 0: Processing proc " << p << endl;
574 size_type theirNumRows = 0;
579 receive (*pComm, p, &theirNumRows);
581 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
582 << theirNumRows <<
" rows" << endl;
584 if (theirNumRows != 0) {
589 ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
590 receive (*pComm, p, as<int> (theirNumRows),
591 theirRows.getRawPtr ());
600 const global_size_t numRows = pRowMap->getGlobalNumElements ();
601 const GO indexBase = pRowMap->getIndexBase ();
602 bool theirRowsValid =
true;
603 for (size_type k = 0; k < theirNumRows; ++k) {
604 if (theirRows[k] < indexBase ||
605 as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
606 theirRowsValid =
false;
609 if (! theirRowsValid) {
610 TEUCHOS_TEST_FOR_EXCEPTION(
611 ! theirRowsValid, std::logic_error,
612 "Proc " << p <<
" has at least one invalid row index. "
613 "Here are all of them: " <<
614 Teuchos::toString (theirRows ()) <<
". Valid row index "
615 "range (zero-based): [0, " << (numRows - 1) <<
"].");
630 ArrayRCP<size_t> theirNumEntriesPerRow;
631 theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
632 for (size_type k = 0; k < theirNumRows; ++k) {
633 theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
640 send (*pComm, static_cast<int> (theirNumRows),
641 theirNumEntriesPerRow.getRawPtr(), p);
645 std::accumulate (theirNumEntriesPerRow.begin(),
646 theirNumEntriesPerRow.end(), 0);
649 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
650 << theirNumEntries <<
" entries" << endl;
655 if (theirNumEntries == 0) {
664 ArrayRCP<GO> theirColInd (theirNumEntries);
665 ArrayRCP<scalar_type> theirValues (theirNumEntries);
673 for (size_type k = 0; k < theirNumRows;
674 theirCurPos += theirNumEntriesPerRow[k], k++) {
676 const GO theirRow = theirRows[k];
682 if (curNumEntries > 0) {
683 ArrayView<GO> colIndView =
684 colInd (curPos, curNumEntries);
685 ArrayView<GO> theirColIndView =
686 theirColInd (theirCurPos, curNumEntries);
687 std::copy (colIndView.begin(), colIndView.end(),
688 theirColIndView.begin());
690 ArrayView<scalar_type> valuesView =
691 values (curPos, curNumEntries);
692 ArrayView<scalar_type> theirValuesView =
693 theirValues (theirCurPos, curNumEntries);
694 std::copy (valuesView.begin(), valuesView.end(),
695 theirValuesView.begin());
702 send (*pComm, static_cast<int> (theirNumEntries),
703 theirColInd.getRawPtr(), p);
704 send (*pComm, static_cast<int> (theirNumEntries),
705 theirValues.getRawPtr(), p);
708 cerr <<
"-- Proc 0: Finished with proc " << p << endl;
716 numEntriesPerRow = null;
721 if (debug && myRank == 0) {
722 cerr <<
"-- Proc 0: About to fill in myRowPtr" << endl;
730 myRowPtr = arcp<size_t> (myNumRows+1);
732 for (size_type k = 1; k < myNumRows+1; ++k) {
733 myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
735 if (extraDebug && debug) {
736 cerr <<
"Proc " << Teuchos::rank (*(pRowMap->getComm()))
737 <<
": myRowPtr[0.." << myNumRows <<
"] = [";
738 for (size_type k = 0; k < myNumRows+1; ++k) {
744 cerr <<
"]" << endl << endl;
747 if (debug && myRank == 0) {
748 cerr <<
"-- Proc 0: Done with distribute" << endl;
765 static Teuchos::RCP<sparse_matrix_type>
766 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
767 Teuchos::ArrayRCP<size_t>& myRowPtr,
768 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
769 Teuchos::ArrayRCP<scalar_type>& myValues,
770 const Teuchos::RCP<const map_type>& pRowMap,
771 const Teuchos::RCP<const map_type>& pRangeMap,
772 const Teuchos::RCP<const map_type>& pDomainMap,
773 const bool callFillComplete =
true)
775 using Teuchos::ArrayView;
789 TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
790 "makeMatrix: myRowPtr array is null. "
791 "Please report this bug to the Tpetra developers.");
792 TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
793 "makeMatrix: domain map is null. "
794 "Please report this bug to the Tpetra developers.");
795 TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
796 "makeMatrix: range map is null. "
797 "Please report this bug to the Tpetra developers.");
798 TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
799 "makeMatrix: row map is null. "
800 "Please report this bug to the Tpetra developers.");
804 RCP<sparse_matrix_type> A =
810 ArrayView<const GO> myRows = pRowMap->getNodeElementList ();
811 const size_type myNumRows = myRows.size ();
814 const GO indexBase = pRowMap->getIndexBase ();
815 for (size_type i = 0; i < myNumRows; ++i) {
816 const size_type myCurPos = myRowPtr[i];
818 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
819 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
822 for (size_type k = 0; k < curNumEntries; ++k) {
823 curColInd[k] += indexBase;
826 if (curNumEntries > 0) {
827 A->insertGlobalValues (myRows[i], curColInd, curValues);
834 myNumEntriesPerRow = null;
839 if (callFillComplete) {
840 A->fillComplete (pDomainMap, pRangeMap);
850 static Teuchos::RCP<sparse_matrix_type>
851 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
852 Teuchos::ArrayRCP<size_t>& myRowPtr,
853 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
854 Teuchos::ArrayRCP<scalar_type>& myValues,
855 const Teuchos::RCP<const map_type>& pRowMap,
856 const Teuchos::RCP<const map_type>& pRangeMap,
857 const Teuchos::RCP<const map_type>& pDomainMap,
858 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
859 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams)
861 using Teuchos::ArrayView;
875 TEUCHOS_TEST_FOR_EXCEPTION(
876 myRowPtr.is_null(), std::logic_error,
877 "makeMatrix: myRowPtr array is null. "
878 "Please report this bug to the Tpetra developers.");
879 TEUCHOS_TEST_FOR_EXCEPTION(
880 pDomainMap.is_null(), std::logic_error,
881 "makeMatrix: domain map is null. "
882 "Please report this bug to the Tpetra developers.");
883 TEUCHOS_TEST_FOR_EXCEPTION(
884 pRangeMap.is_null(), std::logic_error,
885 "makeMatrix: range map is null. "
886 "Please report this bug to the Tpetra developers.");
887 TEUCHOS_TEST_FOR_EXCEPTION(
888 pRowMap.is_null(), std::logic_error,
889 "makeMatrix: row map is null. "
890 "Please report this bug to the Tpetra developers.");
894 RCP<sparse_matrix_type> A =
896 StaticProfile, constructorParams));
900 ArrayView<const GO> myRows = pRowMap->getNodeElementList();
901 const size_type myNumRows = myRows.size();
904 const GO indexBase = pRowMap->getIndexBase ();
905 for (size_type i = 0; i < myNumRows; ++i) {
906 const size_type myCurPos = myRowPtr[i];
908 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
909 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
912 for (size_type k = 0; k < curNumEntries; ++k) {
913 curColInd[k] += indexBase;
915 if (curNumEntries > 0) {
916 A->insertGlobalValues (myRows[i], curColInd, curValues);
923 myNumEntriesPerRow = null;
928 A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
936 static Teuchos::RCP<sparse_matrix_type>
937 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
938 Teuchos::ArrayRCP<size_t>& myRowPtr,
939 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
940 Teuchos::ArrayRCP<scalar_type>& myValues,
941 const Teuchos::RCP<const map_type>& rowMap,
942 Teuchos::RCP<const map_type>& colMap,
943 const Teuchos::RCP<const map_type>& domainMap,
944 const Teuchos::RCP<const map_type>& rangeMap,
945 const bool callFillComplete =
true)
947 using Teuchos::ArrayView;
956 RCP<sparse_matrix_type> A;
957 if (colMap.is_null ()) {
965 ArrayView<const GO> myRows = rowMap->getNodeElementList ();
966 const size_type myNumRows = myRows.size ();
969 const GO indexBase = rowMap->getIndexBase ();
970 for (size_type i = 0; i < myNumRows; ++i) {
971 const size_type myCurPos = myRowPtr[i];
972 const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
973 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
974 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
977 for (size_type k = 0; k < curNumEntries; ++k) {
978 curColInd[k] += indexBase;
980 if (curNumEntries > 0) {
981 A->insertGlobalValues (myRows[i], curColInd, curValues);
988 myNumEntriesPerRow = null;
993 if (callFillComplete) {
994 A->fillComplete (domainMap, rangeMap);
995 if (colMap.is_null ()) {
996 colMap = A->getColMap ();
1020 static Teuchos::RCP<const Teuchos::MatrixMarket::Banner>
1021 readBanner (std::istream& in,
1023 const bool tolerant=
false,
1025 const bool isGraph=
false)
1027 using Teuchos::MatrixMarket::Banner;
1032 typedef Teuchos::ScalarTraits<scalar_type> STS;
1034 RCP<Banner> pBanner;
1038 const bool readFailed = ! getline(in, line);
1039 TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1040 "Failed to get Matrix Market banner line from input.");
1047 pBanner = rcp (
new Banner (line, tolerant));
1048 }
catch (std::exception& e) {
1049 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1050 "Matrix Market banner line contains syntax error(s): "
1053 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() !=
"matrix",
1054 std::invalid_argument,
"The Matrix Market file does not contain "
1055 "matrix data. Its Banner (first) line says that its object type is \""
1056 << pBanner->matrixType() <<
"\", rather than the required \"matrix\".");
1060 TEUCHOS_TEST_FOR_EXCEPTION(
1061 ! STS::isComplex && pBanner->dataType() ==
"complex",
1062 std::invalid_argument,
1063 "The Matrix Market file contains complex-valued data, but you are "
1064 "trying to read it into a matrix containing entries of the real-"
1065 "valued Scalar type \""
1066 << Teuchos::TypeNameTraits<scalar_type>::name() <<
"\".");
1067 TEUCHOS_TEST_FOR_EXCEPTION(
1069 pBanner->dataType() !=
"real" &&
1070 pBanner->dataType() !=
"complex" &&
1071 pBanner->dataType() !=
"integer",
1072 std::invalid_argument,
1073 "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
1074 "Matrix Market file may not contain a \"pattern\" matrix. A "
1075 "pattern matrix is really just a graph with no weights. It "
1076 "should be stored in a CrsGraph, not a CrsMatrix.");
1078 TEUCHOS_TEST_FOR_EXCEPTION(
1080 pBanner->dataType() !=
"pattern",
1081 std::invalid_argument,
1082 "When reading Matrix Market data into a Tpetra::CrsGraph, the "
1083 "Matrix Market file must contain a \"pattern\" matrix.");
1110 static Teuchos::Tuple<global_ordinal_type, 3>
1111 readCoordDims (std::istream& in,
1113 const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1114 const Teuchos::RCP<const comm_type>& pComm,
1115 const bool tolerant =
false,
1118 using Teuchos::MatrixMarket::readCoordinateDimensions;
1119 using Teuchos::Tuple;
1124 Tuple<global_ordinal_type, 3> dims;
1130 bool success =
false;
1131 if (pComm->getRank() == 0) {
1132 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() !=
"coordinate",
1133 std::invalid_argument,
"The Tpetra::CrsMatrix Matrix Market reader "
1134 "only accepts \"coordinate\" (sparse) matrix data.");
1138 success = readCoordinateDimensions (in, numRows, numCols,
1139 numNonzeros, lineNumber,
1145 dims[2] = numNonzeros;
1153 int the_success = success ? 1 : 0;
1154 Teuchos::broadcast (*pComm, 0, &the_success);
1155 success = (the_success == 1);
1160 Teuchos::broadcast (*pComm, 0, dims);
1168 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1169 "Error reading Matrix Market sparse matrix: failed to read "
1170 "coordinate matrix dimensions.");
1185 typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1187 typedef Teuchos::MatrixMarket::SymmetrizingGraphAdder<Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> > graph_adder_type;
1214 static Teuchos::RCP<adder_type>
1215 makeAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1216 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1217 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1218 const bool tolerant=
false,
1219 const bool debug=
false)
1221 if (pComm->getRank () == 0) {
1222 typedef Teuchos::MatrixMarket::Raw::Adder<
scalar_type,
1225 Teuchos::RCP<raw_adder_type> pRaw =
1226 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1228 return Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
1231 return Teuchos::null;
1260 static Teuchos::RCP<graph_adder_type>
1261 makeGraphAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1262 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1263 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1264 const bool tolerant=
false,
1265 const bool debug=
false)
1267 if (pComm->getRank () == 0) {
1268 typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
1269 Teuchos::RCP<raw_adder_type> pRaw =
1270 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1272 return Teuchos::rcp (
new graph_adder_type (pRaw, pBanner->symmType ()));
1275 return Teuchos::null;
1280 static Teuchos::RCP<sparse_graph_type>
1281 readSparseGraphHelper (std::istream& in,
1282 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1283 const Teuchos::RCP<const map_type>& rowMap,
1284 Teuchos::RCP<const map_type>& colMap,
1285 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1286 const bool tolerant,
1289 using Teuchos::MatrixMarket::Banner;
1292 using Teuchos::Tuple;
1296 const int myRank = pComm->getRank ();
1297 const int rootRank = 0;
1302 size_t lineNumber = 1;
1304 if (debug && myRank == rootRank) {
1305 cerr <<
"Matrix Market reader: readGraph:" << endl
1306 <<
"-- Reading banner line" << endl;
1315 RCP<const Banner> pBanner;
1321 int bannerIsCorrect = 1;
1322 std::ostringstream errMsg;
1324 if (myRank == rootRank) {
1327 pBanner = readBanner (in, lineNumber, tolerant, debug,
true);
1329 catch (std::exception& e) {
1330 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
1331 "threw an exception: " << e.what();
1332 bannerIsCorrect = 0;
1335 if (bannerIsCorrect) {
1340 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
1341 bannerIsCorrect = 0;
1342 errMsg <<
"The Matrix Market input file must contain a "
1343 "\"coordinate\"-format sparse graph in order to create a "
1344 "Tpetra::CrsGraph object from it, but the file's matrix "
1345 "type is \"" << pBanner->matrixType() <<
"\" instead.";
1350 if (tolerant && pBanner->matrixType() ==
"array") {
1351 bannerIsCorrect = 0;
1352 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
1353 "format sparse graph in order to create a Tpetra::CrsGraph "
1354 "object from it, but the file's matrix type is \"array\" "
1355 "instead. That probably means the file contains dense matrix "
1362 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1369 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1370 std::invalid_argument, errMsg.str ());
1372 if (debug && myRank == rootRank) {
1373 cerr <<
"-- Reading dimensions line" << endl;
1381 Tuple<global_ordinal_type, 3> dims =
1382 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1384 if (debug && myRank == rootRank) {
1385 cerr <<
"-- Making Adder for collecting graph data" << endl;
1392 RCP<graph_adder_type> pAdder =
1393 makeGraphAdder (pComm, pBanner, dims, tolerant, debug);
1395 if (debug && myRank == rootRank) {
1396 cerr <<
"-- Reading graph data" << endl;
1406 int readSuccess = 1;
1407 std::ostringstream errMsg;
1408 if (myRank == rootRank) {
1411 typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
1413 reader_type reader (pAdder);
1416 std::pair<bool, std::vector<size_t> > results =
1417 reader.read (in, lineNumber, tolerant, debug);
1418 readSuccess = results.first ? 1 : 0;
1420 catch (std::exception& e) {
1425 broadcast (*pComm, rootRank, ptr (&readSuccess));
1434 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1435 "Failed to read the Matrix Market sparse graph file: "
1439 if (debug && myRank == rootRank) {
1440 cerr <<
"-- Successfully read the Matrix Market data" << endl;
1453 if (debug && myRank == rootRank) {
1454 cerr <<
"-- Tolerant mode: rebroadcasting graph dimensions"
1456 <<
"----- Dimensions before: "
1457 << dims[0] <<
" x " << dims[1]
1461 Tuple<global_ordinal_type, 2> updatedDims;
1462 if (myRank == rootRank) {
1469 std::max (dims[0], pAdder->getAdder()->numRows());
1470 updatedDims[1] = pAdder->getAdder()->numCols();
1472 broadcast (*pComm, rootRank, updatedDims);
1473 dims[0] = updatedDims[0];
1474 dims[1] = updatedDims[1];
1475 if (debug && myRank == rootRank) {
1476 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
1490 if (myRank == rootRank) {
1497 if (dims[0] < pAdder->getAdder ()->numRows ()) {
1501 broadcast (*pComm, 0, ptr (&dimsMatch));
1502 if (dimsMatch == 0) {
1509 Tuple<global_ordinal_type, 2> addersDims;
1510 if (myRank == rootRank) {
1511 addersDims[0] = pAdder->getAdder()->numRows();
1512 addersDims[1] = pAdder->getAdder()->numCols();
1514 broadcast (*pComm, 0, addersDims);
1515 TEUCHOS_TEST_FOR_EXCEPTION(
1516 dimsMatch == 0, std::runtime_error,
1517 "The graph metadata says that the graph is " << dims[0] <<
" x "
1518 << dims[1] <<
", but the actual data says that the graph is "
1519 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
1520 "data includes more rows than reported in the metadata. This "
1521 "is not allowed when parsing in strict mode. Parse the graph in "
1522 "tolerant mode to ignore the metadata when it disagrees with the "
1528 RCP<map_type> proc0Map;
1530 if(Teuchos::is_null(rowMap)) {
1534 indexBase = rowMap->getIndexBase();
1536 if(myRank == rootRank) {
1537 proc0Map = rcp(
new map_type(dims[0],dims[0],indexBase,pComm));
1540 proc0Map = rcp(
new map_type(dims[0],0,indexBase,pComm));
1544 std::map<global_ordinal_type, size_t> numEntriesPerRow_map;
1545 if (myRank == rootRank) {
1546 const auto& entries = pAdder()->getAdder()->getEntries();
1551 for (
const auto& entry : entries) {
1553 ++numEntriesPerRow_map[gblRow];
1557 Teuchos::Array<size_t> numEntriesPerRow (proc0Map->getNodeNumElements ());
1558 for (
const auto& ent : numEntriesPerRow_map) {
1560 numEntriesPerRow[lclRow] = ent.second;
1567 std::map<global_ordinal_type, size_t> empty_map;
1568 std::swap (numEntriesPerRow_map, empty_map);
1571 RCP<sparse_graph_type> proc0Graph =
1573 StaticProfile,constructorParams));
1574 if(myRank == rootRank) {
1575 typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
1578 const std::vector<element_type>& entries =
1579 pAdder->getAdder()->getEntries();
1582 for(
size_t curPos=0; curPos<entries.size(); curPos++) {
1583 const element_type& curEntry = entries[curPos];
1586 Teuchos::ArrayView<const global_ordinal_type> colView(&curCol,1);
1587 proc0Graph->insertGlobalIndices(curRow,colView);
1590 proc0Graph->fillComplete();
1592 RCP<sparse_graph_type> distGraph;
1593 if(Teuchos::is_null(rowMap))
1596 RCP<map_type> distMap =
1597 rcp(
new map_type(dims[0],0,pComm,GloballyDistributed));
1600 distGraph = rcp(
new sparse_graph_type(distMap,colMap,0,StaticProfile,constructorParams));
1603 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1604 import_type importer (proc0Map, distMap);
1607 distGraph->doImport(*proc0Graph,importer,
INSERT);
1610 distGraph = rcp(
new sparse_graph_type(rowMap,colMap,0,StaticProfile,constructorParams));
1613 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1614 import_type importer (proc0Map, rowMap);
1617 distGraph->doImport(*proc0Graph,importer,
INSERT);
1647 static Teuchos::RCP<sparse_graph_type>
1649 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1650 const bool callFillComplete=
true,
1651 const bool tolerant=
false,
1652 const bool debug=
false)
1654 using Teuchos::broadcast;
1655 using Teuchos::outArg;
1663 if (comm->getRank () == 0) {
1665 in.open (filename.c_str ());
1672 broadcast<int, int> (*comm, 0, outArg (opened));
1673 TEUCHOS_TEST_FOR_EXCEPTION
1674 (opened == 0, std::runtime_error,
"readSparseGraphFile: "
1675 "Failed to open file \"" << filename <<
"\" on Process 0.");
1684 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
1685 static Teuchos::RCP<sparse_graph_type> TPETRA_DEPRECATED
1693 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1694 const Teuchos::RCP<node_type>& ,
1695 const bool callFillComplete=
true,
1696 const bool tolerant=
false,
1697 const bool debug=
false)
1701 callFillComplete, tolerant, debug);
1703 #endif // TPETRA_ENABLE_DEPRECATED_CODE
1733 static Teuchos::RCP<sparse_graph_type>
1735 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1736 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1737 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1738 const bool tolerant=
false,
1739 const bool debug=
false)
1741 using Teuchos::broadcast;
1742 using Teuchos::outArg;
1750 if (pComm->getRank () == 0) {
1752 in.open (filename.c_str ());
1759 broadcast<int, int> (*pComm, 0, outArg (opened));
1760 TEUCHOS_TEST_FOR_EXCEPTION
1761 (opened == 0, std::runtime_error,
"readSparseGraphFile: "
1762 "Failed to open file \"" << filename <<
"\" on Process 0.");
1763 if (pComm->getRank () == 0) {
1764 in.open (filename.c_str ());
1768 fillCompleteParams, tolerant, debug);
1774 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
1775 static Teuchos::RCP<sparse_graph_type> TPETRA_DEPRECATED
1783 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1784 const Teuchos::RCP<node_type>& ,
1785 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1786 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1787 const bool tolerant=
false,
1788 const bool debug=
false)
1792 constructorParams, fillCompleteParams,
1795 #endif // TPETRA_ENABLE_DEPRECATED_CODE
1834 static Teuchos::RCP<sparse_graph_type>
1836 const Teuchos::RCP<const map_type>& rowMap,
1837 Teuchos::RCP<const map_type>& colMap,
1838 const Teuchos::RCP<const map_type>& domainMap,
1839 const Teuchos::RCP<const map_type>& rangeMap,
1840 const bool callFillComplete=
true,
1841 const bool tolerant=
false,
1842 const bool debug=
false)
1844 using Teuchos::broadcast;
1845 using Teuchos::Comm;
1846 using Teuchos::outArg;
1849 TEUCHOS_TEST_FOR_EXCEPTION
1850 (rowMap.is_null (), std::invalid_argument,
1851 "Input rowMap must be nonnull.");
1852 RCP<const Comm<int> > comm = rowMap->getComm ();
1853 if (comm.is_null ()) {
1856 return Teuchos::null;
1865 if (comm->getRank () == 0) {
1867 in.open (filename.c_str ());
1874 broadcast<int, int> (*comm, 0, outArg (opened));
1875 TEUCHOS_TEST_FOR_EXCEPTION
1876 (opened == 0, std::runtime_error,
"readSparseGraphFile: "
1877 "Failed to open file \"" << filename <<
"\" on Process 0.");
1879 callFillComplete, tolerant, debug);
1907 static Teuchos::RCP<sparse_graph_type>
1909 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1910 const bool callFillComplete=
true,
1911 const bool tolerant=
false,
1912 const bool debug=
false)
1914 Teuchos::RCP<const map_type> fakeRowMap;
1915 Teuchos::RCP<const map_type> fakeColMap;
1916 Teuchos::RCP<Teuchos::ParameterList> fakeCtorParams;
1918 Teuchos::RCP<sparse_graph_type> graph =
1919 readSparseGraphHelper (in, pComm,
1920 fakeRowMap, fakeColMap,
1921 fakeCtorParams, tolerant, debug);
1922 if (callFillComplete) {
1923 graph->fillComplete ();
1928 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
1929 static Teuchos::RCP<sparse_graph_type> TPETRA_DEPRECATED
1932 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1933 const Teuchos::RCP<node_type>& ,
1934 const bool callFillComplete=
true,
1935 const bool tolerant=
false,
1936 const bool debug=
false)
1942 #endif // TPETRA_ENABLE_DEPRECATED_CODE
1973 static Teuchos::RCP<sparse_graph_type>
1975 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1976 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1977 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1978 const bool tolerant=
false,
1979 const bool debug=
false)
1981 Teuchos::RCP<const map_type> fakeRowMap;
1982 Teuchos::RCP<const map_type> fakeColMap;
1983 Teuchos::RCP<sparse_graph_type> graph =
1984 readSparseGraphHelper (in, pComm,
1985 fakeRowMap, fakeColMap,
1986 constructorParams, tolerant, debug);
1987 graph->fillComplete (fillCompleteParams);
1991 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
1992 static Teuchos::RCP<sparse_graph_type> TPETRA_DEPRECATED
1995 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1996 const Teuchos::RCP<node_type>& ,
1997 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1998 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1999 const bool tolerant=
false,
2000 const bool debug=
false)
2004 fillCompleteParams, tolerant, debug);
2006 #endif // TPETRA_ENABLE_DEPRECATED_CODE
2048 static Teuchos::RCP<sparse_graph_type>
2050 const Teuchos::RCP<const map_type>& rowMap,
2051 Teuchos::RCP<const map_type>& colMap,
2052 const Teuchos::RCP<const map_type>& domainMap,
2053 const Teuchos::RCP<const map_type>& rangeMap,
2054 const bool callFillComplete=
true,
2055 const bool tolerant=
false,
2056 const bool debug=
false)
2058 Teuchos::RCP<sparse_graph_type> graph =
2059 readSparseGraphHelper (in, rowMap->getComm (),
2060 rowMap, colMap, Teuchos::null, tolerant,
2062 if (callFillComplete) {
2063 graph->fillComplete (domainMap, rangeMap);
2091 static Teuchos::RCP<sparse_matrix_type>
2093 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2094 const bool callFillComplete=
true,
2095 const bool tolerant=
false,
2096 const bool debug=
false)
2098 const int myRank = pComm->getRank ();
2103 in.open (filename.c_str ());
2111 return readSparse (in, pComm, callFillComplete, tolerant, debug);
2117 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
2118 static Teuchos::RCP<sparse_matrix_type> TPETRA_DEPRECATED
2121 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2122 const Teuchos::RCP<node_type>& ,
2123 const bool callFillComplete=
true,
2124 const bool tolerant=
false,
2125 const bool debug=
false)
2127 return readSparseFile (filename, pComm, callFillComplete, tolerant, debug);
2129 #endif // TPETRA_ENABLE_DEPRECATED_CODE
2159 static Teuchos::RCP<sparse_matrix_type>
2161 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2162 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2163 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2164 const bool tolerant=
false,
2165 const bool debug=
false)
2168 if (pComm->getRank () == 0) {
2169 in.open (filename.c_str ());
2171 return readSparse (in, pComm, constructorParams,
2172 fillCompleteParams, tolerant, debug);
2175 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
2176 static Teuchos::RCP<sparse_matrix_type> TPETRA_DEPRECATED
2179 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2180 const Teuchos::RCP<node_type>& ,
2181 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2182 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2183 const bool tolerant=
false,
2184 const bool debug=
false)
2187 constructorParams, fillCompleteParams,
2190 #endif // TPETRA_ENABLE_DEPRECATED_CODE
2229 static Teuchos::RCP<sparse_matrix_type>
2231 const Teuchos::RCP<const map_type>& rowMap,
2232 Teuchos::RCP<const map_type>& colMap,
2233 const Teuchos::RCP<const map_type>& domainMap,
2234 const Teuchos::RCP<const map_type>& rangeMap,
2235 const bool callFillComplete=
true,
2236 const bool tolerant=
false,
2237 const bool debug=
false)
2239 using Teuchos::broadcast;
2240 using Teuchos::Comm;
2241 using Teuchos::outArg;
2244 TEUCHOS_TEST_FOR_EXCEPTION(
2245 rowMap.is_null (), std::invalid_argument,
2246 "Row Map must be nonnull.");
2248 RCP<const Comm<int> > comm = rowMap->getComm ();
2249 const int myRank = comm->getRank ();
2259 in.open (filename.c_str ());
2266 broadcast<int, int> (*comm, 0, outArg (opened));
2267 TEUCHOS_TEST_FOR_EXCEPTION(
2268 opened == 0, std::runtime_error,
2269 "readSparseFile: Failed to open file \"" << filename <<
"\" on "
2271 return readSparse (in, rowMap, colMap, domainMap, rangeMap,
2272 callFillComplete, tolerant, debug);
2300 static Teuchos::RCP<sparse_matrix_type>
2302 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2303 const bool callFillComplete=
true,
2304 const bool tolerant=
false,
2305 const bool debug=
false)
2307 using Teuchos::MatrixMarket::Banner;
2308 using Teuchos::arcp;
2309 using Teuchos::ArrayRCP;
2310 using Teuchos::broadcast;
2311 using Teuchos::null;
2314 using Teuchos::REDUCE_MAX;
2315 using Teuchos::reduceAll;
2316 using Teuchos::Tuple;
2319 typedef Teuchos::ScalarTraits<scalar_type> STS;
2321 const bool extraDebug =
false;
2322 const int myRank = pComm->getRank ();
2323 const int rootRank = 0;
2328 size_t lineNumber = 1;
2330 if (debug && myRank == rootRank) {
2331 cerr <<
"Matrix Market reader: readSparse:" << endl
2332 <<
"-- Reading banner line" << endl;
2341 RCP<const Banner> pBanner;
2347 int bannerIsCorrect = 1;
2348 std::ostringstream errMsg;
2350 if (myRank == rootRank) {
2353 pBanner = readBanner (in, lineNumber, tolerant, debug);
2355 catch (std::exception& e) {
2356 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2357 "threw an exception: " << e.what();
2358 bannerIsCorrect = 0;
2361 if (bannerIsCorrect) {
2366 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2367 bannerIsCorrect = 0;
2368 errMsg <<
"The Matrix Market input file must contain a "
2369 "\"coordinate\"-format sparse matrix in order to create a "
2370 "Tpetra::CrsMatrix object from it, but the file's matrix "
2371 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2376 if (tolerant && pBanner->matrixType() ==
"array") {
2377 bannerIsCorrect = 0;
2378 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2379 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2380 "object from it, but the file's matrix type is \"array\" "
2381 "instead. That probably means the file contains dense matrix "
2388 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2395 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2396 std::invalid_argument, errMsg.str ());
2398 if (debug && myRank == rootRank) {
2399 cerr <<
"-- Reading dimensions line" << endl;
2407 Tuple<global_ordinal_type, 3> dims =
2408 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2410 if (debug && myRank == rootRank) {
2411 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2416 RCP<adder_type> pAdder =
2417 makeAdder (pComm, pBanner, dims, tolerant, debug);
2419 if (debug && myRank == rootRank) {
2420 cerr <<
"-- Reading matrix data" << endl;
2430 int readSuccess = 1;
2431 std::ostringstream errMsg;
2432 if (myRank == rootRank) {
2435 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2437 reader_type reader (pAdder);
2440 std::pair<bool, std::vector<size_t> > results =
2441 reader.read (in, lineNumber, tolerant, debug);
2442 readSuccess = results.first ? 1 : 0;
2444 catch (std::exception& e) {
2449 broadcast (*pComm, rootRank, ptr (&readSuccess));
2458 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2459 "Failed to read the Matrix Market sparse matrix file: "
2463 if (debug && myRank == rootRank) {
2464 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2477 if (debug && myRank == rootRank) {
2478 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2480 <<
"----- Dimensions before: "
2481 << dims[0] <<
" x " << dims[1]
2485 Tuple<global_ordinal_type, 2> updatedDims;
2486 if (myRank == rootRank) {
2493 std::max (dims[0], pAdder->getAdder()->numRows());
2494 updatedDims[1] = pAdder->getAdder()->numCols();
2496 broadcast (*pComm, rootRank, updatedDims);
2497 dims[0] = updatedDims[0];
2498 dims[1] = updatedDims[1];
2499 if (debug && myRank == rootRank) {
2500 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
2514 if (myRank == rootRank) {
2521 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2525 broadcast (*pComm, 0, ptr (&dimsMatch));
2526 if (dimsMatch == 0) {
2533 Tuple<global_ordinal_type, 2> addersDims;
2534 if (myRank == rootRank) {
2535 addersDims[0] = pAdder->getAdder()->numRows();
2536 addersDims[1] = pAdder->getAdder()->numCols();
2538 broadcast (*pComm, 0, addersDims);
2539 TEUCHOS_TEST_FOR_EXCEPTION(
2540 dimsMatch == 0, std::runtime_error,
2541 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
2542 << dims[1] <<
", but the actual data says that the matrix is "
2543 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
2544 "data includes more rows than reported in the metadata. This "
2545 "is not allowed when parsing in strict mode. Parse the matrix in "
2546 "tolerant mode to ignore the metadata when it disagrees with the "
2551 if (debug && myRank == rootRank) {
2552 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
2564 ArrayRCP<size_t> numEntriesPerRow;
2565 ArrayRCP<size_t> rowPtr;
2566 ArrayRCP<global_ordinal_type> colInd;
2567 ArrayRCP<scalar_type> values;
2572 int mergeAndConvertSucceeded = 1;
2573 std::ostringstream errMsg;
2575 if (myRank == rootRank) {
2577 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2587 const size_type numRows = dims[0];
2590 pAdder->getAdder()->merge ();
2593 const std::vector<element_type>& entries =
2594 pAdder->getAdder()->getEntries();
2597 const size_t numEntries = (size_t)entries.size();
2600 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
2601 <<
" rows and numEntries=" << numEntries
2602 <<
" entries." << endl;
2608 numEntriesPerRow = arcp<size_t> (numRows);
2609 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2610 rowPtr = arcp<size_t> (numRows+1);
2611 std::fill (rowPtr.begin(), rowPtr.end(), 0);
2612 colInd = arcp<global_ordinal_type> (numEntries);
2613 values = arcp<scalar_type> (numEntries);
2620 for (curPos = 0; curPos < numEntries; ++curPos) {
2621 const element_type& curEntry = entries[curPos];
2623 TEUCHOS_TEST_FOR_EXCEPTION(
2624 curRow < prvRow, std::logic_error,
2625 "Row indices are out of order, even though they are supposed "
2626 "to be sorted. curRow = " << curRow <<
", prvRow = "
2627 << prvRow <<
", at curPos = " << curPos <<
". Please report "
2628 "this bug to the Tpetra developers.");
2629 if (curRow > prvRow) {
2635 numEntriesPerRow[curRow]++;
2636 colInd[curPos] = curEntry.colIndex();
2637 values[curPos] = curEntry.value();
2642 rowPtr[numRows] = numEntries;
2644 catch (std::exception& e) {
2645 mergeAndConvertSucceeded = 0;
2646 errMsg <<
"Failed to merge sparse matrix entries and convert to "
2647 "CSR format: " << e.what();
2650 if (debug && mergeAndConvertSucceeded) {
2652 const size_type numRows = dims[0];
2653 const size_type maxToDisplay = 100;
2655 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
2656 << (numEntriesPerRow.size()-1) <<
"] ";
2657 if (numRows > maxToDisplay) {
2658 cerr <<
"(only showing first and last few entries) ";
2662 if (numRows > maxToDisplay) {
2663 for (size_type k = 0; k < 2; ++k) {
2664 cerr << numEntriesPerRow[k] <<
" ";
2667 for (size_type k = numRows-2; k < numRows-1; ++k) {
2668 cerr << numEntriesPerRow[k] <<
" ";
2672 for (size_type k = 0; k < numRows-1; ++k) {
2673 cerr << numEntriesPerRow[k] <<
" ";
2676 cerr << numEntriesPerRow[numRows-1];
2678 cerr <<
"]" << endl;
2680 cerr <<
"----- Proc 0: rowPtr ";
2681 if (numRows > maxToDisplay) {
2682 cerr <<
"(only showing first and last few entries) ";
2685 if (numRows > maxToDisplay) {
2686 for (size_type k = 0; k < 2; ++k) {
2687 cerr << rowPtr[k] <<
" ";
2690 for (size_type k = numRows-2; k < numRows; ++k) {
2691 cerr << rowPtr[k] <<
" ";
2695 for (size_type k = 0; k < numRows; ++k) {
2696 cerr << rowPtr[k] <<
" ";
2699 cerr << rowPtr[numRows] <<
"]" << endl;
2710 if (debug && myRank == rootRank) {
2711 cerr <<
"-- Making range, domain, and row maps" << endl;
2718 RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
2719 RCP<const map_type> pDomainMap =
2720 makeDomainMap (pRangeMap, dims[0], dims[1]);
2721 RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
2723 if (debug && myRank == rootRank) {
2724 cerr <<
"-- Distributing the matrix data" << endl;
2738 ArrayRCP<size_t> myNumEntriesPerRow;
2739 ArrayRCP<size_t> myRowPtr;
2740 ArrayRCP<global_ordinal_type> myColInd;
2741 ArrayRCP<scalar_type> myValues;
2743 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
2744 numEntriesPerRow, rowPtr, colInd, values, debug);
2746 if (debug && myRank == rootRank) {
2747 cerr <<
"-- Inserting matrix entries on each processor";
2748 if (callFillComplete) {
2749 cerr <<
" and calling fillComplete()";
2760 RCP<sparse_matrix_type> pMatrix =
2761 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
2762 pRowMap, pRangeMap, pDomainMap, callFillComplete);
2767 int localIsNull = pMatrix.is_null () ? 1 : 0;
2768 int globalIsNull = 0;
2769 reduceAll (*pComm, REDUCE_MAX, localIsNull, ptr (&globalIsNull));
2770 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
2771 "Reader::makeMatrix() returned a null pointer on at least one "
2772 "process. Please report this bug to the Tpetra developers.");
2775 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
2776 "Reader::makeMatrix() returned a null pointer. "
2777 "Please report this bug to the Tpetra developers.");
2789 if (callFillComplete) {
2790 const int numProcs = pComm->getSize ();
2792 if (extraDebug && debug) {
2794 pRangeMap->getGlobalNumElements ();
2796 pDomainMap->getGlobalNumElements ();
2797 if (myRank == rootRank) {
2798 cerr <<
"-- Matrix is "
2799 << globalNumRows <<
" x " << globalNumCols
2800 <<
" with " << pMatrix->getGlobalNumEntries()
2801 <<
" entries, and index base "
2802 << pMatrix->getIndexBase() <<
"." << endl;
2805 for (
int p = 0; p < numProcs; ++p) {
2807 cerr <<
"-- Proc " << p <<
" owns "
2808 << pMatrix->getNodeNumCols() <<
" columns, and "
2809 << pMatrix->getNodeNumEntries() <<
" entries." << endl;
2816 if (debug && myRank == rootRank) {
2817 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
2823 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
2824 static Teuchos::RCP<sparse_matrix_type> TPETRA_DEPRECATED
2827 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2828 const Teuchos::RCP<node_type>& ,
2829 const bool callFillComplete=
true,
2830 const bool tolerant=
false,
2831 const bool debug=
false)
2833 return readSparse (in, pComm, callFillComplete, tolerant, debug);
2835 #endif // TPETRA_ENABLE_DEPRECATED_CODE
2865 static Teuchos::RCP<sparse_matrix_type>
2867 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2868 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2869 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2870 const bool tolerant=
false,
2871 const bool debug=
false)
2873 using Teuchos::MatrixMarket::Banner;
2874 using Teuchos::arcp;
2875 using Teuchos::ArrayRCP;
2876 using Teuchos::broadcast;
2877 using Teuchos::null;
2880 using Teuchos::reduceAll;
2881 using Teuchos::Tuple;
2884 typedef Teuchos::ScalarTraits<scalar_type> STS;
2886 const bool extraDebug =
false;
2887 const int myRank = pComm->getRank ();
2888 const int rootRank = 0;
2893 size_t lineNumber = 1;
2895 if (debug && myRank == rootRank) {
2896 cerr <<
"Matrix Market reader: readSparse:" << endl
2897 <<
"-- Reading banner line" << endl;
2906 RCP<const Banner> pBanner;
2912 int bannerIsCorrect = 1;
2913 std::ostringstream errMsg;
2915 if (myRank == rootRank) {
2918 pBanner = readBanner (in, lineNumber, tolerant, debug);
2920 catch (std::exception& e) {
2921 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2922 "threw an exception: " << e.what();
2923 bannerIsCorrect = 0;
2926 if (bannerIsCorrect) {
2931 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2932 bannerIsCorrect = 0;
2933 errMsg <<
"The Matrix Market input file must contain a "
2934 "\"coordinate\"-format sparse matrix in order to create a "
2935 "Tpetra::CrsMatrix object from it, but the file's matrix "
2936 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2941 if (tolerant && pBanner->matrixType() ==
"array") {
2942 bannerIsCorrect = 0;
2943 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2944 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2945 "object from it, but the file's matrix type is \"array\" "
2946 "instead. That probably means the file contains dense matrix "
2953 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2960 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2961 std::invalid_argument, errMsg.str ());
2963 if (debug && myRank == rootRank) {
2964 cerr <<
"-- Reading dimensions line" << endl;
2972 Tuple<global_ordinal_type, 3> dims =
2973 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2975 if (debug && myRank == rootRank) {
2976 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2981 RCP<adder_type> pAdder =
2982 makeAdder (pComm, pBanner, dims, tolerant, debug);
2984 if (debug && myRank == rootRank) {
2985 cerr <<
"-- Reading matrix data" << endl;
2995 int readSuccess = 1;
2996 std::ostringstream errMsg;
2997 if (myRank == rootRank) {
3000 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3002 reader_type reader (pAdder);
3005 std::pair<bool, std::vector<size_t> > results =
3006 reader.read (in, lineNumber, tolerant, debug);
3007 readSuccess = results.first ? 1 : 0;
3009 catch (std::exception& e) {
3014 broadcast (*pComm, rootRank, ptr (&readSuccess));
3023 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
3024 "Failed to read the Matrix Market sparse matrix file: "
3028 if (debug && myRank == rootRank) {
3029 cerr <<
"-- Successfully read the Matrix Market data" << endl;
3042 if (debug && myRank == rootRank) {
3043 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
3045 <<
"----- Dimensions before: "
3046 << dims[0] <<
" x " << dims[1]
3050 Tuple<global_ordinal_type, 2> updatedDims;
3051 if (myRank == rootRank) {
3058 std::max (dims[0], pAdder->getAdder()->numRows());
3059 updatedDims[1] = pAdder->getAdder()->numCols();
3061 broadcast (*pComm, rootRank, updatedDims);
3062 dims[0] = updatedDims[0];
3063 dims[1] = updatedDims[1];
3064 if (debug && myRank == rootRank) {
3065 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
3079 if (myRank == rootRank) {
3086 if (dims[0] < pAdder->getAdder ()->numRows ()) {
3090 broadcast (*pComm, 0, ptr (&dimsMatch));
3091 if (dimsMatch == 0) {
3098 Tuple<global_ordinal_type, 2> addersDims;
3099 if (myRank == rootRank) {
3100 addersDims[0] = pAdder->getAdder()->numRows();
3101 addersDims[1] = pAdder->getAdder()->numCols();
3103 broadcast (*pComm, 0, addersDims);
3104 TEUCHOS_TEST_FOR_EXCEPTION(
3105 dimsMatch == 0, std::runtime_error,
3106 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
3107 << dims[1] <<
", but the actual data says that the matrix is "
3108 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
3109 "data includes more rows than reported in the metadata. This "
3110 "is not allowed when parsing in strict mode. Parse the matrix in "
3111 "tolerant mode to ignore the metadata when it disagrees with the "
3116 if (debug && myRank == rootRank) {
3117 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
3129 ArrayRCP<size_t> numEntriesPerRow;
3130 ArrayRCP<size_t> rowPtr;
3131 ArrayRCP<global_ordinal_type> colInd;
3132 ArrayRCP<scalar_type> values;
3137 int mergeAndConvertSucceeded = 1;
3138 std::ostringstream errMsg;
3140 if (myRank == rootRank) {
3142 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3152 const size_type numRows = dims[0];
3155 pAdder->getAdder()->merge ();
3158 const std::vector<element_type>& entries =
3159 pAdder->getAdder()->getEntries();
3162 const size_t numEntries = (size_t)entries.size();
3165 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
3166 <<
" rows and numEntries=" << numEntries
3167 <<
" entries." << endl;
3173 numEntriesPerRow = arcp<size_t> (numRows);
3174 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3175 rowPtr = arcp<size_t> (numRows+1);
3176 std::fill (rowPtr.begin(), rowPtr.end(), 0);
3177 colInd = arcp<global_ordinal_type> (numEntries);
3178 values = arcp<scalar_type> (numEntries);
3185 for (curPos = 0; curPos < numEntries; ++curPos) {
3186 const element_type& curEntry = entries[curPos];
3188 TEUCHOS_TEST_FOR_EXCEPTION(
3189 curRow < prvRow, std::logic_error,
3190 "Row indices are out of order, even though they are supposed "
3191 "to be sorted. curRow = " << curRow <<
", prvRow = "
3192 << prvRow <<
", at curPos = " << curPos <<
". Please report "
3193 "this bug to the Tpetra developers.");
3194 if (curRow > prvRow) {
3200 numEntriesPerRow[curRow]++;
3201 colInd[curPos] = curEntry.colIndex();
3202 values[curPos] = curEntry.value();
3207 rowPtr[numRows] = numEntries;
3209 catch (std::exception& e) {
3210 mergeAndConvertSucceeded = 0;
3211 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3212 "CSR format: " << e.what();
3215 if (debug && mergeAndConvertSucceeded) {
3217 const size_type numRows = dims[0];
3218 const size_type maxToDisplay = 100;
3220 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3221 << (numEntriesPerRow.size()-1) <<
"] ";
3222 if (numRows > maxToDisplay) {
3223 cerr <<
"(only showing first and last few entries) ";
3227 if (numRows > maxToDisplay) {
3228 for (size_type k = 0; k < 2; ++k) {
3229 cerr << numEntriesPerRow[k] <<
" ";
3232 for (size_type k = numRows-2; k < numRows-1; ++k) {
3233 cerr << numEntriesPerRow[k] <<
" ";
3237 for (size_type k = 0; k < numRows-1; ++k) {
3238 cerr << numEntriesPerRow[k] <<
" ";
3241 cerr << numEntriesPerRow[numRows-1];
3243 cerr <<
"]" << endl;
3245 cerr <<
"----- Proc 0: rowPtr ";
3246 if (numRows > maxToDisplay) {
3247 cerr <<
"(only showing first and last few entries) ";
3250 if (numRows > maxToDisplay) {
3251 for (size_type k = 0; k < 2; ++k) {
3252 cerr << rowPtr[k] <<
" ";
3255 for (size_type k = numRows-2; k < numRows; ++k) {
3256 cerr << rowPtr[k] <<
" ";
3260 for (size_type k = 0; k < numRows; ++k) {
3261 cerr << rowPtr[k] <<
" ";
3264 cerr << rowPtr[numRows] <<
"]" << endl;
3275 if (debug && myRank == rootRank) {
3276 cerr <<
"-- Making range, domain, and row maps" << endl;
3283 RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
3284 RCP<const map_type> pDomainMap =
3285 makeDomainMap (pRangeMap, dims[0], dims[1]);
3286 RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
3288 if (debug && myRank == rootRank) {
3289 cerr <<
"-- Distributing the matrix data" << endl;
3303 ArrayRCP<size_t> myNumEntriesPerRow;
3304 ArrayRCP<size_t> myRowPtr;
3305 ArrayRCP<global_ordinal_type> myColInd;
3306 ArrayRCP<scalar_type> myValues;
3308 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
3309 numEntriesPerRow, rowPtr, colInd, values, debug);
3311 if (debug && myRank == rootRank) {
3312 cerr <<
"-- Inserting matrix entries on each process "
3313 "and calling fillComplete()" << endl;
3322 Teuchos::RCP<sparse_matrix_type> pMatrix =
3323 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
3324 pRowMap, pRangeMap, pDomainMap, constructorParams,
3325 fillCompleteParams);
3330 int localIsNull = pMatrix.is_null () ? 1 : 0;
3331 int globalIsNull = 0;
3332 reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
3333 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
3334 "Reader::makeMatrix() returned a null pointer on at least one "
3335 "process. Please report this bug to the Tpetra developers.");
3338 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
3339 "Reader::makeMatrix() returned a null pointer. "
3340 "Please report this bug to the Tpetra developers.");
3349 if (extraDebug && debug) {
3350 const int numProcs = pComm->getSize ();
3352 pRangeMap->getGlobalNumElements();
3354 pDomainMap->getGlobalNumElements();
3355 if (myRank == rootRank) {
3356 cerr <<
"-- Matrix is "
3357 << globalNumRows <<
" x " << globalNumCols
3358 <<
" with " << pMatrix->getGlobalNumEntries()
3359 <<
" entries, and index base "
3360 << pMatrix->getIndexBase() <<
"." << endl;
3363 for (
int p = 0; p < numProcs; ++p) {
3365 cerr <<
"-- Proc " << p <<
" owns "
3366 << pMatrix->getNodeNumCols() <<
" columns, and "
3367 << pMatrix->getNodeNumEntries() <<
" entries." << endl;
3373 if (debug && myRank == rootRank) {
3374 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3380 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
3381 static Teuchos::RCP<sparse_matrix_type> TPETRA_DEPRECATED
3384 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
3385 const Teuchos::RCP<node_type>& ,
3386 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
3387 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
3388 const bool tolerant=
false,
3389 const bool debug=
false)
3391 return readSparse (in, pComm, constructorParams,
3392 fillCompleteParams, tolerant, debug);
3394 #endif // TPETRA_ENABLE_DEPRECATED_CODE
3436 static Teuchos::RCP<sparse_matrix_type>
3438 const Teuchos::RCP<const map_type>& rowMap,
3439 Teuchos::RCP<const map_type>& colMap,
3440 const Teuchos::RCP<const map_type>& domainMap,
3441 const Teuchos::RCP<const map_type>& rangeMap,
3442 const bool callFillComplete=
true,
3443 const bool tolerant=
false,
3444 const bool debug=
false)
3446 using Teuchos::MatrixMarket::Banner;
3447 using Teuchos::arcp;
3448 using Teuchos::ArrayRCP;
3449 using Teuchos::ArrayView;
3451 using Teuchos::broadcast;
3452 using Teuchos::Comm;
3453 using Teuchos::null;
3456 using Teuchos::reduceAll;
3457 using Teuchos::Tuple;
3460 typedef Teuchos::ScalarTraits<scalar_type> STS;
3462 RCP<const Comm<int> > pComm = rowMap->getComm ();
3463 const int myRank = pComm->getRank ();
3464 const int rootRank = 0;
3465 const bool extraDebug =
false;
3470 TEUCHOS_TEST_FOR_EXCEPTION(
3471 rowMap.is_null (), std::invalid_argument,
3472 "Row Map must be nonnull.");
3473 TEUCHOS_TEST_FOR_EXCEPTION(
3474 rangeMap.is_null (), std::invalid_argument,
3475 "Range Map must be nonnull.");
3476 TEUCHOS_TEST_FOR_EXCEPTION(
3477 domainMap.is_null (), std::invalid_argument,
3478 "Domain Map must be nonnull.");
3479 TEUCHOS_TEST_FOR_EXCEPTION(
3480 rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
3481 std::invalid_argument,
3482 "The specified row Map's communicator (rowMap->getComm())"
3483 "differs from the given separately supplied communicator pComm.");
3484 TEUCHOS_TEST_FOR_EXCEPTION(
3485 domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
3486 std::invalid_argument,
3487 "The specified domain Map's communicator (domainMap->getComm())"
3488 "differs from the given separately supplied communicator pComm.");
3489 TEUCHOS_TEST_FOR_EXCEPTION(
3490 rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
3491 std::invalid_argument,
3492 "The specified range Map's communicator (rangeMap->getComm())"
3493 "differs from the given separately supplied communicator pComm.");
3498 size_t lineNumber = 1;
3500 if (debug && myRank == rootRank) {
3501 cerr <<
"Matrix Market reader: readSparse:" << endl
3502 <<
"-- Reading banner line" << endl;
3511 RCP<const Banner> pBanner;
3517 int bannerIsCorrect = 1;
3518 std::ostringstream errMsg;
3520 if (myRank == rootRank) {
3523 pBanner = readBanner (in, lineNumber, tolerant, debug);
3525 catch (std::exception& e) {
3526 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
3527 "threw an exception: " << e.what();
3528 bannerIsCorrect = 0;
3531 if (bannerIsCorrect) {
3536 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
3537 bannerIsCorrect = 0;
3538 errMsg <<
"The Matrix Market input file must contain a "
3539 "\"coordinate\"-format sparse matrix in order to create a "
3540 "Tpetra::CrsMatrix object from it, but the file's matrix "
3541 "type is \"" << pBanner->matrixType() <<
"\" instead.";
3546 if (tolerant && pBanner->matrixType() ==
"array") {
3547 bannerIsCorrect = 0;
3548 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
3549 "format sparse matrix in order to create a Tpetra::CrsMatrix "
3550 "object from it, but the file's matrix type is \"array\" "
3551 "instead. That probably means the file contains dense matrix "
3558 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
3565 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
3566 std::invalid_argument, errMsg.str ());
3568 if (debug && myRank == rootRank) {
3569 cerr <<
"-- Reading dimensions line" << endl;
3577 Tuple<global_ordinal_type, 3> dims =
3578 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
3580 if (debug && myRank == rootRank) {
3581 cerr <<
"-- Making Adder for collecting matrix data" << endl;
3588 RCP<adder_type> pAdder =
3589 makeAdder (pComm, pBanner, dims, tolerant, debug);
3591 if (debug && myRank == rootRank) {
3592 cerr <<
"-- Reading matrix data" << endl;
3602 int readSuccess = 1;
3603 std::ostringstream errMsg;
3604 if (myRank == rootRank) {
3607 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3609 reader_type reader (pAdder);
3612 std::pair<bool, std::vector<size_t> > results =
3613 reader.read (in, lineNumber, tolerant, debug);
3614 readSuccess = results.first ? 1 : 0;
3616 catch (std::exception& e) {
3621 broadcast (*pComm, rootRank, ptr (&readSuccess));
3630 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
3631 "Failed to read the Matrix Market sparse matrix file: "
3635 if (debug && myRank == rootRank) {
3636 cerr <<
"-- Successfully read the Matrix Market data" << endl;
3649 if (debug && myRank == rootRank) {
3650 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
3652 <<
"----- Dimensions before: "
3653 << dims[0] <<
" x " << dims[1]
3657 Tuple<global_ordinal_type, 2> updatedDims;
3658 if (myRank == rootRank) {
3665 std::max (dims[0], pAdder->getAdder()->numRows());
3666 updatedDims[1] = pAdder->getAdder()->numCols();
3668 broadcast (*pComm, rootRank, updatedDims);
3669 dims[0] = updatedDims[0];
3670 dims[1] = updatedDims[1];
3671 if (debug && myRank == rootRank) {
3672 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
3686 if (myRank == rootRank) {
3693 if (dims[0] < pAdder->getAdder ()->numRows ()) {
3697 broadcast (*pComm, 0, ptr (&dimsMatch));
3698 if (dimsMatch == 0) {
3705 Tuple<global_ordinal_type, 2> addersDims;
3706 if (myRank == rootRank) {
3707 addersDims[0] = pAdder->getAdder()->numRows();
3708 addersDims[1] = pAdder->getAdder()->numCols();
3710 broadcast (*pComm, 0, addersDims);
3711 TEUCHOS_TEST_FOR_EXCEPTION(
3712 dimsMatch == 0, std::runtime_error,
3713 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
3714 << dims[1] <<
", but the actual data says that the matrix is "
3715 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
3716 "data includes more rows than reported in the metadata. This "
3717 "is not allowed when parsing in strict mode. Parse the matrix in "
3718 "tolerant mode to ignore the metadata when it disagrees with the "
3723 if (debug && myRank == rootRank) {
3724 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
3736 ArrayRCP<size_t> numEntriesPerRow;
3737 ArrayRCP<size_t> rowPtr;
3738 ArrayRCP<global_ordinal_type> colInd;
3739 ArrayRCP<scalar_type> values;
3740 size_t maxNumEntriesPerRow = 0;
3745 int mergeAndConvertSucceeded = 1;
3746 std::ostringstream errMsg;
3748 if (myRank == rootRank) {
3750 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3760 const size_type numRows = dims[0];
3763 pAdder->getAdder()->merge ();
3766 const std::vector<element_type>& entries =
3767 pAdder->getAdder()->getEntries();
3770 const size_t numEntries = (size_t)entries.size();
3773 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
3774 <<
" rows and numEntries=" << numEntries
3775 <<
" entries." << endl;
3781 numEntriesPerRow = arcp<size_t> (numRows);
3782 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3783 rowPtr = arcp<size_t> (numRows+1);
3784 std::fill (rowPtr.begin(), rowPtr.end(), 0);
3785 colInd = arcp<global_ordinal_type> (numEntries);
3786 values = arcp<scalar_type> (numEntries);
3793 for (curPos = 0; curPos < numEntries; ++curPos) {
3794 const element_type& curEntry = entries[curPos];
3796 TEUCHOS_TEST_FOR_EXCEPTION(
3797 curRow < prvRow, std::logic_error,
3798 "Row indices are out of order, even though they are supposed "
3799 "to be sorted. curRow = " << curRow <<
", prvRow = "
3800 << prvRow <<
", at curPos = " << curPos <<
". Please report "
3801 "this bug to the Tpetra developers.");
3802 if (curRow > prvRow) {
3808 numEntriesPerRow[curRow]++;
3809 colInd[curPos] = curEntry.colIndex();
3810 values[curPos] = curEntry.value();
3815 rowPtr[numRows] = numEntries;
3817 catch (std::exception& e) {
3818 mergeAndConvertSucceeded = 0;
3819 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3820 "CSR format: " << e.what();
3823 if (debug && mergeAndConvertSucceeded) {
3825 const size_type numRows = dims[0];
3826 const size_type maxToDisplay = 100;
3828 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3829 << (numEntriesPerRow.size()-1) <<
"] ";
3830 if (numRows > maxToDisplay) {
3831 cerr <<
"(only showing first and last few entries) ";
3835 if (numRows > maxToDisplay) {
3836 for (size_type k = 0; k < 2; ++k) {
3837 cerr << numEntriesPerRow[k] <<
" ";
3840 for (size_type k = numRows-2; k < numRows-1; ++k) {
3841 cerr << numEntriesPerRow[k] <<
" ";
3845 for (size_type k = 0; k < numRows-1; ++k) {
3846 cerr << numEntriesPerRow[k] <<
" ";
3849 cerr << numEntriesPerRow[numRows-1];
3851 cerr <<
"]" << endl;
3853 cerr <<
"----- Proc 0: rowPtr ";
3854 if (numRows > maxToDisplay) {
3855 cerr <<
"(only showing first and last few entries) ";
3858 if (numRows > maxToDisplay) {
3859 for (size_type k = 0; k < 2; ++k) {
3860 cerr << rowPtr[k] <<
" ";
3863 for (size_type k = numRows-2; k < numRows; ++k) {
3864 cerr << rowPtr[k] <<
" ";
3868 for (size_type k = 0; k < numRows; ++k) {
3869 cerr << rowPtr[k] <<
" ";
3872 cerr << rowPtr[numRows] <<
"]" << endl;
3874 cerr <<
"----- Proc 0: colInd = [";
3875 for (
size_t k = 0; k < rowPtr[numRows]; ++k) {
3876 cerr << colInd[k] <<
" ";
3878 cerr <<
"]" << endl;
3892 if (debug && myRank == rootRank) {
3893 cerr <<
"-- Verifying Maps" << endl;
3895 TEUCHOS_TEST_FOR_EXCEPTION(
3896 as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
3897 std::invalid_argument,
3898 "The range Map has " << rangeMap->getGlobalNumElements ()
3899 <<
" entries, but the matrix has a global number of rows " << dims[0]
3901 TEUCHOS_TEST_FOR_EXCEPTION(
3902 as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
3903 std::invalid_argument,
3904 "The domain Map has " << domainMap->getGlobalNumElements ()
3905 <<
" entries, but the matrix has a global number of columns "
3909 RCP<Teuchos::FancyOStream> err = debug ?
3910 Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
3912 RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
3913 ArrayView<const global_ordinal_type> myRows =
3914 gatherRowMap->getNodeElementList ();
3915 const size_type myNumRows = myRows.size ();
3918 ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
3919 for (size_type i_ = 0; i_ < myNumRows; i_++) {
3920 gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_]-indexBase];
3921 if (gatherNumEntriesPerRow[i_] > maxNumEntriesPerRow)
3922 maxNumEntriesPerRow = gatherNumEntriesPerRow[i_];
3928 RCP<sparse_matrix_type> A_proc0 =
3930 Tpetra::StaticProfile));
3931 if (myRank == rootRank) {
3933 cerr <<
"-- Proc 0: Filling gather matrix" << endl;
3936 cerr <<
"---- Rows: " << Teuchos::toString (myRows) << endl;
3940 for (size_type i_ = 0; i_ < myNumRows; ++i_) {
3941 size_type i = myRows[i_] - indexBase;
3943 const size_type curPos = as<size_type> (rowPtr[i]);
3945 ArrayView<global_ordinal_type> curColInd =
3946 colInd.view (curPos, curNumEntries);
3947 ArrayView<scalar_type> curValues =
3948 values.view (curPos, curNumEntries);
3951 for (size_type k = 0; k < curNumEntries; ++k) {
3952 curColInd[k] += indexBase;
3955 cerr <<
"------ Columns: " << Teuchos::toString (curColInd) << endl;
3956 cerr <<
"------ Values: " << Teuchos::toString (curValues) << endl;
3959 if (curNumEntries > 0) {
3960 A_proc0->insertGlobalValues (myRows[i_], curColInd, curValues);
3966 numEntriesPerRow = null;
3972 broadcast<int,size_t> (*pComm, 0, &maxNumEntriesPerRow);
3974 RCP<sparse_matrix_type> A;
3975 if (colMap.is_null ()) {
3981 export_type exp (gatherRowMap, rowMap);
3982 A->doExport (*A_proc0, exp,
INSERT);
3984 if (callFillComplete) {
3985 A->fillComplete (domainMap, rangeMap);
3997 if (callFillComplete) {
3998 const int numProcs = pComm->getSize ();
4000 if (extraDebug && debug) {
4001 const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
4002 const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
4003 if (myRank == rootRank) {
4004 cerr <<
"-- Matrix is "
4005 << globalNumRows <<
" x " << globalNumCols
4006 <<
" with " << A->getGlobalNumEntries()
4007 <<
" entries, and index base "
4008 << A->getIndexBase() <<
"." << endl;
4011 for (
int p = 0; p < numProcs; ++p) {
4013 cerr <<
"-- Proc " << p <<
" owns "
4014 << A->getNodeNumCols() <<
" columns, and "
4015 << A->getNodeNumEntries() <<
" entries." << endl;
4022 if (debug && myRank == rootRank) {
4023 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
4057 static Teuchos::RCP<multivector_type>
4059 const Teuchos::RCP<const comm_type>& comm,
4060 Teuchos::RCP<const map_type>& map,
4061 const bool tolerant=
false,
4062 const bool debug=
false)
4065 if (comm->getRank () == 0) {
4066 in.open (filename.c_str ());
4068 return readDense (in, comm, map, tolerant, debug);
4071 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
4072 static Teuchos::RCP<multivector_type> TPETRA_DEPRECATED
4075 const Teuchos::RCP<const comm_type>& comm,
4076 const Teuchos::RCP<node_type>& ,
4077 Teuchos::RCP<const map_type>& map,
4078 const bool tolerant=
false,
4079 const bool debug=
false)
4083 #endif // TPETRA_ENABLE_DEPRECATED_CODE
4114 static Teuchos::RCP<vector_type>
4116 const Teuchos::RCP<const comm_type>& comm,
4117 Teuchos::RCP<const map_type>& map,
4118 const bool tolerant=
false,
4119 const bool debug=
false)
4122 if (comm->getRank () == 0) {
4123 in.open (filename.c_str ());
4125 return readVector (in, comm, map, tolerant, debug);
4128 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
4129 static Teuchos::RCP<vector_type> TPETRA_DEPRECATED
4133 const Teuchos::RCP<const comm_type>& comm,
4134 const Teuchos::RCP<node_type>& ,
4135 Teuchos::RCP<const map_type>& map,
4136 const bool tolerant=
false,
4137 const bool debug=
false)
4141 #endif // TPETRA_ENABLE_DEPRECATED_CODE
4209 static Teuchos::RCP<multivector_type>
4211 const Teuchos::RCP<const comm_type>& comm,
4212 Teuchos::RCP<const map_type>& map,
4213 const bool tolerant=
false,
4214 const bool debug=
false)
4216 Teuchos::RCP<Teuchos::FancyOStream> err =
4217 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4218 return readDenseImpl<scalar_type> (in, comm, map, err, tolerant, debug);
4221 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
4222 static Teuchos::RCP<multivector_type> TPETRA_DEPRECATED
4225 const Teuchos::RCP<const comm_type>& comm,
4226 const Teuchos::RCP<node_type>& ,
4227 Teuchos::RCP<const map_type>& map,
4228 const bool tolerant=
false,
4229 const bool debug=
false)
4231 return readDense (in, comm, map, tolerant, debug);
4233 #endif // TPETRA_ENABLE_DEPRECATED_CODE
4236 static Teuchos::RCP<vector_type>
4238 const Teuchos::RCP<const comm_type>& comm,
4239 Teuchos::RCP<const map_type>& map,
4240 const bool tolerant=
false,
4241 const bool debug=
false)
4243 Teuchos::RCP<Teuchos::FancyOStream> err =
4244 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4245 return readVectorImpl<scalar_type> (in, comm, map, err, tolerant, debug);
4248 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
4249 static Teuchos::RCP<vector_type> TPETRA_DEPRECATED
4252 const Teuchos::RCP<const comm_type>& comm,
4253 const Teuchos::RCP<node_type>& ,
4254 Teuchos::RCP<const map_type>& map,
4255 const bool tolerant=
false,
4256 const bool debug=
false)
4258 return readVector(in, comm, map, tolerant, debug);
4260 #endif // TPETRA_ENABLE_DEPRECATED_CODE
4282 static Teuchos::RCP<const map_type>
4284 const Teuchos::RCP<const comm_type>& comm,
4285 const bool tolerant=
false,
4286 const bool debug=
false)
4288 using Teuchos::inOutArg;
4289 using Teuchos::broadcast;
4293 if (comm->getRank () == 0) {
4294 in.open (filename.c_str ());
4299 broadcast<int, int> (*comm, 0, inOutArg (success));
4300 TEUCHOS_TEST_FOR_EXCEPTION(
4301 success == 0, std::runtime_error,
4302 "Tpetra::MatrixMarket::Reader::readMapFile: "
4303 "Failed to read file \"" << filename <<
"\" on Process 0.");
4304 return readMap (in, comm, tolerant, debug);
4307 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
4308 static Teuchos::RCP<const map_type> TPETRA_DEPRECATED
4312 const Teuchos::RCP<const comm_type>& comm,
4313 const Teuchos::RCP<node_type>& ,
4314 const bool tolerant=
false,
4315 const bool debug=
false)
4317 return readMapFile (filename, comm, tolerant, debug);
4319 #endif // TPETRA_ENABLE_DEPRECATED_CODE
4322 template<
class MultiVectorScalarType>
4327 readDenseImpl (std::istream& in,
4328 const Teuchos::RCP<const comm_type>& comm,
4329 Teuchos::RCP<const map_type>& map,
4330 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4331 const bool tolerant=
false,
4332 const bool debug=
false)
4334 using Teuchos::MatrixMarket::Banner;
4335 using Teuchos::MatrixMarket::checkCommentLine;
4336 using Teuchos::ArrayRCP;
4338 using Teuchos::broadcast;
4339 using Teuchos::outArg;
4341 using Teuchos::Tuple;
4343 typedef MultiVectorScalarType ST;
4347 typedef Teuchos::ScalarTraits<ST> STS;
4348 typedef typename STS::magnitudeType MT;
4349 typedef Teuchos::ScalarTraits<MT> STM;
4354 const int myRank = comm->getRank ();
4356 if (! err.is_null ()) {
4360 *err << myRank <<
": readDenseImpl" << endl;
4362 if (! err.is_null ()) {
4396 size_t lineNumber = 1;
4399 std::ostringstream exMsg;
4400 int localBannerReadSuccess = 1;
4401 int localDimsReadSuccess = 1;
4406 *err << myRank <<
": readDenseImpl: Reading banner line (dense)" << endl;
4412 RCP<const Banner> pBanner;
4414 pBanner = readBanner (in, lineNumber, tolerant, debug);
4415 }
catch (std::exception& e) {
4417 localBannerReadSuccess = 0;
4420 if (localBannerReadSuccess) {
4421 if (pBanner->matrixType () !=
"array") {
4422 exMsg <<
"The Matrix Market file does not contain dense matrix "
4423 "data. Its banner (first) line says that its matrix type is \""
4424 << pBanner->matrixType () <<
"\", rather that the required "
4426 localBannerReadSuccess = 0;
4427 }
else if (pBanner->dataType() ==
"pattern") {
4428 exMsg <<
"The Matrix Market file's banner (first) "
4429 "line claims that the matrix's data type is \"pattern\". This does "
4430 "not make sense for a dense matrix, yet the file reports the matrix "
4431 "as dense. The only valid data types for a dense matrix are "
4432 "\"real\", \"complex\", and \"integer\".";
4433 localBannerReadSuccess = 0;
4437 dims[2] = encodeDataType (pBanner->dataType ());
4443 if (localBannerReadSuccess) {
4445 *err << myRank <<
": readDenseImpl: Reading dimensions line (dense)" << endl;
4454 bool commentLine =
true;
4456 while (commentLine) {
4459 if (in.eof () || in.fail ()) {
4460 exMsg <<
"Unable to get array dimensions line (at all) from line "
4461 << lineNumber <<
" of input stream. The input stream "
4462 <<
"claims that it is "
4463 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4464 localDimsReadSuccess = 0;
4467 if (getline (in, line)) {
4474 size_t start = 0, size = 0;
4475 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4482 std::istringstream istr (line);
4486 if (istr.eof () || istr.fail ()) {
4487 exMsg <<
"Unable to read any data from line " << lineNumber
4488 <<
" of input; the line should contain the matrix dimensions "
4489 <<
"\"<numRows> <numCols>\".";
4490 localDimsReadSuccess = 0;
4495 exMsg <<
"Failed to get number of rows from line "
4496 << lineNumber <<
" of input; the line should contains the "
4497 <<
"matrix dimensions \"<numRows> <numCols>\".";
4498 localDimsReadSuccess = 0;
4500 dims[0] = theNumRows;
4502 exMsg <<
"No more data after number of rows on line "
4503 << lineNumber <<
" of input; the line should contain the "
4504 <<
"matrix dimensions \"<numRows> <numCols>\".";
4505 localDimsReadSuccess = 0;
4510 exMsg <<
"Failed to get number of columns from line "
4511 << lineNumber <<
" of input; the line should contain "
4512 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4513 localDimsReadSuccess = 0;
4515 dims[1] = theNumCols;
4526 Tuple<GO, 5> bannerDimsReadResult;
4528 bannerDimsReadResult[0] = dims[0];
4529 bannerDimsReadResult[1] = dims[1];
4530 bannerDimsReadResult[2] = dims[2];
4531 bannerDimsReadResult[3] = localBannerReadSuccess;
4532 bannerDimsReadResult[4] = localDimsReadSuccess;
4536 broadcast (*comm, 0, bannerDimsReadResult);
4538 TEUCHOS_TEST_FOR_EXCEPTION(
4539 bannerDimsReadResult[3] == 0, std::runtime_error,
4540 "Failed to read banner line: " << exMsg.str ());
4541 TEUCHOS_TEST_FOR_EXCEPTION(
4542 bannerDimsReadResult[4] == 0, std::runtime_error,
4543 "Failed to read matrix dimensions line: " << exMsg.str ());
4545 dims[0] = bannerDimsReadResult[0];
4546 dims[1] = bannerDimsReadResult[1];
4547 dims[2] = bannerDimsReadResult[2];
4552 const size_t numCols =
static_cast<size_t> (dims[1]);
4557 RCP<const map_type> proc0Map;
4558 if (map.is_null ()) {
4562 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4563 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4565 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4569 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4573 RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
4579 int localReadDataSuccess = 1;
4583 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4588 TEUCHOS_TEST_FOR_EXCEPTION(
4589 ! X->isConstantStride (), std::logic_error,
4590 "Can't get a 1-D view of the entries of the MultiVector X on "
4591 "Process 0, because the stride between the columns of X is not "
4592 "constant. This shouldn't happen because we just created X and "
4593 "haven't filled it in yet. Please report this bug to the Tpetra "
4600 ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4601 TEUCHOS_TEST_FOR_EXCEPTION(
4602 as<global_size_t> (X_view.size ()) < numRows * numCols,
4604 "The view of X has size " << X_view <<
" which is not enough to "
4605 "accommodate the expected number of entries numRows*numCols = "
4606 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
4607 "Please report this bug to the Tpetra developers.");
4608 const size_t stride = X->getStride ();
4615 const bool isComplex = (dims[2] == 1);
4616 size_type count = 0, curRow = 0, curCol = 0;
4619 while (getline (in, line)) {
4623 size_t start = 0, size = 0;
4624 const bool commentLine =
4625 checkCommentLine (line, start, size, lineNumber, tolerant);
4626 if (! commentLine) {
4632 if (count >= X_view.size()) {
4637 TEUCHOS_TEST_FOR_EXCEPTION(
4638 count >= X_view.size(),
4640 "The Matrix Market input stream has more data in it than "
4641 "its metadata reported. Current line number is "
4642 << lineNumber <<
".");
4651 const size_t pos = line.substr (start, size).find (
':');
4652 if (pos != std::string::npos) {
4656 std::istringstream istr (line.substr (start, size));
4660 if (istr.eof() || istr.fail()) {
4667 TEUCHOS_TEST_FOR_EXCEPTION(
4668 ! tolerant && (istr.eof() || istr.fail()),
4670 "Line " << lineNumber <<
" of the Matrix Market file is "
4671 "empty, or we cannot read from it for some other reason.");
4674 ST val = STS::zero();
4677 MT real = STM::zero(), imag = STM::zero();
4694 TEUCHOS_TEST_FOR_EXCEPTION(
4695 ! tolerant && istr.eof(), std::runtime_error,
4696 "Failed to get the real part of a complex-valued matrix "
4697 "entry from line " << lineNumber <<
" of the Matrix Market "
4703 }
else if (istr.eof()) {
4704 TEUCHOS_TEST_FOR_EXCEPTION(
4705 ! tolerant && istr.eof(), std::runtime_error,
4706 "Missing imaginary part of a complex-valued matrix entry "
4707 "on line " << lineNumber <<
" of the Matrix Market file.");
4713 TEUCHOS_TEST_FOR_EXCEPTION(
4714 ! tolerant && istr.fail(), std::runtime_error,
4715 "Failed to get the imaginary part of a complex-valued "
4716 "matrix entry from line " << lineNumber <<
" of the "
4717 "Matrix Market file.");
4724 TEUCHOS_TEST_FOR_EXCEPTION(
4725 ! tolerant && istr.fail(), std::runtime_error,
4726 "Failed to get a real-valued matrix entry from line "
4727 << lineNumber <<
" of the Matrix Market file.");
4730 if (istr.fail() && tolerant) {
4736 TEUCHOS_TEST_FOR_EXCEPTION(
4737 ! tolerant && istr.fail(), std::runtime_error,
4738 "Failed to read matrix data from line " << lineNumber
4739 <<
" of the Matrix Market file.");
4742 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4744 curRow = count % numRows;
4745 curCol = count / numRows;
4746 X_view[curRow + curCol*stride] = val;
4751 TEUCHOS_TEST_FOR_EXCEPTION(
4752 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
4754 "The Matrix Market metadata reports that the dense matrix is "
4755 << numRows <<
" x " << numCols <<
", and thus has "
4756 << numRows*numCols <<
" total entries, but we only found " << count
4757 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
4758 }
catch (std::exception& e) {
4760 localReadDataSuccess = 0;
4765 *err << myRank <<
": readDenseImpl: done reading data" << endl;
4769 int globalReadDataSuccess = localReadDataSuccess;
4770 broadcast (*comm, 0, outArg (globalReadDataSuccess));
4771 TEUCHOS_TEST_FOR_EXCEPTION(
4772 globalReadDataSuccess == 0, std::runtime_error,
4773 "Failed to read the multivector's data: " << exMsg.str ());
4778 if (comm->getSize () == 1 && map.is_null ()) {
4780 if (! err.is_null ()) {
4784 *err << myRank <<
": readDenseImpl: done" << endl;
4786 if (! err.is_null ()) {
4793 *err << myRank <<
": readDenseImpl: Creating target MV" << endl;
4797 RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
4800 *err << myRank <<
": readDenseImpl: Creating Export" << endl;
4806 Export<LO, GO, NT> exporter (proc0Map, map, err);
4809 *err << myRank <<
": readDenseImpl: Exporting" << endl;
4812 Y->doExport (*X, exporter,
INSERT);
4814 if (! err.is_null ()) {
4818 *err << myRank <<
": readDenseImpl: done" << endl;
4820 if (! err.is_null ()) {
4829 template<
class VectorScalarType>
4834 readVectorImpl (std::istream& in,
4835 const Teuchos::RCP<const comm_type>& comm,
4836 Teuchos::RCP<const map_type>& map,
4837 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4838 const bool tolerant=
false,
4839 const bool debug=
false)
4841 using Teuchos::MatrixMarket::Banner;
4842 using Teuchos::MatrixMarket::checkCommentLine;
4844 using Teuchos::broadcast;
4845 using Teuchos::outArg;
4847 using Teuchos::Tuple;
4849 typedef VectorScalarType ST;
4853 typedef Teuchos::ScalarTraits<ST> STS;
4854 typedef typename STS::magnitudeType MT;
4855 typedef Teuchos::ScalarTraits<MT> STM;
4860 const int myRank = comm->getRank ();
4862 if (! err.is_null ()) {
4866 *err << myRank <<
": readVectorImpl" << endl;
4868 if (! err.is_null ()) {
4902 size_t lineNumber = 1;
4905 std::ostringstream exMsg;
4906 int localBannerReadSuccess = 1;
4907 int localDimsReadSuccess = 1;
4912 *err << myRank <<
": readVectorImpl: Reading banner line (dense)" << endl;
4918 RCP<const Banner> pBanner;
4920 pBanner = readBanner (in, lineNumber, tolerant, debug);
4921 }
catch (std::exception& e) {
4923 localBannerReadSuccess = 0;
4926 if (localBannerReadSuccess) {
4927 if (pBanner->matrixType () !=
"array") {
4928 exMsg <<
"The Matrix Market file does not contain dense matrix "
4929 "data. Its banner (first) line says that its matrix type is \""
4930 << pBanner->matrixType () <<
"\", rather that the required "
4932 localBannerReadSuccess = 0;
4933 }
else if (pBanner->dataType() ==
"pattern") {
4934 exMsg <<
"The Matrix Market file's banner (first) "
4935 "line claims that the matrix's data type is \"pattern\". This does "
4936 "not make sense for a dense matrix, yet the file reports the matrix "
4937 "as dense. The only valid data types for a dense matrix are "
4938 "\"real\", \"complex\", and \"integer\".";
4939 localBannerReadSuccess = 0;
4943 dims[2] = encodeDataType (pBanner->dataType ());
4949 if (localBannerReadSuccess) {
4951 *err << myRank <<
": readVectorImpl: Reading dimensions line (dense)" << endl;
4960 bool commentLine =
true;
4962 while (commentLine) {
4965 if (in.eof () || in.fail ()) {
4966 exMsg <<
"Unable to get array dimensions line (at all) from line "
4967 << lineNumber <<
" of input stream. The input stream "
4968 <<
"claims that it is "
4969 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4970 localDimsReadSuccess = 0;
4973 if (getline (in, line)) {
4980 size_t start = 0, size = 0;
4981 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4988 std::istringstream istr (line);
4992 if (istr.eof () || istr.fail ()) {
4993 exMsg <<
"Unable to read any data from line " << lineNumber
4994 <<
" of input; the line should contain the matrix dimensions "
4995 <<
"\"<numRows> <numCols>\".";
4996 localDimsReadSuccess = 0;
5001 exMsg <<
"Failed to get number of rows from line "
5002 << lineNumber <<
" of input; the line should contains the "
5003 <<
"matrix dimensions \"<numRows> <numCols>\".";
5004 localDimsReadSuccess = 0;
5006 dims[0] = theNumRows;
5008 exMsg <<
"No more data after number of rows on line "
5009 << lineNumber <<
" of input; the line should contain the "
5010 <<
"matrix dimensions \"<numRows> <numCols>\".";
5011 localDimsReadSuccess = 0;
5016 exMsg <<
"Failed to get number of columns from line "
5017 << lineNumber <<
" of input; the line should contain "
5018 <<
"the matrix dimensions \"<numRows> <numCols>\".";
5019 localDimsReadSuccess = 0;
5021 dims[1] = theNumCols;
5031 exMsg <<
"File does not contain a 1D Vector";
5032 localDimsReadSuccess = 0;
5038 Tuple<GO, 5> bannerDimsReadResult;
5040 bannerDimsReadResult[0] = dims[0];
5041 bannerDimsReadResult[1] = dims[1];
5042 bannerDimsReadResult[2] = dims[2];
5043 bannerDimsReadResult[3] = localBannerReadSuccess;
5044 bannerDimsReadResult[4] = localDimsReadSuccess;
5049 broadcast (*comm, 0, bannerDimsReadResult);
5051 TEUCHOS_TEST_FOR_EXCEPTION(
5052 bannerDimsReadResult[3] == 0, std::runtime_error,
5053 "Failed to read banner line: " << exMsg.str ());
5054 TEUCHOS_TEST_FOR_EXCEPTION(
5055 bannerDimsReadResult[4] == 0, std::runtime_error,
5056 "Failed to read matrix dimensions line: " << exMsg.str ());
5058 dims[0] = bannerDimsReadResult[0];
5059 dims[1] = bannerDimsReadResult[1];
5060 dims[2] = bannerDimsReadResult[2];
5065 const size_t numCols =
static_cast<size_t> (dims[1]);
5070 RCP<const map_type> proc0Map;
5071 if (map.is_null ()) {
5075 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
5076 const size_t localNumRows = (myRank == 0) ? numRows : 0;
5078 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
5082 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
5086 RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
5092 int localReadDataSuccess = 1;
5096 *err << myRank <<
": readVectorImpl: Reading matrix data (dense)"
5101 TEUCHOS_TEST_FOR_EXCEPTION(
5102 ! X->isConstantStride (), std::logic_error,
5103 "Can't get a 1-D view of the entries of the MultiVector X on "
5104 "Process 0, because the stride between the columns of X is not "
5105 "constant. This shouldn't happen because we just created X and "
5106 "haven't filled it in yet. Please report this bug to the Tpetra "
5113 Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst ();
5114 TEUCHOS_TEST_FOR_EXCEPTION(
5115 as<global_size_t> (X_view.size ()) < numRows * numCols,
5117 "The view of X has size " << X_view <<
" which is not enough to "
5118 "accommodate the expected number of entries numRows*numCols = "
5119 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
5120 "Please report this bug to the Tpetra developers.");
5121 const size_t stride = X->getStride ();
5128 const bool isComplex = (dims[2] == 1);
5129 size_type count = 0, curRow = 0, curCol = 0;
5132 while (getline (in, line)) {
5136 size_t start = 0, size = 0;
5137 const bool commentLine =
5138 checkCommentLine (line, start, size, lineNumber, tolerant);
5139 if (! commentLine) {
5145 if (count >= X_view.size()) {
5150 TEUCHOS_TEST_FOR_EXCEPTION(
5151 count >= X_view.size(),
5153 "The Matrix Market input stream has more data in it than "
5154 "its metadata reported. Current line number is "
5155 << lineNumber <<
".");
5164 const size_t pos = line.substr (start, size).find (
':');
5165 if (pos != std::string::npos) {
5169 std::istringstream istr (line.substr (start, size));
5173 if (istr.eof() || istr.fail()) {
5180 TEUCHOS_TEST_FOR_EXCEPTION(
5181 ! tolerant && (istr.eof() || istr.fail()),
5183 "Line " << lineNumber <<
" of the Matrix Market file is "
5184 "empty, or we cannot read from it for some other reason.");
5187 ST val = STS::zero();
5190 MT real = STM::zero(), imag = STM::zero();
5207 TEUCHOS_TEST_FOR_EXCEPTION(
5208 ! tolerant && istr.eof(), std::runtime_error,
5209 "Failed to get the real part of a complex-valued matrix "
5210 "entry from line " << lineNumber <<
" of the Matrix Market "
5216 }
else if (istr.eof()) {
5217 TEUCHOS_TEST_FOR_EXCEPTION(
5218 ! tolerant && istr.eof(), std::runtime_error,
5219 "Missing imaginary part of a complex-valued matrix entry "
5220 "on line " << lineNumber <<
" of the Matrix Market file.");
5226 TEUCHOS_TEST_FOR_EXCEPTION(
5227 ! tolerant && istr.fail(), std::runtime_error,
5228 "Failed to get the imaginary part of a complex-valued "
5229 "matrix entry from line " << lineNumber <<
" of the "
5230 "Matrix Market file.");
5237 TEUCHOS_TEST_FOR_EXCEPTION(
5238 ! tolerant && istr.fail(), std::runtime_error,
5239 "Failed to get a real-valued matrix entry from line "
5240 << lineNumber <<
" of the Matrix Market file.");
5243 if (istr.fail() && tolerant) {
5249 TEUCHOS_TEST_FOR_EXCEPTION(
5250 ! tolerant && istr.fail(), std::runtime_error,
5251 "Failed to read matrix data from line " << lineNumber
5252 <<
" of the Matrix Market file.");
5255 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
5257 curRow = count % numRows;
5258 curCol = count / numRows;
5259 X_view[curRow + curCol*stride] = val;
5264 TEUCHOS_TEST_FOR_EXCEPTION(
5265 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
5267 "The Matrix Market metadata reports that the dense matrix is "
5268 << numRows <<
" x " << numCols <<
", and thus has "
5269 << numRows*numCols <<
" total entries, but we only found " << count
5270 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
5271 }
catch (std::exception& e) {
5273 localReadDataSuccess = 0;
5278 *err << myRank <<
": readVectorImpl: done reading data" << endl;
5282 int globalReadDataSuccess = localReadDataSuccess;
5283 broadcast (*comm, 0, outArg (globalReadDataSuccess));
5284 TEUCHOS_TEST_FOR_EXCEPTION(
5285 globalReadDataSuccess == 0, std::runtime_error,
5286 "Failed to read the multivector's data: " << exMsg.str ());
5291 if (comm->getSize () == 1 && map.is_null ()) {
5293 if (! err.is_null ()) {
5297 *err << myRank <<
": readVectorImpl: done" << endl;
5299 if (! err.is_null ()) {
5306 *err << myRank <<
": readVectorImpl: Creating target MV" << endl;
5310 RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
5313 *err << myRank <<
": readVectorImpl: Creating Export" << endl;
5319 Export<LO, GO, NT> exporter (proc0Map, map, err);
5322 *err << myRank <<
": readVectorImpl: Exporting" << endl;
5325 Y->doExport (*X, exporter,
INSERT);
5327 if (! err.is_null ()) {
5331 *err << myRank <<
": readVectorImpl: done" << endl;
5333 if (! err.is_null ()) {
5361 static Teuchos::RCP<const map_type>
5363 const Teuchos::RCP<const comm_type>& comm,
5364 const bool tolerant=
false,
5365 const bool debug=
false)
5367 Teuchos::RCP<Teuchos::FancyOStream> err =
5368 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
5369 return readMap (in, comm, err, tolerant, debug);
5372 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
5373 static Teuchos::RCP<const map_type> TPETRA_DEPRECATED
5377 const Teuchos::RCP<const comm_type>& comm,
5378 const Teuchos::RCP<node_type>& ,
5379 const bool tolerant=
false,
5380 const bool debug=
false)
5382 return readMap(in, comm, tolerant, debug);
5384 #endif // TPETRA_ENABLE_DEPRECATED_CODE
5411 static Teuchos::RCP<const map_type>
5413 const Teuchos::RCP<const comm_type>& comm,
5414 const Teuchos::RCP<Teuchos::FancyOStream>& err,
5415 const bool tolerant=
false,
5416 const bool debug=
false)
5418 using Teuchos::arcp;
5419 using Teuchos::Array;
5420 using Teuchos::ArrayRCP;
5422 using Teuchos::broadcast;
5423 using Teuchos::Comm;
5424 using Teuchos::CommRequest;
5425 using Teuchos::inOutArg;
5426 using Teuchos::ireceive;
5427 using Teuchos::outArg;
5429 using Teuchos::receive;
5430 using Teuchos::reduceAll;
5431 using Teuchos::REDUCE_MIN;
5432 using Teuchos::isend;
5433 using Teuchos::SerialComm;
5434 using Teuchos::toString;
5435 using Teuchos::wait;
5438 typedef ptrdiff_t int_type;
5444 const int numProcs = comm->getSize ();
5445 const int myRank = comm->getRank ();
5447 if (err.is_null ()) {
5451 std::ostringstream os;
5452 os << myRank <<
": readMap: " << endl;
5455 if (err.is_null ()) {
5461 const int sizeTag = 1339;
5463 const int dataTag = 1340;
5469 RCP<CommRequest<int> > sizeReq;
5470 RCP<CommRequest<int> > dataReq;
5475 ArrayRCP<int_type> numGidsToRecv (1);
5476 numGidsToRecv[0] = 0;
5478 sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
5481 int readSuccess = 1;
5482 std::ostringstream exMsg;
5491 RCP<const Comm<int> > proc0Comm (
new SerialComm<int> ());
5493 RCP<const map_type> dataMap;
5497 data = readDenseImpl<GO> (in, proc0Comm, dataMap, err, tolerant, debug);
5499 if (data.is_null ()) {
5501 exMsg <<
"readDenseImpl() returned null." << endl;
5503 }
catch (std::exception& e) {
5505 exMsg << e.what () << endl;
5511 std::map<int, Array<GO> > pid2gids;
5516 int_type globalNumGIDs = 0;
5526 if (myRank == 0 && readSuccess == 1) {
5527 if (data->getNumVectors () == 2) {
5528 ArrayRCP<const GO> GIDs = data->getData (0);
5529 ArrayRCP<const GO> PIDs = data->getData (1);
5530 globalNumGIDs = GIDs.size ();
5534 if (globalNumGIDs > 0) {
5535 const int pid =
static_cast<int> (PIDs[0]);
5537 if (pid < 0 || pid >= numProcs) {
5539 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5540 <<
"Encountered invalid PID " << pid <<
"." << endl;
5543 const GO gid = GIDs[0];
5544 pid2gids[pid].push_back (gid);
5548 if (readSuccess == 1) {
5551 for (size_type k = 1; k < globalNumGIDs; ++k) {
5552 const int pid =
static_cast<int> (PIDs[k]);
5553 if (pid < 0 || pid >= numProcs) {
5555 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5556 <<
"Encountered invalid PID " << pid <<
"." << endl;
5559 const int_type gid = GIDs[k];
5560 pid2gids[pid].push_back (gid);
5561 if (gid < indexBase) {
5568 else if (data->getNumVectors () == 1) {
5569 if (data->getGlobalLength () % 2 !=
static_cast<GST
> (0)) {
5571 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data has the "
5572 "wrong format (for Map format 2.0). The global number of rows "
5573 "in the MultiVector must be even (divisible by 2)." << endl;
5576 ArrayRCP<const GO> theData = data->getData (0);
5577 globalNumGIDs =
static_cast<GO
> (data->getGlobalLength ()) /
5578 static_cast<GO> (2);
5582 if (globalNumGIDs > 0) {
5583 const int pid =
static_cast<int> (theData[1]);
5584 if (pid < 0 || pid >= numProcs) {
5586 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5587 <<
"Encountered invalid PID " << pid <<
"." << endl;
5590 const GO gid = theData[0];
5591 pid2gids[pid].push_back (gid);
5597 for (int_type k = 1; k < globalNumGIDs; ++k) {
5598 const int pid =
static_cast<int> (theData[2*k + 1]);
5599 if (pid < 0 || pid >= numProcs) {
5601 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5602 <<
"Encountered invalid PID " << pid <<
"." << endl;
5605 const GO gid = theData[2*k];
5606 pid2gids[pid].push_back (gid);
5607 if (gid < indexBase) {
5616 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data must have "
5617 "either 1 column (for the new Map format 2.0) or 2 columns (for "
5618 "the old Map format 1.0).";
5625 int_type readResults[3];
5626 readResults[0] =
static_cast<int_type
> (indexBase);
5627 readResults[1] =
static_cast<int_type
> (globalNumGIDs);
5628 readResults[2] =
static_cast<int_type
> (readSuccess);
5629 broadcast<int, int_type> (*comm, 0, 3, readResults);
5631 indexBase =
static_cast<GO
> (readResults[0]);
5632 globalNumGIDs =
static_cast<int_type
> (readResults[1]);
5633 readSuccess =
static_cast<int> (readResults[2]);
5639 TEUCHOS_TEST_FOR_EXCEPTION(
5640 readSuccess != 1, std::runtime_error,
5641 "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
5642 "following exception message: " << exMsg.str ());
5646 for (
int p = 1; p < numProcs; ++p) {
5647 ArrayRCP<int_type> numGidsToSend (1);
5649 auto it = pid2gids.find (p);
5650 if (it == pid2gids.end ()) {
5651 numGidsToSend[0] = 0;
5653 numGidsToSend[0] = it->second.size ();
5655 sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
5656 wait<int> (*comm, outArg (sizeReq));
5661 wait<int> (*comm, outArg (sizeReq));
5666 ArrayRCP<GO> myGids;
5667 int_type myNumGids = 0;
5669 GO* myGidsRaw = NULL;
5671 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
5672 if (it != pid2gids.end ()) {
5673 myGidsRaw = it->second.getRawPtr ();
5674 myNumGids = it->second.size ();
5676 myGids = arcp<GO> (myGidsRaw, 0, myNumGids,
false);
5680 myNumGids = numGidsToRecv[0];
5681 myGids = arcp<GO> (myNumGids);
5688 if (myNumGids > 0) {
5689 dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
5693 for (
int p = 1; p < numProcs; ++p) {
5695 ArrayRCP<GO> sendGids;
5696 GO* sendGidsRaw = NULL;
5697 int_type numSendGids = 0;
5699 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
5700 if (it != pid2gids.end ()) {
5701 numSendGids = it->second.size ();
5702 sendGidsRaw = it->second.getRawPtr ();
5703 sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids,
false);
5706 if (numSendGids > 0) {
5707 dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
5709 wait<int> (*comm, outArg (dataReq));
5711 else if (myRank == p) {
5713 wait<int> (*comm, outArg (dataReq));
5718 std::ostringstream os;
5719 os << myRank <<
": readMap: creating Map" << endl;
5722 const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
5723 RCP<const map_type> newMap;
5730 std::ostringstream errStrm;
5732 newMap = rcp (
new map_type (INVALID, myGids (), indexBase, comm));
5734 catch (std::exception& e) {
5736 errStrm <<
"Process " << comm->getRank () <<
" threw an exception: "
5737 << e.what () << std::endl;
5741 errStrm <<
"Process " << comm->getRank () <<
" threw an exception "
5742 "not a subclass of std::exception" << std::endl;
5744 Teuchos::reduceAll<int, int> (*comm, Teuchos::REDUCE_MIN,
5745 lclSuccess, Teuchos::outArg (gblSuccess));
5746 if (gblSuccess != 1) {
5749 TEUCHOS_TEST_FOR_EXCEPTION(gblSuccess != 1, std::runtime_error,
"Map constructor failed!");
5751 if (err.is_null ()) {
5755 std::ostringstream os;
5756 os << myRank <<
": readMap: done" << endl;
5759 if (err.is_null ()) {
5765 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
5766 static Teuchos::RCP<const map_type> TPETRA_DEPRECATED
5770 const Teuchos::RCP<const comm_type>& comm,
5771 const Teuchos::RCP<node_type>& ,
5772 const Teuchos::RCP<Teuchos::FancyOStream>& err,
5773 const bool tolerant=
false,
5774 const bool debug=
false)
5776 return readMap (in, comm, err, tolerant, debug);
5778 #endif // TPETRA_ENABLE_DEPRECATED_CODE
5793 encodeDataType (
const std::string& dataType)
5795 if (dataType ==
"real" || dataType ==
"integer") {
5797 }
else if (dataType ==
"complex") {
5799 }
else if (dataType ==
"pattern") {
5805 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
5806 "Unrecognized Matrix Market data type \"" << dataType
5807 <<
"\". We should never get here. "
5808 "Please report this bug to the Tpetra developers.");
5841 template<
class SparseMatrixType>
5846 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
5908 const std::string& matrixName,
5909 const std::string& matrixDescription,
5910 const bool debug=
false)
5912 Teuchos::RCP<const Teuchos::Comm<int> > comm = matrix.getComm ();
5913 TEUCHOS_TEST_FOR_EXCEPTION
5914 (comm.is_null (), std::invalid_argument,
5915 "The input matrix's communicator (Teuchos::Comm object) is null.");
5916 const int myRank = comm->getRank ();
5921 out.open (filename.c_str ());
5923 writeSparse (out, matrix, matrixName, matrixDescription, debug);
5932 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5933 const std::string& matrixName,
5934 const std::string& matrixDescription,
5935 const bool debug=
false)
5937 TEUCHOS_TEST_FOR_EXCEPTION
5938 (pMatrix.is_null (), std::invalid_argument,
5939 "The input matrix is null.");
5941 matrixDescription, debug);
5966 const bool debug=
false)
5974 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5975 const bool debug=
false)
6013 const std::string& matrixName,
6014 const std::string& matrixDescription,
6015 const bool debug=
false)
6017 using Teuchos::ArrayView;
6018 using Teuchos::Comm;
6019 using Teuchos::FancyOStream;
6020 using Teuchos::getFancyOStream;
6021 using Teuchos::null;
6023 using Teuchos::rcpFromRef;
6029 using STS =
typename Teuchos::ScalarTraits<ST>;
6035 Teuchos::SetScientific<ST> sci (out);
6038 RCP<const Comm<int> > comm = matrix.getComm ();
6039 TEUCHOS_TEST_FOR_EXCEPTION(
6040 comm.is_null (), std::invalid_argument,
6041 "The input matrix's communicator (Teuchos::Comm object) is null.");
6042 const int myRank = comm->getRank ();
6045 RCP<FancyOStream> err =
6046 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6048 std::ostringstream os;
6049 os << myRank <<
": writeSparse" << endl;
6052 os <<
"-- " << myRank <<
": past barrier" << endl;
6057 const bool debugPrint = debug && myRank == 0;
6059 RCP<const map_type> rowMap = matrix.getRowMap ();
6060 RCP<const map_type> colMap = matrix.getColMap ();
6061 RCP<const map_type> domainMap = matrix.getDomainMap ();
6062 RCP<const map_type> rangeMap = matrix.getRangeMap ();
6064 const global_size_t numRows = rangeMap->getGlobalNumElements ();
6065 const global_size_t numCols = domainMap->getGlobalNumElements ();
6067 if (debug && myRank == 0) {
6068 std::ostringstream os;
6069 os <<
"-- Input sparse matrix is:"
6070 <<
"---- " << numRows <<
" x " << numCols << endl
6072 << (matrix.isGloballyIndexed() ?
"Globally" :
"Locally")
6073 <<
" indexed." << endl
6074 <<
"---- Its row map has " << rowMap->getGlobalNumElements ()
6075 <<
" elements." << endl
6076 <<
"---- Its col map has " << colMap->getGlobalNumElements ()
6077 <<
" elements." << endl;
6082 const size_t localNumRows = (myRank == 0) ? numRows : 0;
6084 std::ostringstream os;
6085 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6088 RCP<const map_type> gatherRowMap =
6089 Details::computeGatherMap (rowMap, err, debug);
6097 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6098 RCP<const map_type> gatherColMap =
6099 Details::computeGatherMap (colMap, err, debug);
6103 import_type importer (rowMap, gatherRowMap);
6108 RCP<sparse_matrix_type> newMatrix =
6110 static_cast<size_t> (0)));
6112 newMatrix->doImport (matrix, importer,
INSERT);
6117 RCP<const map_type> gatherDomainMap =
6118 rcp (
new map_type (numCols, localNumCols,
6119 domainMap->getIndexBase (),
6121 RCP<const map_type> gatherRangeMap =
6122 rcp (
new map_type (numRows, localNumRows,
6123 rangeMap->getIndexBase (),
6125 newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
6129 cerr <<
"-- Output sparse matrix is:"
6130 <<
"---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
6132 << newMatrix->getDomainMap ()->getGlobalNumElements ()
6134 << newMatrix->getGlobalNumEntries () <<
" entries;" << endl
6136 << (newMatrix->isGloballyIndexed () ?
"Globally" :
"Locally")
6137 <<
" indexed." << endl
6138 <<
"---- Its row map has "
6139 << newMatrix->getRowMap ()->getGlobalNumElements ()
6140 <<
" elements, with index base "
6141 << newMatrix->getRowMap ()->getIndexBase () <<
"." << endl
6142 <<
"---- Its col map has "
6143 << newMatrix->getColMap ()->getGlobalNumElements ()
6144 <<
" elements, with index base "
6145 << newMatrix->getColMap ()->getIndexBase () <<
"." << endl
6146 <<
"---- Element count of output matrix's column Map may differ "
6147 <<
"from that of the input matrix's column Map, if some columns "
6148 <<
"of the matrix contain no entries." << endl;
6161 out <<
"%%MatrixMarket matrix coordinate "
6162 << (STS::isComplex ?
"complex" :
"real")
6163 <<
" general" << endl;
6166 if (matrixName !=
"") {
6167 printAsComment (out, matrixName);
6169 if (matrixDescription !=
"") {
6170 printAsComment (out, matrixDescription);
6180 out << newMatrix->getRangeMap ()->getGlobalNumElements () <<
" "
6181 << newMatrix->getDomainMap ()->getGlobalNumElements () <<
" "
6182 << newMatrix->getGlobalNumEntries () << endl;
6187 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6188 const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
6196 if (newMatrix->isGloballyIndexed()) {
6199 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6200 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6201 for (GO globalRowIndex = minAllGlobalIndex;
6202 globalRowIndex <= maxAllGlobalIndex;
6204 ArrayView<const GO> ind;
6205 ArrayView<const ST> val;
6206 newMatrix->getGlobalRowView (globalRowIndex, ind, val);
6207 auto indIter = ind.begin ();
6208 auto valIter = val.begin ();
6209 for (; indIter != ind.end() && valIter != val.end();
6210 ++indIter, ++valIter) {
6211 const GO globalColIndex = *indIter;
6214 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6215 << (globalColIndex + 1 - colIndexBase) <<
" ";
6216 if (STS::isComplex) {
6217 out << STS::real (*valIter) <<
" " << STS::imag (*valIter);
6226 using OTG = Teuchos::OrdinalTraits<GO>;
6227 for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
6228 localRowIndex <= gatherRowMap->getMaxLocalIndex();
6231 const GO globalRowIndex =
6232 gatherRowMap->getGlobalElement (localRowIndex);
6233 TEUCHOS_TEST_FOR_EXCEPTION(
6234 globalRowIndex == OTG::invalid(), std::logic_error,
6235 "Failed to convert the supposed local row index "
6236 << localRowIndex <<
" into a global row index. "
6237 "Please report this bug to the Tpetra developers.");
6238 ArrayView<const LO> ind;
6239 ArrayView<const ST> val;
6240 newMatrix->getLocalRowView (localRowIndex, ind, val);
6241 auto indIter = ind.begin ();
6242 auto valIter = val.begin ();
6243 for (; indIter != ind.end() && valIter != val.end();
6244 ++indIter, ++valIter) {
6246 const GO globalColIndex =
6247 newMatrix->getColMap()->getGlobalElement (*indIter);
6248 TEUCHOS_TEST_FOR_EXCEPTION(
6249 globalColIndex == OTG::invalid(), std::logic_error,
6250 "On local row " << localRowIndex <<
" of the sparse matrix: "
6251 "Failed to convert the supposed local column index "
6252 << *indIter <<
" into a global column index. Please report "
6253 "this bug to the Tpetra developers.");
6256 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6257 << (globalColIndex + 1 - colIndexBase) <<
" ";
6258 if (STS::isComplex) {
6259 out << STS::real (*valIter) <<
" " << STS::imag (*valIter);
6273 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6274 const std::string& matrixName,
6275 const std::string& matrixDescription,
6276 const bool debug=
false)
6278 TEUCHOS_TEST_FOR_EXCEPTION
6279 (pMatrix.is_null (), std::invalid_argument,
6280 "The input matrix is null.");
6281 writeSparse (out, *pMatrix, matrixName, matrixDescription, debug);
6317 const std::string& graphName,
6318 const std::string& graphDescription,
6319 const bool debug=
false)
6321 using Teuchos::ArrayView;
6322 using Teuchos::Comm;
6323 using Teuchos::FancyOStream;
6324 using Teuchos::getFancyOStream;
6325 using Teuchos::null;
6327 using Teuchos::rcpFromRef;
6338 if (rowMap.is_null ()) {
6341 auto comm = rowMap->getComm ();
6342 if (comm.is_null ()) {
6345 const int myRank = comm->getRank ();
6348 RCP<FancyOStream> err =
6349 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6351 std::ostringstream os;
6352 os << myRank <<
": writeSparseGraph" << endl;
6355 os <<
"-- " << myRank <<
": past barrier" << endl;
6360 const bool debugPrint = debug && myRank == 0;
6367 const global_size_t numRows = rangeMap->getGlobalNumElements ();
6368 const global_size_t numCols = domainMap->getGlobalNumElements ();
6370 if (debug && myRank == 0) {
6371 std::ostringstream os;
6372 os <<
"-- Input sparse graph is:"
6373 <<
"---- " << numRows <<
" x " << numCols <<
" with "
6377 <<
" indexed." << endl
6378 <<
"---- Its row Map has " << rowMap->getGlobalNumElements ()
6379 <<
" elements." << endl
6380 <<
"---- Its col Map has " << colMap->getGlobalNumElements ()
6381 <<
" elements." << endl;
6386 const size_t localNumRows = (myRank == 0) ? numRows : 0;
6388 std::ostringstream os;
6389 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6392 auto gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
6400 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6401 auto gatherColMap = Details::computeGatherMap (colMap, err, debug);
6410 static_cast<size_t> (0));
6417 RCP<const map_type> gatherDomainMap =
6418 rcp (
new map_type (numCols, localNumCols,
6419 domainMap->getIndexBase (),
6421 RCP<const map_type> gatherRangeMap =
6422 rcp (
new map_type (numRows, localNumRows,
6423 rangeMap->getIndexBase (),
6425 newGraph.
fillComplete (gatherDomainMap, gatherRangeMap);
6429 cerr <<
"-- Output sparse graph is:"
6430 <<
"---- " << newGraph.
getRangeMap ()->getGlobalNumElements ()
6437 <<
" indexed." << endl
6438 <<
"---- Its row map has "
6439 << newGraph.
getRowMap ()->getGlobalNumElements ()
6440 <<
" elements, with index base "
6441 << newGraph.
getRowMap ()->getIndexBase () <<
"." << endl
6442 <<
"---- Its col map has "
6443 << newGraph.
getColMap ()->getGlobalNumElements ()
6444 <<
" elements, with index base "
6445 << newGraph.
getColMap ()->getIndexBase () <<
"." << endl
6446 <<
"---- Element count of output graph's column Map may differ "
6447 <<
"from that of the input matrix's column Map, if some columns "
6448 <<
"of the matrix contain no entries." << endl;
6462 out <<
"%%MatrixMarket matrix coordinate pattern general" << endl;
6465 if (graphName !=
"") {
6466 printAsComment (out, graphName);
6468 if (graphDescription !=
"") {
6469 printAsComment (out, graphDescription);
6480 out << newGraph.
getRangeMap ()->getGlobalNumElements () <<
" "
6481 << newGraph.
getDomainMap ()->getGlobalNumElements () <<
" "
6487 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6488 const GO colIndexBase = newGraph.
getColMap()->getIndexBase ();
6499 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6500 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6501 for (GO globalRowIndex = minAllGlobalIndex;
6502 globalRowIndex <= maxAllGlobalIndex;
6504 ArrayView<const GO> ind;
6506 for (
auto indIter = ind.begin (); indIter != ind.end (); ++indIter) {
6507 const GO globalColIndex = *indIter;
6510 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6511 << (globalColIndex + 1 - colIndexBase) <<
" ";
6517 typedef Teuchos::OrdinalTraits<GO> OTG;
6518 for (LO localRowIndex = gatherRowMap->getMinLocalIndex ();
6519 localRowIndex <= gatherRowMap->getMaxLocalIndex ();
6522 const GO globalRowIndex =
6523 gatherRowMap->getGlobalElement (localRowIndex);
6524 TEUCHOS_TEST_FOR_EXCEPTION
6525 (globalRowIndex == OTG::invalid (), std::logic_error,
"Failed "
6526 "to convert the supposed local row index " << localRowIndex <<
6527 " into a global row index. Please report this bug to the "
6528 "Tpetra developers.");
6529 ArrayView<const LO> ind;
6531 for (
auto indIter = ind.begin (); indIter != ind.end (); ++indIter) {
6533 const GO globalColIndex =
6534 newGraph.
getColMap ()->getGlobalElement (*indIter);
6535 TEUCHOS_TEST_FOR_EXCEPTION(
6536 globalColIndex == OTG::invalid(), std::logic_error,
6537 "On local row " << localRowIndex <<
" of the sparse graph: "
6538 "Failed to convert the supposed local column index "
6539 << *indIter <<
" into a global column index. Please report "
6540 "this bug to the Tpetra developers.");
6543 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6544 << (globalColIndex + 1 - colIndexBase) <<
" ";
6560 const bool debug=
false)
6602 const std::string& graphName,
6603 const std::string& graphDescription,
6604 const bool debug=
false)
6607 if (comm.is_null ()) {
6615 const int myRank = comm->getRank ();
6620 out.open (filename.c_str ());
6635 const bool debug=
false)
6650 const Teuchos::RCP<const crs_graph_type>& pGraph,
6651 const std::string& graphName,
6652 const std::string& graphDescription,
6653 const bool debug=
false)
6669 const Teuchos::RCP<const crs_graph_type>& pGraph,
6670 const bool debug=
false)
6700 const bool debug=
false)
6708 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6709 const bool debug=
false)
6745 const std::string& matrixName,
6746 const std::string& matrixDescription,
6747 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6748 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6750 const int myRank = X.
getMap ().is_null () ? 0 :
6751 (X.
getMap ()->getComm ().is_null () ? 0 :
6752 X.
getMap ()->getComm ()->getRank ());
6756 out.open (filename.c_str());
6759 writeDense (out, X, matrixName, matrixDescription, err, dbg);
6772 const Teuchos::RCP<const multivector_type>& X,
6773 const std::string& matrixName,
6774 const std::string& matrixDescription,
6775 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6776 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6778 TEUCHOS_TEST_FOR_EXCEPTION(
6779 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
6780 "writeDenseFile: The input MultiVector X is null.");
6781 writeDenseFile (filename, *X, matrixName, matrixDescription, err, dbg);
6792 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6793 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6805 const Teuchos::RCP<const multivector_type>& X,
6806 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6807 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6809 TEUCHOS_TEST_FOR_EXCEPTION(
6810 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
6811 "writeDenseFile: The input MultiVector X is null.");
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 using Teuchos::Comm;
6855 using Teuchos::outArg;
6856 using Teuchos::REDUCE_MAX;
6857 using Teuchos::reduceAll;
6861 RCP<const Comm<int> > comm = X.
getMap ().is_null () ?
6862 Teuchos::null : X.
getMap ()->getComm ();
6863 const int myRank = comm.is_null () ? 0 : comm->getRank ();
6868 const bool debug = ! dbg.is_null ();
6871 std::ostringstream os;
6872 os << myRank <<
": writeDense" << endl;
6877 writeDenseHeader (out, X, matrixName, matrixDescription, err, dbg);
6882 for (
size_t j = 0; j < numVecs; ++j) {
6883 writeDenseColumn (out, * (X.
getVector (j)), err, dbg);
6888 std::ostringstream os;
6889 os << myRank <<
": writeDense: Done" << endl;
6923 writeDenseHeader (std::ostream& out,
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;
6933 using Teuchos::REDUCE_MAX;
6934 using Teuchos::reduceAll;
6936 typedef Teuchos::ScalarTraits<scalar_type> STS;
6937 const char prefix[] =
"Tpetra::MatrixMarket::writeDenseHeader: ";
6939 RCP<const Comm<int> > comm = X.getMap ().is_null () ?
6940 Teuchos::null : X.getMap ()->getComm ();
6941 const int myRank = comm.is_null () ? 0 : comm->getRank ();
6948 const bool debug = ! dbg.is_null ();
6952 std::ostringstream os;
6953 os << myRank <<
": writeDenseHeader" << endl;
6972 std::ostringstream hdr;
6974 std::string dataType;
6975 if (STS::isComplex) {
6976 dataType =
"complex";
6977 }
else if (STS::isOrdinal) {
6978 dataType =
"integer";
6982 hdr <<
"%%MatrixMarket matrix array " << dataType <<
" general"
6987 if (matrixName !=
"") {
6988 printAsComment (hdr, matrixName);
6990 if (matrixDescription !=
"") {
6991 printAsComment (hdr, matrixDescription);
6994 hdr << X.getGlobalLength () <<
" " << X.getNumVectors () << endl;
6998 }
catch (std::exception& e) {
6999 if (! err.is_null ()) {
7000 *err << prefix <<
"While writing the Matrix Market header, "
7001 "Process 0 threw an exception: " << e.what () << endl;
7010 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
7011 TEUCHOS_TEST_FOR_EXCEPTION(
7012 gblErr == 1, std::runtime_error, prefix <<
"Some error occurred "
7013 "which prevented this method from completing.");
7017 *dbg << myRank <<
": writeDenseHeader: Done" << endl;
7040 writeDenseColumn (std::ostream& out,
7042 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7043 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7045 using Teuchos::arcp;
7046 using Teuchos::Array;
7047 using Teuchos::ArrayRCP;
7048 using Teuchos::ArrayView;
7049 using Teuchos::Comm;
7050 using Teuchos::CommRequest;
7051 using Teuchos::ireceive;
7052 using Teuchos::isend;
7053 using Teuchos::outArg;
7054 using Teuchos::REDUCE_MAX;
7055 using Teuchos::reduceAll;
7057 using Teuchos::TypeNameTraits;
7058 using Teuchos::wait;
7060 typedef Teuchos::ScalarTraits<scalar_type> STS;
7062 const Comm<int>& comm = * (X.getMap ()->getComm ());
7063 const int myRank = comm.getRank ();
7064 const int numProcs = comm.getSize ();
7071 const bool debug = ! dbg.is_null ();
7075 std::ostringstream os;
7076 os << myRank <<
": writeDenseColumn" << endl;
7085 Teuchos::SetScientific<scalar_type> sci (out);
7087 const size_t myNumRows = X.getLocalLength ();
7088 const size_t numCols = X.getNumVectors ();
7091 const int sizeTag = 1337;
7092 const int dataTag = 1338;
7153 RCP<CommRequest<int> > sendReqSize, sendReqData;
7159 Array<ArrayRCP<size_t> > recvSizeBufs (3);
7160 Array<ArrayRCP<scalar_type> > recvDataBufs (3);
7161 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7162 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7165 ArrayRCP<size_t> sendDataSize (1);
7166 sendDataSize[0] = myNumRows;
7170 std::ostringstream os;
7171 os << myRank <<
": Post receive-size receives from "
7172 "Procs 1 and 2: tag = " << sizeTag << endl;
7176 recvSizeBufs[0].resize (1);
7180 (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7181 recvSizeBufs[1].resize (1);
7182 (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7183 recvSizeBufs[2].resize (1);
7184 (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7187 ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
7191 ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
7194 else if (myRank == 1 || myRank == 2) {
7196 std::ostringstream os;
7197 os << myRank <<
": Post send-size send: size = "
7198 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7205 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7209 std::ostringstream os;
7210 os << myRank <<
": Not posting my send-size send yet" << endl;
7219 std::ostringstream os;
7220 os << myRank <<
": Pack my entries" << endl;
7223 ArrayRCP<scalar_type> sendDataBuf;
7225 sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
7226 X.get1dCopy (sendDataBuf (), myNumRows);
7228 catch (std::exception& e) {
7230 if (! err.is_null ()) {
7231 std::ostringstream os;
7232 os <<
"Process " << myRank <<
": Attempt to pack my MultiVector "
7233 "entries threw an exception: " << e.what () << endl;
7238 std::ostringstream os;
7239 os << myRank <<
": Done packing my entries" << endl;
7248 *dbg << myRank <<
": Post send-data send: tag = " << dataTag
7251 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7259 std::ostringstream os;
7260 os << myRank <<
": Write my entries" << endl;
7266 const size_t printNumRows = myNumRows;
7267 ArrayView<const scalar_type> printData = sendDataBuf ();
7268 const size_t printStride = printNumRows;
7269 if (static_cast<size_t> (printData.size ()) < printStride * numCols) {
7271 if (! err.is_null ()) {
7272 std::ostringstream os;
7273 os <<
"Process " << myRank <<
": My MultiVector data's size "
7274 << printData.size () <<
" does not match my local dimensions "
7275 << printStride <<
" x " << numCols <<
"." << endl;
7283 for (
size_t col = 0; col < numCols; ++col) {
7284 for (
size_t row = 0; row < printNumRows; ++row) {
7285 if (STS::isComplex) {
7286 out << STS::real (printData[row + col * printStride]) <<
" "
7287 << STS::imag (printData[row + col * printStride]) << endl;
7289 out << printData[row + col * printStride] << endl;
7298 const int recvRank = 1;
7299 const int circBufInd = recvRank % 3;
7301 std::ostringstream os;
7302 os << myRank <<
": Wait on receive-size receive from Process "
7303 << recvRank << endl;
7307 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7311 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7312 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7314 if (! err.is_null ()) {
7315 std::ostringstream os;
7316 os << myRank <<
": Result of receive-size receive from Process "
7317 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7318 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7319 "This should never happen, and suggests that the receive never "
7320 "got posted. Please report this bug to the Tpetra developers."
7335 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7337 std::ostringstream os;
7338 os << myRank <<
": Post receive-data receive from Process "
7339 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7340 << recvDataBufs[circBufInd].size () << endl;
7343 if (! recvSizeReqs[circBufInd].is_null ()) {
7345 if (! err.is_null ()) {
7346 std::ostringstream os;
7347 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7348 "null, before posting the receive-data receive from Process "
7349 << recvRank <<
". This should never happen. Please report "
7350 "this bug to the Tpetra developers." << endl;
7354 recvDataReqs[circBufInd] =
7355 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7356 recvRank, dataTag, comm);
7359 else if (myRank == 1) {
7362 std::ostringstream os;
7363 os << myRank <<
": Wait on my send-size send" << endl;
7366 wait<int> (comm, outArg (sendReqSize));
7372 for (
int p = 1; p < numProcs; ++p) {
7374 if (p + 2 < numProcs) {
7376 const int recvRank = p + 2;
7377 const int circBufInd = recvRank % 3;
7379 std::ostringstream os;
7380 os << myRank <<
": Post receive-size receive from Process "
7381 << recvRank <<
": tag = " << sizeTag << endl;
7384 if (! recvSizeReqs[circBufInd].is_null ()) {
7386 if (! err.is_null ()) {
7387 std::ostringstream os;
7388 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7389 <<
"null, for the receive-size receive from Process "
7390 << recvRank <<
"! This may mean that this process never "
7391 <<
"finished waiting for the receive from Process "
7392 << (recvRank - 3) <<
"." << endl;
7396 recvSizeReqs[circBufInd] =
7397 ireceive<int, size_t> (recvSizeBufs[circBufInd],
7398 recvRank, sizeTag, comm);
7401 if (p + 1 < numProcs) {
7402 const int recvRank = p + 1;
7403 const int circBufInd = recvRank % 3;
7407 std::ostringstream os;
7408 os << myRank <<
": Wait on receive-size receive from Process "
7409 << recvRank << endl;
7412 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7416 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7417 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7419 if (! err.is_null ()) {
7420 std::ostringstream os;
7421 os << myRank <<
": Result of receive-size receive from Process "
7422 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7423 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7424 "This should never happen, and suggests that the receive never "
7425 "got posted. Please report this bug to the Tpetra developers."
7439 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7441 std::ostringstream os;
7442 os << myRank <<
": Post receive-data receive from Process "
7443 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7444 << recvDataBufs[circBufInd].size () << endl;
7447 if (! recvDataReqs[circBufInd].is_null ()) {
7449 if (! err.is_null ()) {
7450 std::ostringstream os;
7451 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
7452 <<
"null, for the receive-data receive from Process "
7453 << recvRank <<
"! This may mean that this process never "
7454 <<
"finished waiting for the receive from Process "
7455 << (recvRank - 3) <<
"." << endl;
7459 recvDataReqs[circBufInd] =
7460 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7461 recvRank, dataTag, comm);
7465 const int recvRank = p;
7466 const int circBufInd = recvRank % 3;
7468 std::ostringstream os;
7469 os << myRank <<
": Wait on receive-data receive from Process "
7470 << recvRank << endl;
7473 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
7480 std::ostringstream os;
7481 os << myRank <<
": Write entries from Process " << recvRank
7483 *dbg << os.str () << endl;
7485 size_t printNumRows = (recvSizeBufs[circBufInd])[0];
7486 if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7488 if (! err.is_null ()) {
7489 std::ostringstream os;
7490 os << myRank <<
": Result of receive-size receive from Process "
7491 << recvRank <<
" was Teuchos::OrdinalTraits<size_t>::"
7492 "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
7493 <<
". This should never happen, and suggests that its "
7494 "receive-size receive was never posted. "
7495 "Please report this bug to the Tpetra developers." << endl;
7502 if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
7504 if (! err.is_null ()) {
7505 std::ostringstream os;
7506 os << myRank <<
": Result of receive-size receive from Proc "
7507 << recvRank <<
" was " << printNumRows <<
" > 0, but "
7508 "recvDataBufs[" << circBufInd <<
"] is null. This should "
7509 "never happen. Please report this bug to the Tpetra "
7510 "developers." << endl;
7520 ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
7521 const size_t printStride = printNumRows;
7525 for (
size_t col = 0; col < numCols; ++col) {
7526 for (
size_t row = 0; row < printNumRows; ++row) {
7527 if (STS::isComplex) {
7528 out << STS::real (printData[row + col * printStride]) <<
" "
7529 << STS::imag (printData[row + col * printStride]) << endl;
7531 out << printData[row + col * printStride] << endl;
7536 else if (myRank == p) {
7539 std::ostringstream os;
7540 os << myRank <<
": Wait on my send-data send" << endl;
7543 wait<int> (comm, outArg (sendReqData));
7545 else if (myRank == p + 1) {
7548 std::ostringstream os;
7549 os << myRank <<
": Post send-data send: tag = " << dataTag
7553 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7556 std::ostringstream os;
7557 os << myRank <<
": Wait on my send-size send" << endl;
7560 wait<int> (comm, outArg (sendReqSize));
7562 else if (myRank == p + 2) {
7565 std::ostringstream os;
7566 os << myRank <<
": Post send-size send: size = "
7567 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7570 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7575 reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
7576 TEUCHOS_TEST_FOR_EXCEPTION(
7577 gblErr == 1, std::runtime_error,
"Tpetra::MatrixMarket::writeDense "
7578 "experienced some kind of error and was unable to complete.");
7582 *dbg << myRank <<
": writeDenseColumn: Done" << endl;
7596 const Teuchos::RCP<const multivector_type>& X,
7597 const std::string& matrixName,
7598 const std::string& matrixDescription,
7599 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7600 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7602 TEUCHOS_TEST_FOR_EXCEPTION(
7603 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7604 "writeDense: The input MultiVector X is null.");
7605 writeDense (out, *X, matrixName, matrixDescription, err, dbg);
7616 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7617 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7629 const Teuchos::RCP<const multivector_type>& X,
7630 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7631 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7633 TEUCHOS_TEST_FOR_EXCEPTION(
7634 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7635 "writeDense: The input MultiVector X is null.");
7661 Teuchos::RCP<Teuchos::FancyOStream> err =
7662 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
7677 const Teuchos::RCP<Teuchos::FancyOStream>& err,
7678 const bool debug=
false)
7680 using Teuchos::Array;
7681 using Teuchos::ArrayRCP;
7682 using Teuchos::ArrayView;
7683 using Teuchos::Comm;
7684 using Teuchos::CommRequest;
7685 using Teuchos::ireceive;
7686 using Teuchos::isend;
7688 using Teuchos::TypeNameTraits;
7689 using Teuchos::wait;
7692 typedef int pid_type;
7707 typedef ptrdiff_t int_type;
7708 TEUCHOS_TEST_FOR_EXCEPTION(
7709 sizeof (GO) >
sizeof (int_type), std::logic_error,
7710 "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
7711 <<
" is too big for ptrdiff_t. sizeof(GO) = " <<
sizeof (GO)
7712 <<
" > sizeof(ptrdiff_t) = " <<
sizeof (ptrdiff_t) <<
".");
7713 TEUCHOS_TEST_FOR_EXCEPTION(
7714 sizeof (pid_type) >
sizeof (int_type), std::logic_error,
7715 "The (MPI) process rank type pid_type=" <<
7716 TypeNameTraits<pid_type>::name () <<
" is too big for ptrdiff_t. "
7717 "sizeof(pid_type) = " <<
sizeof (pid_type) <<
" > sizeof(ptrdiff_t)"
7718 " = " <<
sizeof (ptrdiff_t) <<
".");
7720 const Comm<int>& comm = * (map.
getComm ());
7721 const int myRank = comm.getRank ();
7722 const int numProcs = comm.getSize ();
7724 if (! err.is_null ()) {
7728 std::ostringstream os;
7729 os << myRank <<
": writeMap" << endl;
7732 if (! err.is_null ()) {
7739 const int sizeTag = 1337;
7740 const int dataTag = 1338;
7799 RCP<CommRequest<int> > sendReqSize, sendReqData;
7805 Array<ArrayRCP<int_type> > recvSizeBufs (3);
7806 Array<ArrayRCP<int_type> > recvDataBufs (3);
7807 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7808 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7811 ArrayRCP<int_type> sendDataSize (1);
7812 sendDataSize[0] = myNumRows;
7816 std::ostringstream os;
7817 os << myRank <<
": Post receive-size receives from "
7818 "Procs 1 and 2: tag = " << sizeTag << endl;
7822 recvSizeBufs[0].resize (1);
7823 (recvSizeBufs[0])[0] = -1;
7824 recvSizeBufs[1].resize (1);
7825 (recvSizeBufs[1])[0] = -1;
7826 recvSizeBufs[2].resize (1);
7827 (recvSizeBufs[2])[0] = -1;
7830 ireceive<int, int_type> (recvSizeBufs[1], 1, sizeTag, comm);
7834 ireceive<int, int_type> (recvSizeBufs[2], 2, sizeTag, comm);
7837 else if (myRank == 1 || myRank == 2) {
7839 std::ostringstream os;
7840 os << myRank <<
": Post send-size send: size = "
7841 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7848 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
7852 std::ostringstream os;
7853 os << myRank <<
": Not posting my send-size send yet" << endl;
7864 std::ostringstream os;
7865 os << myRank <<
": Pack my GIDs and PIDs" << endl;
7869 ArrayRCP<int_type> sendDataBuf (myNumRows * 2);
7872 const int_type myMinGblIdx =
7874 for (
size_t k = 0; k < myNumRows; ++k) {
7875 const int_type gid = myMinGblIdx +
static_cast<int_type
> (k);
7876 const int_type pid =
static_cast<int_type
> (myRank);
7877 sendDataBuf[2*k] = gid;
7878 sendDataBuf[2*k+1] = pid;
7883 for (
size_t k = 0; k < myNumRows; ++k) {
7884 const int_type gid =
static_cast<int_type
> (myGblInds[k]);
7885 const int_type pid =
static_cast<int_type
> (myRank);
7886 sendDataBuf[2*k] = gid;
7887 sendDataBuf[2*k+1] = pid;
7892 std::ostringstream os;
7893 os << myRank <<
": Done packing my GIDs and PIDs" << endl;
7900 *err << myRank <<
": Post send-data send: tag = " << dataTag
7903 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
7908 *err << myRank <<
": Write MatrixMarket header" << endl;
7913 std::ostringstream hdr;
7917 hdr <<
"%%MatrixMarket matrix array integer general" << endl
7918 <<
"% Format: Version 2.0" << endl
7920 <<
"% This file encodes a Tpetra::Map." << endl
7921 <<
"% It is stored as a dense vector, with twice as many " << endl
7922 <<
"% entries as the global number of GIDs (global indices)." << endl
7923 <<
"% (GID, PID) pairs are stored contiguously, where the PID " << endl
7924 <<
"% is the rank of the process owning that GID." << endl
7929 std::ostringstream os;
7930 os << myRank <<
": Write my GIDs and PIDs" << endl;
7936 const int_type printNumRows = myNumRows;
7937 ArrayView<const int_type> printData = sendDataBuf ();
7938 for (int_type k = 0; k < printNumRows; ++k) {
7939 const int_type gid = printData[2*k];
7940 const int_type pid = printData[2*k+1];
7941 out << gid << endl << pid << endl;
7947 const int recvRank = 1;
7948 const int circBufInd = recvRank % 3;
7950 std::ostringstream os;
7951 os << myRank <<
": Wait on receive-size receive from Process "
7952 << recvRank << endl;
7956 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7960 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
7961 if (debug && recvNumRows == -1) {
7962 std::ostringstream os;
7963 os << myRank <<
": Result of receive-size receive from Process "
7964 << recvRank <<
" is -1. This should never happen, and "
7965 "suggests that the receive never got posted. Please report "
7966 "this bug to the Tpetra developers." << endl;
7971 recvDataBufs[circBufInd].resize (recvNumRows * 2);
7973 std::ostringstream os;
7974 os << myRank <<
": Post receive-data receive from Process "
7975 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7976 << recvDataBufs[circBufInd].size () << endl;
7979 if (! recvSizeReqs[circBufInd].is_null ()) {
7980 std::ostringstream os;
7981 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7982 "null, before posting the receive-data receive from Process "
7983 << recvRank <<
". This should never happen. Please report "
7984 "this bug to the Tpetra developers." << endl;
7987 recvDataReqs[circBufInd] =
7988 ireceive<int, int_type> (recvDataBufs[circBufInd],
7989 recvRank, dataTag, comm);
7992 else if (myRank == 1) {
7995 std::ostringstream os;
7996 os << myRank <<
": Wait on my send-size send" << endl;
7999 wait<int> (comm, outArg (sendReqSize));
8005 for (
int p = 1; p < numProcs; ++p) {
8007 if (p + 2 < numProcs) {
8009 const int recvRank = p + 2;
8010 const int circBufInd = recvRank % 3;
8012 std::ostringstream os;
8013 os << myRank <<
": Post receive-size receive from Process "
8014 << recvRank <<
": tag = " << sizeTag << endl;
8017 if (! recvSizeReqs[circBufInd].is_null ()) {
8018 std::ostringstream os;
8019 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
8020 <<
"null, for the receive-size receive from Process "
8021 << recvRank <<
"! This may mean that this process never "
8022 <<
"finished waiting for the receive from Process "
8023 << (recvRank - 3) <<
"." << endl;
8026 recvSizeReqs[circBufInd] =
8027 ireceive<int, int_type> (recvSizeBufs[circBufInd],
8028 recvRank, sizeTag, comm);
8031 if (p + 1 < numProcs) {
8032 const int recvRank = p + 1;
8033 const int circBufInd = recvRank % 3;
8037 std::ostringstream os;
8038 os << myRank <<
": Wait on receive-size receive from Process "
8039 << recvRank << endl;
8042 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
8046 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
8047 if (debug && recvNumRows == -1) {
8048 std::ostringstream os;
8049 os << myRank <<
": Result of receive-size receive from Process "
8050 << recvRank <<
" is -1. This should never happen, and "
8051 "suggests that the receive never got posted. Please report "
8052 "this bug to the Tpetra developers." << endl;
8057 recvDataBufs[circBufInd].resize (recvNumRows * 2);
8059 std::ostringstream os;
8060 os << myRank <<
": Post receive-data receive from Process "
8061 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
8062 << recvDataBufs[circBufInd].size () << endl;
8065 if (! recvDataReqs[circBufInd].is_null ()) {
8066 std::ostringstream os;
8067 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
8068 <<
"null, for the receive-data receive from Process "
8069 << recvRank <<
"! This may mean that this process never "
8070 <<
"finished waiting for the receive from Process "
8071 << (recvRank - 3) <<
"." << endl;
8074 recvDataReqs[circBufInd] =
8075 ireceive<int, int_type> (recvDataBufs[circBufInd],
8076 recvRank, dataTag, comm);
8080 const int recvRank = p;
8081 const int circBufInd = recvRank % 3;
8083 std::ostringstream os;
8084 os << myRank <<
": Wait on receive-data receive from Process "
8085 << recvRank << endl;
8088 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
8095 std::ostringstream os;
8096 os << myRank <<
": Write GIDs and PIDs from Process "
8097 << recvRank << endl;
8098 *err << os.str () << endl;
8100 const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
8101 if (debug && printNumRows == -1) {
8102 std::ostringstream os;
8103 os << myRank <<
": Result of receive-size receive from Process "
8104 << recvRank <<
" was -1. This should never happen, and "
8105 "suggests that its receive-size receive was never posted. "
8106 "Please report this bug to the Tpetra developers." << endl;
8109 if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
8110 std::ostringstream os;
8111 os << myRank <<
": Result of receive-size receive from Proc "
8112 << recvRank <<
" was " << printNumRows <<
" > 0, but "
8113 "recvDataBufs[" << circBufInd <<
"] is null. This should "
8114 "never happen. Please report this bug to the Tpetra "
8115 "developers." << endl;
8118 ArrayView<const int_type> printData = (recvDataBufs[circBufInd]) ();
8119 for (int_type k = 0; k < printNumRows; ++k) {
8120 const int_type gid = printData[2*k];
8121 const int_type pid = printData[2*k+1];
8122 out << gid << endl << pid << endl;
8125 else if (myRank == p) {
8128 std::ostringstream os;
8129 os << myRank <<
": Wait on my send-data send" << endl;
8132 wait<int> (comm, outArg (sendReqData));
8134 else if (myRank == p + 1) {
8137 std::ostringstream os;
8138 os << myRank <<
": Post send-data send: tag = " << dataTag
8142 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
8145 std::ostringstream os;
8146 os << myRank <<
": Wait on my send-size send" << endl;
8149 wait<int> (comm, outArg (sendReqSize));
8151 else if (myRank == p + 2) {
8154 std::ostringstream os;
8155 os << myRank <<
": Post send-size send: size = "
8156 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
8159 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
8163 if (! err.is_null ()) {
8167 *err << myRank <<
": writeMap: Done" << endl;
8169 if (! err.is_null ()) {
8179 const int myRank = map.
getComm ()->getRank ();
8182 out.open (filename.c_str());
8215 printAsComment (std::ostream& out,
const std::string& str)
8218 std::istringstream inpstream (str);
8221 while (getline (inpstream, line)) {
8222 if (! line.empty()) {
8225 if (line[0] ==
'%') {
8226 out << line << endl;
8229 out <<
"%% " << line << endl;
8257 Teuchos::ParameterList pl;
8283 Teuchos::ParameterList pl;
8326 const Teuchos::ParameterList& params)
8329 std::string tmpFile =
"__TMP__" + fileName;
8330 const int myRank = A.
getDomainMap()->getComm()->getRank();
8331 bool precisionChanged=
false;
8341 if (std::ifstream(tmpFile))
8342 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
8343 "writeOperator: temporary file " << tmpFile <<
" already exists");
8344 out.open(tmpFile.c_str());
8345 if (params.isParameter(
"precision")) {
8346 oldPrecision = out.precision(params.get<
int>(
"precision"));
8347 precisionChanged=
true;
8351 const std::string header = writeOperatorImpl(out, A, params);
8354 if (precisionChanged)
8355 out.precision(oldPrecision);
8357 out.open(fileName.c_str(), std::ios::binary);
8358 bool printMatrixMarketHeader =
true;
8359 if (params.isParameter(
"print MatrixMarket header"))
8360 printMatrixMarketHeader = params.get<
bool>(
"print MatrixMarket header");
8361 if (printMatrixMarketHeader && myRank == 0) {
8366 std::ifstream src(tmpFile, std::ios_base::binary);
8370 remove(tmpFile.c_str());
8415 const Teuchos::ParameterList& params)
8417 const int myRank = A.
getDomainMap ()->getComm ()->getRank ();
8434 std::ostringstream tmpOut;
8436 if (params.isParameter (
"precision") && params.isType<
int> (
"precision")) {
8437 (void) tmpOut.precision (params.get<
int> (
"precision"));
8441 const std::string header = writeOperatorImpl (tmpOut, A, params);
8444 bool printMatrixMarketHeader =
true;
8445 if (params.isParameter (
"print MatrixMarket header") &&
8446 params.isType<
bool> (
"print MatrixMarket header")) {
8447 printMatrixMarketHeader = params.get<
bool> (
"print MatrixMarket header");
8449 if (printMatrixMarketHeader && myRank == 0) {
8465 out << tmpOut.str ();
8479 writeOperatorImpl (std::ostream& os,
8480 const operator_type& A,
8481 const Teuchos::ParameterList& params)
8485 using Teuchos::ArrayRCP;
8486 using Teuchos::Array;
8490 typedef scalar_type Scalar;
8491 typedef Teuchos::OrdinalTraits<LO> TLOT;
8492 typedef Teuchos::OrdinalTraits<GO> TGOT;
8496 const map_type& domainMap = *(A.getDomainMap());
8497 RCP<const map_type> rangeMap = A.getRangeMap();
8498 RCP<const Teuchos::Comm<int> > comm = rangeMap->getComm();
8499 const int myRank = comm->getRank();
8500 const size_t numProcs = comm->getSize();
8503 if (params.isParameter(
"probing size"))
8504 numMVs = params.get<
int>(
"probing size");
8507 GO minColGid = domainMap.getMinAllGlobalIndex();
8508 GO maxColGid = domainMap.getMaxAllGlobalIndex();
8513 GO numGlobElts = maxColGid - minColGid + TGOT::one();
8514 GO numChunks = numGlobElts / numMVs;
8515 GO rem = numGlobElts % numMVs;
8516 GO indexBase = rangeMap->getIndexBase();
8518 int offsetToUseInPrinting = 1 - indexBase;
8519 if (params.isParameter(
"zero-based indexing")) {
8520 if (params.get<
bool>(
"zero-based indexing") ==
true)
8521 offsetToUseInPrinting = -indexBase;
8525 size_t numLocalRangeEntries = rangeMap->getNodeNumElements();
8528 RCP<const map_type> allGidsMap = rcp(
new map_type(TGOT::invalid(), numLocalRangeEntries,
8531 mv_type_go allGids(allGidsMap,1);
8532 Teuchos::ArrayRCP<GO> allGidsData = allGids.getDataNonConst(0);
8534 for (
size_t i=0; i<numLocalRangeEntries; i++)
8535 allGidsData[i] = rangeMap->getGlobalElement(i);
8536 allGidsData = Teuchos::null;
8539 GO numTargetMapEntries=TGOT::zero();
8540 Teuchos::Array<GO> importGidList;
8542 numTargetMapEntries = rangeMap->getGlobalNumElements();
8543 importGidList.reserve(numTargetMapEntries);
8544 for (GO j=0; j<numTargetMapEntries; ++j) importGidList.push_back(j + indexBase);
8546 importGidList.reserve(numTargetMapEntries);
8548 RCP<map_type> importGidMap = rcp(
new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
8551 import_type gidImporter(allGidsMap, importGidMap);
8552 mv_type_go importedGids(importGidMap, 1);
8553 importedGids.doImport(allGids, gidImporter,
INSERT);
8556 ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
8557 RCP<const map_type> importMap = rcp(
new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm) );
8560 import_type importer(rangeMap, importMap);
8562 RCP<mv_type> colsOnPid0 = rcp(
new mv_type(importMap,numMVs));
8564 RCP<mv_type> ei = rcp(
new mv_type(A.getDomainMap(),numMVs));
8565 RCP<mv_type> colsA = rcp(
new mv_type(A.getRangeMap(),numMVs));
8567 Array<GO> globalColsArray, localColsArray;
8568 globalColsArray.reserve(numMVs);
8569 localColsArray.reserve(numMVs);
8571 ArrayRCP<ArrayRCP<Scalar> > eiData(numMVs);
8572 for (
size_t i=0; i<numMVs; ++i)
8573 eiData[i] = ei->getDataNonConst(i);
8578 for (GO k=0; k<numChunks; ++k) {
8579 for (
size_t j=0; j<numMVs; ++j ) {
8581 GO curGlobalCol = minColGid + k*numMVs + j;
8582 globalColsArray.push_back(curGlobalCol);
8584 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8585 if (curLocalCol != TLOT::invalid()) {
8586 eiData[j][curLocalCol] = TGOT::one();
8587 localColsArray.push_back(curLocalCol);
8593 A.apply(*ei,*colsA);
8595 colsOnPid0->doImport(*colsA,importer,
INSERT);
8598 globalNnz += writeColumns(os,*colsOnPid0, numMVs, importedGidsData(),
8599 globalColsArray, offsetToUseInPrinting);
8602 for (
size_t j=0; j<numMVs; ++j ) {
8603 for (
int i=0; i<localColsArray.size(); ++i)
8604 eiData[j][localColsArray[i]] = TGOT::zero();
8606 globalColsArray.clear();
8607 localColsArray.clear();
8615 for (
int j=0; j<rem; ++j ) {
8616 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8617 globalColsArray.push_back(curGlobalCol);
8619 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8620 if (curLocalCol != TLOT::invalid()) {
8621 eiData[j][curLocalCol] = TGOT::one();
8622 localColsArray.push_back(curLocalCol);
8628 A.apply(*ei,*colsA);
8630 colsOnPid0->doImport(*colsA,importer,
INSERT);
8632 globalNnz += writeColumns(os,*colsOnPid0, rem, importedGidsData(),
8633 globalColsArray, offsetToUseInPrinting);
8636 for (
int j=0; j<rem; ++j ) {
8637 for (
int i=0; i<localColsArray.size(); ++i)
8638 eiData[j][localColsArray[i]] = TGOT::zero();
8640 globalColsArray.clear();
8641 localColsArray.clear();
8650 std::ostringstream oss;
8652 oss <<
"%%MatrixMarket matrix coordinate ";
8653 if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
8658 oss <<
" general" << std::endl;
8659 oss <<
"% Tpetra::Operator" << std::endl;
8660 std::time_t now = std::time(NULL);
8661 oss <<
"% time stamp: " << ctime(&now);
8662 oss <<
"% written from " << numProcs <<
" processes" << std::endl;
8663 size_t numRows = rangeMap->getGlobalNumElements();
8664 size_t numCols = domainMap.getGlobalNumElements();
8665 oss << numRows <<
" " << numCols <<
" " << globalNnz << std::endl;
8672 writeColumns(std::ostream& os, mv_type
const &colsA,
size_t const &numCols,
8673 Teuchos::ArrayView<const global_ordinal_type>
const &rowGids,
8674 Teuchos::Array<global_ordinal_type>
const &colsArray,
8678 typedef scalar_type Scalar;
8679 typedef Teuchos::ScalarTraits<Scalar> STS;
8682 const Scalar zero = STS::zero();
8683 const size_t numRows = colsA.getGlobalLength();
8684 for (
size_t j=0; j<numCols; ++j) {
8685 Teuchos::ArrayRCP<const Scalar>
const curCol = colsA.getData(j);
8686 const GO J = colsArray[j];
8687 for (
size_t i=0; i<numRows; ++i) {
8688 const Scalar val = curCol[i];
8690 os << rowGids[i]+indexBase <<
" " << J+indexBase <<
" " << val << std::endl;
8707 #endif // __MatrixMarket_Tpetra_hpp
static Teuchos::RCP< vector_type > readVectorFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read a Vector from the given Matrix Market file.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Teuchos::RCP< const map_type > getRowMap() const override
Returns the Map that describes the row distribution in this graph.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market input stream.
static void writeSparseGraphFile(const std::string &filename, const Teuchos::RCP< const crs_graph_type > &pGraph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), with no comments...
static void writeOperator(std::ostream &out, const operator_type &A)
Write a Tpetra::Operator to an output stream.
GlobalOrdinal getMinGlobalIndex() const
The minimum global index owned by the calling process.
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.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
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.
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.
void getLocalRowView(const local_ordinal_type lclRow, Teuchos::ArrayView< const local_ordinal_type > &lclColInds) const override
Get a const, non-persisting view of the given local row's local column indices, as a Teuchos::ArrayVi...
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
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< const map_type > readMapFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given Matrix Market file.
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 void writeDenseFile(const std::string &filename, const Teuchos::RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market file.
SparseMatrixType::global_ordinal_type global_ordinal_type
Type of indices as read from the Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream, with provided Maps. ...
static Teuchos::RCP< multivector_type > readDense(std::istream &in, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market input stream.
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 void writeSparseGraph(std::ostream &out, const crs_graph_type &graph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given output stream, with no comments...
size_t global_size_t
Global size_t object.
static Teuchos::RCP< vector_type > readVector(std::istream &in, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream.
SparseMatrixType::node_type node_type
The Kokkos Node type; fourth template parameter of Tpetra::CrsMatrix.
static void writeOperator(const std::string &fileName, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to a file, with options.
static Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream.
Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Specialization of Tpetra::Map that matches SparseMatrixType.
Insert new values that don't currently exist.
SparseMatrixType::scalar_type scalar_type
Type of the entries of the sparse matrix.
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...
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 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.
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.
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 Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream, with optional debugging output stream...
static Teuchos::RCP< multivector_type > readDenseFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market file.
static void writeSparseGraphFile(const std::string &filename, const crs_graph_type &graph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), with no comments...
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.
void getGlobalRowView(const global_ordinal_type gblRow, Teuchos::ArrayView< const global_ordinal_type > &gblColInds) const override
Get a const, non-persisting view of the given global row's global column indices, as a Teuchos::Array...
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.
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).
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a NONOWNING view of the global indices owned by this 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.
size_t getNodeNumElements() const
The number of elements belonging to the calling process.
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.
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.
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.