42 #ifndef __MatrixMarket_Tpetra_hpp 43 #define __MatrixMarket_Tpetra_hpp 56 #include "Tpetra_ConfigDefs.hpp" 57 #include "Tpetra_CrsMatrix.hpp" 58 #include "Tpetra_Operator.hpp" 59 #include "Tpetra_Vector.hpp" 61 #include "Teuchos_MatrixMarket_Raw_Adder.hpp" 62 #include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp" 63 #include "Teuchos_MatrixMarket_assignScalar.hpp" 64 #include "Teuchos_MatrixMarket_Banner.hpp" 65 #include "Teuchos_MatrixMarket_CoordDataReader.hpp" 66 #include "Teuchos_MatrixMarket_SetScientific.hpp" 161 template<
class SparseMatrixType>
166 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
181 typedef typename SparseMatrixType::global_ordinal_type
198 typedef Teuchos::Comm<int> comm_type;
202 typedef Teuchos::RCP<const comm_type> comm_ptr;
203 typedef Teuchos::RCP<const map_type> map_ptr;
204 typedef Teuchos::RCP<node_type> node_ptr;
212 typedef Teuchos::ArrayRCP<int>::size_type size_type;
225 static Teuchos::RCP<const map_type>
226 makeRangeMap (
const Teuchos::RCP<const comm_type>& pComm,
227 const Teuchos::RCP<node_type>& pNode,
228 const global_ordinal_type numRows)
231 if (pNode.is_null ()) {
232 return rcp (
new map_type (static_cast<global_size_t> (numRows),
233 static_cast<global_ordinal_type> (0),
234 pComm, GloballyDistributed));
237 return rcp (
new map_type (static_cast<global_size_t> (numRows),
238 static_cast<global_ordinal_type> (0),
239 pComm, GloballyDistributed, pNode));
271 static RCP<const map_type>
272 makeRowMap (
const RCP<const map_type>& pRowMap,
273 const RCP<const comm_type>& pComm,
274 const RCP<node_type>& pNode,
275 const global_ordinal_type numRows)
279 if (pRowMap.is_null ()) {
280 if (pNode.is_null ()) {
281 return rcp (
new map_type (static_cast<global_size_t> (numRows),
282 static_cast<global_ordinal_type> (0),
283 pComm, GloballyDistributed));
286 return rcp (
new map_type (static_cast<global_size_t> (numRows),
287 static_cast<global_ordinal_type> (0),
288 pComm, GloballyDistributed, pNode));
292 TEUCHOS_TEST_FOR_EXCEPTION
293 (! pRowMap->isDistributed () && pComm->getSize () > 1,
294 std::invalid_argument,
"The specified row map is not distributed, " 295 "but the given communicator includes more than one process (in " 296 "fact, there are " << pComm->getSize () <<
" processes).");
297 TEUCHOS_TEST_FOR_EXCEPTION
298 (pRowMap->getComm () != pComm, std::invalid_argument,
299 "The specified row Map's communicator (pRowMap->getComm()) " 300 "differs from the given separately supplied communicator pComm.");
301 TEUCHOS_TEST_FOR_EXCEPTION
302 (pRowMap->getNode () != pNode, std::invalid_argument,
303 "The specified row Map's node (pRowMap->getNode()) differs from " 304 "the given separately supplied node pNode.");
323 static Teuchos::RCP<const map_type>
324 makeDomainMap (
const Teuchos::RCP<const map_type>& pRangeMap,
325 const global_ordinal_type numRows,
326 const global_ordinal_type numCols)
329 typedef local_ordinal_type LO;
330 typedef global_ordinal_type GO;
331 typedef node_type NT;
333 if (numRows == numCols) {
336 return createUniformContigMapWithNode<LO,GO,NT> (numCols,
337 pRangeMap->getComm (),
338 pRangeMap->getNode ());
415 distribute (ArrayRCP<size_t>& myNumEntriesPerRow,
416 ArrayRCP<size_t>& myRowPtr,
417 ArrayRCP<global_ordinal_type>& myColInd,
418 ArrayRCP<scalar_type>& myValues,
419 const RCP<const map_type>& pRowMap,
420 ArrayRCP<size_t>& numEntriesPerRow,
421 ArrayRCP<size_t>& rowPtr,
422 ArrayRCP<global_ordinal_type>& colInd,
423 ArrayRCP<scalar_type>& values,
424 const bool debug=
false)
428 using Teuchos::CommRequest;
430 using Teuchos::receive;
435 const bool extraDebug =
false;
436 RCP<const comm_type> pComm = pRowMap->getComm ();
437 const int numProcs = pComm->getSize ();
438 const int myRank = pComm->getRank ();
439 const int rootRank = 0;
442 typedef global_ordinal_type GO;
446 ArrayView<const GO> myRows = pRowMap->getNodeElementList();
447 const size_type myNumRows = myRows.size();
448 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
449 pRowMap->getNodeNumElements(),
451 "pRowMap->getNodeElementList().size() = " 453 <<
" != pRowMap->getNodeNumElements() = " 454 << pRowMap->getNodeNumElements() <<
". " 455 "Please report this bug to the Tpetra developers.");
456 TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
458 "On Proc 0: numEntriesPerRow.size() = " 459 << numEntriesPerRow.size()
460 <<
" != pRowMap->getNodeElementList().size() = " 461 << myNumRows <<
". Please report this bug to the " 462 "Tpetra developers.");
466 myNumEntriesPerRow = arcp<size_t> (myNumRows);
468 if (myRank != rootRank) {
472 send (*pComm, myNumRows, rootRank);
473 if (myNumRows != 0) {
477 send (*pComm, static_cast<int> (myNumRows),
478 myRows.getRawPtr(), rootRank);
488 receive (*pComm, rootRank,
489 static_cast<int> (myNumRows),
490 myNumEntriesPerRow.getRawPtr());
494 const local_ordinal_type myNumEntries =
495 std::accumulate (myNumEntriesPerRow.begin(),
496 myNumEntriesPerRow.end(), 0);
502 myColInd = arcp<GO> (myNumEntries);
503 myValues = arcp<scalar_type> (myNumEntries);
504 if (myNumEntries > 0) {
507 receive (*pComm, rootRank,
508 static_cast<int> (myNumEntries),
509 myColInd.getRawPtr());
510 receive (*pComm, rootRank,
511 static_cast<int> (myNumEntries),
512 myValues.getRawPtr());
518 cerr <<
"-- Proc 0: Copying my data from global arrays" << endl;
522 for (size_type k = 0; k < myNumRows; ++k) {
523 const GO myCurRow = myRows[k];
524 const local_ordinal_type numEntriesInThisRow = numEntriesPerRow[myCurRow];
525 myNumEntriesPerRow[k] = numEntriesInThisRow;
527 if (extraDebug && debug) {
528 cerr <<
"Proc " << pRowMap->getComm ()->getRank ()
529 <<
": myNumEntriesPerRow[0.." << (myNumRows-1) <<
"] = [";
530 for (size_type k = 0; k < myNumRows; ++k) {
531 cerr << myNumEntriesPerRow[k];
532 if (k < myNumRows-1) {
539 const local_ordinal_type myNumEntries =
540 std::accumulate (myNumEntriesPerRow.begin(),
541 myNumEntriesPerRow.end(), 0);
543 cerr <<
"-- Proc 0: I own " << myNumRows <<
" rows and " 544 << myNumEntries <<
" entries" << endl;
546 myColInd = arcp<GO> (myNumEntries);
547 myValues = arcp<scalar_type> (myNumEntries);
554 local_ordinal_type myCurPos = 0;
555 for (size_type k = 0; k < myNumRows;
556 myCurPos += myNumEntriesPerRow[k], ++k) {
557 const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
558 const GO myRow = myRows[k];
559 const size_t curPos = rowPtr[myRow];
562 if (curNumEntries > 0) {
563 ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
564 ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
565 std::copy (colIndView.begin(), colIndView.end(),
566 myColIndView.begin());
568 ArrayView<scalar_type> valuesView =
569 values (curPos, curNumEntries);
570 ArrayView<scalar_type> myValuesView =
571 myValues (myCurPos, curNumEntries);
572 std::copy (valuesView.begin(), valuesView.end(),
573 myValuesView.begin());
578 for (
int p = 1; p < numProcs; ++p) {
580 cerr <<
"-- Proc 0: Processing proc " << p << endl;
583 size_type theirNumRows = 0;
588 receive (*pComm, p, &theirNumRows);
590 cerr <<
"-- Proc 0: Proc " << p <<
" owns " 591 << theirNumRows <<
" rows" << endl;
593 if (theirNumRows != 0) {
598 ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
599 receive (*pComm, p, as<int> (theirNumRows),
600 theirRows.getRawPtr ());
609 const global_size_t numRows = pRowMap->getGlobalNumElements ();
610 const GO indexBase = pRowMap->getIndexBase ();
611 bool theirRowsValid =
true;
612 for (size_type k = 0; k < theirNumRows; ++k) {
613 if (theirRows[k] < indexBase ||
614 as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
615 theirRowsValid =
false;
618 if (! theirRowsValid) {
619 TEUCHOS_TEST_FOR_EXCEPTION(
620 ! theirRowsValid, std::logic_error,
621 "Proc " << p <<
" has at least one invalid row index. " 622 "Here are all of them: " <<
623 Teuchos::toString (theirRows ()) <<
". Valid row index " 624 "range (zero-based): [0, " << (numRows - 1) <<
"].");
639 ArrayRCP<size_t> theirNumEntriesPerRow;
640 theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
641 for (size_type k = 0; k < theirNumRows; ++k) {
642 theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
649 send (*pComm, static_cast<int> (theirNumRows),
650 theirNumEntriesPerRow.getRawPtr(), p);
653 const local_ordinal_type theirNumEntries =
654 std::accumulate (theirNumEntriesPerRow.begin(),
655 theirNumEntriesPerRow.end(), 0);
658 cerr <<
"-- Proc 0: Proc " << p <<
" owns " 659 << theirNumEntries <<
" entries" << endl;
664 if (theirNumEntries == 0) {
673 ArrayRCP<GO> theirColInd (theirNumEntries);
674 ArrayRCP<scalar_type> theirValues (theirNumEntries);
681 local_ordinal_type theirCurPos = 0;
682 for (size_type k = 0; k < theirNumRows;
683 theirCurPos += theirNumEntriesPerRow[k], k++) {
684 const local_ordinal_type curNumEntries = theirNumEntriesPerRow[k];
685 const GO theirRow = theirRows[k];
686 const local_ordinal_type curPos = rowPtr[theirRow];
691 if (curNumEntries > 0) {
692 ArrayView<GO> colIndView =
693 colInd (curPos, curNumEntries);
694 ArrayView<GO> theirColIndView =
695 theirColInd (theirCurPos, curNumEntries);
696 std::copy (colIndView.begin(), colIndView.end(),
697 theirColIndView.begin());
699 ArrayView<scalar_type> valuesView =
700 values (curPos, curNumEntries);
701 ArrayView<scalar_type> theirValuesView =
702 theirValues (theirCurPos, curNumEntries);
703 std::copy (valuesView.begin(), valuesView.end(),
704 theirValuesView.begin());
711 send (*pComm, static_cast<int> (theirNumEntries),
712 theirColInd.getRawPtr(), p);
713 send (*pComm, static_cast<int> (theirNumEntries),
714 theirValues.getRawPtr(), p);
717 cerr <<
"-- Proc 0: Finished with proc " << p << endl;
725 numEntriesPerRow = null;
730 if (debug && myRank == 0) {
731 cerr <<
"-- Proc 0: About to fill in myRowPtr" << endl;
739 myRowPtr = arcp<size_t> (myNumRows+1);
741 for (size_type k = 1; k < myNumRows+1; ++k) {
742 myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
744 if (extraDebug && debug) {
745 cerr <<
"Proc " << Teuchos::rank (*(pRowMap->getComm()))
746 <<
": myRowPtr[0.." << myNumRows <<
"] = [";
747 for (size_type k = 0; k < myNumRows+1; ++k) {
753 cerr <<
"]" << endl << endl;
756 if (debug && myRank == 0) {
757 cerr <<
"-- Proc 0: Done with distribute" << endl;
774 static Teuchos::RCP<sparse_matrix_type>
775 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
776 Teuchos::ArrayRCP<size_t>& myRowPtr,
777 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
778 Teuchos::ArrayRCP<scalar_type>& myValues,
779 const Teuchos::RCP<const map_type>& pRowMap,
780 const Teuchos::RCP<const map_type>& pRangeMap,
781 const Teuchos::RCP<const map_type>& pDomainMap,
782 const bool callFillComplete =
true)
784 using Teuchos::ArrayView;
791 typedef global_ordinal_type GO;
798 TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
799 "makeMatrix: myRowPtr array is null. " 800 "Please report this bug to the Tpetra developers.");
801 TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
802 "makeMatrix: domain map is null. " 803 "Please report this bug to the Tpetra developers.");
804 TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
805 "makeMatrix: range map is null. " 806 "Please report this bug to the Tpetra developers.");
807 TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
808 "makeMatrix: row map is null. " 809 "Please report this bug to the Tpetra developers.");
816 RCP<sparse_matrix_type> A =
817 rcp (
new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
822 ArrayView<const GO> myRows = pRowMap->getNodeElementList ();
823 const size_type myNumRows = myRows.size ();
826 const GO indexBase = pRowMap->getIndexBase ();
827 for (size_type i = 0; i < myNumRows; ++i) {
828 const size_type myCurPos = myRowPtr[i];
829 const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
830 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
831 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
834 for (size_type k = 0; k < curNumEntries; ++k) {
835 curColInd[k] += indexBase;
838 if (curNumEntries > 0) {
839 A->insertGlobalValues (myRows[i], curColInd, curValues);
846 myNumEntriesPerRow = null;
851 if (callFillComplete) {
852 A->fillComplete (pDomainMap, pRangeMap);
862 static Teuchos::RCP<sparse_matrix_type>
863 makeMatrix (ArrayRCP<size_t>& myNumEntriesPerRow,
864 ArrayRCP<size_t>& myRowPtr,
865 ArrayRCP<global_ordinal_type>& myColInd,
866 ArrayRCP<scalar_type>& myValues,
867 const Teuchos::RCP<const map_type>& pRowMap,
868 const Teuchos::RCP<const map_type>& pRangeMap,
869 const Teuchos::RCP<const map_type>& pDomainMap,
870 const RCP<Teuchos::ParameterList>& constructorParams,
871 const RCP<Teuchos::ParameterList>& fillCompleteParams)
878 typedef global_ordinal_type GO;
885 TEUCHOS_TEST_FOR_EXCEPTION(
886 myRowPtr.is_null(), std::logic_error,
887 "makeMatrix: myRowPtr array is null. " 888 "Please report this bug to the Tpetra developers.");
889 TEUCHOS_TEST_FOR_EXCEPTION(
890 pDomainMap.is_null(), std::logic_error,
891 "makeMatrix: domain map is null. " 892 "Please report this bug to the Tpetra developers.");
893 TEUCHOS_TEST_FOR_EXCEPTION(
894 pRangeMap.is_null(), std::logic_error,
895 "makeMatrix: range map is null. " 896 "Please report this bug to the Tpetra developers.");
897 TEUCHOS_TEST_FOR_EXCEPTION(
898 pRowMap.is_null(), std::logic_error,
899 "makeMatrix: row map is null. " 900 "Please report this bug to the Tpetra developers.");
907 RCP<sparse_matrix_type> A =
908 rcp (
new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
913 ArrayView<const GO> myRows = pRowMap->getNodeElementList();
914 const size_type myNumRows = myRows.size();
917 const GO indexBase = pRowMap->getIndexBase ();
918 for (size_type i = 0; i < myNumRows; ++i) {
919 const size_type myCurPos = myRowPtr[i];
920 const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
921 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
922 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
925 for (size_type k = 0; k < curNumEntries; ++k) {
926 curColInd[k] += indexBase;
928 if (curNumEntries > 0) {
929 A->insertGlobalValues (myRows[i], curColInd, curValues);
936 myNumEntriesPerRow = null;
941 A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
949 static Teuchos::RCP<sparse_matrix_type>
950 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
951 Teuchos::ArrayRCP<size_t>& myRowPtr,
952 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
953 Teuchos::ArrayRCP<scalar_type>& myValues,
954 const Teuchos::RCP<const map_type>& rowMap,
955 Teuchos::RCP<const map_type>& colMap,
956 const Teuchos::RCP<const map_type>& domainMap,
957 const Teuchos::RCP<const map_type>& rangeMap,
958 const bool callFillComplete =
true)
960 using Teuchos::ArrayView;
965 typedef global_ordinal_type GO;
966 typedef typename ArrayView<const GO>::size_type size_type;
972 RCP<sparse_matrix_type> A;
973 if (colMap.is_null ()) {
974 A = rcp (
new sparse_matrix_type (rowMap, myNumEntriesPerRow,
DynamicProfile));
976 A = rcp (
new sparse_matrix_type (rowMap, colMap, myNumEntriesPerRow,
DynamicProfile));
981 ArrayView<const GO> myRows = rowMap->getNodeElementList ();
982 const size_type myNumRows = myRows.size ();
985 const GO indexBase = rowMap->getIndexBase ();
986 for (size_type i = 0; i < myNumRows; ++i) {
987 const size_type myCurPos = myRowPtr[i];
988 const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
989 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
990 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
993 for (size_type k = 0; k < curNumEntries; ++k) {
994 curColInd[k] += indexBase;
996 if (curNumEntries > 0) {
997 A->insertGlobalValues (myRows[i], curColInd, curValues);
1004 myNumEntriesPerRow = null;
1009 if (callFillComplete) {
1010 A->fillComplete (domainMap, rangeMap);
1011 if (colMap.is_null ()) {
1012 colMap = A->getColMap ();
1036 static RCP<const Teuchos::MatrixMarket::Banner>
1037 readBanner (std::istream& in,
1039 const bool tolerant=
false,
1040 const bool debug=
false)
1042 using Teuchos::MatrixMarket::Banner;
1045 typedef Teuchos::ScalarTraits<scalar_type> STS;
1047 RCP<Banner> pBanner;
1051 const bool readFailed = ! getline(in, line);
1052 TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1053 "Failed to get Matrix Market banner line from input.");
1060 pBanner = rcp (
new Banner (line, tolerant));
1061 }
catch (std::exception& e) {
1062 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1063 "Matrix Market banner line contains syntax error(s): " 1066 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() !=
"matrix",
1067 std::invalid_argument,
"The Matrix Market file does not contain " 1068 "matrix data. Its Banner (first) line says that its object type is \"" 1069 << pBanner->matrixType() <<
"\", rather than the required \"matrix\".");
1073 TEUCHOS_TEST_FOR_EXCEPTION(
1074 ! STS::isComplex && pBanner->dataType() ==
"complex",
1075 std::invalid_argument,
1076 "The Matrix Market file contains complex-valued data, but you are " 1077 "trying to read it into a matrix containing entries of the real-" 1078 "valued Scalar type \"" 1079 << Teuchos::TypeNameTraits<scalar_type>::name() <<
"\".");
1080 TEUCHOS_TEST_FOR_EXCEPTION(
1081 pBanner->dataType() !=
"real" &&
1082 pBanner->dataType() !=
"complex" &&
1083 pBanner->dataType() !=
"integer",
1084 std::invalid_argument,
1085 "When reading Matrix Market data into a Tpetra::CrsMatrix, the " 1086 "Matrix Market file may not contain a \"pattern\" matrix. A " 1087 "pattern matrix is really just a graph with no weights. It " 1088 "should be stored in a CrsGraph, not a CrsMatrix.");
1115 static Tuple<global_ordinal_type, 3>
1116 readCoordDims (std::istream& in,
1118 const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1119 const Teuchos::RCP<const comm_type>& pComm,
1120 const bool tolerant =
false,
1121 const bool debug =
false)
1123 using Teuchos::MatrixMarket::readCoordinateDimensions;
1124 using Teuchos::Tuple;
1129 Tuple<global_ordinal_type, 3> dims;
1135 bool success =
false;
1136 if (pComm->getRank() == 0) {
1137 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() !=
"coordinate",
1138 std::invalid_argument,
"The Tpetra::CrsMatrix Matrix Market reader " 1139 "only accepts \"coordinate\" (sparse) matrix data.");
1141 global_ordinal_type numRows, numCols, numNonzeros;
1143 success = readCoordinateDimensions (in, numRows, numCols,
1144 numNonzeros, lineNumber,
1150 dims[2] = numNonzeros;
1158 int the_success = success ? 1 : 0;
1159 Teuchos::broadcast (*pComm, 0, &the_success);
1160 success = (the_success == 1);
1165 Teuchos::broadcast (*pComm, 0, dims);
1173 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1174 "Error reading Matrix Market sparse matrix: failed to read " 1175 "coordinate matrix dimensions.");
1190 typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1217 static RCP<adder_type>
1218 makeAdder (
const Teuchos::RCP<
const Comm<int> >& pComm,
1219 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1220 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1221 const bool tolerant=
false,
1222 const bool debug=
false)
1224 if (pComm->getRank () == 0) {
1225 typedef Teuchos::MatrixMarket::Raw::Adder<
scalar_type,
1226 global_ordinal_type>
1228 Teuchos::RCP<raw_adder_type> pRaw =
1229 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1231 return Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
1234 return Teuchos::null;
1262 static Teuchos::RCP<sparse_matrix_type>
1264 const RCP<
const Comm<int> >& pComm,
1265 const bool callFillComplete=
true,
1266 const bool tolerant=
false,
1267 const bool debug=
false)
1269 return readSparseFile (filename, pComm, Teuchos::null, callFillComplete, tolerant, debug);
1273 static Teuchos::RCP<sparse_matrix_type>
1275 const RCP<
const Comm<int> >& pComm,
1276 const RCP<node_type>& pNode,
1277 const bool callFillComplete=
true,
1278 const bool tolerant=
false,
1279 const bool debug=
false)
1281 const int myRank = pComm->getRank ();
1286 in.open (filename.c_str ());
1294 return readSparse (in, pComm, pNode, callFillComplete, tolerant, debug);
1328 static Teuchos::RCP<sparse_matrix_type>
1330 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1331 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1332 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1333 const bool tolerant=
false,
1334 const bool debug=
false)
1337 constructorParams, fillCompleteParams,
1342 static Teuchos::RCP<sparse_matrix_type>
1344 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1345 const Teuchos::RCP<node_type>& pNode,
1346 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1347 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1348 const bool tolerant=
false,
1349 const bool debug=
false)
1352 if (pComm->getRank () == 0) {
1353 in.open (filename.c_str ());
1355 return readSparse (in, pComm, pNode, constructorParams,
1356 fillCompleteParams, tolerant, debug);
1396 static Teuchos::RCP<sparse_matrix_type>
1398 const RCP<const map_type>& rowMap,
1399 RCP<const map_type>& colMap,
1400 const RCP<const map_type>& domainMap,
1401 const RCP<const map_type>& rangeMap,
1402 const bool callFillComplete=
true,
1403 const bool tolerant=
false,
1404 const bool debug=
false)
1406 using Teuchos::broadcast;
1407 using Teuchos::outArg;
1408 TEUCHOS_TEST_FOR_EXCEPTION(
1409 rowMap.is_null (), std::invalid_argument,
1410 "Row Map must be nonnull.");
1412 RCP<const Comm<int> > comm = rowMap->getComm ();
1413 const int myRank = comm->getRank ();
1423 in.open (filename.c_str ());
1430 broadcast<int, int> (*comm, 0, outArg (opened));
1431 TEUCHOS_TEST_FOR_EXCEPTION(
1432 opened == 0, std::runtime_error,
1433 "readSparseFile: Failed to open file \"" << filename <<
"\" on " 1435 return readSparse (in, rowMap, colMap, domainMap, rangeMap,
1436 callFillComplete, tolerant, debug);
1464 static Teuchos::RCP<sparse_matrix_type>
1466 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1467 const bool callFillComplete=
true,
1468 const bool tolerant=
false,
1469 const bool debug=
false)
1471 return readSparse (in, pComm, Teuchos::null, callFillComplete, tolerant, debug);
1475 static Teuchos::RCP<sparse_matrix_type>
1477 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1478 const Teuchos::RCP<node_type>& pNode,
1479 const bool callFillComplete=
true,
1480 const bool tolerant=
false,
1481 const bool debug=
false)
1483 using Teuchos::MatrixMarket::Banner;
1484 using Teuchos::arcp;
1485 using Teuchos::ArrayRCP;
1486 using Teuchos::broadcast;
1487 using Teuchos::null;
1490 using Teuchos::REDUCE_MAX;
1491 using Teuchos::reduceAll;
1492 using Teuchos::Tuple;
1495 typedef Teuchos::ScalarTraits<scalar_type> STS;
1497 const bool extraDebug =
false;
1498 const int myRank = pComm->getRank ();
1499 const int rootRank = 0;
1504 size_t lineNumber = 1;
1506 if (debug && myRank == rootRank) {
1507 cerr <<
"Matrix Market reader: readSparse:" << endl
1508 <<
"-- Reading banner line" << endl;
1517 RCP<const Banner> pBanner;
1523 int bannerIsCorrect = 1;
1524 std::ostringstream errMsg;
1526 if (myRank == rootRank) {
1529 pBanner = readBanner (in, lineNumber, tolerant, debug);
1531 catch (std::exception& e) {
1532 errMsg <<
"Attempt to read the Matrix Market file's Banner line " 1533 "threw an exception: " << e.what();
1534 bannerIsCorrect = 0;
1537 if (bannerIsCorrect) {
1542 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
1543 bannerIsCorrect = 0;
1544 errMsg <<
"The Matrix Market input file must contain a " 1545 "\"coordinate\"-format sparse matrix in order to create a " 1546 "Tpetra::CrsMatrix object from it, but the file's matrix " 1547 "type is \"" << pBanner->matrixType() <<
"\" instead.";
1552 if (tolerant && pBanner->matrixType() ==
"array") {
1553 bannerIsCorrect = 0;
1554 errMsg <<
"Matrix Market file must contain a \"coordinate\"-" 1555 "format sparse matrix in order to create a Tpetra::CrsMatrix " 1556 "object from it, but the file's matrix type is \"array\" " 1557 "instead. That probably means the file contains dense matrix " 1564 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1571 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1572 std::invalid_argument, errMsg.str ());
1574 if (debug && myRank == rootRank) {
1575 cerr <<
"-- Reading dimensions line" << endl;
1583 Tuple<global_ordinal_type, 3> dims =
1584 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1586 if (debug && myRank == rootRank) {
1587 cerr <<
"-- Making Adder for collecting matrix data" << endl;
1592 RCP<adder_type> pAdder =
1593 makeAdder (pComm, pBanner, dims, tolerant, debug);
1595 if (debug && myRank == rootRank) {
1596 cerr <<
"-- Reading matrix data" << endl;
1606 int readSuccess = 1;
1607 std::ostringstream errMsg;
1608 if (myRank == rootRank) {
1611 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
1613 reader_type reader (pAdder);
1616 std::pair<bool, std::vector<size_t> > results =
1617 reader.read (in, lineNumber, tolerant, debug);
1618 readSuccess = results.first ? 1 : 0;
1620 catch (std::exception& e) {
1625 broadcast (*pComm, rootRank, ptr (&readSuccess));
1634 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1635 "Failed to read the Matrix Market sparse matrix file: " 1639 if (debug && myRank == rootRank) {
1640 cerr <<
"-- Successfully read the Matrix Market data" << endl;
1653 if (debug && myRank == rootRank) {
1654 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions" 1656 <<
"----- Dimensions before: " 1657 << dims[0] <<
" x " << dims[1]
1661 Tuple<global_ordinal_type, 2> updatedDims;
1662 if (myRank == rootRank) {
1669 std::max (dims[0], pAdder->getAdder()->numRows());
1670 updatedDims[1] = pAdder->getAdder()->numCols();
1672 broadcast (*pComm, rootRank, updatedDims);
1673 dims[0] = updatedDims[0];
1674 dims[1] = updatedDims[1];
1675 if (debug && myRank == rootRank) {
1676 cerr <<
"----- Dimensions after: " << dims[0] <<
" x " 1690 if (myRank == rootRank) {
1697 if (dims[0] < pAdder->getAdder ()->numRows ()) {
1701 broadcast (*pComm, 0, ptr (&dimsMatch));
1702 if (dimsMatch == 0) {
1709 Tuple<global_ordinal_type, 2> addersDims;
1710 if (myRank == rootRank) {
1711 addersDims[0] = pAdder->getAdder()->numRows();
1712 addersDims[1] = pAdder->getAdder()->numCols();
1714 broadcast (*pComm, 0, addersDims);
1715 TEUCHOS_TEST_FOR_EXCEPTION(
1716 dimsMatch == 0, std::runtime_error,
1717 "The matrix metadata says that the matrix is " << dims[0] <<
" x " 1718 << dims[1] <<
", but the actual data says that the matrix is " 1719 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the " 1720 "data includes more rows than reported in the metadata. This " 1721 "is not allowed when parsing in strict mode. Parse the matrix in " 1722 "tolerant mode to ignore the metadata when it disagrees with the " 1727 if (debug && myRank == rootRank) {
1728 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
1740 ArrayRCP<size_t> numEntriesPerRow;
1741 ArrayRCP<size_t> rowPtr;
1742 ArrayRCP<global_ordinal_type> colInd;
1743 ArrayRCP<scalar_type> values;
1748 int mergeAndConvertSucceeded = 1;
1749 std::ostringstream errMsg;
1751 if (myRank == rootRank) {
1753 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
1754 global_ordinal_type> element_type;
1763 const size_type numRows = dims[0];
1766 pAdder->getAdder()->merge ();
1769 const std::vector<element_type>& entries =
1770 pAdder->getAdder()->getEntries();
1773 const size_t numEntries = (size_t)entries.size();
1776 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
1777 <<
" rows and numEntries=" << numEntries
1778 <<
" entries." << endl;
1784 numEntriesPerRow = arcp<size_t> (numRows);
1785 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
1786 rowPtr = arcp<size_t> (numRows+1);
1787 std::fill (rowPtr.begin(), rowPtr.end(), 0);
1788 colInd = arcp<global_ordinal_type> (numEntries);
1789 values = arcp<scalar_type> (numEntries);
1793 global_ordinal_type prvRow = 0;
1796 for (curPos = 0; curPos < numEntries; ++curPos) {
1797 const element_type& curEntry = entries[curPos];
1798 const global_ordinal_type curRow = curEntry.rowIndex();
1799 TEUCHOS_TEST_FOR_EXCEPTION(
1800 curRow < prvRow, std::logic_error,
1801 "Row indices are out of order, even though they are supposed " 1802 "to be sorted. curRow = " << curRow <<
", prvRow = " 1803 << prvRow <<
", at curPos = " << curPos <<
". Please report " 1804 "this bug to the Tpetra developers.");
1805 if (curRow > prvRow) {
1806 for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
1811 numEntriesPerRow[curRow]++;
1812 colInd[curPos] = curEntry.colIndex();
1813 values[curPos] = curEntry.value();
1818 rowPtr[numRows] = numEntries;
1820 catch (std::exception& e) {
1821 mergeAndConvertSucceeded = 0;
1822 errMsg <<
"Failed to merge sparse matrix entries and convert to " 1823 "CSR format: " << e.what();
1826 if (debug && mergeAndConvertSucceeded) {
1828 const size_type numRows = dims[0];
1829 const size_type maxToDisplay = 100;
1831 cerr <<
"----- Proc 0: numEntriesPerRow[0.." 1832 << (numEntriesPerRow.size()-1) <<
"] ";
1833 if (numRows > maxToDisplay) {
1834 cerr <<
"(only showing first and last few entries) ";
1838 if (numRows > maxToDisplay) {
1839 for (size_type k = 0; k < 2; ++k) {
1840 cerr << numEntriesPerRow[k] <<
" ";
1843 for (size_type k = numRows-2; k < numRows-1; ++k) {
1844 cerr << numEntriesPerRow[k] <<
" ";
1848 for (size_type k = 0; k < numRows-1; ++k) {
1849 cerr << numEntriesPerRow[k] <<
" ";
1852 cerr << numEntriesPerRow[numRows-1];
1854 cerr <<
"]" << endl;
1856 cerr <<
"----- Proc 0: rowPtr ";
1857 if (numRows > maxToDisplay) {
1858 cerr <<
"(only showing first and last few entries) ";
1861 if (numRows > maxToDisplay) {
1862 for (size_type k = 0; k < 2; ++k) {
1863 cerr << rowPtr[k] <<
" ";
1866 for (size_type k = numRows-2; k < numRows; ++k) {
1867 cerr << rowPtr[k] <<
" ";
1871 for (size_type k = 0; k < numRows; ++k) {
1872 cerr << rowPtr[k] <<
" ";
1875 cerr << rowPtr[numRows] <<
"]" << endl;
1886 if (debug && myRank == rootRank) {
1887 cerr <<
"-- Making range, domain, and row maps" << endl;
1894 RCP<const map_type> pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
1895 RCP<const map_type> pDomainMap =
1896 makeDomainMap (pRangeMap, dims[0], dims[1]);
1897 RCP<const map_type> pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
1899 if (debug && myRank == rootRank) {
1900 cerr <<
"-- Distributing the matrix data" << endl;
1914 ArrayRCP<size_t> myNumEntriesPerRow;
1915 ArrayRCP<size_t> myRowPtr;
1916 ArrayRCP<global_ordinal_type> myColInd;
1917 ArrayRCP<scalar_type> myValues;
1919 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
1920 numEntriesPerRow, rowPtr, colInd, values, debug);
1922 if (debug && myRank == rootRank) {
1923 cerr <<
"-- Inserting matrix entries on each processor";
1924 if (callFillComplete) {
1925 cerr <<
" and calling fillComplete()";
1936 RCP<sparse_matrix_type> pMatrix =
1937 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
1938 pRowMap, pRangeMap, pDomainMap, callFillComplete);
1943 int localIsNull = pMatrix.is_null () ? 1 : 0;
1944 int globalIsNull = 0;
1945 reduceAll (*pComm, REDUCE_MAX, localIsNull, ptr (&globalIsNull));
1946 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
1947 "Reader::makeMatrix() returned a null pointer on at least one " 1948 "process. Please report this bug to the Tpetra developers.");
1951 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
1952 "Reader::makeMatrix() returned a null pointer. " 1953 "Please report this bug to the Tpetra developers.");
1965 if (callFillComplete) {
1966 const int numProcs = pComm->getSize ();
1968 if (extraDebug && debug) {
1970 pRangeMap->getGlobalNumElements ();
1972 pDomainMap->getGlobalNumElements ();
1973 if (myRank == rootRank) {
1974 cerr <<
"-- Matrix is " 1975 << globalNumRows <<
" x " << globalNumCols
1976 <<
" with " << pMatrix->getGlobalNumEntries()
1977 <<
" entries, and index base " 1978 << pMatrix->getIndexBase() <<
"." << endl;
1981 for (
int p = 0; p < numProcs; ++p) {
1983 cerr <<
"-- Proc " << p <<
" owns " 1984 << pMatrix->getNodeNumCols() <<
" columns, and " 1985 << pMatrix->getNodeNumEntries() <<
" entries." << endl;
1992 if (debug && myRank == rootRank) {
1993 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data" 2028 static Teuchos::RCP<sparse_matrix_type>
2030 const RCP<
const Comm<int> >& pComm,
2031 const RCP<Teuchos::ParameterList>& constructorParams,
2032 const RCP<Teuchos::ParameterList>& fillCompleteParams,
2033 const bool tolerant=
false,
2034 const bool debug=
false)
2036 return readSparse (in, pComm, Teuchos::null, constructorParams,
2037 fillCompleteParams, tolerant, debug);
2041 static Teuchos::RCP<sparse_matrix_type>
2043 const RCP<
const Comm<int> >& pComm,
2044 const RCP<node_type>& pNode,
2045 const RCP<Teuchos::ParameterList>& constructorParams,
2046 const RCP<Teuchos::ParameterList>& fillCompleteParams,
2047 const bool tolerant=
false,
2048 const bool debug=
false)
2050 using Teuchos::MatrixMarket::Banner;
2051 using Teuchos::broadcast;
2053 using Teuchos::reduceAll;
2056 typedef Teuchos::ScalarTraits<scalar_type> STS;
2058 const bool extraDebug =
false;
2059 const int myRank = pComm->getRank ();
2060 const int rootRank = 0;
2065 size_t lineNumber = 1;
2067 if (debug && myRank == rootRank) {
2068 cerr <<
"Matrix Market reader: readSparse:" << endl
2069 <<
"-- Reading banner line" << endl;
2078 RCP<const Banner> pBanner;
2084 int bannerIsCorrect = 1;
2085 std::ostringstream errMsg;
2087 if (myRank == rootRank) {
2090 pBanner = readBanner (in, lineNumber, tolerant, debug);
2092 catch (std::exception& e) {
2093 errMsg <<
"Attempt to read the Matrix Market file's Banner line " 2094 "threw an exception: " << e.what();
2095 bannerIsCorrect = 0;
2098 if (bannerIsCorrect) {
2103 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2104 bannerIsCorrect = 0;
2105 errMsg <<
"The Matrix Market input file must contain a " 2106 "\"coordinate\"-format sparse matrix in order to create a " 2107 "Tpetra::CrsMatrix object from it, but the file's matrix " 2108 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2113 if (tolerant && pBanner->matrixType() ==
"array") {
2114 bannerIsCorrect = 0;
2115 errMsg <<
"Matrix Market file must contain a \"coordinate\"-" 2116 "format sparse matrix in order to create a Tpetra::CrsMatrix " 2117 "object from it, but the file's matrix type is \"array\" " 2118 "instead. That probably means the file contains dense matrix " 2125 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2132 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2133 std::invalid_argument, errMsg.str ());
2135 if (debug && myRank == rootRank) {
2136 cerr <<
"-- Reading dimensions line" << endl;
2144 Tuple<global_ordinal_type, 3> dims =
2145 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2147 if (debug && myRank == rootRank) {
2148 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2153 RCP<adder_type> pAdder =
2154 makeAdder (pComm, pBanner, dims, tolerant, debug);
2156 if (debug && myRank == rootRank) {
2157 cerr <<
"-- Reading matrix data" << endl;
2167 int readSuccess = 1;
2168 std::ostringstream errMsg;
2169 if (myRank == rootRank) {
2172 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2174 reader_type reader (pAdder);
2177 std::pair<bool, std::vector<size_t> > results =
2178 reader.read (in, lineNumber, tolerant, debug);
2179 readSuccess = results.first ? 1 : 0;
2181 catch (std::exception& e) {
2186 broadcast (*pComm, rootRank, ptr (&readSuccess));
2195 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2196 "Failed to read the Matrix Market sparse matrix file: " 2200 if (debug && myRank == rootRank) {
2201 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2214 if (debug && myRank == rootRank) {
2215 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions" 2217 <<
"----- Dimensions before: " 2218 << dims[0] <<
" x " << dims[1]
2222 Tuple<global_ordinal_type, 2> updatedDims;
2223 if (myRank == rootRank) {
2230 std::max (dims[0], pAdder->getAdder()->numRows());
2231 updatedDims[1] = pAdder->getAdder()->numCols();
2233 broadcast (*pComm, rootRank, updatedDims);
2234 dims[0] = updatedDims[0];
2235 dims[1] = updatedDims[1];
2236 if (debug && myRank == rootRank) {
2237 cerr <<
"----- Dimensions after: " << dims[0] <<
" x " 2251 if (myRank == rootRank) {
2258 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2262 broadcast (*pComm, 0, ptr (&dimsMatch));
2263 if (dimsMatch == 0) {
2270 Tuple<global_ordinal_type, 2> addersDims;
2271 if (myRank == rootRank) {
2272 addersDims[0] = pAdder->getAdder()->numRows();
2273 addersDims[1] = pAdder->getAdder()->numCols();
2275 broadcast (*pComm, 0, addersDims);
2276 TEUCHOS_TEST_FOR_EXCEPTION(
2277 dimsMatch == 0, std::runtime_error,
2278 "The matrix metadata says that the matrix is " << dims[0] <<
" x " 2279 << dims[1] <<
", but the actual data says that the matrix is " 2280 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the " 2281 "data includes more rows than reported in the metadata. This " 2282 "is not allowed when parsing in strict mode. Parse the matrix in " 2283 "tolerant mode to ignore the metadata when it disagrees with the " 2288 if (debug && myRank == rootRank) {
2289 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
2301 ArrayRCP<size_t> numEntriesPerRow;
2302 ArrayRCP<size_t> rowPtr;
2303 ArrayRCP<global_ordinal_type> colInd;
2304 ArrayRCP<scalar_type> values;
2309 int mergeAndConvertSucceeded = 1;
2310 std::ostringstream errMsg;
2312 if (myRank == rootRank) {
2314 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2315 global_ordinal_type> element_type;
2324 const size_type numRows = dims[0];
2327 pAdder->getAdder()->merge ();
2330 const std::vector<element_type>& entries =
2331 pAdder->getAdder()->getEntries();
2334 const size_t numEntries = (size_t)entries.size();
2337 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
2338 <<
" rows and numEntries=" << numEntries
2339 <<
" entries." << endl;
2345 numEntriesPerRow = arcp<size_t> (numRows);
2346 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2347 rowPtr = arcp<size_t> (numRows+1);
2348 std::fill (rowPtr.begin(), rowPtr.end(), 0);
2349 colInd = arcp<global_ordinal_type> (numEntries);
2350 values = arcp<scalar_type> (numEntries);
2354 global_ordinal_type prvRow = 0;
2357 for (curPos = 0; curPos < numEntries; ++curPos) {
2358 const element_type& curEntry = entries[curPos];
2359 const global_ordinal_type curRow = curEntry.rowIndex();
2360 TEUCHOS_TEST_FOR_EXCEPTION(
2361 curRow < prvRow, std::logic_error,
2362 "Row indices are out of order, even though they are supposed " 2363 "to be sorted. curRow = " << curRow <<
", prvRow = " 2364 << prvRow <<
", at curPos = " << curPos <<
". Please report " 2365 "this bug to the Tpetra developers.");
2366 if (curRow > prvRow) {
2367 for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
2372 numEntriesPerRow[curRow]++;
2373 colInd[curPos] = curEntry.colIndex();
2374 values[curPos] = curEntry.value();
2379 rowPtr[numRows] = numEntries;
2381 catch (std::exception& e) {
2382 mergeAndConvertSucceeded = 0;
2383 errMsg <<
"Failed to merge sparse matrix entries and convert to " 2384 "CSR format: " << e.what();
2387 if (debug && mergeAndConvertSucceeded) {
2389 const size_type numRows = dims[0];
2390 const size_type maxToDisplay = 100;
2392 cerr <<
"----- Proc 0: numEntriesPerRow[0.." 2393 << (numEntriesPerRow.size()-1) <<
"] ";
2394 if (numRows > maxToDisplay) {
2395 cerr <<
"(only showing first and last few entries) ";
2399 if (numRows > maxToDisplay) {
2400 for (size_type k = 0; k < 2; ++k) {
2401 cerr << numEntriesPerRow[k] <<
" ";
2404 for (size_type k = numRows-2; k < numRows-1; ++k) {
2405 cerr << numEntriesPerRow[k] <<
" ";
2409 for (size_type k = 0; k < numRows-1; ++k) {
2410 cerr << numEntriesPerRow[k] <<
" ";
2413 cerr << numEntriesPerRow[numRows-1];
2415 cerr <<
"]" << endl;
2417 cerr <<
"----- Proc 0: rowPtr ";
2418 if (numRows > maxToDisplay) {
2419 cerr <<
"(only showing first and last few entries) ";
2422 if (numRows > maxToDisplay) {
2423 for (size_type k = 0; k < 2; ++k) {
2424 cerr << rowPtr[k] <<
" ";
2427 for (size_type k = numRows-2; k < numRows; ++k) {
2428 cerr << rowPtr[k] <<
" ";
2432 for (size_type k = 0; k < numRows; ++k) {
2433 cerr << rowPtr[k] <<
" ";
2436 cerr << rowPtr[numRows] <<
"]" << endl;
2447 if (debug && myRank == rootRank) {
2448 cerr <<
"-- Making range, domain, and row maps" << endl;
2455 RCP<const map_type> pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
2456 RCP<const map_type> pDomainMap =
2457 makeDomainMap (pRangeMap, dims[0], dims[1]);
2458 RCP<const map_type> pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
2460 if (debug && myRank == rootRank) {
2461 cerr <<
"-- Distributing the matrix data" << endl;
2475 ArrayRCP<size_t> myNumEntriesPerRow;
2476 ArrayRCP<size_t> myRowPtr;
2477 ArrayRCP<global_ordinal_type> myColInd;
2478 ArrayRCP<scalar_type> myValues;
2480 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
2481 numEntriesPerRow, rowPtr, colInd, values, debug);
2483 if (debug && myRank == rootRank) {
2484 cerr <<
"-- Inserting matrix entries on each process " 2485 "and calling fillComplete()" << endl;
2494 Teuchos::RCP<sparse_matrix_type> pMatrix =
2495 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
2496 pRowMap, pRangeMap, pDomainMap, constructorParams,
2497 fillCompleteParams);
2502 int localIsNull = pMatrix.is_null () ? 1 : 0;
2503 int globalIsNull = 0;
2504 reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
2505 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
2506 "Reader::makeMatrix() returned a null pointer on at least one " 2507 "process. Please report this bug to the Tpetra developers.");
2510 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
2511 "Reader::makeMatrix() returned a null pointer. " 2512 "Please report this bug to the Tpetra developers.");
2521 if (extraDebug && debug) {
2522 const int numProcs = pComm->getSize ();
2524 pRangeMap->getGlobalNumElements();
2526 pDomainMap->getGlobalNumElements();
2527 if (myRank == rootRank) {
2528 cerr <<
"-- Matrix is " 2529 << globalNumRows <<
" x " << globalNumCols
2530 <<
" with " << pMatrix->getGlobalNumEntries()
2531 <<
" entries, and index base " 2532 << pMatrix->getIndexBase() <<
"." << endl;
2535 for (
int p = 0; p < numProcs; ++p) {
2537 cerr <<
"-- Proc " << p <<
" owns " 2538 << pMatrix->getNodeNumCols() <<
" columns, and " 2539 << pMatrix->getNodeNumEntries() <<
" entries." << endl;
2545 if (debug && myRank == rootRank) {
2546 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data" 2592 static Teuchos::RCP<sparse_matrix_type>
2594 const Teuchos::RCP<const map_type>& rowMap,
2595 Teuchos::RCP<const map_type>& colMap,
2596 const Teuchos::RCP<const map_type>& domainMap,
2597 const Teuchos::RCP<const map_type>& rangeMap,
2598 const bool callFillComplete=
true,
2599 const bool tolerant=
false,
2600 const bool debug=
false)
2602 using Teuchos::MatrixMarket::Banner;
2604 using Teuchos::broadcast;
2605 using Teuchos::Comm;
2608 using Teuchos::reduceAll;
2609 using Teuchos::Tuple;
2612 typedef Teuchos::ScalarTraits<scalar_type> STS;
2614 RCP<const Comm<int> > pComm = rowMap->getComm ();
2615 const int myRank = pComm->getRank ();
2616 const int rootRank = 0;
2617 const bool extraDebug =
false;
2622 TEUCHOS_TEST_FOR_EXCEPTION(
2623 rowMap.is_null (), std::invalid_argument,
2624 "Row Map must be nonnull.");
2625 TEUCHOS_TEST_FOR_EXCEPTION(
2626 rangeMap.is_null (), std::invalid_argument,
2627 "Range Map must be nonnull.");
2628 TEUCHOS_TEST_FOR_EXCEPTION(
2629 domainMap.is_null (), std::invalid_argument,
2630 "Domain Map must be nonnull.");
2631 TEUCHOS_TEST_FOR_EXCEPTION(
2632 rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
2633 std::invalid_argument,
2634 "The specified row Map's communicator (rowMap->getComm())" 2635 "differs from the given separately supplied communicator pComm.");
2636 TEUCHOS_TEST_FOR_EXCEPTION(
2637 domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
2638 std::invalid_argument,
2639 "The specified domain Map's communicator (domainMap->getComm())" 2640 "differs from the given separately supplied communicator pComm.");
2641 TEUCHOS_TEST_FOR_EXCEPTION(
2642 rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
2643 std::invalid_argument,
2644 "The specified range Map's communicator (rangeMap->getComm())" 2645 "differs from the given separately supplied communicator pComm.");
2650 size_t lineNumber = 1;
2652 if (debug && myRank == rootRank) {
2653 cerr <<
"Matrix Market reader: readSparse:" << endl
2654 <<
"-- Reading banner line" << endl;
2663 RCP<const Banner> pBanner;
2669 int bannerIsCorrect = 1;
2670 std::ostringstream errMsg;
2672 if (myRank == rootRank) {
2675 pBanner = readBanner (in, lineNumber, tolerant, debug);
2677 catch (std::exception& e) {
2678 errMsg <<
"Attempt to read the Matrix Market file's Banner line " 2679 "threw an exception: " << e.what();
2680 bannerIsCorrect = 0;
2683 if (bannerIsCorrect) {
2688 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2689 bannerIsCorrect = 0;
2690 errMsg <<
"The Matrix Market input file must contain a " 2691 "\"coordinate\"-format sparse matrix in order to create a " 2692 "Tpetra::CrsMatrix object from it, but the file's matrix " 2693 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2698 if (tolerant && pBanner->matrixType() ==
"array") {
2699 bannerIsCorrect = 0;
2700 errMsg <<
"Matrix Market file must contain a \"coordinate\"-" 2701 "format sparse matrix in order to create a Tpetra::CrsMatrix " 2702 "object from it, but the file's matrix type is \"array\" " 2703 "instead. That probably means the file contains dense matrix " 2710 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2717 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2718 std::invalid_argument, errMsg.str ());
2720 if (debug && myRank == rootRank) {
2721 cerr <<
"-- Reading dimensions line" << endl;
2729 Tuple<global_ordinal_type, 3> dims =
2730 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2732 if (debug && myRank == rootRank) {
2733 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2740 RCP<adder_type> pAdder =
2741 makeAdder (pComm, pBanner, dims, tolerant, debug);
2743 if (debug && myRank == rootRank) {
2744 cerr <<
"-- Reading matrix data" << endl;
2754 int readSuccess = 1;
2755 std::ostringstream errMsg;
2756 if (myRank == rootRank) {
2759 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2761 reader_type reader (pAdder);
2764 std::pair<bool, std::vector<size_t> > results =
2765 reader.read (in, lineNumber, tolerant, debug);
2766 readSuccess = results.first ? 1 : 0;
2768 catch (std::exception& e) {
2773 broadcast (*pComm, rootRank, ptr (&readSuccess));
2782 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2783 "Failed to read the Matrix Market sparse matrix file: " 2787 if (debug && myRank == rootRank) {
2788 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2801 if (debug && myRank == rootRank) {
2802 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions" 2804 <<
"----- Dimensions before: " 2805 << dims[0] <<
" x " << dims[1]
2809 Tuple<global_ordinal_type, 2> updatedDims;
2810 if (myRank == rootRank) {
2817 std::max (dims[0], pAdder->getAdder()->numRows());
2818 updatedDims[1] = pAdder->getAdder()->numCols();
2820 broadcast (*pComm, rootRank, updatedDims);
2821 dims[0] = updatedDims[0];
2822 dims[1] = updatedDims[1];
2823 if (debug && myRank == rootRank) {
2824 cerr <<
"----- Dimensions after: " << dims[0] <<
" x " 2838 if (myRank == rootRank) {
2845 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2849 broadcast (*pComm, 0, ptr (&dimsMatch));
2850 if (dimsMatch == 0) {
2857 Tuple<global_ordinal_type, 2> addersDims;
2858 if (myRank == rootRank) {
2859 addersDims[0] = pAdder->getAdder()->numRows();
2860 addersDims[1] = pAdder->getAdder()->numCols();
2862 broadcast (*pComm, 0, addersDims);
2863 TEUCHOS_TEST_FOR_EXCEPTION(
2864 dimsMatch == 0, std::runtime_error,
2865 "The matrix metadata says that the matrix is " << dims[0] <<
" x " 2866 << dims[1] <<
", but the actual data says that the matrix is " 2867 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the " 2868 "data includes more rows than reported in the metadata. This " 2869 "is not allowed when parsing in strict mode. Parse the matrix in " 2870 "tolerant mode to ignore the metadata when it disagrees with the " 2875 if (debug && myRank == rootRank) {
2876 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
2888 ArrayRCP<size_t> numEntriesPerRow;
2889 ArrayRCP<size_t> rowPtr;
2890 ArrayRCP<global_ordinal_type> colInd;
2891 ArrayRCP<scalar_type> values;
2896 int mergeAndConvertSucceeded = 1;
2897 std::ostringstream errMsg;
2899 if (myRank == rootRank) {
2901 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2902 global_ordinal_type> element_type;
2911 const size_type numRows = dims[0];
2914 pAdder->getAdder()->merge ();
2917 const std::vector<element_type>& entries =
2918 pAdder->getAdder()->getEntries();
2921 const size_t numEntries = (size_t)entries.size();
2924 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
2925 <<
" rows and numEntries=" << numEntries
2926 <<
" entries." << endl;
2932 numEntriesPerRow = arcp<size_t> (numRows);
2933 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2934 rowPtr = arcp<size_t> (numRows+1);
2935 std::fill (rowPtr.begin(), rowPtr.end(), 0);
2936 colInd = arcp<global_ordinal_type> (numEntries);
2937 values = arcp<scalar_type> (numEntries);
2941 global_ordinal_type prvRow = 0;
2944 for (curPos = 0; curPos < numEntries; ++curPos) {
2945 const element_type& curEntry = entries[curPos];
2946 const global_ordinal_type curRow = curEntry.rowIndex();
2947 TEUCHOS_TEST_FOR_EXCEPTION(
2948 curRow < prvRow, std::logic_error,
2949 "Row indices are out of order, even though they are supposed " 2950 "to be sorted. curRow = " << curRow <<
", prvRow = " 2951 << prvRow <<
", at curPos = " << curPos <<
". Please report " 2952 "this bug to the Tpetra developers.");
2953 if (curRow > prvRow) {
2954 for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
2959 numEntriesPerRow[curRow]++;
2960 colInd[curPos] = curEntry.colIndex();
2961 values[curPos] = curEntry.value();
2966 rowPtr[numRows] = numEntries;
2968 catch (std::exception& e) {
2969 mergeAndConvertSucceeded = 0;
2970 errMsg <<
"Failed to merge sparse matrix entries and convert to " 2971 "CSR format: " << e.what();
2974 if (debug && mergeAndConvertSucceeded) {
2976 const size_type numRows = dims[0];
2977 const size_type maxToDisplay = 100;
2979 cerr <<
"----- Proc 0: numEntriesPerRow[0.." 2980 << (numEntriesPerRow.size()-1) <<
"] ";
2981 if (numRows > maxToDisplay) {
2982 cerr <<
"(only showing first and last few entries) ";
2986 if (numRows > maxToDisplay) {
2987 for (size_type k = 0; k < 2; ++k) {
2988 cerr << numEntriesPerRow[k] <<
" ";
2991 for (size_type k = numRows-2; k < numRows-1; ++k) {
2992 cerr << numEntriesPerRow[k] <<
" ";
2996 for (size_type k = 0; k < numRows-1; ++k) {
2997 cerr << numEntriesPerRow[k] <<
" ";
3000 cerr << numEntriesPerRow[numRows-1];
3002 cerr <<
"]" << endl;
3004 cerr <<
"----- Proc 0: rowPtr ";
3005 if (numRows > maxToDisplay) {
3006 cerr <<
"(only showing first and last few entries) ";
3009 if (numRows > maxToDisplay) {
3010 for (size_type k = 0; k < 2; ++k) {
3011 cerr << rowPtr[k] <<
" ";
3014 for (size_type k = numRows-2; k < numRows; ++k) {
3015 cerr << rowPtr[k] <<
" ";
3019 for (size_type k = 0; k < numRows; ++k) {
3020 cerr << rowPtr[k] <<
" ";
3023 cerr << rowPtr[numRows] <<
"]" << endl;
3025 cerr <<
"----- Proc 0: colInd = [";
3026 for (
size_t k = 0; k < rowPtr[numRows]; ++k) {
3027 cerr << colInd[k] <<
" ";
3029 cerr <<
"]" << endl;
3043 if (debug && myRank == rootRank) {
3044 cerr <<
"-- Verifying Maps" << endl;
3046 TEUCHOS_TEST_FOR_EXCEPTION(
3047 as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
3048 std::invalid_argument,
3049 "The range Map has " << rangeMap->getGlobalNumElements ()
3050 <<
" entries, but the matrix has a global number of rows " << dims[0]
3052 TEUCHOS_TEST_FOR_EXCEPTION(
3053 as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
3054 std::invalid_argument,
3055 "The domain Map has " << domainMap->getGlobalNumElements ()
3056 <<
" entries, but the matrix has a global number of columns " 3060 RCP<Teuchos::FancyOStream> err = debug ?
3061 Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
3063 RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
3064 ArrayView<const global_ordinal_type> myRows =
3065 gatherRowMap->getNodeElementList ();
3066 const size_type myNumRows = myRows.size ();
3067 const global_ordinal_type indexBase = gatherRowMap->getIndexBase ();
3069 ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
3070 for (size_type i_ = 0; i_ < myNumRows; i_++) {
3071 gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_]-indexBase];
3077 RCP<sparse_matrix_type> A_proc0 =
3078 rcp (
new sparse_matrix_type (gatherRowMap, gatherNumEntriesPerRow,
3080 if (myRank == rootRank) {
3082 cerr <<
"-- Proc 0: Filling gather matrix" << endl;
3085 cerr <<
"---- Rows: " << Teuchos::toString (myRows) << endl;
3089 for (size_type i_ = 0; i_ < myNumRows; ++i_) {
3090 size_type i = myRows[i_] - indexBase;
3092 const size_type curPos = as<size_type> (rowPtr[i]);
3093 const local_ordinal_type curNumEntries = numEntriesPerRow[i];
3094 ArrayView<global_ordinal_type> curColInd =
3095 colInd.view (curPos, curNumEntries);
3096 ArrayView<scalar_type> curValues =
3097 values.view (curPos, curNumEntries);
3100 for (size_type k = 0; k < curNumEntries; ++k) {
3101 curColInd[k] += indexBase;
3104 cerr <<
"------ Columns: " << Teuchos::toString (curColInd) << endl;
3105 cerr <<
"------ Values: " << Teuchos::toString (curValues) << endl;
3108 if (curNumEntries > 0) {
3109 A_proc0->insertGlobalValues (myRows[i_], curColInd, curValues);
3115 numEntriesPerRow = null;
3121 RCP<sparse_matrix_type> A;
3122 if (colMap.is_null ()) {
3123 A = rcp (
new sparse_matrix_type (rowMap, 0));
3125 A = rcp (
new sparse_matrix_type (rowMap, colMap, 0));
3128 export_type exp (gatherRowMap, rowMap);
3129 A->doExport (*A_proc0, exp,
INSERT);
3131 if (callFillComplete) {
3132 A->fillComplete (domainMap, rangeMap);
3144 if (callFillComplete) {
3145 const int numProcs = pComm->getSize ();
3147 if (extraDebug && debug) {
3148 const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
3149 const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
3150 if (myRank == rootRank) {
3151 cerr <<
"-- Matrix is " 3152 << globalNumRows <<
" x " << globalNumCols
3153 <<
" with " << A->getGlobalNumEntries()
3154 <<
" entries, and index base " 3155 << A->getIndexBase() <<
"." << endl;
3158 for (
int p = 0; p < numProcs; ++p) {
3160 cerr <<
"-- Proc " << p <<
" owns " 3161 << A->getNodeNumCols() <<
" columns, and " 3162 << A->getNodeNumEntries() <<
" entries." << endl;
3169 if (debug && myRank == rootRank) {
3170 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data" 3205 static RCP<multivector_type>
3207 const RCP<const comm_type>& comm,
3208 RCP<const map_type>& map,
3209 const bool tolerant=
false,
3210 const bool debug=
false)
3213 if (comm->getRank () == 0) {
3214 in.open (filename.c_str ());
3216 return readDense (in, comm, map, tolerant, debug);
3220 static RCP<multivector_type>
3222 const RCP<const comm_type>& comm,
3223 const RCP<node_type>& node,
3224 RCP<const map_type>& map,
3225 const bool tolerant=
false,
3226 const bool debug=
false)
3229 if (comm->getRank () == 0) {
3230 in.open (filename.c_str ());
3232 return readDense (in, comm, node, map, tolerant, debug);
3264 static RCP<vector_type>
3266 const RCP<const comm_type>& comm,
3267 RCP<const map_type>& map,
3268 const bool tolerant=
false,
3269 const bool debug=
false)
3272 if (comm->getRank () == 0) {
3273 in.open (filename.c_str ());
3275 return readVector (in, comm, map, tolerant, debug);
3280 static RCP<vector_type>
3282 const RCP<const comm_type>& comm,
3283 const RCP<node_type>& node,
3284 RCP<const map_type>& map,
3285 const bool tolerant=
false,
3286 const bool debug=
false)
3289 if (comm->getRank () == 0) {
3290 in.open (filename.c_str ());
3292 return readVector (in, comm, node, map, tolerant, debug);
3362 static RCP<multivector_type>
3364 const RCP<const comm_type>& comm,
3365 RCP<const map_type>& map,
3366 const bool tolerant=
false,
3367 const bool debug=
false)
3369 return readDense (in, comm, Teuchos::null, map, tolerant, debug);
3373 static RCP<multivector_type>
3375 const RCP<const comm_type>& comm,
3376 const RCP<node_type>& node,
3377 RCP<const map_type>& map,
3378 const bool tolerant=
false,
3379 const bool debug=
false)
3381 Teuchos::RCP<Teuchos::FancyOStream> err =
3382 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
3383 return readDenseImpl<scalar_type> (in, comm, node, map,
3384 err, tolerant, debug);
3388 static RCP<vector_type>
3390 const RCP<const comm_type>& comm,
3391 RCP<const map_type>& map,
3392 const bool tolerant=
false,
3393 const bool debug=
false)
3395 Teuchos::RCP<Teuchos::FancyOStream> err =
3396 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
3397 return readVectorImpl<scalar_type> (in, comm, Teuchos::null, map,
3398 err, tolerant, debug);
3402 static RCP<vector_type>
3404 const RCP<const comm_type>& comm,
3405 const RCP<node_type>& node,
3406 RCP<const map_type>& map,
3407 const bool tolerant=
false,
3408 const bool debug=
false)
3410 Teuchos::RCP<Teuchos::FancyOStream> err =
3411 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
3412 return readVectorImpl<scalar_type> (in, comm, node, map,
3413 err, tolerant, debug);
3437 static RCP<const map_type>
3439 const RCP<const comm_type>& comm,
3440 const RCP<node_type>& node,
3441 const bool tolerant=
false,
3442 const bool debug=
false)
3444 using Teuchos::inOutArg;
3445 using Teuchos::broadcast;
3449 if (comm->getRank () == 0) {
3450 in.open (filename.c_str ());
3455 broadcast<int, int> (*comm, 0, inOutArg (success));
3456 TEUCHOS_TEST_FOR_EXCEPTION(
3457 success == 0, std::runtime_error,
3458 "Tpetra::MatrixMarket::Reader::readMapFile: " 3459 "Failed to read file \"" << filename <<
"\" on Process 0.");
3460 return readMap (in, comm, node, tolerant, debug);
3464 template<
class MultiVectorScalarType>
3469 readDenseImpl (std::istream& in,
3470 const RCP<const comm_type>& comm,
3471 const RCP<node_type>& node,
3472 RCP<const map_type>& map,
3473 const Teuchos::RCP<Teuchos::FancyOStream>& err,
3474 const bool tolerant=
false,
3475 const bool debug=
false)
3477 using Teuchos::MatrixMarket::Banner;
3478 using Teuchos::MatrixMarket::checkCommentLine;
3480 using Teuchos::broadcast;
3481 using Teuchos::outArg;
3483 typedef MultiVectorScalarType ST;
3484 typedef local_ordinal_type LO;
3485 typedef global_ordinal_type GO;
3486 typedef node_type NT;
3487 typedef Teuchos::ScalarTraits<ST> STS;
3488 typedef typename STS::magnitudeType MT;
3489 typedef Teuchos::ScalarTraits<MT> STM;
3494 const int myRank = comm->getRank ();
3496 if (! err.is_null ()) {
3500 *err << myRank <<
": readDenseImpl" << endl;
3502 if (! err.is_null ()) {
3536 size_t lineNumber = 1;
3539 std::ostringstream exMsg;
3540 int localBannerReadSuccess = 1;
3541 int localDimsReadSuccess = 1;
3546 *err << myRank <<
": readDenseImpl: Reading banner line (dense)" << endl;
3552 RCP<const Banner> pBanner;
3554 pBanner = readBanner (in, lineNumber, tolerant, debug);
3555 }
catch (std::exception& e) {
3557 localBannerReadSuccess = 0;
3560 if (localBannerReadSuccess) {
3561 if (pBanner->matrixType () !=
"array") {
3562 exMsg <<
"The Matrix Market file does not contain dense matrix " 3563 "data. Its banner (first) line says that its matrix type is \"" 3564 << pBanner->matrixType () <<
"\", rather that the required " 3566 localBannerReadSuccess = 0;
3567 }
else if (pBanner->dataType() ==
"pattern") {
3568 exMsg <<
"The Matrix Market file's banner (first) " 3569 "line claims that the matrix's data type is \"pattern\". This does " 3570 "not make sense for a dense matrix, yet the file reports the matrix " 3571 "as dense. The only valid data types for a dense matrix are " 3572 "\"real\", \"complex\", and \"integer\".";
3573 localBannerReadSuccess = 0;
3577 dims[2] = encodeDataType (pBanner->dataType ());
3583 if (localBannerReadSuccess) {
3585 *err << myRank <<
": readDenseImpl: Reading dimensions line (dense)" << endl;
3594 bool commentLine =
true;
3596 while (commentLine) {
3599 if (in.eof () || in.fail ()) {
3600 exMsg <<
"Unable to get array dimensions line (at all) from line " 3601 << lineNumber <<
" of input stream. The input stream " 3602 <<
"claims that it is " 3603 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
3604 localDimsReadSuccess = 0;
3607 if (getline (in, line)) {
3614 size_t start = 0, size = 0;
3615 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
3622 std::istringstream istr (line);
3626 if (istr.eof () || istr.fail ()) {
3627 exMsg <<
"Unable to read any data from line " << lineNumber
3628 <<
" of input; the line should contain the matrix dimensions " 3629 <<
"\"<numRows> <numCols>\".";
3630 localDimsReadSuccess = 0;
3635 exMsg <<
"Failed to get number of rows from line " 3636 << lineNumber <<
" of input; the line should contains the " 3637 <<
"matrix dimensions \"<numRows> <numCols>\".";
3638 localDimsReadSuccess = 0;
3640 dims[0] = theNumRows;
3642 exMsg <<
"No more data after number of rows on line " 3643 << lineNumber <<
" of input; the line should contain the " 3644 <<
"matrix dimensions \"<numRows> <numCols>\".";
3645 localDimsReadSuccess = 0;
3650 exMsg <<
"Failed to get number of columns from line " 3651 << lineNumber <<
" of input; the line should contain " 3652 <<
"the matrix dimensions \"<numRows> <numCols>\".";
3653 localDimsReadSuccess = 0;
3655 dims[1] = theNumCols;
3666 Tuple<GO, 5> bannerDimsReadResult;
3668 bannerDimsReadResult[0] = dims[0];
3669 bannerDimsReadResult[1] = dims[1];
3670 bannerDimsReadResult[2] = dims[2];
3671 bannerDimsReadResult[3] = localBannerReadSuccess;
3672 bannerDimsReadResult[4] = localDimsReadSuccess;
3676 broadcast (*comm, 0, bannerDimsReadResult);
3678 TEUCHOS_TEST_FOR_EXCEPTION(
3679 bannerDimsReadResult[3] == 0, std::runtime_error,
3680 "Failed to read banner line: " << exMsg.str ());
3681 TEUCHOS_TEST_FOR_EXCEPTION(
3682 bannerDimsReadResult[4] == 0, std::runtime_error,
3683 "Failed to read matrix dimensions line: " << exMsg.str ());
3685 dims[0] = bannerDimsReadResult[0];
3686 dims[1] = bannerDimsReadResult[1];
3687 dims[2] = bannerDimsReadResult[2];
3692 const size_t numCols =
static_cast<size_t> (dims[1]);
3697 RCP<const map_type> proc0Map;
3698 if (map.is_null ()) {
3702 if (node.is_null ()) {
3703 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
3705 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
3707 const size_t localNumRows = (myRank == 0) ? numRows : 0;
3709 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
3710 comm, map->getNode ());
3713 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
3717 RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
3723 int localReadDataSuccess = 1;
3727 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)" 3732 TEUCHOS_TEST_FOR_EXCEPTION(
3733 ! X->isConstantStride (), std::logic_error,
3734 "Can't get a 1-D view of the entries of the MultiVector X on " 3735 "Process 0, because the stride between the columns of X is not " 3736 "constant. This shouldn't happen because we just created X and " 3737 "haven't filled it in yet. Please report this bug to the Tpetra " 3744 ArrayRCP<ST> X_view = X->get1dViewNonConst ();
3745 TEUCHOS_TEST_FOR_EXCEPTION(
3746 as<global_size_t> (X_view.size ()) < numRows * numCols,
3748 "The view of X has size " << X_view <<
" which is not enough to " 3749 "accommodate the expected number of entries numRows*numCols = " 3750 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". " 3751 "Please report this bug to the Tpetra developers.");
3752 const size_t stride = X->getStride ();
3759 const bool isComplex = (dims[2] == 1);
3760 size_type count = 0, curRow = 0, curCol = 0;
3763 while (getline (in, line)) {
3767 size_t start = 0, size = 0;
3768 const bool commentLine =
3769 checkCommentLine (line, start, size, lineNumber, tolerant);
3770 if (! commentLine) {
3776 if (count >= X_view.size()) {
3781 TEUCHOS_TEST_FOR_EXCEPTION(
3782 count >= X_view.size(),
3784 "The Matrix Market input stream has more data in it than " 3785 "its metadata reported. Current line number is " 3786 << lineNumber <<
".");
3795 const size_t pos = line.substr (start, size).find (
':');
3796 if (pos != std::string::npos) {
3800 std::istringstream istr (line.substr (start, size));
3804 if (istr.eof() || istr.fail()) {
3811 TEUCHOS_TEST_FOR_EXCEPTION(
3812 ! tolerant && (istr.eof() || istr.fail()),
3814 "Line " << lineNumber <<
" of the Matrix Market file is " 3815 "empty, or we cannot read from it for some other reason.");
3818 ST val = STS::zero();
3821 MT real = STM::zero(), imag = STM::zero();
3838 TEUCHOS_TEST_FOR_EXCEPTION(
3839 ! tolerant && istr.eof(), std::runtime_error,
3840 "Failed to get the real part of a complex-valued matrix " 3841 "entry from line " << lineNumber <<
" of the Matrix Market " 3847 }
else if (istr.eof()) {
3848 TEUCHOS_TEST_FOR_EXCEPTION(
3849 ! tolerant && istr.eof(), std::runtime_error,
3850 "Missing imaginary part of a complex-valued matrix entry " 3851 "on line " << lineNumber <<
" of the Matrix Market file.");
3857 TEUCHOS_TEST_FOR_EXCEPTION(
3858 ! tolerant && istr.fail(), std::runtime_error,
3859 "Failed to get the imaginary part of a complex-valued " 3860 "matrix entry from line " << lineNumber <<
" of the " 3861 "Matrix Market file.");
3868 TEUCHOS_TEST_FOR_EXCEPTION(
3869 ! tolerant && istr.fail(), std::runtime_error,
3870 "Failed to get a real-valued matrix entry from line " 3871 << lineNumber <<
" of the Matrix Market file.");
3874 if (istr.fail() && tolerant) {
3880 TEUCHOS_TEST_FOR_EXCEPTION(
3881 ! tolerant && istr.fail(), std::runtime_error,
3882 "Failed to read matrix data from line " << lineNumber
3883 <<
" of the Matrix Market file.");
3886 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
3888 curRow = count % numRows;
3889 curCol = count / numRows;
3890 X_view[curRow + curCol*stride] = val;
3895 TEUCHOS_TEST_FOR_EXCEPTION(
3896 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
3898 "The Matrix Market metadata reports that the dense matrix is " 3899 << numRows <<
" x " << numCols <<
", and thus has " 3900 << numRows*numCols <<
" total entries, but we only found " << count
3901 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
3902 }
catch (std::exception& e) {
3904 localReadDataSuccess = 0;
3909 *err << myRank <<
": readDenseImpl: done reading data" << endl;
3913 int globalReadDataSuccess = localReadDataSuccess;
3914 broadcast (*comm, 0, outArg (globalReadDataSuccess));
3915 TEUCHOS_TEST_FOR_EXCEPTION(
3916 globalReadDataSuccess == 0, std::runtime_error,
3917 "Failed to read the multivector's data: " << exMsg.str ());
3922 if (comm->getSize () == 1 && map.is_null ()) {
3924 if (! err.is_null ()) {
3928 *err << myRank <<
": readDenseImpl: done" << endl;
3930 if (! err.is_null ()) {
3937 *err << myRank <<
": readDenseImpl: Creating target MV" << endl;
3941 RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
3944 *err << myRank <<
": readDenseImpl: Creating Export" << endl;
3953 *err << myRank <<
": readDenseImpl: Exporting" << endl;
3956 Y->doExport (*X, exporter,
INSERT);
3958 if (! err.is_null ()) {
3962 *err << myRank <<
": readDenseImpl: done" << endl;
3964 if (! err.is_null ()) {
3973 template<
class VectorScalarType>
3978 readVectorImpl (std::istream& in,
3979 const RCP<const comm_type>& comm,
3980 const RCP<node_type>& node,
3981 RCP<const map_type>& map,
3982 const Teuchos::RCP<Teuchos::FancyOStream>& err,
3983 const bool tolerant=
false,
3984 const bool debug=
false)
3986 using Teuchos::MatrixMarket::Banner;
3987 using Teuchos::MatrixMarket::checkCommentLine;
3989 using Teuchos::broadcast;
3990 using Teuchos::outArg;
3992 typedef VectorScalarType ST;
3993 typedef local_ordinal_type LO;
3994 typedef global_ordinal_type GO;
3995 typedef node_type NT;
3996 typedef Teuchos::ScalarTraits<ST> STS;
3997 typedef typename STS::magnitudeType MT;
3998 typedef Teuchos::ScalarTraits<MT> STM;
4003 const int myRank = comm->getRank ();
4005 if (! err.is_null ()) {
4009 *err << myRank <<
": readVectorImpl" << endl;
4011 if (! err.is_null ()) {
4045 size_t lineNumber = 1;
4048 std::ostringstream exMsg;
4049 int localBannerReadSuccess = 1;
4050 int localDimsReadSuccess = 1;
4055 *err << myRank <<
": readVectorImpl: Reading banner line (dense)" << endl;
4061 RCP<const Banner> pBanner;
4063 pBanner = readBanner (in, lineNumber, tolerant, debug);
4064 }
catch (std::exception& e) {
4066 localBannerReadSuccess = 0;
4069 if (localBannerReadSuccess) {
4070 if (pBanner->matrixType () !=
"array") {
4071 exMsg <<
"The Matrix Market file does not contain dense matrix " 4072 "data. Its banner (first) line says that its matrix type is \"" 4073 << pBanner->matrixType () <<
"\", rather that the required " 4075 localBannerReadSuccess = 0;
4076 }
else if (pBanner->dataType() ==
"pattern") {
4077 exMsg <<
"The Matrix Market file's banner (first) " 4078 "line claims that the matrix's data type is \"pattern\". This does " 4079 "not make sense for a dense matrix, yet the file reports the matrix " 4080 "as dense. The only valid data types for a dense matrix are " 4081 "\"real\", \"complex\", and \"integer\".";
4082 localBannerReadSuccess = 0;
4086 dims[2] = encodeDataType (pBanner->dataType ());
4092 if (localBannerReadSuccess) {
4094 *err << myRank <<
": readVectorImpl: Reading dimensions line (dense)" << endl;
4103 bool commentLine =
true;
4105 while (commentLine) {
4108 if (in.eof () || in.fail ()) {
4109 exMsg <<
"Unable to get array dimensions line (at all) from line " 4110 << lineNumber <<
" of input stream. The input stream " 4111 <<
"claims that it is " 4112 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4113 localDimsReadSuccess = 0;
4116 if (getline (in, line)) {
4123 size_t start = 0, size = 0;
4124 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4131 std::istringstream istr (line);
4135 if (istr.eof () || istr.fail ()) {
4136 exMsg <<
"Unable to read any data from line " << lineNumber
4137 <<
" of input; the line should contain the matrix dimensions " 4138 <<
"\"<numRows> <numCols>\".";
4139 localDimsReadSuccess = 0;
4144 exMsg <<
"Failed to get number of rows from line " 4145 << lineNumber <<
" of input; the line should contains the " 4146 <<
"matrix dimensions \"<numRows> <numCols>\".";
4147 localDimsReadSuccess = 0;
4149 dims[0] = theNumRows;
4151 exMsg <<
"No more data after number of rows on line " 4152 << lineNumber <<
" of input; the line should contain the " 4153 <<
"matrix dimensions \"<numRows> <numCols>\".";
4154 localDimsReadSuccess = 0;
4159 exMsg <<
"Failed to get number of columns from line " 4160 << lineNumber <<
" of input; the line should contain " 4161 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4162 localDimsReadSuccess = 0;
4164 dims[1] = theNumCols;
4174 exMsg <<
"File does not contain a 1D Vector";
4175 localDimsReadSuccess = 0;
4181 Tuple<GO, 5> bannerDimsReadResult;
4183 bannerDimsReadResult[0] = dims[0];
4184 bannerDimsReadResult[1] = dims[1];
4185 bannerDimsReadResult[2] = dims[2];
4186 bannerDimsReadResult[3] = localBannerReadSuccess;
4187 bannerDimsReadResult[4] = localDimsReadSuccess;
4192 broadcast (*comm, 0, bannerDimsReadResult);
4194 TEUCHOS_TEST_FOR_EXCEPTION(
4195 bannerDimsReadResult[3] == 0, std::runtime_error,
4196 "Failed to read banner line: " << exMsg.str ());
4197 TEUCHOS_TEST_FOR_EXCEPTION(
4198 bannerDimsReadResult[4] == 0, std::runtime_error,
4199 "Failed to read matrix dimensions line: " << exMsg.str ());
4201 dims[0] = bannerDimsReadResult[0];
4202 dims[1] = bannerDimsReadResult[1];
4203 dims[2] = bannerDimsReadResult[2];
4208 const size_t numCols =
static_cast<size_t> (dims[1]);
4213 RCP<const map_type> proc0Map;
4214 if (map.is_null ()) {
4218 if (node.is_null ()) {
4219 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4221 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
4223 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4225 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4226 comm, map->getNode ());
4229 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4233 RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
4239 int localReadDataSuccess = 1;
4243 *err << myRank <<
": readVectorImpl: Reading matrix data (dense)" 4248 TEUCHOS_TEST_FOR_EXCEPTION(
4249 ! X->isConstantStride (), std::logic_error,
4250 "Can't get a 1-D view of the entries of the MultiVector X on " 4251 "Process 0, because the stride between the columns of X is not " 4252 "constant. This shouldn't happen because we just created X and " 4253 "haven't filled it in yet. Please report this bug to the Tpetra " 4260 ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4261 TEUCHOS_TEST_FOR_EXCEPTION(
4262 as<global_size_t> (X_view.size ()) < numRows * numCols,
4264 "The view of X has size " << X_view <<
" which is not enough to " 4265 "accommodate the expected number of entries numRows*numCols = " 4266 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". " 4267 "Please report this bug to the Tpetra developers.");
4268 const size_t stride = X->getStride ();
4275 const bool isComplex = (dims[2] == 1);
4276 size_type count = 0, curRow = 0, curCol = 0;
4279 while (getline (in, line)) {
4283 size_t start = 0, size = 0;
4284 const bool commentLine =
4285 checkCommentLine (line, start, size, lineNumber, tolerant);
4286 if (! commentLine) {
4292 if (count >= X_view.size()) {
4297 TEUCHOS_TEST_FOR_EXCEPTION(
4298 count >= X_view.size(),
4300 "The Matrix Market input stream has more data in it than " 4301 "its metadata reported. Current line number is " 4302 << lineNumber <<
".");
4311 const size_t pos = line.substr (start, size).find (
':');
4312 if (pos != std::string::npos) {
4316 std::istringstream istr (line.substr (start, size));
4320 if (istr.eof() || istr.fail()) {
4327 TEUCHOS_TEST_FOR_EXCEPTION(
4328 ! tolerant && (istr.eof() || istr.fail()),
4330 "Line " << lineNumber <<
" of the Matrix Market file is " 4331 "empty, or we cannot read from it for some other reason.");
4334 ST val = STS::zero();
4337 MT real = STM::zero(), imag = STM::zero();
4354 TEUCHOS_TEST_FOR_EXCEPTION(
4355 ! tolerant && istr.eof(), std::runtime_error,
4356 "Failed to get the real part of a complex-valued matrix " 4357 "entry from line " << lineNumber <<
" of the Matrix Market " 4363 }
else if (istr.eof()) {
4364 TEUCHOS_TEST_FOR_EXCEPTION(
4365 ! tolerant && istr.eof(), std::runtime_error,
4366 "Missing imaginary part of a complex-valued matrix entry " 4367 "on line " << lineNumber <<
" of the Matrix Market file.");
4373 TEUCHOS_TEST_FOR_EXCEPTION(
4374 ! tolerant && istr.fail(), std::runtime_error,
4375 "Failed to get the imaginary part of a complex-valued " 4376 "matrix entry from line " << lineNumber <<
" of the " 4377 "Matrix Market file.");
4384 TEUCHOS_TEST_FOR_EXCEPTION(
4385 ! tolerant && istr.fail(), std::runtime_error,
4386 "Failed to get a real-valued matrix entry from line " 4387 << lineNumber <<
" of the Matrix Market file.");
4390 if (istr.fail() && tolerant) {
4396 TEUCHOS_TEST_FOR_EXCEPTION(
4397 ! tolerant && istr.fail(), std::runtime_error,
4398 "Failed to read matrix data from line " << lineNumber
4399 <<
" of the Matrix Market file.");
4402 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4404 curRow = count % numRows;
4405 curCol = count / numRows;
4406 X_view[curRow + curCol*stride] = val;
4411 TEUCHOS_TEST_FOR_EXCEPTION(
4412 ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
4414 "The Matrix Market metadata reports that the dense matrix is " 4415 << numRows <<
" x " << numCols <<
", and thus has " 4416 << numRows*numCols <<
" total entries, but we only found " << count
4417 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
4418 }
catch (std::exception& e) {
4420 localReadDataSuccess = 0;
4425 *err << myRank <<
": readVectorImpl: done reading data" << endl;
4429 int globalReadDataSuccess = localReadDataSuccess;
4430 broadcast (*comm, 0, outArg (globalReadDataSuccess));
4431 TEUCHOS_TEST_FOR_EXCEPTION(
4432 globalReadDataSuccess == 0, std::runtime_error,
4433 "Failed to read the multivector's data: " << exMsg.str ());
4438 if (comm->getSize () == 1 && map.is_null ()) {
4440 if (! err.is_null ()) {
4444 *err << myRank <<
": readVectorImpl: done" << endl;
4446 if (! err.is_null ()) {
4453 *err << myRank <<
": readVectorImpl: Creating target MV" << endl;
4457 RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
4460 *err << myRank <<
": readVectorImpl: Creating Export" << endl;
4469 *err << myRank <<
": readVectorImpl: Exporting" << endl;
4472 Y->doExport (*X, exporter,
INSERT);
4474 if (! err.is_null ()) {
4478 *err << myRank <<
": readVectorImpl: done" << endl;
4480 if (! err.is_null ()) {
4509 static RCP<const map_type>
4511 const RCP<const comm_type>& comm,
4512 const RCP<node_type>& node,
4513 const bool tolerant=
false,
4514 const bool debug=
false)
4516 Teuchos::RCP<Teuchos::FancyOStream> err =
4517 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4518 return readMap (in, comm, node, err, tolerant, debug);
4547 static RCP<const map_type>
4549 const RCP<const comm_type>& comm,
4550 const RCP<node_type>& node,
4551 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4552 const bool tolerant=
false,
4553 const bool debug=
false)
4555 using Teuchos::arcp;
4556 using Teuchos::Array;
4557 using Teuchos::ArrayRCP;
4559 using Teuchos::broadcast;
4560 using Teuchos::Comm;
4561 using Teuchos::CommRequest;
4562 using Teuchos::inOutArg;
4563 using Teuchos::ireceive;
4564 using Teuchos::outArg;
4565 using Teuchos::receive;
4566 using Teuchos::reduceAll;
4567 using Teuchos::REDUCE_MIN;
4568 using Teuchos::isend;
4569 using Teuchos::SerialComm;
4570 using Teuchos::toString;
4571 using Teuchos::wait;
4574 typedef ptrdiff_t int_type;
4575 typedef local_ordinal_type LO;
4576 typedef global_ordinal_type GO;
4577 typedef node_type NT;
4580 const int numProcs = comm->getSize ();
4581 const int myRank = comm->getRank ();
4583 if (err.is_null ()) {
4587 std::ostringstream os;
4588 os << myRank <<
": readMap: " << endl;
4591 if (err.is_null ()) {
4597 const int sizeTag = 1339;
4599 const int dataTag = 1340;
4605 RCP<CommRequest<int> > sizeReq;
4606 RCP<CommRequest<int> > dataReq;
4611 ArrayRCP<int_type> numGidsToRecv (1);
4612 numGidsToRecv[0] = 0;
4614 sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
4617 int readSuccess = 1;
4618 std::ostringstream exMsg;
4627 RCP<const Comm<int> > proc0Comm (
new SerialComm<int> ());
4629 RCP<const map_type> dataMap;
4633 data = readDenseImpl<GO> (in, proc0Comm, node, dataMap,
4634 err, tolerant, debug);
4636 if (data.is_null ()) {
4638 exMsg <<
"readDenseImpl() returned null." << endl;
4640 }
catch (std::exception& e) {
4642 exMsg << e.what () << endl;
4648 std::map<int, Array<GO> > pid2gids;
4653 int_type globalNumGIDs = 0;
4663 if (myRank == 0 && readSuccess == 1) {
4664 if (data->getNumVectors () == 2) {
4665 ArrayRCP<const GO> GIDs = data->getData (0);
4666 ArrayRCP<const GO> PIDs = data->getData (1);
4667 globalNumGIDs = GIDs.size ();
4671 if (globalNumGIDs > 0) {
4672 const int pid =
static_cast<int> (PIDs[0]);
4674 if (pid < 0 || pid >= numProcs) {
4676 exMsg <<
"Tpetra::MatrixMarket::readMap: " 4677 <<
"Encountered invalid PID " << pid <<
"." << endl;
4680 const GO gid = GIDs[0];
4681 pid2gids[pid].push_back (gid);
4685 if (readSuccess == 1) {
4688 for (size_type k = 1; k < globalNumGIDs; ++k) {
4689 const int pid =
static_cast<int> (PIDs[k]);
4690 if (pid < 0 || pid >= numProcs) {
4692 exMsg <<
"Tpetra::MatrixMarket::readMap: " 4693 <<
"Encountered invalid PID " << pid <<
"." << endl;
4696 const int_type gid = GIDs[k];
4697 pid2gids[pid].push_back (gid);
4698 if (gid < indexBase) {
4705 else if (data->getNumVectors () == 1) {
4706 if (data->getGlobalLength () % 2 !=
static_cast<GST
> (0)) {
4708 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data has the " 4709 "wrong format (for Map format 2.0). The global number of rows " 4710 "in the MultiVector must be even (divisible by 2)." << endl;
4713 ArrayRCP<const GO> theData = data->getData (0);
4714 globalNumGIDs =
static_cast<GO
> (data->getGlobalLength ()) /
4715 static_cast<GO> (2);
4719 if (globalNumGIDs > 0) {
4720 const int pid =
static_cast<int> (theData[1]);
4721 if (pid < 0 || pid >= numProcs) {
4723 exMsg <<
"Tpetra::MatrixMarket::readMap: " 4724 <<
"Encountered invalid PID " << pid <<
"." << endl;
4727 const GO gid = theData[0];
4728 pid2gids[pid].push_back (gid);
4734 for (int_type k = 1; k < globalNumGIDs; ++k) {
4735 const int pid =
static_cast<int> (theData[2*k + 1]);
4736 if (pid < 0 || pid >= numProcs) {
4738 exMsg <<
"Tpetra::MatrixMarket::readMap: " 4739 <<
"Encountered invalid PID " << pid <<
"." << endl;
4742 const GO gid = theData[2*k];
4743 pid2gids[pid].push_back (gid);
4744 if (gid < indexBase) {
4753 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data must have " 4754 "either 1 column (for the new Map format 2.0) or 2 columns (for " 4755 "the old Map format 1.0).";
4762 int_type readResults[3];
4763 readResults[0] =
static_cast<int_type
> (indexBase);
4764 readResults[1] =
static_cast<int_type
> (globalNumGIDs);
4765 readResults[2] =
static_cast<int_type
> (readSuccess);
4766 broadcast<int, int_type> (*comm, 0, 3, readResults);
4768 indexBase =
static_cast<GO
> (readResults[0]);
4769 globalNumGIDs =
static_cast<int_type
> (readResults[1]);
4770 readSuccess =
static_cast<int> (readResults[2]);
4776 TEUCHOS_TEST_FOR_EXCEPTION(
4777 readSuccess != 1, std::runtime_error,
4778 "Tpetra::MatrixMarket::readMap: Reading the Map failed with the " 4779 "following exception message: " << exMsg.str ());
4783 for (
int p = 1; p < numProcs; ++p) {
4784 ArrayRCP<int_type> numGidsToSend (1);
4786 typename std::map<int, Array<GO> >::const_iterator it = pid2gids.find (p);
4787 if (it == pid2gids.end ()) {
4788 numGidsToSend[0] = 0;
4790 numGidsToSend[0] = it->second.size ();
4792 sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
4793 wait<int> (*comm, outArg (sizeReq));
4798 wait<int> (*comm, outArg (sizeReq));
4803 ArrayRCP<GO> myGids;
4804 int_type myNumGids = 0;
4806 GO* myGidsRaw = NULL;
4808 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
4809 if (it != pid2gids.end ()) {
4810 myGidsRaw = it->second.getRawPtr ();
4811 myNumGids = it->second.size ();
4813 myGids = arcp<GO> (myGidsRaw, 0, myNumGids,
false);
4817 myNumGids = numGidsToRecv[0];
4818 myGids = arcp<GO> (myNumGids);
4825 if (myNumGids > 0) {
4826 dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
4830 for (
int p = 1; p < numProcs; ++p) {
4832 ArrayRCP<GO> sendGids;
4833 GO* sendGidsRaw = NULL;
4834 int_type numSendGids = 0;
4836 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
4837 if (it != pid2gids.end ()) {
4838 numSendGids = it->second.size ();
4839 sendGidsRaw = it->second.getRawPtr ();
4840 sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids,
false);
4843 if (numSendGids > 0) {
4844 dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
4846 wait<int> (*comm, outArg (dataReq));
4848 else if (myRank == p) {
4850 wait<int> (*comm, outArg (dataReq));
4855 std::ostringstream os;
4856 os << myRank <<
": readMap: creating Map" << endl;
4859 const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
4860 RCP<const map_type> newMap =
4861 rcp (
new map_type (INVALID, myGids (), indexBase, comm, node));
4863 if (err.is_null ()) {
4867 std::ostringstream os;
4868 os << myRank <<
": readMap: done" << endl;
4871 if (err.is_null ()) {
4890 encodeDataType (
const std::string& dataType)
4892 if (dataType ==
"real" || dataType ==
"integer") {
4894 }
else if (dataType ==
"complex") {
4896 }
else if (dataType ==
"pattern") {
4902 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
4903 "Unrecognized Matrix Market data type \"" << dataType
4904 <<
"\". We should never get here. " 4905 "Please report this bug to the Tpetra developers.");
4938 template<
class SparseMatrixType>
4943 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
5002 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5003 const std::string& matrixName,
5004 const std::string& matrixDescription,
5005 const bool debug=
false)
5007 TEUCHOS_TEST_FOR_EXCEPTION(
5008 pMatrix.is_null (), std::invalid_argument,
5009 "The input matrix is null.");
5010 Teuchos::RCP<const Teuchos::Comm<int> > comm = pMatrix->getComm ();
5011 TEUCHOS_TEST_FOR_EXCEPTION(
5012 comm.is_null (), std::invalid_argument,
5013 "The input matrix's communicator (Teuchos::Comm object) is null.");
5014 const int myRank = comm->getRank ();
5019 out.open (filename.c_str ());
5021 writeSparse (out, pMatrix, matrixName, matrixDescription, debug);
5049 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5050 const bool debug=
false)
5052 writeSparseFile (filename, pMatrix,
"",
"", debug);
5087 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5088 const std::string& matrixName,
5089 const std::string& matrixDescription,
5090 const bool debug=
false)
5092 using Teuchos::Comm;
5093 using Teuchos::FancyOStream;
5094 using Teuchos::getFancyOStream;
5095 using Teuchos::null;
5097 using Teuchos::rcpFromRef;
5100 typedef scalar_type ST;
5101 typedef local_ordinal_type LO;
5102 typedef global_ordinal_type GO;
5103 typedef typename Teuchos::ScalarTraits<ST> STS;
5104 typedef typename ArrayView<const LO>::const_iterator lo_iter;
5105 typedef typename ArrayView<const GO>::const_iterator go_iter;
5106 typedef typename ArrayView<const ST>::const_iterator st_iter;
5108 TEUCHOS_TEST_FOR_EXCEPTION(
5109 pMatrix.is_null (), std::invalid_argument,
5110 "The input matrix is null.");
5116 Teuchos::MatrixMarket::details::SetScientific<ST> sci (out);
5119 RCP<const Comm<int> > comm = pMatrix->getComm ();
5120 TEUCHOS_TEST_FOR_EXCEPTION(
5121 comm.is_null (), std::invalid_argument,
5122 "The input matrix's communicator (Teuchos::Comm object) is null.");
5123 const int myRank = comm->getRank ();
5126 RCP<FancyOStream> err =
5127 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
5129 std::ostringstream os;
5130 os << myRank <<
": writeSparse" << endl;
5133 os <<
"-- " << myRank <<
": past barrier" << endl;
5138 const bool debugPrint = debug && myRank == 0;
5140 RCP<const map_type> rowMap = pMatrix->getRowMap ();
5141 RCP<const map_type> colMap = pMatrix->getColMap ();
5142 RCP<const map_type> domainMap = pMatrix->getDomainMap ();
5143 RCP<const map_type> rangeMap = pMatrix->getRangeMap ();
5145 const global_size_t numRows = rangeMap->getGlobalNumElements ();
5146 const global_size_t numCols = domainMap->getGlobalNumElements ();
5148 if (debug && myRank == 0) {
5149 std::ostringstream os;
5150 os <<
"-- Input sparse matrix is:" 5151 <<
"---- " << numRows <<
" x " << numCols <<
" with " 5152 << pMatrix->getGlobalNumEntries() <<
" entries;" << endl
5154 << (pMatrix->isGloballyIndexed() ?
"Globally" :
"Locally")
5155 <<
" indexed." << endl
5156 <<
"---- Its row map has " << rowMap->getGlobalNumElements ()
5157 <<
" elements." << endl
5158 <<
"---- Its col map has " << colMap->getGlobalNumElements ()
5159 <<
" elements." << endl;
5164 const size_t localNumRows = (myRank == 0) ? numRows : 0;
5165 RCP<node_type> node = rowMap->getNode();
5167 std::ostringstream os;
5168 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
5171 RCP<const map_type> gatherRowMap =
5172 Details::computeGatherMap (rowMap, err, debug);
5180 const size_t localNumCols = (myRank == 0) ? numCols : 0;
5181 RCP<const map_type> gatherColMap =
5182 Details::computeGatherMap (colMap, err, debug);
5186 import_type importer (rowMap, gatherRowMap);
5191 RCP<sparse_matrix_type> newMatrix =
5192 rcp (
new sparse_matrix_type (gatherRowMap, gatherColMap,
5193 static_cast<size_t> (0)));
5195 newMatrix->doImport (*pMatrix, importer,
INSERT);
5200 RCP<const map_type> gatherDomainMap =
5201 rcp (
new map_type (numCols, localNumCols,
5202 domainMap->getIndexBase (),
5204 RCP<const map_type> gatherRangeMap =
5205 rcp (
new map_type (numRows, localNumRows,
5206 rangeMap->getIndexBase (),
5208 newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
5212 cerr <<
"-- Output sparse matrix is:" 5213 <<
"---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
5215 << newMatrix->getDomainMap ()->getGlobalNumElements ()
5217 << newMatrix->getGlobalNumEntries () <<
" entries;" << endl
5219 << (newMatrix->isGloballyIndexed () ?
"Globally" :
"Locally")
5220 <<
" indexed." << endl
5221 <<
"---- Its row map has " 5222 << newMatrix->getRowMap ()->getGlobalNumElements ()
5223 <<
" elements, with index base " 5224 << newMatrix->getRowMap ()->getIndexBase () <<
"." << endl
5225 <<
"---- Its col map has " 5226 << newMatrix->getColMap ()->getGlobalNumElements ()
5227 <<
" elements, with index base " 5228 << newMatrix->getColMap ()->getIndexBase () <<
"." << endl
5229 <<
"---- Element count of output matrix's column Map may differ " 5230 <<
"from that of the input matrix's column Map, if some columns " 5231 <<
"of the matrix contain no entries." << endl;
5244 out <<
"%%MatrixMarket matrix coordinate " 5245 << (STS::isComplex ?
"complex" :
"real")
5246 <<
" general" << endl;
5249 if (matrixName !=
"") {
5250 printAsComment (out, matrixName);
5252 if (matrixDescription !=
"") {
5253 printAsComment (out, matrixDescription);
5263 out << newMatrix->getRangeMap ()->getGlobalNumElements () <<
" " 5264 << newMatrix->getDomainMap ()->getGlobalNumElements () <<
" " 5265 << newMatrix->getGlobalNumEntries () << endl;
5270 const GO rowIndexBase = gatherRowMap->getIndexBase ();
5271 const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
5279 if (newMatrix->isGloballyIndexed()) {
5282 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
5283 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
5284 for (GO globalRowIndex = minAllGlobalIndex;
5285 globalRowIndex <= maxAllGlobalIndex;
5287 ArrayView<const GO> ind;
5288 ArrayView<const ST> val;
5289 newMatrix->getGlobalRowView (globalRowIndex, ind, val);
5290 go_iter indIter = ind.begin ();
5291 st_iter valIter = val.begin ();
5292 for (; indIter != ind.end() && valIter != val.end();
5293 ++indIter, ++valIter) {
5294 const GO globalColIndex = *indIter;
5297 out << (globalRowIndex + 1 - rowIndexBase) <<
" " 5298 << (globalColIndex + 1 - colIndexBase) <<
" ";
5299 if (STS::isComplex) {
5300 out << STS::real (*valIter) <<
" " << STS::imag (*valIter);
5308 typedef OrdinalTraits<GO> OTG;
5309 for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
5310 localRowIndex <= gatherRowMap->getMaxLocalIndex();
5313 const GO globalRowIndex =
5314 gatherRowMap->getGlobalElement (localRowIndex);
5315 TEUCHOS_TEST_FOR_EXCEPTION(
5316 globalRowIndex == OTG::invalid(), std::logic_error,
5317 "Failed to convert the supposed local row index " 5318 << localRowIndex <<
" into a global row index. " 5319 "Please report this bug to the Tpetra developers.");
5320 ArrayView<const LO> ind;
5321 ArrayView<const ST> val;
5322 newMatrix->getLocalRowView (localRowIndex, ind, val);
5323 lo_iter indIter = ind.begin ();
5324 st_iter valIter = val.begin ();
5325 for (; indIter != ind.end() && valIter != val.end();
5326 ++indIter, ++valIter) {
5328 const GO globalColIndex =
5329 newMatrix->getColMap()->getGlobalElement (*indIter);
5330 TEUCHOS_TEST_FOR_EXCEPTION(
5331 globalColIndex == OTG::invalid(), std::logic_error,
5332 "On local row " << localRowIndex <<
" of the sparse matrix: " 5333 "Failed to convert the supposed local column index " 5334 << *indIter <<
" into a global column index. Please report " 5335 "this bug to the Tpetra developers.");
5338 out << (globalRowIndex + 1 - rowIndexBase) <<
" " 5339 << (globalColIndex + 1 - colIndexBase) <<
" ";
5340 if (STS::isComplex) {
5341 out << STS::real (*valIter) <<
" " << STS::imag (*valIter);
5376 const RCP<const sparse_matrix_type>& pMatrix,
5377 const bool debug=
false)
5379 writeSparse (out, pMatrix,
"",
"", debug);
5412 const multivector_type& X,
5413 const std::string& matrixName,
5414 const std::string& matrixDescription,
5415 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5416 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5418 const int myRank = X.
getMap ().is_null () ? 0 :
5419 (X.
getMap ()->getComm ().is_null () ? 0 :
5420 X.
getMap ()->getComm ()->getRank ());
5424 out.open (filename.c_str());
5427 writeDense (out, X, matrixName, matrixDescription, err, dbg);
5440 const Teuchos::RCP<const multivector_type>& X,
5441 const std::string& matrixName,
5442 const std::string& matrixDescription,
5443 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5444 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5446 TEUCHOS_TEST_FOR_EXCEPTION(
5447 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 5448 "writeDenseFile: The input MultiVector X is null.");
5449 writeDenseFile (filename, *X, matrixName, matrixDescription, err, dbg);
5459 const multivector_type& X,
5460 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5461 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5463 writeDenseFile (filename, X,
"",
"", err, dbg);
5473 const RCP<const multivector_type>& X,
5474 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5475 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5477 TEUCHOS_TEST_FOR_EXCEPTION(
5478 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 5479 "writeDenseFile: The input MultiVector X is null.");
5480 writeDenseFile (filename, *X, err, dbg);
5516 const multivector_type& X,
5517 const std::string& matrixName,
5518 const std::string& matrixDescription,
5519 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5520 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5522 using Teuchos::Comm;
5523 using Teuchos::outArg;
5524 using Teuchos::REDUCE_MAX;
5525 using Teuchos::reduceAll;
5529 RCP<const Comm<int> > comm = X.
getMap ().is_null () ?
5530 Teuchos::null : X.
getMap ()->getComm ();
5531 const int myRank = comm.is_null () ? 0 : comm->getRank ();
5536 const bool debug = ! dbg.is_null ();
5539 std::ostringstream os;
5540 os << myRank <<
": writeDense" << endl;
5545 writeDenseHeader (out, X, matrixName, matrixDescription, err, dbg);
5550 for (
size_t j = 0; j < numVecs; ++j) {
5551 writeDenseColumn (out, * (X.
getVector (j)), err, dbg);
5556 std::ostringstream os;
5557 os << myRank <<
": writeDense: Done" << endl;
5591 writeDenseHeader (std::ostream& out,
5592 const multivector_type& X,
5593 const std::string& matrixName,
5594 const std::string& matrixDescription,
5595 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5596 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5598 using Teuchos::Comm;
5599 using Teuchos::outArg;
5601 using Teuchos::REDUCE_MAX;
5602 using Teuchos::reduceAll;
5605 typedef Teuchos::ScalarTraits<scalar_type> STS;
5606 const char prefix[] =
"Tpetra::MatrixMarket::writeDenseHeader: ";
5608 RCP<const Comm<int> > comm = X.
getMap ().is_null () ?
5609 Teuchos::null : X.
getMap ()->getComm ();
5610 const int myRank = comm.is_null () ? 0 : comm->getRank ();
5617 const bool debug = ! dbg.is_null ();
5621 std::ostringstream os;
5622 os << myRank <<
": writeDenseHeader" << endl;
5641 std::ostringstream hdr;
5643 std::string dataType;
5644 if (STS::isComplex) {
5645 dataType =
"complex";
5646 }
else if (STS::isOrdinal) {
5647 dataType =
"integer";
5651 hdr <<
"%%MatrixMarket matrix array " << dataType <<
" general" 5656 if (matrixName !=
"") {
5657 printAsComment (hdr, matrixName);
5659 if (matrixDescription !=
"") {
5660 printAsComment (hdr, matrixDescription);
5667 }
catch (std::exception& e) {
5668 if (! err.is_null ()) {
5669 *err << prefix <<
"While writing the Matrix Market header, " 5670 "Process 0 threw an exception: " << e.what () << endl;
5679 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
5680 TEUCHOS_TEST_FOR_EXCEPTION(
5681 gblErr == 1, std::runtime_error, prefix <<
"Some error occurred " 5682 "which prevented this method from completing.");
5686 *dbg << myRank <<
": writeDenseHeader: Done" << endl;
5709 writeDenseColumn (std::ostream& out,
5710 const multivector_type& X,
5711 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5712 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5714 using Teuchos::arcp;
5715 using Teuchos::ArrayRCP;
5716 using Teuchos::ArrayView;
5717 using Teuchos::CommRequest;
5718 using Teuchos::ireceive;
5719 using Teuchos::isend;
5720 using Teuchos::outArg;
5721 using Teuchos::REDUCE_MAX;
5722 using Teuchos::reduceAll;
5724 using Teuchos::TypeNameTraits;
5725 using Teuchos::wait;
5728 typedef Teuchos::ScalarTraits<scalar_type> STS;
5730 const Comm<int>& comm = * (X.
getMap ()->getComm ());
5731 const int myRank = comm.getRank ();
5732 const int numProcs = comm.getSize ();
5739 const bool debug = ! dbg.is_null ();
5743 std::ostringstream os;
5744 os << myRank <<
": writeDenseColumn" << endl;
5753 Teuchos::MatrixMarket::details::SetScientific<scalar_type> sci (out);
5759 const int sizeTag = 1337;
5760 const int dataTag = 1338;
5821 RCP<CommRequest<int> > sendReqSize, sendReqData;
5827 Array<ArrayRCP<size_t> > recvSizeBufs (3);
5828 Array<ArrayRCP<scalar_type> > recvDataBufs (3);
5829 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
5830 Array<RCP<CommRequest<int> > > recvDataReqs (3);
5833 ArrayRCP<size_t> sendDataSize (1);
5834 sendDataSize[0] = myNumRows;
5838 std::ostringstream os;
5839 os << myRank <<
": Post receive-size receives from " 5840 "Procs 1 and 2: tag = " << sizeTag << endl;
5844 recvSizeBufs[0].resize (1);
5848 (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
5849 recvSizeBufs[1].resize (1);
5850 (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
5851 recvSizeBufs[2].resize (1);
5852 (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
5855 ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
5859 ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
5862 else if (myRank == 1 || myRank == 2) {
5864 std::ostringstream os;
5865 os << myRank <<
": Post send-size send: size = " 5866 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
5873 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
5877 std::ostringstream os;
5878 os << myRank <<
": Not posting my send-size send yet" << endl;
5887 std::ostringstream os;
5888 os << myRank <<
": Pack my entries" << endl;
5891 ArrayRCP<scalar_type> sendDataBuf;
5893 sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
5894 X.
get1dCopy (sendDataBuf (), myNumRows);
5896 catch (std::exception& e) {
5898 if (! err.is_null ()) {
5899 std::ostringstream os;
5900 os <<
"Process " << myRank <<
": Attempt to pack my MultiVector " 5901 "entries threw an exception: " << e.what () << endl;
5906 std::ostringstream os;
5907 os << myRank <<
": Done packing my entries" << endl;
5916 *dbg << myRank <<
": Post send-data send: tag = " << dataTag
5919 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
5927 std::ostringstream os;
5928 os << myRank <<
": Write my entries" << endl;
5934 const size_t printNumRows = myNumRows;
5935 ArrayView<const scalar_type> printData = sendDataBuf ();
5936 const size_t printStride = printNumRows;
5937 if (static_cast<size_t> (printData.size ()) < printStride * numCols) {
5939 if (! err.is_null ()) {
5940 std::ostringstream os;
5941 os <<
"Process " << myRank <<
": My MultiVector data's size " 5942 << printData.size () <<
" does not match my local dimensions " 5943 << printStride <<
" x " << numCols <<
"." << endl;
5951 for (
size_t col = 0; col < numCols; ++col) {
5952 for (
size_t row = 0; row < printNumRows; ++row) {
5953 if (STS::isComplex) {
5954 out << STS::real (printData[row + col * printStride]) <<
" " 5955 << STS::imag (printData[row + col * printStride]) << endl;
5957 out << printData[row + col * printStride] << endl;
5966 const int recvRank = 1;
5967 const int circBufInd = recvRank % 3;
5969 std::ostringstream os;
5970 os << myRank <<
": Wait on receive-size receive from Process " 5971 << recvRank << endl;
5975 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
5979 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
5980 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
5982 if (! err.is_null ()) {
5983 std::ostringstream os;
5984 os << myRank <<
": Result of receive-size receive from Process " 5985 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() " 5986 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". " 5987 "This should never happen, and suggests that the receive never " 5988 "got posted. Please report this bug to the Tpetra developers." 6003 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
6005 std::ostringstream os;
6006 os << myRank <<
": Post receive-data receive from Process " 6007 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 6008 << recvDataBufs[circBufInd].size () << endl;
6011 if (! recvSizeReqs[circBufInd].is_null ()) {
6013 if (! err.is_null ()) {
6014 std::ostringstream os;
6015 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 6016 "null, before posting the receive-data receive from Process " 6017 << recvRank <<
". This should never happen. Please report " 6018 "this bug to the Tpetra developers." << endl;
6022 recvDataReqs[circBufInd] =
6023 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
6024 recvRank, dataTag, comm);
6027 else if (myRank == 1) {
6030 std::ostringstream os;
6031 os << myRank <<
": Wait on my send-size send" << endl;
6034 wait<int> (comm, outArg (sendReqSize));
6040 for (
int p = 1; p < numProcs; ++p) {
6042 if (p + 2 < numProcs) {
6044 const int recvRank = p + 2;
6045 const int circBufInd = recvRank % 3;
6047 std::ostringstream os;
6048 os << myRank <<
": Post receive-size receive from Process " 6049 << recvRank <<
": tag = " << sizeTag << endl;
6052 if (! recvSizeReqs[circBufInd].is_null ()) {
6054 if (! err.is_null ()) {
6055 std::ostringstream os;
6056 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 6057 <<
"null, for the receive-size receive from Process " 6058 << recvRank <<
"! This may mean that this process never " 6059 <<
"finished waiting for the receive from Process " 6060 << (recvRank - 3) <<
"." << endl;
6064 recvSizeReqs[circBufInd] =
6065 ireceive<int, size_t> (recvSizeBufs[circBufInd],
6066 recvRank, sizeTag, comm);
6069 if (p + 1 < numProcs) {
6070 const int recvRank = p + 1;
6071 const int circBufInd = recvRank % 3;
6075 std::ostringstream os;
6076 os << myRank <<
": Wait on receive-size receive from Process " 6077 << recvRank << endl;
6080 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
6084 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
6085 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
6087 if (! err.is_null ()) {
6088 std::ostringstream os;
6089 os << myRank <<
": Result of receive-size receive from Process " 6090 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() " 6091 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". " 6092 "This should never happen, and suggests that the receive never " 6093 "got posted. Please report this bug to the Tpetra developers." 6107 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
6109 std::ostringstream os;
6110 os << myRank <<
": Post receive-data receive from Process " 6111 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 6112 << recvDataBufs[circBufInd].size () << endl;
6115 if (! recvDataReqs[circBufInd].is_null ()) {
6117 if (! err.is_null ()) {
6118 std::ostringstream os;
6119 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not " 6120 <<
"null, for the receive-data receive from Process " 6121 << recvRank <<
"! This may mean that this process never " 6122 <<
"finished waiting for the receive from Process " 6123 << (recvRank - 3) <<
"." << endl;
6127 recvDataReqs[circBufInd] =
6128 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
6129 recvRank, dataTag, comm);
6133 const int recvRank = p;
6134 const int circBufInd = recvRank % 3;
6136 std::ostringstream os;
6137 os << myRank <<
": Wait on receive-data receive from Process " 6138 << recvRank << endl;
6141 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
6148 std::ostringstream os;
6149 os << myRank <<
": Write entries from Process " << recvRank
6151 *dbg << os.str () << endl;
6153 size_t printNumRows = (recvSizeBufs[circBufInd])[0];
6154 if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
6156 if (! err.is_null ()) {
6157 std::ostringstream os;
6158 os << myRank <<
": Result of receive-size receive from Process " 6159 << recvRank <<
" was Teuchos::OrdinalTraits<size_t>::" 6160 "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
6161 <<
". This should never happen, and suggests that its " 6162 "receive-size receive was never posted. " 6163 "Please report this bug to the Tpetra developers." << endl;
6170 if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
6172 if (! err.is_null ()) {
6173 std::ostringstream os;
6174 os << myRank <<
": Result of receive-size receive from Proc " 6175 << recvRank <<
" was " << printNumRows <<
" > 0, but " 6176 "recvDataBufs[" << circBufInd <<
"] is null. This should " 6177 "never happen. Please report this bug to the Tpetra " 6178 "developers." << endl;
6188 ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
6189 const size_t printStride = printNumRows;
6193 for (
size_t col = 0; col < numCols; ++col) {
6194 for (
size_t row = 0; row < printNumRows; ++row) {
6195 if (STS::isComplex) {
6196 out << STS::real (printData[row + col * printStride]) <<
" " 6197 << STS::imag (printData[row + col * printStride]) << endl;
6199 out << printData[row + col * printStride] << endl;
6204 else if (myRank == p) {
6207 std::ostringstream os;
6208 os << myRank <<
": Wait on my send-data send" << endl;
6211 wait<int> (comm, outArg (sendReqData));
6213 else if (myRank == p + 1) {
6216 std::ostringstream os;
6217 os << myRank <<
": Post send-data send: tag = " << dataTag
6221 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
6224 std::ostringstream os;
6225 os << myRank <<
": Wait on my send-size send" << endl;
6228 wait<int> (comm, outArg (sendReqSize));
6230 else if (myRank == p + 2) {
6233 std::ostringstream os;
6234 os << myRank <<
": Post send-size send: size = " 6235 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
6238 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
6243 reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
6244 TEUCHOS_TEST_FOR_EXCEPTION(
6245 gblErr == 1, std::runtime_error,
"Tpetra::MatrixMarket::writeDense " 6246 "experienced some kind of error and was unable to complete.");
6250 *dbg << myRank <<
": writeDenseColumn: Done" << endl;
6264 const RCP<const multivector_type>& X,
6265 const std::string& matrixName,
6266 const std::string& matrixDescription,
6267 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6268 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6270 TEUCHOS_TEST_FOR_EXCEPTION(
6271 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 6272 "writeDense: The input MultiVector X is null.");
6273 writeDense (out, *X, matrixName, matrixDescription, err, dbg);
6283 const multivector_type& X,
6284 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6285 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6287 writeDense (out, X,
"",
"", err, dbg);
6297 const RCP<const multivector_type>& X,
6298 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6299 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6301 TEUCHOS_TEST_FOR_EXCEPTION(
6302 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::" 6303 "writeDense: The input MultiVector X is null.");
6304 writeDense (out, *X,
"",
"", err, dbg);
6327 writeMap (std::ostream& out,
const map_type& map,
const bool debug=
false)
6329 Teuchos::RCP<Teuchos::FancyOStream> err =
6330 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
6331 writeMap (out, map, err, debug);
6344 const map_type& map,
6345 const Teuchos::RCP<Teuchos::FancyOStream>& err,
6346 const bool debug=
false)
6348 using Teuchos::ArrayRCP;
6349 using Teuchos::ArrayView;
6350 using Teuchos::CommRequest;
6351 using Teuchos::ireceive;
6352 using Teuchos::isend;
6354 using Teuchos::TypeNameTraits;
6355 using Teuchos::wait;
6357 typedef global_ordinal_type GO;
6358 typedef int pid_type;
6373 typedef ptrdiff_t int_type;
6374 TEUCHOS_TEST_FOR_EXCEPTION(
6375 sizeof (GO) >
sizeof (int_type), std::logic_error,
6376 "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
6377 <<
" is too big for ptrdiff_t. sizeof(GO) = " <<
sizeof (GO)
6378 <<
" > sizeof(ptrdiff_t) = " <<
sizeof (ptrdiff_t) <<
".");
6379 TEUCHOS_TEST_FOR_EXCEPTION(
6380 sizeof (pid_type) >
sizeof (int_type), std::logic_error,
6381 "The (MPI) process rank type pid_type=" <<
6382 TypeNameTraits<pid_type>::name () <<
" is too big for ptrdiff_t. " 6383 "sizeof(pid_type) = " <<
sizeof (pid_type) <<
" > sizeof(ptrdiff_t)" 6384 " = " <<
sizeof (ptrdiff_t) <<
".");
6386 const Comm<int>& comm = * (map.
getComm ());
6387 const int myRank = comm.getRank ();
6388 const int numProcs = comm.getSize ();
6390 if (! err.is_null ()) {
6394 std::ostringstream os;
6395 os << myRank <<
": writeMap" << endl;
6398 if (! err.is_null ()) {
6405 const int sizeTag = 1337;
6406 const int dataTag = 1338;
6465 RCP<CommRequest<int> > sendReqSize, sendReqData;
6471 Array<ArrayRCP<int_type> > recvSizeBufs (3);
6472 Array<ArrayRCP<int_type> > recvDataBufs (3);
6473 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
6474 Array<RCP<CommRequest<int> > > recvDataReqs (3);
6477 ArrayRCP<int_type> sendDataSize (1);
6478 sendDataSize[0] = myNumRows;
6482 std::ostringstream os;
6483 os << myRank <<
": Post receive-size receives from " 6484 "Procs 1 and 2: tag = " << sizeTag << endl;
6488 recvSizeBufs[0].resize (1);
6489 (recvSizeBufs[0])[0] = -1;
6490 recvSizeBufs[1].resize (1);
6491 (recvSizeBufs[1])[0] = -1;
6492 recvSizeBufs[2].resize (1);
6493 (recvSizeBufs[2])[0] = -1;
6496 ireceive<int, int_type> (recvSizeBufs[1], 1, sizeTag, comm);
6500 ireceive<int, int_type> (recvSizeBufs[2], 2, sizeTag, comm);
6503 else if (myRank == 1 || myRank == 2) {
6505 std::ostringstream os;
6506 os << myRank <<
": Post send-size send: size = " 6507 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
6514 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
6518 std::ostringstream os;
6519 os << myRank <<
": Not posting my send-size send yet" << endl;
6530 std::ostringstream os;
6531 os << myRank <<
": Pack my GIDs and PIDs" << endl;
6535 ArrayRCP<int_type> sendDataBuf (myNumRows * 2);
6538 const int_type myMinGblIdx =
6540 for (
size_t k = 0; k < myNumRows; ++k) {
6541 const int_type gid = myMinGblIdx +
static_cast<int_type
> (k);
6542 const int_type pid =
static_cast<int_type
> (myRank);
6543 sendDataBuf[2*k] = gid;
6544 sendDataBuf[2*k+1] = pid;
6549 for (
size_t k = 0; k < myNumRows; ++k) {
6550 const int_type gid =
static_cast<int_type
> (myGblInds[k]);
6551 const int_type pid =
static_cast<int_type
> (myRank);
6552 sendDataBuf[2*k] = gid;
6553 sendDataBuf[2*k+1] = pid;
6558 std::ostringstream os;
6559 os << myRank <<
": Done packing my GIDs and PIDs" << endl;
6566 *err << myRank <<
": Post send-data send: tag = " << dataTag
6569 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
6574 *err << myRank <<
": Write MatrixMarket header" << endl;
6579 std::ostringstream hdr;
6583 hdr <<
"%%MatrixMarket matrix array integer general" << endl
6584 <<
"% Format: Version 2.0" << endl
6586 <<
"% This file encodes a Tpetra::Map." << endl
6587 <<
"% It is stored as a dense vector, with twice as many " << endl
6588 <<
"% entries as the global number of GIDs (global indices)." << endl
6589 <<
"% (GID, PID) pairs are stored contiguously, where the PID " << endl
6590 <<
"% is the rank of the process owning that GID." << endl
6595 std::ostringstream os;
6596 os << myRank <<
": Write my GIDs and PIDs" << endl;
6602 const int_type printNumRows = myNumRows;
6603 ArrayView<const int_type> printData = sendDataBuf ();
6604 for (int_type k = 0; k < printNumRows; ++k) {
6605 const int_type gid = printData[2*k];
6606 const int_type pid = printData[2*k+1];
6607 out << gid << endl << pid << endl;
6613 const int recvRank = 1;
6614 const int circBufInd = recvRank % 3;
6616 std::ostringstream os;
6617 os << myRank <<
": Wait on receive-size receive from Process " 6618 << recvRank << endl;
6622 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
6626 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
6627 if (debug && recvNumRows == -1) {
6628 std::ostringstream os;
6629 os << myRank <<
": Result of receive-size receive from Process " 6630 << recvRank <<
" is -1. This should never happen, and " 6631 "suggests that the receive never got posted. Please report " 6632 "this bug to the Tpetra developers." << endl;
6637 recvDataBufs[circBufInd].resize (recvNumRows * 2);
6639 std::ostringstream os;
6640 os << myRank <<
": Post receive-data receive from Process " 6641 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 6642 << recvDataBufs[circBufInd].size () << endl;
6645 if (! recvSizeReqs[circBufInd].is_null ()) {
6646 std::ostringstream os;
6647 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 6648 "null, before posting the receive-data receive from Process " 6649 << recvRank <<
". This should never happen. Please report " 6650 "this bug to the Tpetra developers." << endl;
6653 recvDataReqs[circBufInd] =
6654 ireceive<int, int_type> (recvDataBufs[circBufInd],
6655 recvRank, dataTag, comm);
6658 else if (myRank == 1) {
6661 std::ostringstream os;
6662 os << myRank <<
": Wait on my send-size send" << endl;
6665 wait<int> (comm, outArg (sendReqSize));
6671 for (
int p = 1; p < numProcs; ++p) {
6673 if (p + 2 < numProcs) {
6675 const int recvRank = p + 2;
6676 const int circBufInd = recvRank % 3;
6678 std::ostringstream os;
6679 os << myRank <<
": Post receive-size receive from Process " 6680 << recvRank <<
": tag = " << sizeTag << endl;
6683 if (! recvSizeReqs[circBufInd].is_null ()) {
6684 std::ostringstream os;
6685 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not " 6686 <<
"null, for the receive-size receive from Process " 6687 << recvRank <<
"! This may mean that this process never " 6688 <<
"finished waiting for the receive from Process " 6689 << (recvRank - 3) <<
"." << endl;
6692 recvSizeReqs[circBufInd] =
6693 ireceive<int, int_type> (recvSizeBufs[circBufInd],
6694 recvRank, sizeTag, comm);
6697 if (p + 1 < numProcs) {
6698 const int recvRank = p + 1;
6699 const int circBufInd = recvRank % 3;
6703 std::ostringstream os;
6704 os << myRank <<
": Wait on receive-size receive from Process " 6705 << recvRank << endl;
6708 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
6712 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
6713 if (debug && recvNumRows == -1) {
6714 std::ostringstream os;
6715 os << myRank <<
": Result of receive-size receive from Process " 6716 << recvRank <<
" is -1. This should never happen, and " 6717 "suggests that the receive never got posted. Please report " 6718 "this bug to the Tpetra developers." << endl;
6723 recvDataBufs[circBufInd].resize (recvNumRows * 2);
6725 std::ostringstream os;
6726 os << myRank <<
": Post receive-data receive from Process " 6727 << recvRank <<
": tag = " << dataTag <<
", buffer size = " 6728 << recvDataBufs[circBufInd].size () << endl;
6731 if (! recvDataReqs[circBufInd].is_null ()) {
6732 std::ostringstream os;
6733 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not " 6734 <<
"null, for the receive-data receive from Process " 6735 << recvRank <<
"! This may mean that this process never " 6736 <<
"finished waiting for the receive from Process " 6737 << (recvRank - 3) <<
"." << endl;
6740 recvDataReqs[circBufInd] =
6741 ireceive<int, int_type> (recvDataBufs[circBufInd],
6742 recvRank, dataTag, comm);
6746 const int recvRank = p;
6747 const int circBufInd = recvRank % 3;
6749 std::ostringstream os;
6750 os << myRank <<
": Wait on receive-data receive from Process " 6751 << recvRank << endl;
6754 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
6761 std::ostringstream os;
6762 os << myRank <<
": Write GIDs and PIDs from Process " 6763 << recvRank << endl;
6764 *err << os.str () << endl;
6766 const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
6767 if (debug && printNumRows == -1) {
6768 std::ostringstream os;
6769 os << myRank <<
": Result of receive-size receive from Process " 6770 << recvRank <<
" was -1. This should never happen, and " 6771 "suggests that its receive-size receive was never posted. " 6772 "Please report this bug to the Tpetra developers." << endl;
6775 if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
6776 std::ostringstream os;
6777 os << myRank <<
": Result of receive-size receive from Proc " 6778 << recvRank <<
" was " << printNumRows <<
" > 0, but " 6779 "recvDataBufs[" << circBufInd <<
"] is null. This should " 6780 "never happen. Please report this bug to the Tpetra " 6781 "developers." << endl;
6784 ArrayView<const int_type> printData = (recvDataBufs[circBufInd]) ();
6785 for (int_type k = 0; k < printNumRows; ++k) {
6786 const int_type gid = printData[2*k];
6787 const int_type pid = printData[2*k+1];
6788 out << gid << endl << pid << endl;
6791 else if (myRank == p) {
6794 std::ostringstream os;
6795 os << myRank <<
": Wait on my send-data send" << endl;
6798 wait<int> (comm, outArg (sendReqData));
6800 else if (myRank == p + 1) {
6803 std::ostringstream os;
6804 os << myRank <<
": Post send-data send: tag = " << dataTag
6808 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
6811 std::ostringstream os;
6812 os << myRank <<
": Wait on my send-size send" << endl;
6815 wait<int> (comm, outArg (sendReqSize));
6817 else if (myRank == p + 2) {
6820 std::ostringstream os;
6821 os << myRank <<
": Post send-size send: size = " 6822 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
6825 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
6829 if (! err.is_null ()) {
6833 *err << myRank <<
": writeMap: Done" << endl;
6835 if (! err.is_null ()) {
6843 const map_type& map)
6845 const int myRank = map.
getComm ()->getRank ();
6848 out.open (filename.c_str());
6850 writeMap (out, map);
6881 printAsComment (std::ostream& out,
const std::string& str)
6884 std::istringstream inpstream (str);
6887 while (getline (inpstream, line)) {
6888 if (! line.empty()) {
6891 if (line[0] ==
'%') {
6892 out << line << endl;
6895 out <<
"%% " << line << endl;
6923 Teuchos::ParameterList pl;
6924 writeOperator(fileName, A, pl);
6949 Teuchos::ParameterList pl;
6950 writeOperator (out, A, pl);
6991 const operator_type& A,
6992 const Teuchos::ParameterList& params)
6995 std::string tmpFile =
"__TMP__" + fileName;
6996 const int myRank = A.
getDomainMap()->getComm()->getRank();
6997 bool precisionChanged=
false;
7007 if (std::ifstream(tmpFile))
7008 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
7009 "writeOperator: temporary file " << tmpFile <<
" already exists");
7010 out.open(tmpFile.c_str());
7011 if (params.isParameter(
"precision")) {
7012 oldPrecision = out.precision(params.get<
int>(
"precision"));
7013 precisionChanged=
true;
7017 const std::string header = writeOperatorImpl(out, A, params);
7020 if (precisionChanged)
7021 out.precision(oldPrecision);
7023 out.open(fileName.c_str(), std::ios::binary);
7024 bool printMatrixMarketHeader =
true;
7025 if (params.isParameter(
"print MatrixMarket header"))
7026 printMatrixMarketHeader = params.get<
bool>(
"print MatrixMarket header");
7027 if (printMatrixMarketHeader && myRank == 0) {
7032 std::ifstream src(tmpFile, std::ios_base::binary);
7036 remove(tmpFile.c_str());
7080 const operator_type& A,
7081 const Teuchos::ParameterList& params)
7083 const int myRank = A.
getDomainMap ()->getComm ()->getRank ();
7100 std::ostringstream tmpOut;
7102 if (params.isParameter (
"precision") && params.isType<
int> (
"precision")) {
7103 (void) tmpOut.precision (params.get<
int> (
"precision"));
7107 const std::string header = writeOperatorImpl (tmpOut, A, params);
7110 bool printMatrixMarketHeader =
true;
7111 if (params.isParameter (
"print MatrixMarket header") &&
7112 params.isType<
bool> (
"print MatrixMarket header")) {
7113 printMatrixMarketHeader = params.get<
bool> (
"print MatrixMarket header");
7115 if (printMatrixMarketHeader && myRank == 0) {
7131 out << tmpOut.str ();
7145 writeOperatorImpl (std::ostream& os,
7146 const operator_type& A,
7147 const Teuchos::ParameterList& params)
7151 using Teuchos::ArrayRCP;
7152 using Teuchos::Array;
7154 typedef local_ordinal_type LO;
7155 typedef global_ordinal_type GO;
7156 typedef scalar_type Scalar;
7157 typedef Teuchos::OrdinalTraits<LO> TLOT;
7158 typedef Teuchos::OrdinalTraits<GO> TGOT;
7164 RCP<const Teuchos::Comm<int> > comm = rangeMap->getComm();
7165 const int myRank = comm->getRank();
7166 const size_t numProcs = comm->getSize();
7169 if (params.isParameter(
"probing size"))
7170 numMVs = params.get<
int>(
"probing size");
7179 GO numGlobElts = maxColGid - minColGid + TGOT::one();
7180 GO numChunks = numGlobElts / numMVs;
7181 GO rem = numGlobElts % numMVs;
7182 GO indexBase = rangeMap->getIndexBase();
7184 int offsetToUseInPrinting = 1 - indexBase;
7185 if (params.isParameter(
"zero-based indexing")) {
7186 if (params.get<
bool>(
"zero-based indexing") ==
true)
7187 offsetToUseInPrinting = -indexBase;
7191 size_t numLocalRangeEntries = rangeMap->getNodeNumElements();
7194 RCP<const map_type> allGidsMap = rcp(
new map_type(TGOT::invalid(), numLocalRangeEntries,
7197 mv_type_go allGids(allGidsMap,1);
7198 Teuchos::ArrayRCP<GO> allGidsData = allGids.getDataNonConst(0);
7200 for (
size_t i=0; i<numLocalRangeEntries; i++)
7201 allGidsData[i] = rangeMap->getGlobalElement(i);
7202 allGidsData = Teuchos::null;
7205 GO numTargetMapEntries=TGOT::zero();
7206 Teuchos::Array<GO> importGidList;
7208 numTargetMapEntries = rangeMap->getGlobalNumElements();
7209 importGidList.reserve(numTargetMapEntries);
7210 for (GO j=0; j<numTargetMapEntries; ++j) importGidList.push_back(j + indexBase);
7212 importGidList.reserve(numTargetMapEntries);
7214 RCP<map_type> importGidMap = rcp(
new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
7217 import_type gidImporter(allGidsMap, importGidMap);
7218 mv_type_go importedGids(importGidMap, 1);
7219 importedGids.doImport(allGids, gidImporter,
INSERT);
7222 ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
7223 RCP<const map_type> importMap = rcp(
new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm) );
7226 import_type importer(rangeMap, importMap);
7228 RCP<mv_type> colsOnPid0 = rcp(
new mv_type(importMap,numMVs));
7230 RCP<mv_type> ei = rcp(
new mv_type(A.
getDomainMap(),numMVs));
7231 RCP<mv_type> colsA = rcp(
new mv_type(A.
getRangeMap(),numMVs));
7233 Array<GO> globalColsArray, localColsArray;
7234 globalColsArray.reserve(numMVs);
7235 localColsArray.reserve(numMVs);
7237 ArrayRCP<ArrayRCP<Scalar> > eiData(numMVs);
7238 for (
size_t i=0; i<numMVs; ++i)
7239 eiData[i] = ei->getDataNonConst(i);
7244 for (GO k=0; k<numChunks; ++k) {
7245 for (
size_t j=0; j<numMVs; ++j ) {
7247 GO curGlobalCol = minColGid + k*numMVs + j;
7248 globalColsArray.push_back(curGlobalCol);
7251 if (curLocalCol != TLOT::invalid()) {
7252 eiData[j][curLocalCol] = TGOT::one();
7253 localColsArray.push_back(curLocalCol);
7259 A.
apply(*ei,*colsA);
7261 colsOnPid0->doImport(*colsA,importer,
INSERT);
7264 globalNnz += writeColumns(os,*colsOnPid0, numMVs, importedGidsData(),
7265 globalColsArray, offsetToUseInPrinting);
7268 for (
size_t j=0; j<numMVs; ++j ) {
7269 for (
int i=0; i<localColsArray.size(); ++i)
7270 eiData[j][localColsArray[i]] = TGOT::zero();
7272 globalColsArray.clear();
7273 localColsArray.clear();
7281 for (
int j=0; j<rem; ++j ) {
7282 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
7283 globalColsArray.push_back(curGlobalCol);
7286 if (curLocalCol != TLOT::invalid()) {
7287 eiData[j][curLocalCol] = TGOT::one();
7288 localColsArray.push_back(curLocalCol);
7294 A.
apply(*ei,*colsA);
7296 colsOnPid0->doImport(*colsA,importer,
INSERT);
7298 globalNnz += writeColumns(os,*colsOnPid0, rem, importedGidsData(),
7299 globalColsArray, offsetToUseInPrinting);
7302 for (
int j=0; j<rem; ++j ) {
7303 for (
int i=0; i<localColsArray.size(); ++i)
7304 eiData[j][localColsArray[i]] = TGOT::zero();
7306 globalColsArray.clear();
7307 localColsArray.clear();
7316 std::ostringstream oss;
7318 oss <<
"%%MatrixMarket matrix coordinate ";
7319 if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
7324 oss <<
" general" << std::endl;
7325 oss <<
"% Tpetra::Operator" << std::endl;
7326 std::time_t now = std::time(NULL);
7327 oss <<
"% time stamp: " << ctime(&now);
7328 oss <<
"% written from " << numProcs <<
" processes" << std::endl;
7329 size_t numRows = rangeMap->getGlobalNumElements();
7331 oss << numRows <<
" " << numCols <<
" " << globalNnz << std::endl;
7337 static global_ordinal_type
7338 writeColumns(std::ostream& os, mv_type
const &colsA,
size_t const &numCols,
7339 Teuchos::ArrayView<const global_ordinal_type>
const &rowGids,
7340 Teuchos::Array<global_ordinal_type>
const &colsArray,
7341 global_ordinal_type
const & indexBase) {
7343 typedef global_ordinal_type GO;
7344 typedef scalar_type Scalar;
7345 typedef Teuchos::ScalarTraits<Scalar> STS;
7348 const Scalar zero = STS::zero();
7350 for (
size_t j=0; j<numCols; ++j) {
7351 ArrayRCP<const Scalar>
const curCol = colsA.
getData(j);
7352 const GO J = colsArray[j];
7353 for (
size_t i=0; i<numRows; ++i) {
7354 const Scalar val = curCol[i];
7356 os << rowGids[i]+indexBase <<
" " << J+indexBase <<
" " << val << std::endl;
7373 #endif // __MatrixMarket_Tpetra_hpp size_t getLocalLength() const
Local number of rows on the calling process.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
static RCP< multivector_type > readDense(std::istream &in, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market input stream.
void get1dCopy(const Teuchos::ArrayView< Scalar > &A, const size_t LDA) const
Fill the given array with a copy of this multivector's local values.
static void writeOperator(std::ostream &out, const operator_type &A)
Write a Tpetra::Operator to an output stream.
GlobalOrdinal getMinGlobalIndex() const
The minimum global index owned by the calling process.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
SparseMatrixType::node_type node_type
The fourth template parameter of CrsMatrix and MultiVector.
static void writeDense(std::ostream &out, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
bool isContiguous() const
True if this Map is distributed contiguously, else false.
static RCP< multivector_type > readDense(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Variant of readDense (see above) that takes a Node.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const RCP< const Comm< int > > &pComm, const RCP< node_type > &pNode, const RCP< Teuchos::ParameterList > &constructorParams, const RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of the above readSparse() method that takes a Kokkos Node.
static void writeDenseFile(const std::string &filename, const RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of readSparseFile above that takes a Node object.
SparseMatrixType::local_ordinal_type local_ordinal_type
SparseMatrixType::global_ordinal_type global_ordinal_type
One or more distributed dense vectors.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream.
SparseMatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the sparse matrix.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
static void writeDense(std::ostream &out, const RCP< const multivector_type > &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and or description.
virtual void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const =0
Computes the operator-multivector application.
static void writeDenseFile(const std::string &filename, const Teuchos::RCP< const multivector_type > &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
SparseMatrixType::scalar_type scalar_type
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
Specialization of Tpetra::MultiVector that matches SparseMatrixType.
static RCP< const map_type > readMapFile(const std::string &filename, const RCP< const comm_type > &comm, const RCP< node_type > &node, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given Matrix Market file.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
GlobalOrdinal getMinAllGlobalIndex() const
The minimum global index over all processes in the communicator.
static void writeMap(std::ostream &out, const map_type &map, const bool debug=false)
Print the Map to the given output stream.
static void writeDenseFile(const std::string &filename, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > getVector(const size_t j) const
Return a Vector which is a const view of column j.
SparseMatrixType::global_ordinal_type global_ordinal_type
Type of indices as read from the Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream, with provided Maps. ...
GlobalOrdinal getMaxAllGlobalIndex() const
The maximum global index over all processes in the communicator.
size_t global_size_t
Global size_t object.
SparseMatrixType::node_type node_type
The Kokkos Node type; fourth template parameter of Tpetra::CrsMatrix.
static void writeOperator(const std::string &fileName, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to a file, with options.
Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Specialization of Tpetra::Map that matches SparseMatrixType.
Insert new values that don't currently exist.
SparseMatrixType::scalar_type scalar_type
Type of the entries of the sparse matrix.
Scalar scalar_type
This class' first template parameter; the type of each entry in the MultiVector.
static void writeMapFile(const std::string &filename, const map_type &map)
Write the Map to the given file.
Abstract interface for operators (e.g., matrices and preconditioners).
static RCP< vector_type > readVector(std::istream &in, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
static RCP< multivector_type > readDenseFile(const std::string &filename, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Variant of readDenseMatrix (see above) that takes a Node.
From a distributed map build a map with all GIDs on the root node.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
static RCP< multivector_type > readDenseFile(const std::string &filename, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market file.
static void writeMap(std::ostream &out, const map_type &map, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool debug=false)
Print the Map to the given output stream out.
Teuchos::ArrayRCP< const Scalar > getData(size_t j) const
Const view of the local values in a particular vector of this multivector.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const RCP< const Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const RCP< const Comm< int > > &pComm, const RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparseFile that takes a Node object.
static void writeDense(std::ostream &out, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static void writeSparse(std::ostream &out, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
global_size_t getGlobalLength() const
Global number of rows in the multivector.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
The MultiVector specialization associated with SparseMatrixType.
LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const
The local index corresponding to the given global index.
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a view of the global indices owned by this process.
Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > vector_type
The Vector specialization associated with SparseMatrixType.
size_t getNodeNumElements() const
The number of elements belonging to the calling process.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const RCP< const Comm< int > > &pComm, const RCP< Teuchos::ParameterList > &constructorParams, const RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream.
Describes a parallel distribution of objects over processes.
static void writeDense(std::ostream &out, const RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static RCP< vector_type > readVector(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream, with a supplied Node.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparse() above that takes a Node object.
A distributed dense vector.
static RCP< const map_type > readMap(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream.
size_t getNumVectors() const
Number of columns in the multivector.
SparseMatrixType sparse_matrix_type
This class' template parameter; a specialization of CrsMatrix.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
SparseMatrixType sparse_matrix_type
Template parameter of this class; specialization of CrsMatrix.
static void writeSparse(std::ostream &out, const RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
static void writeDenseFile(const std::string &filename, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
static RCP< vector_type > readVectorFile(const std::string &filename, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Like readVectorFile() (see above), but with a supplied Node object.
static RCP< vector_type > readVectorFile(const std::string &filename, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read a Vector from the given Matrix Market file.
static void writeOperator(std::ostream &out, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to an output stream, with options.
Matrix Market file reader for CrsMatrix and MultiVector.
static RCP< const map_type > readMap(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream, with optional debugging output stream...
Matrix Market file writer for CrsMatrix and MultiVector.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const RCP< const map_type > &rowMap, RCP< const map_type > &colMap, const RCP< const map_type > &domainMap, const RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file, with provided Maps.
Matrix Market file readers and writers for sparse and dense matrices (as CrsMatrix resp...
global_size_t getGlobalNumElements() const
The number of elements in this Map.
static void writeOperator(const std::string &fileName, operator_type const &A)
Write a Tpetra::Operator to a file.