176 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
191 typedef typename SparseMatrixType::global_ordinal_type
213 typedef Teuchos::Comm<int> comm_type;
223 typedef Teuchos::ArrayRCP<int>::size_type size_type;
235 static Teuchos::RCP<const map_type>
236 makeRangeMap (
const Teuchos::RCP<const comm_type>&
pComm,
242 pComm, GloballyDistributed));
272 static Teuchos::RCP<const map_type>
273 makeRowMap (
const Teuchos::RCP<const map_type>&
pRowMap,
274 const Teuchos::RCP<const comm_type>&
pComm,
282 pComm, GloballyDistributed));
287 std::invalid_argument,
"The specified row map is not distributed, "
288 "but the given communicator includes more than one process (in "
289 "fact, there are " <<
pComm->getSize () <<
" processes).");
292 "The specified row Map's communicator (pRowMap->getComm()) "
293 "differs from the given separately supplied communicator pComm.");
312 static Teuchos::RCP<const map_type>
313 makeDomainMap (
const Teuchos::RCP<const map_type>& pRangeMap,
322 if (numRows == numCols) {
325 return createUniformContigMapWithNode<LO,GO,NT> (numCols,
326 pRangeMap->getComm ());
403 distribute (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
404 Teuchos::ArrayRCP<size_t>& myRowPtr,
405 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
406 Teuchos::ArrayRCP<scalar_type>& myValues,
407 const Teuchos::RCP<const map_type>& pRowMap,
408 Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
409 Teuchos::ArrayRCP<size_t>& rowPtr,
410 Teuchos::ArrayRCP<global_ordinal_type>& colInd,
411 Teuchos::ArrayRCP<scalar_type>& values,
412 const bool debug=
false)
415 using Teuchos::ArrayRCP;
416 using Teuchos::ArrayView;
419 using Teuchos::CommRequest;
422 using Teuchos::receive;
427 const bool extraDebug =
false;
428 RCP<const comm_type> pComm = pRowMap->getComm ();
429 const int numProcs = pComm->getSize ();
430 const int myRank = pComm->getRank ();
431 const int rootRank = 0;
438 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
439 const size_type myNumRows = myRows.size();
440 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<size_t>(myNumRows) !=
441 pRowMap->getLocalNumElements(),
443 "pRowMap->getLocalElementList().size() = "
445 <<
" != pRowMap->getLocalNumElements() = "
446 << pRowMap->getLocalNumElements() <<
". "
447 "Please report this bug to the Tpetra developers.");
448 TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
450 "On Proc 0: numEntriesPerRow.size() = "
451 << numEntriesPerRow.size()
452 <<
" != pRowMap->getLocalElementList().size() = "
453 << myNumRows <<
". Please report this bug to the "
454 "Tpetra developers.");
458 myNumEntriesPerRow = arcp<size_t> (myNumRows);
460 if (myRank != rootRank) {
464 send (*pComm, myNumRows, rootRank);
465 if (myNumRows != 0) {
469 send (*pComm,
static_cast<int> (myNumRows),
470 myRows.getRawPtr(), rootRank);
480 receive (*pComm, rootRank,
481 static_cast<int> (myNumRows),
482 myNumEntriesPerRow.getRawPtr());
487 std::accumulate (myNumEntriesPerRow.begin(),
488 myNumEntriesPerRow.end(), 0);
494 myColInd = arcp<GO> (myNumEntries);
495 myValues = arcp<scalar_type> (myNumEntries);
496 if (myNumEntries > 0) {
499 receive (*pComm, rootRank,
500 static_cast<int> (myNumEntries),
501 myColInd.getRawPtr());
502 receive (*pComm, rootRank,
503 static_cast<int> (myNumEntries),
504 myValues.getRawPtr());
510 cerr <<
"-- Proc 0: Copying my data from global arrays" << endl;
514 for (size_type k = 0; k < myNumRows; ++k) {
515 const GO myCurRow = myRows[k];
517 myNumEntriesPerRow[k] = numEntriesInThisRow;
519 if (extraDebug && debug) {
520 cerr <<
"Proc " << pRowMap->getComm ()->getRank ()
521 <<
": myNumEntriesPerRow[0.." << (myNumRows-1) <<
"] = [";
522 for (size_type k = 0; k < myNumRows; ++k) {
523 cerr << myNumEntriesPerRow[k];
524 if (k < myNumRows-1) {
532 std::accumulate (myNumEntriesPerRow.begin(),
533 myNumEntriesPerRow.end(), 0);
535 cerr <<
"-- Proc 0: I own " << myNumRows <<
" rows and "
536 << myNumEntries <<
" entries" << endl;
538 myColInd = arcp<GO> (myNumEntries);
539 myValues = arcp<scalar_type> (myNumEntries);
547 for (size_type k = 0; k < myNumRows;
548 myCurPos += myNumEntriesPerRow[k], ++k) {
550 const GO myRow = myRows[k];
551 const size_t curPos = rowPtr[myRow];
554 if (curNumEntries > 0) {
555 ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
556 ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
557 std::copy (colIndView.begin(), colIndView.end(),
558 myColIndView.begin());
560 ArrayView<scalar_type> valuesView =
561 values (curPos, curNumEntries);
562 ArrayView<scalar_type> myValuesView =
563 myValues (myCurPos, curNumEntries);
564 std::copy (valuesView.begin(), valuesView.end(),
565 myValuesView.begin());
570 for (
int p = 1; p < numProcs; ++p) {
572 cerr <<
"-- Proc 0: Processing proc " << p << endl;
575 size_type theirNumRows = 0;
580 receive (*pComm, p, &theirNumRows);
582 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
583 << theirNumRows <<
" rows" << endl;
585 if (theirNumRows != 0) {
590 ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
591 receive (*pComm, p, as<int> (theirNumRows),
592 theirRows.getRawPtr ());
601 const global_size_t numRows = pRowMap->getGlobalNumElements ();
602 const GO indexBase = pRowMap->getIndexBase ();
603 bool theirRowsValid =
true;
604 for (size_type k = 0; k < theirNumRows; ++k) {
605 if (theirRows[k] < indexBase ||
606 as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
607 theirRowsValid =
false;
610 if (! theirRowsValid) {
611 TEUCHOS_TEST_FOR_EXCEPTION(
612 ! theirRowsValid, std::logic_error,
613 "Proc " << p <<
" has at least one invalid row index. "
614 "Here are all of them: " <<
615 Teuchos::toString (theirRows ()) <<
". Valid row index "
616 "range (zero-based): [0, " << (numRows - 1) <<
"].");
631 ArrayRCP<size_t> theirNumEntriesPerRow;
632 theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
633 for (size_type k = 0; k < theirNumRows; ++k) {
634 theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
641 send (*pComm,
static_cast<int> (theirNumRows),
642 theirNumEntriesPerRow.getRawPtr(), p);
646 std::accumulate (theirNumEntriesPerRow.begin(),
647 theirNumEntriesPerRow.end(), 0);
650 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
651 << theirNumEntries <<
" entries" << endl;
656 if (theirNumEntries == 0) {
665 ArrayRCP<GO> theirColInd (theirNumEntries);
666 ArrayRCP<scalar_type> theirValues (theirNumEntries);
674 for (size_type k = 0; k < theirNumRows;
675 theirCurPos += theirNumEntriesPerRow[k], k++) {
677 const GO theirRow = theirRows[k];
683 if (curNumEntries > 0) {
684 ArrayView<GO> colIndView =
685 colInd (curPos, curNumEntries);
686 ArrayView<GO> theirColIndView =
687 theirColInd (theirCurPos, curNumEntries);
688 std::copy (colIndView.begin(), colIndView.end(),
689 theirColIndView.begin());
691 ArrayView<scalar_type> valuesView =
692 values (curPos, curNumEntries);
693 ArrayView<scalar_type> theirValuesView =
694 theirValues (theirCurPos, curNumEntries);
695 std::copy (valuesView.begin(), valuesView.end(),
696 theirValuesView.begin());
703 send (*pComm,
static_cast<int> (theirNumEntries),
704 theirColInd.getRawPtr(), p);
705 send (*pComm,
static_cast<int> (theirNumEntries),
706 theirValues.getRawPtr(), p);
709 cerr <<
"-- Proc 0: Finished with proc " << p << endl;
717 numEntriesPerRow = null;
722 if (debug && myRank == 0) {
723 cerr <<
"-- Proc 0: About to fill in myRowPtr" << endl;
731 myRowPtr = arcp<size_t> (myNumRows+1);
733 for (size_type k = 1; k < myNumRows+1; ++k) {
734 myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
736 if (extraDebug && debug) {
737 cerr <<
"Proc " << Teuchos::rank (*(pRowMap->getComm()))
738 <<
": myRowPtr[0.." << myNumRows <<
"] = [";
739 for (size_type k = 0; k < myNumRows+1; ++k) {
745 cerr <<
"]" << endl << endl;
748 if (debug && myRank == 0) {
749 cerr <<
"-- Proc 0: Done with distribute" << endl;
766 static Teuchos::RCP<sparse_matrix_type>
767 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
768 Teuchos::ArrayRCP<size_t>& myRowPtr,
769 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
770 Teuchos::ArrayRCP<scalar_type>& myValues,
771 const Teuchos::RCP<const map_type>& pRowMap,
772 const Teuchos::RCP<const map_type>& pRangeMap,
773 const Teuchos::RCP<const map_type>& pDomainMap,
774 const bool callFillComplete =
true)
776 using Teuchos::ArrayView;
790 TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
791 "makeMatrix: myRowPtr array is null. "
792 "Please report this bug to the Tpetra developers.");
793 TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
794 "makeMatrix: domain map is null. "
795 "Please report this bug to the Tpetra developers.");
796 TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
797 "makeMatrix: range map is null. "
798 "Please report this bug to the Tpetra developers.");
799 TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
800 "makeMatrix: row map is null. "
801 "Please report this bug to the Tpetra developers.");
805 RCP<sparse_matrix_type> A =
810 ArrayView<const GO> myRows = pRowMap->getLocalElementList ();
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 =
900 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
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->getLocalElementList ();
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) {
1050 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1051 "Matrix Market banner line contains syntax error(s): "
1055 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() !=
"matrix",
1056 std::invalid_argument,
"The Matrix Market file does not contain "
1057 "matrix data. Its Banner (first) line says that its object type is \""
1058 << pBanner->matrixType() <<
"\", rather than the required \"matrix\".");
1062 TEUCHOS_TEST_FOR_EXCEPTION(
1063 ! STS::isComplex && pBanner->dataType() ==
"complex",
1064 std::invalid_argument,
1065 "The Matrix Market file contains complex-valued data, but you are "
1066 "trying to read it into a matrix containing entries of the real-"
1067 "valued Scalar type \""
1068 << Teuchos::TypeNameTraits<scalar_type>::name() <<
"\".");
1069 TEUCHOS_TEST_FOR_EXCEPTION(
1071 pBanner->dataType() !=
"real" &&
1072 pBanner->dataType() !=
"complex" &&
1073 pBanner->dataType() !=
"integer",
1074 std::invalid_argument,
1075 "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
1076 "Matrix Market file may not contain a \"pattern\" matrix. A "
1077 "pattern matrix is really just a graph with no weights. It "
1078 "should be stored in a CrsGraph, not a CrsMatrix.");
1080 TEUCHOS_TEST_FOR_EXCEPTION(
1082 pBanner->dataType() !=
"pattern",
1083 std::invalid_argument,
1084 "When reading Matrix Market data into a Tpetra::CrsGraph, the "
1085 "Matrix Market file must contain a \"pattern\" matrix.");
1112 static Teuchos::Tuple<global_ordinal_type, 3>
1113 readCoordDims (std::istream& in,
1115 const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1116 const Teuchos::RCP<const comm_type>& pComm,
1117 const bool tolerant =
false,
1120 using Teuchos::MatrixMarket::readCoordinateDimensions;
1121 using Teuchos::Tuple;
1126 Tuple<global_ordinal_type, 3> dims;
1132 bool success =
false;
1133 if (pComm->getRank() == 0) {
1134 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() !=
"coordinate",
1135 std::invalid_argument,
"The Tpetra::CrsMatrix Matrix Market reader "
1136 "only accepts \"coordinate\" (sparse) matrix data.");
1140 success = readCoordinateDimensions (in, numRows, numCols,
1141 numNonzeros, lineNumber,
1147 dims[2] = numNonzeros;
1155 int the_success = success ? 1 : 0;
1156 Teuchos::broadcast (*pComm, 0, &the_success);
1157 success = (the_success == 1);
1162 Teuchos::broadcast (*pComm, 0, dims);
1170 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1171 "Error reading Matrix Market sparse matrix: failed to read "
1172 "coordinate matrix dimensions.");
1187 typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1189 typedef Teuchos::MatrixMarket::SymmetrizingGraphAdder<Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> > graph_adder_type;
1216 static Teuchos::RCP<adder_type>
1217 makeAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1218 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1219 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1220 const bool tolerant=
false,
1221 const bool debug=
false)
1223 if (pComm->getRank () == 0) {
1224 typedef Teuchos::MatrixMarket::Raw::Adder<
scalar_type,
1227 Teuchos::RCP<raw_adder_type> pRaw =
1228 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1230 return Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
1233 return Teuchos::null;
1262 static Teuchos::RCP<graph_adder_type>
1263 makeGraphAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1264 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1265 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1266 const bool tolerant=
false,
1267 const bool debug=
false)
1269 if (pComm->getRank () == 0) {
1270 typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
1271 Teuchos::RCP<raw_adder_type> pRaw =
1272 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1274 return Teuchos::rcp (
new graph_adder_type (pRaw, pBanner->symmType ()));
1277 return Teuchos::null;
1282 static Teuchos::RCP<sparse_graph_type>
1283 readSparseGraphHelper (std::istream& in,
1284 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1285 const Teuchos::RCP<const map_type>& rowMap,
1286 Teuchos::RCP<const map_type>& colMap,
1287 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1288 const bool tolerant,
1291 using Teuchos::MatrixMarket::Banner;
1294 using Teuchos::Tuple;
1298 const int myRank = pComm->getRank ();
1299 const int rootRank = 0;
1304 size_t lineNumber = 1;
1306 if (debug && myRank == rootRank) {
1307 cerr <<
"Matrix Market reader: readGraph:" << endl
1308 <<
"-- Reading banner line" << endl;
1317 RCP<const Banner> pBanner;
1323 int bannerIsCorrect = 1;
1324 std::ostringstream errMsg;
1326 if (myRank == rootRank) {
1329 pBanner = readBanner (in, lineNumber, tolerant, debug,
true);
1331 catch (std::exception& e) {
1332 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
1333 "threw an exception: " << e.what();
1334 bannerIsCorrect = 0;
1337 if (bannerIsCorrect) {
1342 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
1343 bannerIsCorrect = 0;
1344 errMsg <<
"The Matrix Market input file must contain a "
1345 "\"coordinate\"-format sparse graph in order to create a "
1346 "Tpetra::CrsGraph object from it, but the file's matrix "
1347 "type is \"" << pBanner->matrixType() <<
"\" instead.";
1352 if (tolerant && pBanner->matrixType() ==
"array") {
1353 bannerIsCorrect = 0;
1354 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
1355 "format sparse graph in order to create a Tpetra::CrsGraph "
1356 "object from it, but the file's matrix type is \"array\" "
1357 "instead. That probably means the file contains dense matrix "
1364 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1371 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1372 std::invalid_argument, errMsg.str ());
1374 if (debug && myRank == rootRank) {
1375 cerr <<
"-- Reading dimensions line" << endl;
1383 Tuple<global_ordinal_type, 3> dims =
1384 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1386 if (debug && myRank == rootRank) {
1387 cerr <<
"-- Making Adder for collecting graph data" << endl;
1394 RCP<graph_adder_type> pAdder =
1395 makeGraphAdder (pComm, pBanner, dims, tolerant, debug);
1397 if (debug && myRank == rootRank) {
1398 cerr <<
"-- Reading graph data" << endl;
1408 int readSuccess = 1;
1409 std::ostringstream errMsg;
1410 if (myRank == rootRank) {
1413 typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
1415 reader_type reader (pAdder);
1418 std::pair<bool, std::vector<size_t> > results =
1419 reader.read (in, lineNumber, tolerant, debug);
1420 readSuccess = results.first ? 1 : 0;
1422 catch (std::exception& e) {
1427 broadcast (*pComm, rootRank, ptr (&readSuccess));
1436 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1437 "Failed to read the Matrix Market sparse graph file: "
1441 if (debug && myRank == rootRank) {
1442 cerr <<
"-- Successfully read the Matrix Market data" << endl;
1455 if (debug && myRank == rootRank) {
1456 cerr <<
"-- Tolerant mode: rebroadcasting graph dimensions"
1458 <<
"----- Dimensions before: "
1459 << dims[0] <<
" x " << dims[1]
1463 Tuple<global_ordinal_type, 2> updatedDims;
1464 if (myRank == rootRank) {
1471 std::max (dims[0], pAdder->getAdder()->numRows());
1472 updatedDims[1] = pAdder->getAdder()->numCols();
1474 broadcast (*pComm, rootRank, updatedDims);
1475 dims[0] = updatedDims[0];
1476 dims[1] = updatedDims[1];
1477 if (debug && myRank == rootRank) {
1478 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
1492 if (myRank == rootRank) {
1499 if (dims[0] < pAdder->getAdder ()->numRows ()) {
1503 broadcast (*pComm, 0, ptr (&dimsMatch));
1504 if (dimsMatch == 0) {
1511 Tuple<global_ordinal_type, 2> addersDims;
1512 if (myRank == rootRank) {
1513 addersDims[0] = pAdder->getAdder()->numRows();
1514 addersDims[1] = pAdder->getAdder()->numCols();
1516 broadcast (*pComm, 0, addersDims);
1517 TEUCHOS_TEST_FOR_EXCEPTION(
1518 dimsMatch == 0, std::runtime_error,
1519 "The graph metadata says that the graph is " << dims[0] <<
" x "
1520 << dims[1] <<
", but the actual data says that the graph is "
1521 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
1522 "data includes more rows than reported in the metadata. This "
1523 "is not allowed when parsing in strict mode. Parse the graph in "
1524 "tolerant mode to ignore the metadata when it disagrees with the "
1530 RCP<map_type> proc0Map;
1532 if(Teuchos::is_null(rowMap)) {
1536 indexBase = rowMap->getIndexBase();
1538 if(myRank == rootRank) {
1539 proc0Map = rcp(
new map_type(dims[0],dims[0],indexBase,pComm));
1542 proc0Map = rcp(
new map_type(dims[0],0,indexBase,pComm));
1546 std::map<global_ordinal_type, size_t> numEntriesPerRow_map;
1547 if (myRank == rootRank) {
1548 const auto& entries = pAdder()->getAdder()->getEntries();
1553 for (
const auto& entry : entries) {
1555 ++numEntriesPerRow_map[gblRow];
1559 Teuchos::Array<size_t> numEntriesPerRow (proc0Map->getLocalNumElements ());
1560 for (
const auto& ent : numEntriesPerRow_map) {
1562 numEntriesPerRow[lclRow] = ent.second;
1569 std::map<global_ordinal_type, size_t> empty_map;
1570 std::swap (numEntriesPerRow_map, empty_map);
1573 RCP<sparse_graph_type> proc0Graph =
1575 constructorParams));
1576 if(myRank == rootRank) {
1577 typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
1580 const std::vector<element_type>& entries =
1581 pAdder->getAdder()->getEntries();
1584 for(
size_t curPos=0; curPos<entries.size(); curPos++) {
1585 const element_type& curEntry = entries[curPos];
1588 Teuchos::ArrayView<const global_ordinal_type> colView(&curCol,1);
1589 proc0Graph->insertGlobalIndices(curRow,colView);
1592 proc0Graph->fillComplete();
1594 RCP<sparse_graph_type> distGraph;
1595 if(Teuchos::is_null(rowMap))
1598 RCP<map_type> distMap =
1599 rcp(
new map_type(dims[0],0,pComm,GloballyDistributed));
1605 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1606 import_type importer (proc0Map, distMap);
1609 distGraph->doImport(*proc0Graph,importer,
INSERT);
1615 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1616 import_type importer (proc0Map, rowMap);
1619 distGraph->doImport(*proc0Graph,importer,
INSERT);
1649 static Teuchos::RCP<sparse_graph_type>
1651 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1654 const bool debug=
false)
1656 using Teuchos::broadcast;
1657 using Teuchos::outArg;
1665 if (comm->getRank () == 0) {
1676 (
opened == 0, std::runtime_error,
"readSparseGraphFile: "
1677 "Failed to open file \"" <<
filename <<
"\" on Process 0.");
1715 static Teuchos::RCP<sparse_graph_type>
1717 const Teuchos::RCP<
const Teuchos::Comm<int> >&
pComm,
1721 const bool debug=
false)
1723 using Teuchos::broadcast;
1724 using Teuchos::outArg;
1732 if (
pComm->getRank () == 0) {
1743 (
opened == 0, std::runtime_error,
"readSparseGraphFile: "
1744 "Failed to open file \"" <<
filename <<
"\" on Process 0.");
1745 if (
pComm->getRank () == 0) {
1794 static Teuchos::RCP<sparse_graph_type>
1796 const Teuchos::RCP<const map_type>& rowMap,
1797 Teuchos::RCP<const map_type>& colMap,
1798 const Teuchos::RCP<const map_type>& domainMap,
1799 const Teuchos::RCP<const map_type>&
rangeMap,
1802 const bool debug=
false)
1804 using Teuchos::broadcast;
1805 using Teuchos::Comm;
1806 using Teuchos::outArg;
1810 (rowMap.is_null (), std::invalid_argument,
1811 "Input rowMap must be nonnull.");
1813 if (comm.is_null ()) {
1816 return Teuchos::null;
1825 if (comm->getRank () == 0) {
1836 (
opened == 0, std::runtime_error,
"readSparseGraphFile: "
1837 "Failed to open file \"" <<
filename <<
"\" on Process 0.");
1867 static Teuchos::RCP<sparse_graph_type>
1869 const Teuchos::RCP<
const Teuchos::Comm<int> >&
pComm,
1872 const bool debug=
false)
1878 Teuchos::RCP<sparse_graph_type>
graph =
1879 readSparseGraphHelper (
in,
pComm,
1883 graph->fillComplete ();
1918 static Teuchos::RCP<sparse_graph_type>
1920 const Teuchos::RCP<
const Teuchos::Comm<int> >&
pComm,
1924 const bool debug=
false)
1928 Teuchos::RCP<sparse_graph_type>
graph =
1929 readSparseGraphHelper (
in,
pComm,
1977 static Teuchos::RCP<sparse_graph_type>
1979 const Teuchos::RCP<const map_type>& rowMap,
1980 Teuchos::RCP<const map_type>& colMap,
1981 const Teuchos::RCP<const map_type>& domainMap,
1982 const Teuchos::RCP<const map_type>&
rangeMap,
1985 const bool debug=
false)
1987 Teuchos::RCP<sparse_graph_type>
graph =
1988 readSparseGraphHelper (
in, rowMap->getComm (),
1989 rowMap, colMap, Teuchos::null,
tolerant,
1997#include "MatrixMarket_TpetraNew.hpp"
2022 static Teuchos::RCP<sparse_matrix_type>
2024 const Teuchos::RCP<
const Teuchos::Comm<int> >&
pComm,
2027 const bool debug=
false)
2077 static Teuchos::RCP<sparse_matrix_type>
2079 const Teuchos::RCP<
const Teuchos::Comm<int> >&
pComm,
2083 const bool debug=
false)
2086 if (
pComm->getRank () == 0) {
2131 static Teuchos::RCP<sparse_matrix_type>
2133 const Teuchos::RCP<const map_type>& rowMap,
2134 Teuchos::RCP<const map_type>& colMap,
2135 const Teuchos::RCP<const map_type>& domainMap,
2136 const Teuchos::RCP<const map_type>&
rangeMap,
2139 const bool debug=
false)
2141 using Teuchos::broadcast;
2142 using Teuchos::Comm;
2143 using Teuchos::outArg;
2147 rowMap.is_null (), std::invalid_argument,
2148 "Row Map must be nonnull.");
2151 const int myRank = comm->getRank ();
2170 opened == 0, std::runtime_error,
2171 "readSparseFile: Failed to open file \"" <<
filename <<
"\" on "
2202 static Teuchos::RCP<sparse_matrix_type>
2204 const Teuchos::RCP<
const Teuchos::Comm<int> >&
pComm,
2207 const bool debug=
false)
2209 using Teuchos::MatrixMarket::Banner;
2210 using Teuchos::arcp;
2211 using Teuchos::ArrayRCP;
2212 using Teuchos::broadcast;
2213 using Teuchos::null;
2216 using Teuchos::REDUCE_MAX;
2217 using Teuchos::reduceAll;
2218 using Teuchos::Tuple;
2221 typedef Teuchos::ScalarTraits<scalar_type> STS;
2233 cerr <<
"Matrix Market reader: readSparse:" <<
endl
2234 <<
"-- Reading banner line" <<
endl;
2250 std::ostringstream
errMsg;
2257 catch (std::exception&
e) {
2258 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2259 "threw an exception: " <<
e.what();
2270 errMsg <<
"The Matrix Market input file must contain a "
2271 "\"coordinate\"-format sparse matrix in order to create a "
2272 "Tpetra::CrsMatrix object from it, but the file's matrix "
2273 "type is \"" <<
pBanner->matrixType() <<
"\" instead.";
2280 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2281 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2282 "object from it, but the file's matrix type is \"array\" "
2283 "instead. That probably means the file contains dense matrix "
2298 std::invalid_argument,
errMsg.str ());
2301 cerr <<
"-- Reading dimensions line" <<
endl;
2313 cerr <<
"-- Making Adder for collecting matrix data" <<
endl;
2322 cerr <<
"-- Reading matrix data" <<
endl;
2333 std::ostringstream
errMsg;
2337 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2342 std::pair<bool, std::vector<size_t> >
results =
2346 catch (std::exception&
e) {
2361 "Failed to read the Matrix Market sparse matrix file: "
2366 cerr <<
"-- Successfully read the Matrix Market data" <<
endl;
2380 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2382 <<
"----- Dimensions before: "
2395 std::max (
dims[0],
pAdder->getAdder()->numRows());
2402 cerr <<
"----- Dimensions after: " <<
dims[0] <<
" x "
2423 if (
dims[0] <
pAdder->getAdder ()->numRows ()) {
2443 "The matrix metadata says that the matrix is " <<
dims[0] <<
" x "
2444 <<
dims[1] <<
", but the actual data says that the matrix is "
2446 "data includes more rows than reported in the metadata. This "
2447 "is not allowed when parsing in strict mode. Parse the matrix in "
2448 "tolerant mode to ignore the metadata when it disagrees with the "
2454 cerr <<
"-- Converting matrix data into CSR format on Proc 0" <<
endl;
2475 std::ostringstream
errMsg;
2479 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2492 pAdder->getAdder()->merge ();
2495 const std::vector<element_type>& entries =
2496 pAdder->getAdder()->getEntries();
2499 const size_t numEntries = (
size_t)entries.size();
2502 cerr <<
"----- Proc 0: Matrix has numRows=" <<
numRows
2503 <<
" rows and numEntries=" << numEntries
2504 <<
" entries." <<
endl;
2527 "Row indices are out of order, even though they are supposed "
2528 "to be sorted. curRow = " <<
curRow <<
", prvRow = "
2529 <<
prvRow <<
", at curPos = " <<
curPos <<
". Please report "
2530 "this bug to the Tpetra developers.");
2546 catch (std::exception&
e) {
2548 errMsg <<
"Failed to merge sparse matrix entries and convert to "
2549 "CSR format: " <<
e.what();
2557 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
2560 cerr <<
"(only showing first and last few entries) ";
2565 for (size_type
k = 0;
k < 2; ++
k) {
2582 cerr <<
"----- Proc 0: rowPtr ";
2584 cerr <<
"(only showing first and last few entries) ";
2588 for (size_type
k = 0;
k < 2; ++
k) {
2613 cerr <<
"-- Making range, domain, and row maps" <<
endl;
2626 cerr <<
"-- Distributing the matrix data" <<
endl;
2649 cerr <<
"-- Inserting matrix entries on each processor";
2651 cerr <<
" and calling fillComplete()";
2673 "Reader::makeMatrix() returned a null pointer on at least one "
2674 "process. Please report this bug to the Tpetra developers.");
2678 "Reader::makeMatrix() returned a null pointer. "
2679 "Please report this bug to the Tpetra developers.");
2700 cerr <<
"-- Matrix is "
2702 <<
" with " <<
pMatrix->getGlobalNumEntries()
2703 <<
" entries, and index base "
2709 cerr <<
"-- Proc " <<
p <<
" owns "
2710 <<
pMatrix->getLocalNumCols() <<
" columns, and "
2711 <<
pMatrix->getLocalNumEntries() <<
" entries." <<
endl;
2719 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
2754 static Teuchos::RCP<sparse_matrix_type>
2756 const Teuchos::RCP<
const Teuchos::Comm<int> >&
pComm,
2760 const bool debug=
false)
2762 using Teuchos::MatrixMarket::Banner;
2763 using Teuchos::arcp;
2764 using Teuchos::ArrayRCP;
2765 using Teuchos::broadcast;
2766 using Teuchos::null;
2769 using Teuchos::reduceAll;
2770 using Teuchos::Tuple;
2773 typedef Teuchos::ScalarTraits<scalar_type> STS;
2785 cerr <<
"Matrix Market reader: readSparse:" <<
endl
2786 <<
"-- Reading banner line" <<
endl;
2802 std::ostringstream
errMsg;
2809 catch (std::exception&
e) {
2810 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2811 "threw an exception: " <<
e.what();
2822 errMsg <<
"The Matrix Market input file must contain a "
2823 "\"coordinate\"-format sparse matrix in order to create a "
2824 "Tpetra::CrsMatrix object from it, but the file's matrix "
2825 "type is \"" <<
pBanner->matrixType() <<
"\" instead.";
2832 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2833 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2834 "object from it, but the file's matrix type is \"array\" "
2835 "instead. That probably means the file contains dense matrix "
2850 std::invalid_argument,
errMsg.str ());
2853 cerr <<
"-- Reading dimensions line" <<
endl;
2865 cerr <<
"-- Making Adder for collecting matrix data" <<
endl;
2874 cerr <<
"-- Reading matrix data" <<
endl;
2885 std::ostringstream
errMsg;
2889 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2894 std::pair<bool, std::vector<size_t> >
results =
2898 catch (std::exception&
e) {
2913 "Failed to read the Matrix Market sparse matrix file: "
2918 cerr <<
"-- Successfully read the Matrix Market data" <<
endl;
2932 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2934 <<
"----- Dimensions before: "
2947 std::max (
dims[0],
pAdder->getAdder()->numRows());
2954 cerr <<
"----- Dimensions after: " <<
dims[0] <<
" x "
2975 if (
dims[0] <
pAdder->getAdder ()->numRows ()) {
2995 "The matrix metadata says that the matrix is " <<
dims[0] <<
" x "
2996 <<
dims[1] <<
", but the actual data says that the matrix is "
2998 "data includes more rows than reported in the metadata. This "
2999 "is not allowed when parsing in strict mode. Parse the matrix in "
3000 "tolerant mode to ignore the metadata when it disagrees with the "
3006 cerr <<
"-- Converting matrix data into CSR format on Proc 0" <<
endl;
3027 std::ostringstream
errMsg;
3031 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3044 pAdder->getAdder()->merge ();
3047 const std::vector<element_type>& entries =
3048 pAdder->getAdder()->getEntries();
3051 const size_t numEntries = (
size_t)entries.size();
3054 cerr <<
"----- Proc 0: Matrix has numRows=" <<
numRows
3055 <<
" rows and numEntries=" << numEntries
3056 <<
" entries." <<
endl;
3079 "Row indices are out of order, even though they are supposed "
3080 "to be sorted. curRow = " <<
curRow <<
", prvRow = "
3081 <<
prvRow <<
", at curPos = " <<
curPos <<
". Please report "
3082 "this bug to the Tpetra developers.");
3098 catch (std::exception&
e) {
3100 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3101 "CSR format: " <<
e.what();
3109 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3112 cerr <<
"(only showing first and last few entries) ";
3117 for (size_type
k = 0;
k < 2; ++
k) {
3134 cerr <<
"----- Proc 0: rowPtr ";
3136 cerr <<
"(only showing first and last few entries) ";
3140 for (size_type
k = 0;
k < 2; ++
k) {
3165 cerr <<
"-- Making range, domain, and row maps" <<
endl;
3178 cerr <<
"-- Distributing the matrix data" <<
endl;
3201 cerr <<
"-- Inserting matrix entries on each process "
3202 "and calling fillComplete()" <<
endl;
3211 Teuchos::RCP<sparse_matrix_type>
pMatrix =
3223 "Reader::makeMatrix() returned a null pointer on at least one "
3224 "process. Please report this bug to the Tpetra developers.");
3228 "Reader::makeMatrix() returned a null pointer. "
3229 "Please report this bug to the Tpetra developers.");
3245 cerr <<
"-- Matrix is "
3247 <<
" with " <<
pMatrix->getGlobalNumEntries()
3248 <<
" entries, and index base "
3254 cerr <<
"-- Proc " <<
p <<
" owns "
3255 <<
pMatrix->getLocalNumCols() <<
" columns, and "
3256 <<
pMatrix->getLocalNumEntries() <<
" entries." <<
endl;
3263 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3310 static Teuchos::RCP<sparse_matrix_type>
3312 const Teuchos::RCP<const map_type>& rowMap,
3313 Teuchos::RCP<const map_type>& colMap,
3314 const Teuchos::RCP<const map_type>& domainMap,
3315 const Teuchos::RCP<const map_type>&
rangeMap,
3318 const bool debug=
false)
3320 using Teuchos::MatrixMarket::Banner;
3321 using Teuchos::arcp;
3322 using Teuchos::ArrayRCP;
3323 using Teuchos::ArrayView;
3325 using Teuchos::broadcast;
3326 using Teuchos::Comm;
3327 using Teuchos::null;
3330 using Teuchos::reduceAll;
3331 using Teuchos::Tuple;
3334 typedef Teuchos::ScalarTraits<scalar_type> STS;
3345 rowMap.is_null (), std::invalid_argument,
3346 "Row Map must be nonnull.");
3348 rangeMap.is_null (), std::invalid_argument,
3349 "Range Map must be nonnull.");
3351 domainMap.is_null (), std::invalid_argument,
3352 "Domain Map must be nonnull.");
3354 rowMap->getComm().getRawPtr() !=
pComm.getRawPtr(),
3355 std::invalid_argument,
3356 "The specified row Map's communicator (rowMap->getComm())"
3357 "differs from the given separately supplied communicator pComm.");
3359 domainMap->getComm().getRawPtr() !=
pComm.getRawPtr(),
3360 std::invalid_argument,
3361 "The specified domain Map's communicator (domainMap->getComm())"
3362 "differs from the given separately supplied communicator pComm.");
3365 std::invalid_argument,
3366 "The specified range Map's communicator (rangeMap->getComm())"
3367 "differs from the given separately supplied communicator pComm.");
3375 cerr <<
"Matrix Market reader: readSparse:" <<
endl
3376 <<
"-- Reading banner line" <<
endl;
3392 std::ostringstream
errMsg;
3399 catch (std::exception&
e) {
3400 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
3401 "threw an exception: " <<
e.what();
3412 errMsg <<
"The Matrix Market input file must contain a "
3413 "\"coordinate\"-format sparse matrix in order to create a "
3414 "Tpetra::CrsMatrix object from it, but the file's matrix "
3415 "type is \"" <<
pBanner->matrixType() <<
"\" instead.";
3422 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
3423 "format sparse matrix in order to create a Tpetra::CrsMatrix "
3424 "object from it, but the file's matrix type is \"array\" "
3425 "instead. That probably means the file contains dense matrix "
3440 std::invalid_argument,
errMsg.str ());
3443 cerr <<
"-- Reading dimensions line" <<
endl;
3455 cerr <<
"-- Making Adder for collecting matrix data" <<
endl;
3466 cerr <<
"-- Reading matrix data" <<
endl;
3477 std::ostringstream
errMsg;
3481 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3486 std::pair<bool, std::vector<size_t> >
results =
3490 catch (std::exception&
e) {
3505 "Failed to read the Matrix Market sparse matrix file: "
3510 cerr <<
"-- Successfully read the Matrix Market data" <<
endl;
3524 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
3526 <<
"----- Dimensions before: "
3539 std::max (
dims[0],
pAdder->getAdder()->numRows());
3546 cerr <<
"----- Dimensions after: " <<
dims[0] <<
" x "
3567 if (
dims[0] <
pAdder->getAdder ()->numRows ()) {
3587 "The matrix metadata says that the matrix is " <<
dims[0] <<
" x "
3588 <<
dims[1] <<
", but the actual data says that the matrix is "
3590 "data includes more rows than reported in the metadata. This "
3591 "is not allowed when parsing in strict mode. Parse the matrix in "
3592 "tolerant mode to ignore the metadata when it disagrees with the "
3598 cerr <<
"-- Converting matrix data into CSR format on Proc 0" <<
endl;
3619 std::ostringstream
errMsg;
3623 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3636 pAdder->getAdder()->merge ();
3639 const std::vector<element_type>& entries =
3640 pAdder->getAdder()->getEntries();
3643 const size_t numEntries = (
size_t)entries.size();
3646 cerr <<
"----- Proc 0: Matrix has numRows=" <<
numRows
3647 <<
" rows and numEntries=" << numEntries
3648 <<
" entries." <<
endl;
3671 "Row indices are out of order, even though they are supposed "
3672 "to be sorted. curRow = " <<
curRow <<
", prvRow = "
3673 <<
prvRow <<
", at curPos = " <<
curPos <<
". Please report "
3674 "this bug to the Tpetra developers.");
3690 catch (std::exception&
e) {
3692 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3693 "CSR format: " <<
e.what();
3701 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3704 cerr <<
"(only showing first and last few entries) ";
3709 for (size_type
k = 0;
k < 2; ++
k) {
3726 cerr <<
"----- Proc 0: rowPtr ";
3728 cerr <<
"(only showing first and last few entries) ";
3732 for (size_type
k = 0;
k < 2; ++
k) {
3747 cerr <<
"----- Proc 0: colInd = [";
3766 cerr <<
"-- Verifying Maps" <<
endl;
3770 std::invalid_argument,
3771 "The range Map has " <<
rangeMap->getGlobalNumElements ()
3772 <<
" entries, but the matrix has a global number of rows " <<
dims[0]
3776 std::invalid_argument,
3777 "The domain Map has " << domainMap->getGlobalNumElements ()
3778 <<
" entries, but the matrix has a global number of columns "
3783 Teuchos::getFancyOStream (Teuchos::rcpFromRef (
cerr)) :
null;
3803 cerr <<
"-- Proc 0: Filling gather matrix" <<
endl;
3863 if (colMap.is_null ()) {
3889 cerr <<
"-- Matrix is "
3891 <<
" with " <<
A->getGlobalNumEntries()
3892 <<
" entries, and index base "
3893 <<
A->getIndexBase() <<
"." <<
endl;
3898 cerr <<
"-- Proc " <<
p <<
" owns "
3899 <<
A->getLocalNumCols() <<
" columns, and "
3900 <<
A->getLocalNumEntries() <<
" entries." <<
endl;
3908 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3943 static Teuchos::RCP<multivector_type>
3945 const Teuchos::RCP<const comm_type>& comm,
3946 Teuchos::RCP<const map_type>&
map,
3948 const bool debug=
false,
3951 using Teuchos::broadcast;
3952 using Teuchos::outArg;
3956 if (comm->getRank() == 0) {
3961 in.open (
filename.c_str (), std::ios::binary);
3970 opened == 0, std::runtime_error,
3971 "readDenseFile: Failed to open file \"" <<
filename <<
"\" on "
4006 static Teuchos::RCP<vector_type>
4008 const Teuchos::RCP<const comm_type>& comm,
4009 Teuchos::RCP<const map_type>&
map,
4011 const bool debug=
false)
4013 using Teuchos::broadcast;
4014 using Teuchos::outArg;
4018 if (comm->getRank() == 0) {
4029 opened == 0, std::runtime_error,
4030 "readVectorFile: Failed to open file \"" <<
filename <<
"\" on "
4104 static Teuchos::RCP<multivector_type>
4106 const Teuchos::RCP<const comm_type>& comm,
4107 Teuchos::RCP<const map_type>&
map,
4109 const bool debug=
false,
4112 Teuchos::RCP<Teuchos::FancyOStream>
err =
4113 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4119 static Teuchos::RCP<vector_type>
4121 const Teuchos::RCP<const comm_type>& comm,
4122 Teuchos::RCP<const map_type>&
map,
4124 const bool debug=
false)
4126 Teuchos::RCP<Teuchos::FancyOStream>
err =
4127 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4153 static Teuchos::RCP<const map_type>
4155 const Teuchos::RCP<const comm_type>& comm,
4157 const bool debug=
false,
4160 using Teuchos::inOutArg;
4161 using Teuchos::broadcast;
4165 if (comm->getRank () == 0) {
4167 in.open (
filename.c_str (), std::ios::binary);
4176 success == 0, std::runtime_error,
4177 "Tpetra::MatrixMarket::Reader::readMapFile: "
4178 "Failed to read file \"" <<
filename <<
"\" on Process 0.");
4184 template<
class MultiVectorScalarType>
4189 readDenseImpl (std::istream&
in,
4190 const Teuchos::RCP<const comm_type>& comm,
4191 Teuchos::RCP<const map_type>&
map,
4192 const Teuchos::RCP<Teuchos::FancyOStream>&
err,
4194 const bool debug=
false,
4197 using Teuchos::MatrixMarket::Banner;
4198 using Teuchos::MatrixMarket::checkCommentLine;
4199 using Teuchos::ArrayRCP;
4201 using Teuchos::broadcast;
4202 using Teuchos::outArg;
4204 using Teuchos::Tuple;
4210 typedef Teuchos::ScalarTraits<ST> STS;
4211 typedef typename STS::magnitudeType
MT;
4212 typedef Teuchos::ScalarTraits<MT> STM;
4217 const int myRank = comm->getRank ();
4219 if (!
err.is_null ()) {
4225 if (! err.is_null ()) {
4260 size_t lineNumber = 1;
4263 std::ostringstream exMsg;
4264 int localBannerReadSuccess = 1;
4265 int localDimsReadSuccess = 1;
4271 *err << myRank <<
": readDenseImpl: Reading banner line (dense)" << endl;
4277 RCP<const Banner> pBanner;
4279 pBanner = readBanner (in, lineNumber, tolerant, debug);
4280 }
catch (std::exception& e) {
4282 localBannerReadSuccess = 0;
4285 if (localBannerReadSuccess) {
4286 if (pBanner->matrixType () !=
"array") {
4287 exMsg <<
"The Matrix Market file does not contain dense matrix "
4288 "data. Its banner (first) line says that its matrix type is \""
4289 << pBanner->matrixType () <<
"\", rather that the required "
4291 localBannerReadSuccess = 0;
4292 }
else if (pBanner->dataType() ==
"pattern") {
4293 exMsg <<
"The Matrix Market file's banner (first) "
4294 "line claims that the matrix's data type is \"pattern\". This does "
4295 "not make sense for a dense matrix, yet the file reports the matrix "
4296 "as dense. The only valid data types for a dense matrix are "
4297 "\"real\", \"complex\", and \"integer\".";
4298 localBannerReadSuccess = 0;
4302 dims[2] = encodeDataType (pBanner->dataType ());
4308 if (localBannerReadSuccess) {
4310 *err << myRank <<
": readDenseImpl: Reading dimensions line (dense)" << endl;
4319 bool commentLine =
true;
4321 while (commentLine) {
4324 if (in.eof () || in.fail ()) {
4325 exMsg <<
"Unable to get array dimensions line (at all) from line "
4326 << lineNumber <<
" of input stream. The input stream "
4327 <<
"claims that it is "
4328 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4329 localDimsReadSuccess = 0;
4332 if (getline (in, line)) {
4339 size_t start = 0, size = 0;
4340 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4347 std::istringstream istr (line);
4351 if (istr.eof () || istr.fail ()) {
4352 exMsg <<
"Unable to read any data from line " << lineNumber
4353 <<
" of input; the line should contain the matrix dimensions "
4354 <<
"\"<numRows> <numCols>\".";
4355 localDimsReadSuccess = 0;
4360 exMsg <<
"Failed to get number of rows from line "
4361 << lineNumber <<
" of input; the line should contains the "
4362 <<
"matrix dimensions \"<numRows> <numCols>\".";
4363 localDimsReadSuccess = 0;
4365 dims[0] = theNumRows;
4367 exMsg <<
"No more data after number of rows on line "
4368 << lineNumber <<
" of input; the line should contain the "
4369 <<
"matrix dimensions \"<numRows> <numCols>\".";
4370 localDimsReadSuccess = 0;
4375 exMsg <<
"Failed to get number of columns from line "
4376 << lineNumber <<
" of input; the line should contain "
4377 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4378 localDimsReadSuccess = 0;
4380 dims[1] = theNumCols;
4389 in.read(
reinterpret_cast<char*
>(&numRows),
sizeof(numRows));
4390 in.read(
reinterpret_cast<char*
>(&numCols),
sizeof(numCols));
4391 dims[0] = Teuchos::as<GO>(numRows);
4392 dims[1] = Teuchos::as<GO>(numCols);
4393 if ((
typeid(ST) ==
typeid(
double)) || Teuchos::ScalarTraits<ST>::isOrdinal) {
4395 }
else if (Teuchos::ScalarTraits<ST>::isComplex) {
4398 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
4399 "Unrecognized Matrix Market data type. "
4400 "We should never get here. "
4401 "Please report this bug to the Tpetra developers.");
4403 localDimsReadSuccess =
true;
4404 localDimsReadSuccess =
true;
4411 Tuple<GO, 5> bannerDimsReadResult;
4413 bannerDimsReadResult[0] = dims[0];
4414 bannerDimsReadResult[1] = dims[1];
4415 bannerDimsReadResult[2] = dims[2];
4416 bannerDimsReadResult[3] = localBannerReadSuccess;
4417 bannerDimsReadResult[4] = localDimsReadSuccess;
4421 broadcast (*comm, 0, bannerDimsReadResult);
4423 TEUCHOS_TEST_FOR_EXCEPTION(
4424 bannerDimsReadResult[3] == 0, std::runtime_error,
4425 "Failed to read banner line: " << exMsg.str ());
4426 TEUCHOS_TEST_FOR_EXCEPTION(
4427 bannerDimsReadResult[4] == 0, std::runtime_error,
4428 "Failed to read matrix dimensions line: " << exMsg.str ());
4430 dims[0] = bannerDimsReadResult[0];
4431 dims[1] = bannerDimsReadResult[1];
4432 dims[2] = bannerDimsReadResult[2];
4437 const size_t numCols =
static_cast<size_t> (dims[1]);
4442 RCP<const map_type> proc0Map;
4443 if (map.is_null ()) {
4447 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4448 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4450 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4454 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4458 RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
4464 int localReadDataSuccess = 1;
4469 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4474 TEUCHOS_TEST_FOR_EXCEPTION(
4475 ! X->isConstantStride (), std::logic_error,
4476 "Can't get a 1-D view of the entries of the MultiVector X on "
4477 "Process 0, because the stride between the columns of X is not "
4478 "constant. This shouldn't happen because we just created X and "
4479 "haven't filled it in yet. Please report this bug to the Tpetra "
4486 ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4487 TEUCHOS_TEST_FOR_EXCEPTION(
4488 as<global_size_t> (X_view.size ()) < numRows * numCols,
4490 "The view of X has size " << X_view.size() <<
" which is not enough to "
4491 "accommodate the expected number of entries numRows*numCols = "
4492 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
4493 "Please report this bug to the Tpetra developers.");
4494 const size_t stride = X->getStride ();
4501 const bool isComplex = (dims[2] == 1);
4502 size_type count = 0, curRow = 0, curCol = 0;
4505 while (getline (in, line)) {
4509 size_t start = 0, size = 0;
4510 const bool commentLine =
4511 checkCommentLine (line, start, size, lineNumber, tolerant);
4512 if (! commentLine) {
4518 if (count >= X_view.size()) {
4523 TEUCHOS_TEST_FOR_EXCEPTION(
4524 count >= X_view.size(),
4526 "The Matrix Market input stream has more data in it than "
4527 "its metadata reported. Current line number is "
4528 << lineNumber <<
".");
4537 const size_t pos = line.substr (start, size).find (
':');
4538 if (pos != std::string::npos) {
4542 std::istringstream istr (line.substr (start, size));
4546 if (istr.eof() || istr.fail()) {
4553 TEUCHOS_TEST_FOR_EXCEPTION(
4554 ! tolerant && (istr.eof() || istr.fail()),
4556 "Line " << lineNumber <<
" of the Matrix Market file is "
4557 "empty, or we cannot read from it for some other reason.");
4560 ST val = STS::zero();
4563 MT real = STM::zero(), imag = STM::zero();
4580 TEUCHOS_TEST_FOR_EXCEPTION(
4581 ! tolerant && istr.eof(), std::runtime_error,
4582 "Failed to get the real part of a complex-valued matrix "
4583 "entry from line " << lineNumber <<
" of the Matrix Market "
4589 }
else if (istr.eof()) {
4590 TEUCHOS_TEST_FOR_EXCEPTION(
4591 ! tolerant && istr.eof(), std::runtime_error,
4592 "Missing imaginary part of a complex-valued matrix entry "
4593 "on line " << lineNumber <<
" of the Matrix Market file.");
4599 TEUCHOS_TEST_FOR_EXCEPTION(
4600 ! tolerant && istr.fail(), std::runtime_error,
4601 "Failed to get the imaginary part of a complex-valued "
4602 "matrix entry from line " << lineNumber <<
" of the "
4603 "Matrix Market file.");
4610 TEUCHOS_TEST_FOR_EXCEPTION(
4611 ! tolerant && istr.fail(), std::runtime_error,
4612 "Failed to get a real-valued matrix entry from line "
4613 << lineNumber <<
" of the Matrix Market file.");
4616 if (istr.fail() && tolerant) {
4622 TEUCHOS_TEST_FOR_EXCEPTION(
4623 ! tolerant && istr.fail(), std::runtime_error,
4624 "Failed to read matrix data from line " << lineNumber
4625 <<
" of the Matrix Market file.");
4628 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4630 curRow = count % numRows;
4631 curCol = count / numRows;
4632 X_view[curRow + curCol*stride] = val;
4637 TEUCHOS_TEST_FOR_EXCEPTION(
4638 ! tolerant &&
static_cast<global_size_t> (count) < numRows * numCols,
4640 "The Matrix Market metadata reports that the dense matrix is "
4641 << numRows <<
" x " << numCols <<
", and thus has "
4642 << numRows*numCols <<
" total entries, but we only found " << count
4643 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
4644 }
catch (std::exception& e) {
4646 localReadDataSuccess = 0;
4651 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4656 TEUCHOS_TEST_FOR_EXCEPTION(
4657 ! X->isConstantStride (), std::logic_error,
4658 "Can't get a 1-D view of the entries of the MultiVector X on "
4659 "Process 0, because the stride between the columns of X is not "
4660 "constant. This shouldn't happen because we just created X and "
4661 "haven't filled it in yet. Please report this bug to the Tpetra "
4668 auto X_view = X->getLocalViewHost (Access::OverwriteAll);
4670 TEUCHOS_TEST_FOR_EXCEPTION(
4671 as<global_size_t> (X_view.extent(0)) < numRows,
4673 "The view of X has " << X_view.extent(0) <<
" rows which is not enough to "
4674 "accommodate the expected number of entries numRows = "
4676 "Please report this bug to the Tpetra developers.");
4677 TEUCHOS_TEST_FOR_EXCEPTION(
4678 as<global_size_t> (X_view.extent(1)) < numCols,
4680 "The view of X has " << X_view.extent(1) <<
" colums which is not enough to "
4681 "accommodate the expected number of entries numRows = "
4683 "Please report this bug to the Tpetra developers.");
4685 for (
size_t curRow = 0; curRow < numRows; ++curRow) {
4686 for (
size_t curCol = 0; curCol < numCols; ++curCol) {
4687 if (Teuchos::ScalarTraits<ST>::isOrdinal){
4689 in.read(
reinterpret_cast<char*
>(&val),
sizeof(val));
4690 X_view(curRow, curCol) = val;
4693 in.read(
reinterpret_cast<char*
>(&val),
sizeof(val));
4694 X_view(curRow, curCol) = val;
4702 *err << myRank <<
": readDenseImpl: done reading data" << endl;
4706 int globalReadDataSuccess = localReadDataSuccess;
4707 broadcast (*comm, 0, outArg (globalReadDataSuccess));
4708 TEUCHOS_TEST_FOR_EXCEPTION(
4709 globalReadDataSuccess == 0, std::runtime_error,
4710 "Failed to read the multivector's data: " << exMsg.str ());
4715 if (comm->getSize () == 1 && map.is_null ()) {
4717 if (! err.is_null ()) {
4721 *err << myRank <<
": readDenseImpl: done" << endl;
4723 if (! err.is_null ()) {
4730 *err << myRank <<
": readDenseImpl: Creating target MV" << endl;
4734 RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
4737 *err << myRank <<
": readDenseImpl: Creating Export" << endl;
4743 Export<LO, GO, NT> exporter (proc0Map, map, err);
4746 *err << myRank <<
": readDenseImpl: Exporting" << endl;
4749 Y->doExport (*X, exporter,
INSERT);
4751 if (! err.is_null ()) {
4755 *err << myRank <<
": readDenseImpl: done" << endl;
4757 if (! err.is_null ()) {
4766 template<
class VectorScalarType>
4771 readVectorImpl (std::istream& in,
4772 const Teuchos::RCP<const comm_type>& comm,
4773 Teuchos::RCP<const map_type>& map,
4774 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4775 const bool tolerant=
false,
4776 const bool debug=
false)
4778 using Teuchos::MatrixMarket::Banner;
4779 using Teuchos::MatrixMarket::checkCommentLine;
4781 using Teuchos::broadcast;
4782 using Teuchos::outArg;
4784 using Teuchos::Tuple;
4786 typedef VectorScalarType ST;
4790 typedef Teuchos::ScalarTraits<ST> STS;
4791 typedef typename STS::magnitudeType MT;
4792 typedef Teuchos::ScalarTraits<MT> STM;
4797 const int myRank = comm->getRank ();
4799 if (! err.is_null ()) {
4803 *err << myRank <<
": readVectorImpl" << endl;
4805 if (! err.is_null ()) {
4839 size_t lineNumber = 1;
4842 std::ostringstream exMsg;
4843 int localBannerReadSuccess = 1;
4844 int localDimsReadSuccess = 1;
4849 *err << myRank <<
": readVectorImpl: Reading banner line (dense)" << endl;
4855 RCP<const Banner> pBanner;
4857 pBanner = readBanner (in, lineNumber, tolerant, debug);
4858 }
catch (std::exception& e) {
4860 localBannerReadSuccess = 0;
4863 if (localBannerReadSuccess) {
4864 if (pBanner->matrixType () !=
"array") {
4865 exMsg <<
"The Matrix Market file does not contain dense matrix "
4866 "data. Its banner (first) line says that its matrix type is \""
4867 << pBanner->matrixType () <<
"\", rather that the required "
4869 localBannerReadSuccess = 0;
4870 }
else if (pBanner->dataType() ==
"pattern") {
4871 exMsg <<
"The Matrix Market file's banner (first) "
4872 "line claims that the matrix's data type is \"pattern\". This does "
4873 "not make sense for a dense matrix, yet the file reports the matrix "
4874 "as dense. The only valid data types for a dense matrix are "
4875 "\"real\", \"complex\", and \"integer\".";
4876 localBannerReadSuccess = 0;
4880 dims[2] = encodeDataType (pBanner->dataType ());
4886 if (localBannerReadSuccess) {
4888 *err << myRank <<
": readVectorImpl: Reading dimensions line (dense)" << endl;
4897 bool commentLine =
true;
4899 while (commentLine) {
4902 if (in.eof () || in.fail ()) {
4903 exMsg <<
"Unable to get array dimensions line (at all) from line "
4904 << lineNumber <<
" of input stream. The input stream "
4905 <<
"claims that it is "
4906 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4907 localDimsReadSuccess = 0;
4910 if (getline (in, line)) {
4917 size_t start = 0, size = 0;
4918 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4925 std::istringstream istr (line);
4929 if (istr.eof () || istr.fail ()) {
4930 exMsg <<
"Unable to read any data from line " << lineNumber
4931 <<
" of input; the line should contain the matrix dimensions "
4932 <<
"\"<numRows> <numCols>\".";
4933 localDimsReadSuccess = 0;
4938 exMsg <<
"Failed to get number of rows from line "
4939 << lineNumber <<
" of input; the line should contains the "
4940 <<
"matrix dimensions \"<numRows> <numCols>\".";
4941 localDimsReadSuccess = 0;
4943 dims[0] = theNumRows;
4945 exMsg <<
"No more data after number of rows on line "
4946 << lineNumber <<
" of input; the line should contain the "
4947 <<
"matrix dimensions \"<numRows> <numCols>\".";
4948 localDimsReadSuccess = 0;
4953 exMsg <<
"Failed to get number of columns from line "
4954 << lineNumber <<
" of input; the line should contain "
4955 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4956 localDimsReadSuccess = 0;
4958 dims[1] = theNumCols;
4968 exMsg <<
"File does not contain a 1D Vector";
4969 localDimsReadSuccess = 0;
4975 Tuple<GO, 5> bannerDimsReadResult;
4977 bannerDimsReadResult[0] = dims[0];
4978 bannerDimsReadResult[1] = dims[1];
4979 bannerDimsReadResult[2] = dims[2];
4980 bannerDimsReadResult[3] = localBannerReadSuccess;
4981 bannerDimsReadResult[4] = localDimsReadSuccess;
4986 broadcast (*comm, 0, bannerDimsReadResult);
4988 TEUCHOS_TEST_FOR_EXCEPTION(
4989 bannerDimsReadResult[3] == 0, std::runtime_error,
4990 "Failed to read banner line: " << exMsg.str ());
4991 TEUCHOS_TEST_FOR_EXCEPTION(
4992 bannerDimsReadResult[4] == 0, std::runtime_error,
4993 "Failed to read matrix dimensions line: " << exMsg.str ());
4995 dims[0] = bannerDimsReadResult[0];
4996 dims[1] = bannerDimsReadResult[1];
4997 dims[2] = bannerDimsReadResult[2];
5002 const size_t numCols =
static_cast<size_t> (dims[1]);
5007 RCP<const map_type> proc0Map;
5008 if (map.is_null ()) {
5012 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
5013 const size_t localNumRows = (myRank == 0) ? numRows : 0;
5015 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
5019 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
5023 RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
5029 int localReadDataSuccess = 1;
5033 *err << myRank <<
": readVectorImpl: Reading matrix data (dense)"
5038 TEUCHOS_TEST_FOR_EXCEPTION(
5039 ! X->isConstantStride (), std::logic_error,
5040 "Can't get a 1-D view of the entries of the MultiVector X on "
5041 "Process 0, because the stride between the columns of X is not "
5042 "constant. This shouldn't happen because we just created X and "
5043 "haven't filled it in yet. Please report this bug to the Tpetra "
5050 Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst ();
5051 TEUCHOS_TEST_FOR_EXCEPTION(
5052 as<global_size_t> (X_view.size ()) < numRows * numCols,
5054 "The view of X has size " << X_view <<
" which is not enough to "
5055 "accommodate the expected number of entries numRows*numCols = "
5056 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
5057 "Please report this bug to the Tpetra developers.");
5058 const size_t stride = X->getStride ();
5065 const bool isComplex = (dims[2] == 1);
5066 size_type count = 0, curRow = 0, curCol = 0;
5069 while (getline (in, line)) {
5073 size_t start = 0, size = 0;
5074 const bool commentLine =
5075 checkCommentLine (line, start, size, lineNumber, tolerant);
5076 if (! commentLine) {
5082 if (count >= X_view.size()) {
5087 TEUCHOS_TEST_FOR_EXCEPTION(
5088 count >= X_view.size(),
5090 "The Matrix Market input stream has more data in it than "
5091 "its metadata reported. Current line number is "
5092 << lineNumber <<
".");
5101 const size_t pos = line.substr (start, size).find (
':');
5102 if (pos != std::string::npos) {
5106 std::istringstream istr (line.substr (start, size));
5110 if (istr.eof() || istr.fail()) {
5117 TEUCHOS_TEST_FOR_EXCEPTION(
5118 ! tolerant && (istr.eof() || istr.fail()),
5120 "Line " << lineNumber <<
" of the Matrix Market file is "
5121 "empty, or we cannot read from it for some other reason.");
5124 ST val = STS::zero();
5127 MT real = STM::zero(), imag = STM::zero();
5144 TEUCHOS_TEST_FOR_EXCEPTION(
5145 ! tolerant && istr.eof(), std::runtime_error,
5146 "Failed to get the real part of a complex-valued matrix "
5147 "entry from line " << lineNumber <<
" of the Matrix Market "
5153 }
else if (istr.eof()) {
5154 TEUCHOS_TEST_FOR_EXCEPTION(
5155 ! tolerant && istr.eof(), std::runtime_error,
5156 "Missing imaginary part of a complex-valued matrix entry "
5157 "on line " << lineNumber <<
" of the Matrix Market file.");
5163 TEUCHOS_TEST_FOR_EXCEPTION(
5164 ! tolerant && istr.fail(), std::runtime_error,
5165 "Failed to get the imaginary part of a complex-valued "
5166 "matrix entry from line " << lineNumber <<
" of the "
5167 "Matrix Market file.");
5174 TEUCHOS_TEST_FOR_EXCEPTION(
5175 ! tolerant && istr.fail(), std::runtime_error,
5176 "Failed to get a real-valued matrix entry from line "
5177 << lineNumber <<
" of the Matrix Market file.");
5180 if (istr.fail() && tolerant) {
5186 TEUCHOS_TEST_FOR_EXCEPTION(
5187 ! tolerant && istr.fail(), std::runtime_error,
5188 "Failed to read matrix data from line " << lineNumber
5189 <<
" of the Matrix Market file.");
5192 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
5194 curRow = count % numRows;
5195 curCol = count / numRows;
5196 X_view[curRow + curCol*stride] = val;
5201 TEUCHOS_TEST_FOR_EXCEPTION(
5202 ! tolerant &&
static_cast<global_size_t> (count) < numRows * numCols,
5204 "The Matrix Market metadata reports that the dense matrix is "
5205 << numRows <<
" x " << numCols <<
", and thus has "
5206 << numRows*numCols <<
" total entries, but we only found " << count
5207 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
5208 }
catch (std::exception& e) {
5210 localReadDataSuccess = 0;
5215 *err << myRank <<
": readVectorImpl: done reading data" << endl;
5219 int globalReadDataSuccess = localReadDataSuccess;
5220 broadcast (*comm, 0, outArg (globalReadDataSuccess));
5221 TEUCHOS_TEST_FOR_EXCEPTION(
5222 globalReadDataSuccess == 0, std::runtime_error,
5223 "Failed to read the multivector's data: " << exMsg.str ());
5228 if (comm->getSize () == 1 && map.is_null ()) {
5230 if (! err.is_null ()) {
5234 *err << myRank <<
": readVectorImpl: done" << endl;
5236 if (! err.is_null ()) {
5243 *err << myRank <<
": readVectorImpl: Creating target MV" << endl;
5247 RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
5250 *err << myRank <<
": readVectorImpl: Creating Export" << endl;
5256 Export<LO, GO, NT> exporter (proc0Map, map, err);
5259 *err << myRank <<
": readVectorImpl: Exporting" << endl;
5262 Y->doExport (*X, exporter,
INSERT);
5264 if (! err.is_null ()) {
5268 *err << myRank <<
": readVectorImpl: done" << endl;
5270 if (! err.is_null ()) {
5299 static Teuchos::RCP<const map_type>
5301 const Teuchos::RCP<const comm_type>& comm,
5303 const bool debug=
false,
5306 Teuchos::RCP<Teuchos::FancyOStream>
err =
5307 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
5338 static Teuchos::RCP<const map_type>
5340 const Teuchos::RCP<const comm_type>& comm,
5341 const Teuchos::RCP<Teuchos::FancyOStream>&
err,
5343 const bool debug=
false,
5346 using Teuchos::arcp;
5347 using Teuchos::Array;
5348 using Teuchos::ArrayRCP;
5350 using Teuchos::broadcast;
5351 using Teuchos::Comm;
5352 using Teuchos::CommRequest;
5353 using Teuchos::inOutArg;
5354 using Teuchos::ireceive;
5355 using Teuchos::outArg;
5357 using Teuchos::receive;
5358 using Teuchos::reduceAll;
5359 using Teuchos::REDUCE_MIN;
5360 using Teuchos::isend;
5361 using Teuchos::SerialComm;
5362 using Teuchos::toString;
5363 using Teuchos::wait;
5372 const int numProcs = comm->getSize ();
5373 const int myRank = comm->getRank ();
5375 if (
err.is_null ()) {
5379 std::ostringstream
os;
5383 if (
err.is_null ()) {
5410 std::ostringstream
exMsg;
5427 if (
data.is_null ()) {
5429 exMsg <<
"readDenseImpl() returned null." <<
endl;
5431 }
catch (std::exception&
e) {
5439 std::map<int, Array<GO> >
pid2gids;
5455 if (
data->getNumVectors () == 2) {
5463 const int pid =
static_cast<int> (
PIDs[0]);
5467 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5468 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5480 const int pid =
static_cast<int> (
PIDs[
k]);
5483 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5484 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5496 else if (
data->getNumVectors () == 1) {
5497 if (
data->getGlobalLength () % 2 !=
static_cast<GST> (0)) {
5499 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data has the "
5500 "wrong format (for Map format 2.0). The global number of rows "
5501 "in the MultiVector must be even (divisible by 2)." <<
endl;
5506 static_cast<GO
> (2);
5511 const int pid =
static_cast<int> (
theData[1]);
5514 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5515 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5526 const int pid =
static_cast<int> (
theData[2*
k + 1]);
5529 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5530 <<
"Encountered invalid PID " <<
pid <<
"." <<
endl;
5544 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data must have "
5545 "either 1 column (for the new Map format 2.0) or 2 columns (for "
5546 "the old Map format 1.0).";
5569 "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
5570 "following exception message: " <<
exMsg.str ());
5599 typename std::map<int, Array<GO> >::iterator
it =
pid2gids.find (0);
5627 typename std::map<int, Array<GO> >::iterator
it =
pid2gids.find (
p);
5646 std::ostringstream
os;
5650 const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
5662 catch (std::exception&
e) {
5664 errStrm <<
"Process " << comm->getRank () <<
" threw an exception: "
5665 <<
e.what () << std::endl;
5669 errStrm <<
"Process " << comm->getRank () <<
" threw an exception "
5670 "not a subclass of std::exception" << std::endl;
5672 Teuchos::reduceAll<int, int> (*comm, Teuchos::REDUCE_MIN,
5679 if (
err.is_null ()) {
5683 std::ostringstream
os;
5687 if (
err.is_null ()) {
5707 encodeDataType (
const std::string&
dataType)
5711 }
else if (
dataType ==
"complex") {
5713 }
else if (dataType ==
"pattern") {
5719 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
5720 "Unrecognized Matrix Market data type \"" << dataType
5721 <<
"\". We should never get here. "
5722 "Please report this bug to the Tpetra developers.");
5757 static Teuchos::RCP<sparse_matrix_type>
5760 const Teuchos::RCP<const map_type>& rowMap,
5761 Teuchos::RCP<const map_type>& colMap,
5762 const Teuchos::RCP<const map_type>& domainMap,
5763 const Teuchos::RCP<const map_type>&
rangeMap,
5767 const bool debug=
false)
5772 using STS =
typename Teuchos::ScalarTraits<ST>;
5774 using Teuchos::ArrayRCP;
5775 using Teuchos::arcp;
5783 rowMap.is_null (), std::invalid_argument,
5784 "Row Map must be nonnull.");
5785 Teuchos::RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
5787 (comm.is_null (), std::invalid_argument,
5788 "The input row map's communicator (Teuchos::Comm object) is null.");
5790 rangeMap.is_null (), std::invalid_argument,
5791 "Range Map must be nonnull.");
5793 domainMap.is_null (), std::invalid_argument,
5794 "Domain Map must be nonnull.");
5796 domainMap->getComm().getRawPtr() != comm.getRawPtr(),
5797 std::invalid_argument,
5798 "The specified domain Map's communicator (domainMap->getComm())"
5799 "differs from row Map's communicator");
5801 rangeMap->getComm().getRawPtr() != comm.getRawPtr(),
5802 std::invalid_argument,
5803 "The specified range Map's communicator (rangeMap->getComm())"
5804 "differs from row Map's communicator");
5807 const int myRank = comm->getRank();
5808 const int numProc = comm->getSize();
5819 std::ostringstream
errMsg;
5836 using Teuchos::MatrixMarket::Banner;
5842 catch (std::exception&
e) {
5843 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
5844 "threw an exception: " <<
e.what();
5854 errMsg <<
"The Matrix Market input file must contain a "
5855 "\"coordinate\"-format sparse matrix in order to create a "
5856 "Tpetra::CrsMatrix object from it, but the file's matrix "
5857 "type is \"" <<
pBanner->matrixType() <<
"\" instead.";
5864 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
5865 "format sparse matrix in order to create a Tpetra::CrsMatrix "
5866 "object from it, but the file's matrix type is \"array\" "
5867 "instead. That probably means the file contains dense matrix "
5873 using Teuchos::MatrixMarket::readCoordinateDimensions;
5880 "# rows in file does not match rowmap.");
5882 "# rows in file does not match colmap.");
5886 typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type,global_ordinal_type>
raw_adder_type;
5888 Teuchos::RCP<raw_adder_type>
pRaw =
5893 std::cerr <<
"-- Reading matrix data" << std::endl;
5898 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
5908 catch (std::exception&
e) {
5915 typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,global_ordinal_type>
element_type;
5918 pAdder->getAdder()->merge ();
5921 const std::vector<element_type>& entries =
pAdder->getAdder()->getEntries();
5924 const size_t numEntries = (
size_t)entries.size();
5927 std::cerr <<
"----- Proc "<<
myRank<<
": Matrix has numRows=" <<
numRows
5928 <<
" rows and numEntries=" << numEntries
5929 <<
" entries." << std::endl;
5947 LO
INVALID = Teuchos::OrdinalTraits<LO>::invalid();
5957 "Current global row "<<
curRow <<
" is invalid.");
5960 "Row indices are out of order, even though they are supposed "
5961 "to be sorted. curRow = " <<
l_curRow <<
", prvRow = "
5963 "this bug to the Tpetra developers.");
5990 if(colMap.is_null()) {
5992 for(
size_t i=0;
i<rowMap->getLocalNumElements();
i++) {
5993 GO
g_row = rowMap->getGlobalElement(
i);
6002 throw std::runtime_error(
"Reading with a column map is not yet implemented");
6053 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
6117 const bool debug=
false)
6119 Teuchos::RCP<const Teuchos::Comm<int> > comm =
matrix.getComm ();
6121 (comm.is_null (), std::invalid_argument,
6122 "The input matrix's communicator (Teuchos::Comm object) is null.");
6123 const int myRank = comm->getRank ();
6139 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
6142 const bool debug=
false)
6145 (
pMatrix.is_null (), std::invalid_argument,
6146 "The input matrix is null.");
6173 const bool debug=
false)
6181 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
6182 const bool debug=
false)
6222 const bool debug=
false)
6224 using Teuchos::ArrayView;
6225 using Teuchos::Comm;
6226 using Teuchos::FancyOStream;
6227 using Teuchos::getFancyOStream;
6228 using Teuchos::null;
6230 using Teuchos::rcpFromRef;
6236 using STS =
typename Teuchos::ScalarTraits<ST>;
6242 Teuchos::SetScientific<ST>
sci (
out);
6247 comm.is_null (), std::invalid_argument,
6248 "The input matrix's communicator (Teuchos::Comm object) is null.");
6249 const int myRank = comm->getRank ();
6255 std::ostringstream
os;
6272 const global_size_t numCols = domainMap->getGlobalNumElements ();
6274 if (debug &&
myRank == 0) {
6275 std::ostringstream
os;
6276 os <<
"-- Input sparse matrix is:"
6279 << (
matrix.isGloballyIndexed() ?
"Globally" :
"Locally")
6280 <<
" indexed." <<
endl
6281 <<
"---- Its row map has " << rowMap->getGlobalNumElements ()
6282 <<
" elements." <<
endl
6283 <<
"---- Its col map has " << colMap->getGlobalNumElements ()
6284 <<
" elements." <<
endl;
6291 std::ostringstream
os;
6292 os <<
"-- " <<
myRank <<
": making gatherRowMap" <<
endl;
6296 Details::computeGatherMap (rowMap,
err, debug);
6306 Details::computeGatherMap (colMap,
err, debug);
6317 static_cast<size_t> (0)));
6326 domainMap->getIndexBase (),
6336 cerr <<
"-- Output sparse matrix is:"
6337 <<
"---- " <<
newMatrix->getRangeMap ()->getGlobalNumElements ()
6339 <<
newMatrix->getDomainMap ()->getGlobalNumElements ()
6343 << (
newMatrix->isGloballyIndexed () ?
"Globally" :
"Locally")
6344 <<
" indexed." <<
endl
6345 <<
"---- Its row map has "
6346 <<
newMatrix->getRowMap ()->getGlobalNumElements ()
6347 <<
" elements, with index base "
6349 <<
"---- Its col map has "
6350 <<
newMatrix->getColMap ()->getGlobalNumElements ()
6351 <<
" elements, with index base "
6353 <<
"---- Element count of output matrix's column Map may differ "
6354 <<
"from that of the input matrix's column Map, if some columns "
6355 <<
"of the matrix contain no entries." <<
endl;
6368 out <<
"%%MatrixMarket matrix coordinate "
6369 << (STS::isComplex ?
"complex" :
"real")
6370 <<
" general" <<
endl;
6387 out <<
newMatrix->getRangeMap ()->getGlobalNumElements () <<
" "
6388 <<
newMatrix->getDomainMap ()->getGlobalNumElements () <<
" "
6411 typename sparse_matrix_type::global_inds_host_view_type
ind;
6412 typename sparse_matrix_type::values_host_view_type
val;
6414 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6420 if (STS::isComplex) {
6430 using OTG = Teuchos::OrdinalTraits<GO>;
6439 "Failed to convert the supposed local row index "
6441 "Please report this bug to the Tpetra developers.");
6442 typename sparse_matrix_type::local_inds_host_view_type
ind;
6443 typename sparse_matrix_type::values_host_view_type
val;
6445 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6451 "On local row " <<
localRowIndex <<
" of the sparse matrix: "
6452 "Failed to convert the supposed local column index "
6453 <<
ind(
ii) <<
" into a global column index. Please report "
6454 "this bug to the Tpetra developers.");
6459 if (STS::isComplex) {
6474 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
6477 const bool debug=
false)
6480 (
pMatrix.is_null (), std::invalid_argument,
6481 "The input matrix is null.");
6520 const bool debug=
false)
6522 using Teuchos::ArrayView;
6523 using Teuchos::Comm;
6524 using Teuchos::FancyOStream;
6525 using Teuchos::getFancyOStream;
6526 using Teuchos::null;
6528 using Teuchos::rcpFromRef;
6538 auto rowMap =
graph.getRowMap ();
6539 if (rowMap.is_null ()) {
6542 auto comm = rowMap->getComm ();
6543 if (comm.is_null ()) {
6546 const int myRank = comm->getRank ();
6552 std::ostringstream
os;
6564 auto colMap =
graph.getColMap ();
6565 auto domainMap =
graph.getDomainMap ();
6569 const global_size_t numCols = domainMap->getGlobalNumElements ();
6571 if (debug &&
myRank == 0) {
6572 std::ostringstream
os;
6573 os <<
"-- Input sparse graph is:"
6574 <<
"---- " <<
numRows <<
" x " << numCols <<
" with "
6575 <<
graph.getGlobalNumEntries () <<
" entries;" <<
endl
6577 << (
graph.isGloballyIndexed () ?
"Globally" :
"Locally")
6578 <<
" indexed." <<
endl
6579 <<
"---- Its row Map has " << rowMap->getGlobalNumElements ()
6580 <<
" elements." <<
endl
6581 <<
"---- Its col Map has " << colMap->getGlobalNumElements ()
6582 <<
" elements." <<
endl;
6589 std::ostringstream
os;
6590 os <<
"-- " <<
myRank <<
": making gatherRowMap" <<
endl;
6611 static_cast<size_t> (0));
6620 domainMap->getIndexBase (),
6630 cerr <<
"-- Output sparse graph is:"
6631 <<
"---- " <<
newGraph.getRangeMap ()->getGlobalNumElements ()
6633 <<
newGraph.getDomainMap ()->getGlobalNumElements ()
6635 <<
newGraph.getGlobalNumEntries () <<
" entries;" <<
endl
6637 << (
newGraph.isGloballyIndexed () ?
"Globally" :
"Locally")
6638 <<
" indexed." <<
endl
6639 <<
"---- Its row map has "
6640 <<
newGraph.getRowMap ()->getGlobalNumElements ()
6641 <<
" elements, with index base "
6642 <<
newGraph.getRowMap ()->getIndexBase () <<
"." <<
endl
6643 <<
"---- Its col map has "
6644 <<
newGraph.getColMap ()->getGlobalNumElements ()
6645 <<
" elements, with index base "
6646 <<
newGraph.getColMap ()->getIndexBase () <<
"." <<
endl
6647 <<
"---- Element count of output graph's column Map may differ "
6648 <<
"from that of the input matrix's column Map, if some columns "
6649 <<
"of the matrix contain no entries." <<
endl;
6663 out <<
"%%MatrixMarket matrix coordinate pattern general" <<
endl;
6681 out <<
newGraph.getRangeMap ()->getGlobalNumElements () <<
" "
6682 <<
newGraph.getDomainMap ()->getGlobalNumElements () <<
" "
6697 if (
newGraph.isGloballyIndexed ()) {
6705 typename crs_graph_type::global_inds_host_view_type
ind;
6707 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6718 typedef Teuchos::OrdinalTraits<GO>
OTG;
6727 "to convert the supposed local row index " <<
localRowIndex <<
6728 " into a global row index. Please report this bug to the "
6729 "Tpetra developers.");
6730 typename crs_graph_type::local_inds_host_view_type
ind;
6732 for (
size_t ii = 0;
ii <
ind.extent(0);
ii++) {
6738 "On local row " <<
localRowIndex <<
" of the sparse graph: "
6739 "Failed to convert the supposed local column index "
6740 <<
ind(
ii) <<
" into a global column index. Please report "
6741 "this bug to the Tpetra developers.");
6761 const bool debug=
false)
6805 const bool debug=
false)
6807 auto comm =
graph.getComm ();
6808 if (comm.is_null ()) {
6816 const int myRank = comm->getRank ();
6836 const bool debug=
false)
6851 const Teuchos::RCP<const crs_graph_type>&
pGraph,
6854 const bool debug=
false)
6870 const Teuchos::RCP<const crs_graph_type>&
pGraph,
6871 const bool debug=
false)
6901 const bool debug=
false)
6909 const Teuchos::RCP<const sparse_matrix_type>&
pMatrix,
6910 const bool debug=
false)
6948 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6949 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
6951 const int myRank =
X.getMap ().is_null () ? 0 :
6952 (
X.getMap ()->getComm ().is_null () ? 0 :
6953 X.getMap ()->getComm ()->getRank ());
6973 const Teuchos::RCP<const multivector_type>&
X,
6976 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6977 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
6980 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
6981 "writeDenseFile: The input MultiVector X is null.");
6993 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
6994 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
7006 const Teuchos::RCP<const multivector_type>&
X,
7007 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7008 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
7011 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7012 "writeDenseFile: The input MultiVector X is null.");
7052 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7053 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
7055 using Teuchos::Comm;
7056 using Teuchos::outArg;
7057 using Teuchos::REDUCE_MAX;
7058 using Teuchos::reduceAll;
7063 Teuchos::null :
X.getMap ()->getComm ();
7064 const int myRank = comm.is_null () ? 0 : comm->getRank ();
7069 const bool debug = !
dbg.is_null ();
7072 std::ostringstream
os;
7082 const size_t numVecs =
X.getNumVectors ();
7084 writeDenseColumn (
out, * (
X.getVector (
j)),
err,
dbg);
7089 std::ostringstream
os;
7124 writeDenseHeader (std::ostream&
out,
7128 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7129 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
7131 using Teuchos::Comm;
7132 using Teuchos::outArg;
7134 using Teuchos::REDUCE_MAX;
7135 using Teuchos::reduceAll;
7137 typedef Teuchos::ScalarTraits<scalar_type> STS;
7138 const char prefix[] =
"Tpetra::MatrixMarket::writeDenseHeader: ";
7141 Teuchos::null :
X.getMap ()->getComm ();
7142 const int myRank = comm.is_null () ? 0 : comm->getRank ();
7149 const bool debug = !
dbg.is_null ();
7153 std::ostringstream
os;
7173 std::ostringstream
hdr;
7176 if (STS::isComplex) {
7178 }
else if (STS::isOrdinal) {
7179 dataType =
"integer";
7183 hdr <<
"%%MatrixMarket matrix array " << dataType <<
" general"
7188 if (matrixName !=
"") {
7189 printAsComment (hdr, matrixName);
7191 if (matrixDescription !=
"") {
7192 printAsComment (hdr, matrixDescription);
7195 hdr << X.getGlobalLength () <<
" " << X.getNumVectors () << endl;
7199 }
catch (std::exception& e) {
7200 if (! err.is_null ()) {
7201 *err << prefix <<
"While writing the Matrix Market header, "
7202 "Process 0 threw an exception: " << e.what () << endl;
7211 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
7212 TEUCHOS_TEST_FOR_EXCEPTION(
7213 gblErr == 1, std::runtime_error, prefix <<
"Some error occurred "
7214 "which prevented this method from completing.");
7218 *dbg << myRank <<
": writeDenseHeader: Done" << endl;
7241 writeDenseColumn (std::ostream& out,
7243 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7244 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7246 using Teuchos::arcp;
7247 using Teuchos::Array;
7248 using Teuchos::ArrayRCP;
7249 using Teuchos::ArrayView;
7250 using Teuchos::Comm;
7251 using Teuchos::CommRequest;
7252 using Teuchos::ireceive;
7253 using Teuchos::isend;
7254 using Teuchos::outArg;
7255 using Teuchos::REDUCE_MAX;
7256 using Teuchos::reduceAll;
7258 using Teuchos::TypeNameTraits;
7259 using Teuchos::wait;
7261 typedef Teuchos::ScalarTraits<scalar_type> STS;
7263 const Comm<int>& comm = * (X.getMap ()->getComm ());
7264 const int myRank = comm.getRank ();
7265 const int numProcs = comm.getSize ();
7272 const bool debug = ! dbg.is_null ();
7276 std::ostringstream os;
7277 os << myRank <<
": writeDenseColumn" << endl;
7286 Teuchos::SetScientific<scalar_type> sci (out);
7288 const size_t myNumRows = X.getLocalLength ();
7289 const size_t numCols = X.getNumVectors ();
7292 const int sizeTag = 1337;
7293 const int dataTag = 1338;
7354 RCP<CommRequest<int> > sendReqSize, sendReqData;
7360 Array<ArrayRCP<size_t> > recvSizeBufs (3);
7361 Array<ArrayRCP<scalar_type> > recvDataBufs (3);
7362 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7363 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7366 ArrayRCP<size_t> sendDataSize (1);
7367 sendDataSize[0] = myNumRows;
7371 std::ostringstream os;
7372 os << myRank <<
": Post receive-size receives from "
7373 "Procs 1 and 2: tag = " << sizeTag << endl;
7377 recvSizeBufs[0].resize (1);
7381 (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7382 recvSizeBufs[1].resize (1);
7383 (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7384 recvSizeBufs[2].resize (1);
7385 (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7388 ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
7392 ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
7395 else if (myRank == 1 || myRank == 2) {
7397 std::ostringstream os;
7398 os << myRank <<
": Post send-size send: size = "
7399 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7406 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7410 std::ostringstream os;
7411 os << myRank <<
": Not posting my send-size send yet" << endl;
7420 std::ostringstream os;
7421 os << myRank <<
": Pack my entries" << endl;
7424 ArrayRCP<scalar_type> sendDataBuf;
7426 sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
7427 X.get1dCopy (sendDataBuf (), myNumRows);
7429 catch (std::exception& e) {
7431 if (! err.is_null ()) {
7432 std::ostringstream os;
7433 os <<
"Process " << myRank <<
": Attempt to pack my MultiVector "
7434 "entries threw an exception: " << e.what () << endl;
7439 std::ostringstream os;
7440 os << myRank <<
": Done packing my entries" << endl;
7449 *dbg << myRank <<
": Post send-data send: tag = " << dataTag
7452 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7460 std::ostringstream os;
7461 os << myRank <<
": Write my entries" << endl;
7467 const size_t printNumRows = myNumRows;
7468 ArrayView<const scalar_type> printData = sendDataBuf ();
7469 const size_t printStride = printNumRows;
7470 if (
static_cast<size_t> (printData.size ()) < printStride * numCols) {
7472 if (! err.is_null ()) {
7473 std::ostringstream os;
7474 os <<
"Process " << myRank <<
": My MultiVector data's size "
7475 << printData.size () <<
" does not match my local dimensions "
7476 << printStride <<
" x " << numCols <<
"." << endl;
7484 for (
size_t col = 0; col < numCols; ++col) {
7485 for (
size_t row = 0; row < printNumRows; ++row) {
7486 if (STS::isComplex) {
7487 out << STS::real (printData[row + col * printStride]) <<
" "
7488 << STS::imag (printData[row + col * printStride]) << endl;
7490 out << printData[row + col * printStride] << endl;
7499 const int recvRank = 1;
7500 const int circBufInd = recvRank % 3;
7502 std::ostringstream os;
7503 os << myRank <<
": Wait on receive-size receive from Process "
7504 << recvRank << endl;
7508 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7512 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7513 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7515 if (! err.is_null ()) {
7516 std::ostringstream os;
7517 os << myRank <<
": Result of receive-size receive from Process "
7518 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7519 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7520 "This should never happen, and suggests that the receive never "
7521 "got posted. Please report this bug to the Tpetra developers."
7536 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7538 std::ostringstream os;
7539 os << myRank <<
": Post receive-data receive from Process "
7540 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7541 << recvDataBufs[circBufInd].size () << endl;
7544 if (! recvSizeReqs[circBufInd].is_null ()) {
7546 if (! err.is_null ()) {
7547 std::ostringstream os;
7548 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7549 "null, before posting the receive-data receive from Process "
7550 << recvRank <<
". This should never happen. Please report "
7551 "this bug to the Tpetra developers." << endl;
7555 recvDataReqs[circBufInd] =
7556 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7557 recvRank, dataTag, comm);
7560 else if (myRank == 1) {
7563 std::ostringstream os;
7564 os << myRank <<
": Wait on my send-size send" << endl;
7567 wait<int> (comm, outArg (sendReqSize));
7573 for (
int p = 1; p < numProcs; ++p) {
7575 if (p + 2 < numProcs) {
7577 const int recvRank = p + 2;
7578 const int circBufInd = recvRank % 3;
7580 std::ostringstream os;
7581 os << myRank <<
": Post receive-size receive from Process "
7582 << recvRank <<
": tag = " << sizeTag << endl;
7585 if (! recvSizeReqs[circBufInd].is_null ()) {
7587 if (! err.is_null ()) {
7588 std::ostringstream os;
7589 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7590 <<
"null, for the receive-size receive from Process "
7591 << recvRank <<
"! This may mean that this process never "
7592 <<
"finished waiting for the receive from Process "
7593 << (recvRank - 3) <<
"." << endl;
7597 recvSizeReqs[circBufInd] =
7598 ireceive<int, size_t> (recvSizeBufs[circBufInd],
7599 recvRank, sizeTag, comm);
7602 if (p + 1 < numProcs) {
7603 const int recvRank = p + 1;
7604 const int circBufInd = recvRank % 3;
7608 std::ostringstream os;
7609 os << myRank <<
": Wait on receive-size receive from Process "
7610 << recvRank << endl;
7613 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7617 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7618 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7620 if (! err.is_null ()) {
7621 std::ostringstream os;
7622 os << myRank <<
": Result of receive-size receive from Process "
7623 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7624 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7625 "This should never happen, and suggests that the receive never "
7626 "got posted. Please report this bug to the Tpetra developers."
7640 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7642 std::ostringstream os;
7643 os << myRank <<
": Post receive-data receive from Process "
7644 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7645 << recvDataBufs[circBufInd].size () << endl;
7648 if (! recvDataReqs[circBufInd].is_null ()) {
7650 if (! err.is_null ()) {
7651 std::ostringstream os;
7652 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
7653 <<
"null, for the receive-data receive from Process "
7654 << recvRank <<
"! This may mean that this process never "
7655 <<
"finished waiting for the receive from Process "
7656 << (recvRank - 3) <<
"." << endl;
7660 recvDataReqs[circBufInd] =
7661 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7662 recvRank, dataTag, comm);
7666 const int recvRank = p;
7667 const int circBufInd = recvRank % 3;
7669 std::ostringstream os;
7670 os << myRank <<
": Wait on receive-data receive from Process "
7671 << recvRank << endl;
7674 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
7681 std::ostringstream os;
7682 os << myRank <<
": Write entries from Process " << recvRank
7684 *dbg << os.str () << endl;
7686 size_t printNumRows = (recvSizeBufs[circBufInd])[0];
7687 if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7689 if (! err.is_null ()) {
7690 std::ostringstream os;
7691 os << myRank <<
": Result of receive-size receive from Process "
7692 << recvRank <<
" was Teuchos::OrdinalTraits<size_t>::"
7693 "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
7694 <<
". This should never happen, and suggests that its "
7695 "receive-size receive was never posted. "
7696 "Please report this bug to the Tpetra developers." << endl;
7703 if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
7705 if (! err.is_null ()) {
7706 std::ostringstream os;
7707 os << myRank <<
": Result of receive-size receive from Proc "
7708 << recvRank <<
" was " << printNumRows <<
" > 0, but "
7709 "recvDataBufs[" << circBufInd <<
"] is null. This should "
7710 "never happen. Please report this bug to the Tpetra "
7711 "developers." << endl;
7721 ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
7722 const size_t printStride = printNumRows;
7726 for (
size_t col = 0; col < numCols; ++col) {
7727 for (
size_t row = 0; row < printNumRows; ++row) {
7728 if (STS::isComplex) {
7729 out << STS::real (printData[row + col * printStride]) <<
" "
7730 << STS::imag (printData[row + col * printStride]) << endl;
7732 out << printData[row + col * printStride] << endl;
7737 else if (myRank == p) {
7740 std::ostringstream os;
7741 os << myRank <<
": Wait on my send-data send" << endl;
7744 wait<int> (comm, outArg (sendReqData));
7746 else if (myRank == p + 1) {
7749 std::ostringstream os;
7750 os << myRank <<
": Post send-data send: tag = " << dataTag
7754 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7757 std::ostringstream os;
7758 os << myRank <<
": Wait on my send-size send" << endl;
7761 wait<int> (comm, outArg (sendReqSize));
7763 else if (myRank == p + 2) {
7766 std::ostringstream os;
7767 os << myRank <<
": Post send-size send: size = "
7768 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7771 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7776 reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
7777 TEUCHOS_TEST_FOR_EXCEPTION(
7778 gblErr == 1, std::runtime_error,
"Tpetra::MatrixMarket::writeDense "
7779 "experienced some kind of error and was unable to complete.");
7783 *dbg << myRank <<
": writeDenseColumn: Done" << endl;
7797 const Teuchos::RCP<const multivector_type>&
X,
7800 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7801 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
7804 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7805 "writeDense: The input MultiVector X is null.");
7817 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7818 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
7830 const Teuchos::RCP<const multivector_type>&
X,
7831 const Teuchos::RCP<Teuchos::FancyOStream>&
err = Teuchos::null,
7832 const Teuchos::RCP<Teuchos::FancyOStream>&
dbg = Teuchos::null)
7835 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7836 "writeDense: The input MultiVector X is null.");
7862 Teuchos::RCP<Teuchos::FancyOStream>
err =
7863 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
7878 const Teuchos::RCP<Teuchos::FancyOStream>&
err,
7879 const bool debug=
false)
7881 using Teuchos::Array;
7882 using Teuchos::ArrayRCP;
7883 using Teuchos::ArrayView;
7884 using Teuchos::Comm;
7885 using Teuchos::CommRequest;
7886 using Teuchos::ireceive;
7887 using Teuchos::isend;
7889 using Teuchos::TypeNameTraits;
7890 using Teuchos::wait;
7910 sizeof (GO) >
sizeof (
int_type), std::logic_error,
7912 <<
" is too big for ptrdiff_t. sizeof(GO) = " <<
sizeof (GO)
7913 <<
" > sizeof(ptrdiff_t) = " <<
sizeof (
ptrdiff_t) <<
".");
7916 "The (MPI) process rank type pid_type=" <<
7918 "sizeof(pid_type) = " <<
sizeof (
pid_type) <<
" > sizeof(ptrdiff_t)"
7922 const int myRank = comm.getRank ();
7923 const int numProcs = comm.getSize ();
7925 if (!
err.is_null ()) {
7929 std::ostringstream
os;
7933 if (!
err.is_null ()) {
8017 std::ostringstream
os;
8018 os <<
myRank <<
": Post receive-size receives from "
8040 std::ostringstream
os;
8041 os <<
myRank <<
": Post send-size send: size = "
8053 std::ostringstream
os;
8054 os <<
myRank <<
": Not posting my send-size send yet" <<
endl;
8065 std::ostringstream
os;
8072 if (
map.isContiguous ()) {
8093 std::ostringstream
os;
8094 os <<
myRank <<
": Done packing my GIDs and PIDs" <<
endl;
8114 std::ostringstream
hdr;
8118 hdr <<
"%%MatrixMarket matrix array integer general" <<
endl
8119 <<
"% Format: Version 2.0" <<
endl
8121 <<
"% This file encodes a Tpetra::Map." <<
endl
8122 <<
"% It is stored as a dense vector, with twice as many " <<
endl
8123 <<
"% entries as the global number of GIDs (global indices)." <<
endl
8124 <<
"% (GID, PID) pairs are stored contiguously, where the PID " <<
endl
8125 <<
"% is the rank of the process owning that GID." <<
endl
8126 << (2 *
map.getGlobalNumElements ()) <<
" " << 1 <<
endl;
8130 std::ostringstream
os;
8151 std::ostringstream
os;
8152 os <<
myRank <<
": Wait on receive-size receive from Process "
8163 std::ostringstream
os;
8164 os <<
myRank <<
": Result of receive-size receive from Process "
8165 <<
recvRank <<
" is -1. This should never happen, and "
8166 "suggests that the receive never got posted. Please report "
8167 "this bug to the Tpetra developers." <<
endl;
8174 std::ostringstream
os;
8175 os <<
myRank <<
": Post receive-data receive from Process "
8181 std::ostringstream
os;
8183 "null, before posting the receive-data receive from Process "
8184 <<
recvRank <<
". This should never happen. Please report "
8185 "this bug to the Tpetra developers." <<
endl;
8196 std::ostringstream
os;
8213 std::ostringstream
os;
8214 os <<
myRank <<
": Post receive-size receive from Process "
8219 std::ostringstream
os;
8221 <<
"null, for the receive-size receive from Process "
8222 <<
recvRank <<
"! This may mean that this process never "
8223 <<
"finished waiting for the receive from Process "
8238 std::ostringstream
os;
8239 os <<
myRank <<
": Wait on receive-size receive from Process "
8249 std::ostringstream
os;
8250 os <<
myRank <<
": Result of receive-size receive from Process "
8251 <<
recvRank <<
" is -1. This should never happen, and "
8252 "suggests that the receive never got posted. Please report "
8253 "this bug to the Tpetra developers." <<
endl;
8260 std::ostringstream
os;
8261 os <<
myRank <<
": Post receive-data receive from Process "
8267 std::ostringstream
os;
8269 <<
"null, for the receive-data receive from Process "
8270 <<
recvRank <<
"! This may mean that this process never "
8271 <<
"finished waiting for the receive from Process "
8284 std::ostringstream
os;
8285 os <<
myRank <<
": Wait on receive-data receive from Process "
8296 std::ostringstream
os;
8297 os <<
myRank <<
": Write GIDs and PIDs from Process "
8303 std::ostringstream
os;
8304 os <<
myRank <<
": Result of receive-size receive from Process "
8305 <<
recvRank <<
" was -1. This should never happen, and "
8306 "suggests that its receive-size receive was never posted. "
8307 "Please report this bug to the Tpetra developers." <<
endl;
8311 std::ostringstream
os;
8312 os <<
myRank <<
": Result of receive-size receive from Proc "
8314 "recvDataBufs[" <<
circBufInd <<
"] is null. This should "
8315 "never happen. Please report this bug to the Tpetra "
8316 "developers." <<
endl;
8329 std::ostringstream
os;
8338 std::ostringstream
os;
8346 std::ostringstream
os;
8355 std::ostringstream
os;
8356 os <<
myRank <<
": Post send-size send: size = "
8364 if (!
err.is_null ()) {
8370 if (!
err.is_null ()) {
8380 const int myRank =
map.getComm ()->getRank ();
8416 printAsComment (std::ostream&
out,
const std::string&
str)
8423 if (!
line.empty()) {
8426 if (
line[0] ==
'%') {
8458 Teuchos::ParameterList
pl;
8484 Teuchos::ParameterList
pl;
8527 const Teuchos::ParameterList&
params)
8531 const int myRank =
A.getDomainMap()->getComm()->getRank();
8544 "writeOperator: temporary file " <<
tmpFile <<
" already exists");
8546 if (
params.isParameter(
"precision")) {
8560 if (
params.isParameter(
"print MatrixMarket header"))
8567 std::ifstream src(
tmpFile, std::ios_base::binary);
8616 const Teuchos::ParameterList&
params)
8618 const int myRank =
A.getDomainMap ()->getComm ()->getRank ();
8635 std::ostringstream
tmpOut;
8637 if (
params.isParameter (
"precision") &&
params.isType<
int> (
"precision")) {
8646 if (
params.isParameter (
"print MatrixMarket header") &&
8647 params.isType<
bool> (
"print MatrixMarket header")) {
8680 writeOperatorImpl (std::ostream&
os,
8681 const operator_type&
A,
8682 const Teuchos::ParameterList&
params)
8686 using Teuchos::ArrayRCP;
8687 using Teuchos::Array;
8692 typedef Teuchos::OrdinalTraits<LO>
TLOT;
8693 typedef Teuchos::OrdinalTraits<GO>
TGOT;
8697 const map_type& domainMap = *(
A.getDomainMap());
8700 const int myRank = comm->getRank();
8701 const size_t numProcs = comm->getSize();
8704 if (
params.isParameter(
"probing size"))
8708 GO
minColGid = domainMap.getMinAllGlobalIndex();
8709 GO
maxColGid = domainMap.getMaxAllGlobalIndex();
8720 if (
params.isParameter(
"zero-based indexing")) {
8721 if (
params.get<
bool>(
"zero-based indexing") ==
true)
8747 importGidList.reserve(numTargetMapEntries);
8749 RCP<map_type> importGidMap = rcp(
new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
8752 import_type gidImporter(allGidsMap, importGidMap);
8753 mv_type_go importedGids(importGidMap, 1);
8754 importedGids.doImport(allGids, gidImporter,
INSERT);
8757 ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
8758 RCP<const map_type> importMap = rcp(
new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm) );
8761 import_type importer(rangeMap, importMap);
8763 RCP<mv_type> colsOnPid0 = rcp(
new mv_type(importMap,numMVs));
8765 RCP<mv_type> ei = rcp(
new mv_type(A.getDomainMap(),numMVs));
8766 RCP<mv_type> colsA = rcp(
new mv_type(A.getRangeMap(),numMVs));
8768 Array<GO> globalColsArray, localColsArray;
8769 globalColsArray.reserve(numMVs);
8770 localColsArray.reserve(numMVs);
8772 ArrayRCP<ArrayRCP<Scalar> > eiData(numMVs);
8773 for (
size_t i=0; i<numMVs; ++i)
8774 eiData[i] = ei->getDataNonConst(i);
8779 for (GO k=0; k<numChunks; ++k) {
8780 for (
size_t j=0; j<numMVs; ++j ) {
8782 GO curGlobalCol = minColGid + k*numMVs + j;
8783 globalColsArray.push_back(curGlobalCol);
8785 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8786 if (curLocalCol != TLOT::invalid()) {
8787 eiData[j][curLocalCol] = TGOT::one();
8788 localColsArray.push_back(curLocalCol);
8794 A.apply(*ei,*colsA);
8796 colsOnPid0->doImport(*colsA,importer,
INSERT);
8799 globalNnz += writeColumns(os,*colsOnPid0, numMVs, importedGidsData(),
8800 globalColsArray, offsetToUseInPrinting);
8803 for (
size_t j=0; j<numMVs; ++j ) {
8804 for (
int i=0; i<localColsArray.size(); ++i)
8805 eiData[j][localColsArray[i]] = TGOT::zero();
8807 globalColsArray.clear();
8808 localColsArray.clear();
8816 for (
int j=0; j<rem; ++j ) {
8817 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8818 globalColsArray.push_back(curGlobalCol);
8820 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8821 if (curLocalCol != TLOT::invalid()) {
8822 eiData[j][curLocalCol] = TGOT::one();
8823 localColsArray.push_back(curLocalCol);
8829 A.apply(*ei,*colsA);
8831 colsOnPid0->doImport(*colsA,importer,
INSERT);
8833 globalNnz += writeColumns(os,*colsOnPid0, rem, importedGidsData(),
8834 globalColsArray, offsetToUseInPrinting);
8837 for (
int j=0; j<rem; ++j ) {
8838 for (
int i=0; i<localColsArray.size(); ++i)
8839 eiData[j][localColsArray[i]] = TGOT::zero();
8841 globalColsArray.clear();
8842 localColsArray.clear();
8851 std::ostringstream oss;
8853 oss <<
"%%MatrixMarket matrix coordinate ";
8854 if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
8859 oss <<
" general" << std::endl;
8860 oss <<
"% Tpetra::Operator" << std::endl;
8861 std::time_t now = std::time(NULL);
8862 oss <<
"% time stamp: " << ctime(&now);
8863 oss <<
"% written from " << numProcs <<
" processes" << std::endl;
8864 size_t numRows = rangeMap->getGlobalNumElements();
8865 size_t numCols = domainMap.getGlobalNumElements();
8866 oss << numRows <<
" " << numCols <<
" " << globalNnz << std::endl;
8873 writeColumns(std::ostream& os, mv_type
const &colsA,
size_t const &numCols,
8874 Teuchos::ArrayView<const global_ordinal_type>
const &rowGids,
8875 Teuchos::Array<global_ordinal_type>
const &colsArray,
8880 typedef Teuchos::ScalarTraits<Scalar> STS;
8883 const Scalar zero = STS::zero();
8884 const size_t numRows = colsA.getGlobalLength();
8885 for (
size_t j=0; j<numCols; ++j) {
8886 Teuchos::ArrayRCP<const Scalar>
const curCol = colsA.getData(j);
8887 const GO J = colsArray[j];
8888 for (
size_t i=0; i<numRows; ++i) {
8889 const Scalar val = curCol[i];
8891 os << rowGids[i]+indexBase <<
" " << J+indexBase <<
" " << val << std::endl;
8919 const bool debug=
false) {
8924 using STS =
typename Teuchos::ScalarTraits<ST>;
8928 Teuchos::RCP<const Teuchos::Comm<int> > comm =
matrix.getComm ();
8930 (comm.is_null (), std::invalid_argument,
8931 "The input matrix's communicator (Teuchos::Comm object) is null.");
8933 (
matrix.isGloballyIndexed() || !
matrix.isFillComplete(), std::invalid_argument,
8934 "The input matrix must not be GloballyIndexed and must be fillComplete.");
8937 const int myRank = comm->getRank();
8938 const int numProc = comm->getSize();
8960 out <<
"%%MatrixMarket matrix coordinate "
8961 << (STS::isComplex ?
"complex" :
"real")
8962 <<
" general" << std::endl;
8982 Teuchos::SetScientific<ST>
sci (
out);
8987 typename sparse_matrix_type::local_inds_host_view_type indices;
8988 typename sparse_matrix_type::values_host_view_type
values;
8990 for (
size_t ii = 0;
ii < indices.extent(0);
ii++) {
8991 const GO
g_col = colMap->getGlobalElement(indices(
ii));
8996 if (STS::isComplex) {