42#ifndef __MatrixMarket_Tpetra_hpp
43#define __MatrixMarket_Tpetra_hpp
57#include "Tpetra_CrsMatrix.hpp"
58#include "Tpetra_Operator.hpp"
59#include "Tpetra_Vector.hpp"
61#include "Teuchos_MatrixMarket_Raw_Adder.hpp"
62#include "Teuchos_MatrixMarket_Raw_Graph_Adder.hpp"
63#include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp"
64#include "Teuchos_MatrixMarket_SymmetrizingGraphAdder.hpp"
65#include "Teuchos_MatrixMarket_assignScalar.hpp"
66#include "Teuchos_MatrixMarket_Banner.hpp"
67#include "Teuchos_MatrixMarket_CoordDataReader.hpp"
68#include "Teuchos_SetScientific.hpp"
69#include "Teuchos_TimeMonitor.hpp"
72#include "mmio_Tpetra.h"
74#include "Tpetra_Distribution.hpp"
171 template<
class SparseMatrixType>
176 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
191 typedef typename SparseMatrixType::global_ordinal_type
213 typedef Teuchos::Comm<int> comm_type;
223 typedef Teuchos::ArrayRCP<int>::size_type size_type;
235 static Teuchos::RCP<const map_type>
236 makeRangeMap (
const Teuchos::RCP<const comm_type>& pComm,
242 pComm, GloballyDistributed));
272 static Teuchos::RCP<const map_type>
273 makeRowMap (
const Teuchos::RCP<const map_type>& pRowMap,
274 const Teuchos::RCP<const comm_type>& pComm,
279 if (pRowMap.is_null ()) {
280 return rcp (
new map_type (
static_cast<global_size_t> (numRows),
282 pComm, GloballyDistributed));
285 TEUCHOS_TEST_FOR_EXCEPTION
286 (! pRowMap->isDistributed () && pComm->getSize () > 1,
287 std::invalid_argument,
"The specified row map is not distributed, "
288 "but the given communicator includes more than one process (in "
289 "fact, there are " << pComm->getSize () <<
" processes).");
290 TEUCHOS_TEST_FOR_EXCEPTION
291 (pRowMap->getComm () != pComm, std::invalid_argument,
292 "The specified row Map's communicator (pRowMap->getComm()) "
293 "differs from the given separately supplied communicator pComm.");
312 static Teuchos::RCP<const map_type>
313 makeDomainMap (
const Teuchos::RCP<const map_type>& pRangeMap,
322 if (numRows == numCols) {
325 return createUniformContigMapWithNode<LO,GO,NT> (numCols,
326 pRangeMap->getComm ());
403 distribute (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
404 Teuchos::ArrayRCP<size_t>& myRowPtr,
405 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
406 Teuchos::ArrayRCP<scalar_type>& myValues,
407 const Teuchos::RCP<const map_type>& pRowMap,
408 Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
409 Teuchos::ArrayRCP<size_t>& rowPtr,
410 Teuchos::ArrayRCP<global_ordinal_type>& colInd,
411 Teuchos::ArrayRCP<scalar_type>& values,
412 const bool debug=
false)
415 using Teuchos::ArrayRCP;
416 using Teuchos::ArrayView;
419 using Teuchos::CommRequest;
422 using Teuchos::receive;
427 const bool extraDebug =
false;
428 RCP<const comm_type> pComm = pRowMap->getComm ();
429 const int numProcs = pComm->getSize ();
430 const int myRank = pComm->getRank ();
431 const int rootRank = 0;
438 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
439 const size_type myNumRows = myRows.size();
440 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<size_t>(myNumRows) !=
441 pRowMap->getLocalNumElements(),
443 "pRowMap->getLocalElementList().size() = "
445 <<
" != pRowMap->getLocalNumElements() = "
446 << pRowMap->getLocalNumElements() <<
". "
447 "Please report this bug to the Tpetra developers.");
448 TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
450 "On Proc 0: numEntriesPerRow.size() = "
451 << numEntriesPerRow.size()
452 <<
" != pRowMap->getLocalElementList().size() = "
453 << myNumRows <<
". Please report this bug to the "
454 "Tpetra developers.");
458 myNumEntriesPerRow = arcp<size_t> (myNumRows);
460 if (myRank != rootRank) {
464 send (*pComm, myNumRows, rootRank);
465 if (myNumRows != 0) {
469 send (*pComm,
static_cast<int> (myNumRows),
470 myRows.getRawPtr(), rootRank);
480 receive (*pComm, rootRank,
481 static_cast<int> (myNumRows),
482 myNumEntriesPerRow.getRawPtr());
487 std::accumulate (myNumEntriesPerRow.begin(),
488 myNumEntriesPerRow.end(), 0);
494 myColInd = arcp<GO> (myNumEntries);
495 myValues = arcp<scalar_type> (myNumEntries);
496 if (myNumEntries > 0) {
499 receive (*pComm, rootRank,
500 static_cast<int> (myNumEntries),
501 myColInd.getRawPtr());
502 receive (*pComm, rootRank,
503 static_cast<int> (myNumEntries),
504 myValues.getRawPtr());
510 cerr <<
"-- Proc 0: Copying my data from global arrays" << endl;
514 for (size_type k = 0; k < myNumRows; ++k) {
515 const GO myCurRow = myRows[k];
517 myNumEntriesPerRow[k] = numEntriesInThisRow;
519 if (extraDebug && debug) {
520 cerr <<
"Proc " << pRowMap->getComm ()->getRank ()
521 <<
": myNumEntriesPerRow[0.." << (myNumRows-1) <<
"] = [";
522 for (size_type k = 0; k < myNumRows; ++k) {
523 cerr << myNumEntriesPerRow[k];
524 if (k < myNumRows-1) {
532 std::accumulate (myNumEntriesPerRow.begin(),
533 myNumEntriesPerRow.end(), 0);
535 cerr <<
"-- Proc 0: I own " << myNumRows <<
" rows and "
536 << myNumEntries <<
" entries" << endl;
538 myColInd = arcp<GO> (myNumEntries);
539 myValues = arcp<scalar_type> (myNumEntries);
547 for (size_type k = 0; k < myNumRows;
548 myCurPos += myNumEntriesPerRow[k], ++k) {
550 const GO myRow = myRows[k];
551 const size_t curPos = rowPtr[myRow];
554 if (curNumEntries > 0) {
555 ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
556 ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
557 std::copy (colIndView.begin(), colIndView.end(),
558 myColIndView.begin());
560 ArrayView<scalar_type> valuesView =
561 values (curPos, curNumEntries);
562 ArrayView<scalar_type> myValuesView =
563 myValues (myCurPos, curNumEntries);
564 std::copy (valuesView.begin(), valuesView.end(),
565 myValuesView.begin());
570 for (
int p = 1; p < numProcs; ++p) {
572 cerr <<
"-- Proc 0: Processing proc " << p << endl;
575 size_type theirNumRows = 0;
580 receive (*pComm, p, &theirNumRows);
582 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
583 << theirNumRows <<
" rows" << endl;
585 if (theirNumRows != 0) {
590 ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
591 receive (*pComm, p, as<int> (theirNumRows),
592 theirRows.getRawPtr ());
601 const global_size_t numRows = pRowMap->getGlobalNumElements ();
602 const GO indexBase = pRowMap->getIndexBase ();
603 bool theirRowsValid =
true;
604 for (size_type k = 0; k < theirNumRows; ++k) {
605 if (theirRows[k] < indexBase ||
606 as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
607 theirRowsValid =
false;
610 if (! theirRowsValid) {
611 TEUCHOS_TEST_FOR_EXCEPTION(
612 ! theirRowsValid, std::logic_error,
613 "Proc " << p <<
" has at least one invalid row index. "
614 "Here are all of them: " <<
615 Teuchos::toString (theirRows ()) <<
". Valid row index "
616 "range (zero-based): [0, " << (numRows - 1) <<
"].");
631 ArrayRCP<size_t> theirNumEntriesPerRow;
632 theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
633 for (size_type k = 0; k < theirNumRows; ++k) {
634 theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
641 send (*pComm,
static_cast<int> (theirNumRows),
642 theirNumEntriesPerRow.getRawPtr(), p);
646 std::accumulate (theirNumEntriesPerRow.begin(),
647 theirNumEntriesPerRow.end(), 0);
650 cerr <<
"-- Proc 0: Proc " << p <<
" owns "
651 << theirNumEntries <<
" entries" << endl;
656 if (theirNumEntries == 0) {
665 ArrayRCP<GO> theirColInd (theirNumEntries);
666 ArrayRCP<scalar_type> theirValues (theirNumEntries);
674 for (size_type k = 0; k < theirNumRows;
675 theirCurPos += theirNumEntriesPerRow[k], k++) {
677 const GO theirRow = theirRows[k];
683 if (curNumEntries > 0) {
684 ArrayView<GO> colIndView =
685 colInd (curPos, curNumEntries);
686 ArrayView<GO> theirColIndView =
687 theirColInd (theirCurPos, curNumEntries);
688 std::copy (colIndView.begin(), colIndView.end(),
689 theirColIndView.begin());
691 ArrayView<scalar_type> valuesView =
692 values (curPos, curNumEntries);
693 ArrayView<scalar_type> theirValuesView =
694 theirValues (theirCurPos, curNumEntries);
695 std::copy (valuesView.begin(), valuesView.end(),
696 theirValuesView.begin());
703 send (*pComm,
static_cast<int> (theirNumEntries),
704 theirColInd.getRawPtr(), p);
705 send (*pComm,
static_cast<int> (theirNumEntries),
706 theirValues.getRawPtr(), p);
709 cerr <<
"-- Proc 0: Finished with proc " << p << endl;
717 numEntriesPerRow = null;
722 if (debug && myRank == 0) {
723 cerr <<
"-- Proc 0: About to fill in myRowPtr" << endl;
731 myRowPtr = arcp<size_t> (myNumRows+1);
733 for (size_type k = 1; k < myNumRows+1; ++k) {
734 myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
736 if (extraDebug && debug) {
737 cerr <<
"Proc " << Teuchos::rank (*(pRowMap->getComm()))
738 <<
": myRowPtr[0.." << myNumRows <<
"] = [";
739 for (size_type k = 0; k < myNumRows+1; ++k) {
745 cerr <<
"]" << endl << endl;
748 if (debug && myRank == 0) {
749 cerr <<
"-- Proc 0: Done with distribute" << endl;
766 static Teuchos::RCP<sparse_matrix_type>
767 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
768 Teuchos::ArrayRCP<size_t>& myRowPtr,
769 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
770 Teuchos::ArrayRCP<scalar_type>& myValues,
771 const Teuchos::RCP<const map_type>& pRowMap,
772 const Teuchos::RCP<const map_type>& pRangeMap,
773 const Teuchos::RCP<const map_type>& pDomainMap,
774 const bool callFillComplete =
true)
776 using Teuchos::ArrayView;
790 TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
791 "makeMatrix: myRowPtr array is null. "
792 "Please report this bug to the Tpetra developers.");
793 TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
794 "makeMatrix: domain map is null. "
795 "Please report this bug to the Tpetra developers.");
796 TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
797 "makeMatrix: range map is null. "
798 "Please report this bug to the Tpetra developers.");
799 TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
800 "makeMatrix: row map is null. "
801 "Please report this bug to the Tpetra developers.");
805 RCP<sparse_matrix_type> A =
810 ArrayView<const GO> myRows = pRowMap->getLocalElementList ();
811 const size_type myNumRows = myRows.size ();
814 const GO indexBase = pRowMap->getIndexBase ();
815 for (size_type i = 0; i < myNumRows; ++i) {
816 const size_type myCurPos = myRowPtr[i];
818 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
819 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
822 for (size_type k = 0; k < curNumEntries; ++k) {
823 curColInd[k] += indexBase;
826 if (curNumEntries > 0) {
827 A->insertGlobalValues (myRows[i], curColInd, curValues);
834 myNumEntriesPerRow = null;
839 if (callFillComplete) {
840 A->fillComplete (pDomainMap, pRangeMap);
850 static Teuchos::RCP<sparse_matrix_type>
851 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
852 Teuchos::ArrayRCP<size_t>& myRowPtr,
853 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
854 Teuchos::ArrayRCP<scalar_type>& myValues,
855 const Teuchos::RCP<const map_type>& pRowMap,
856 const Teuchos::RCP<const map_type>& pRangeMap,
857 const Teuchos::RCP<const map_type>& pDomainMap,
858 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
859 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams)
861 using Teuchos::ArrayView;
875 TEUCHOS_TEST_FOR_EXCEPTION(
876 myRowPtr.is_null(), std::logic_error,
877 "makeMatrix: myRowPtr array is null. "
878 "Please report this bug to the Tpetra developers.");
879 TEUCHOS_TEST_FOR_EXCEPTION(
880 pDomainMap.is_null(), std::logic_error,
881 "makeMatrix: domain map is null. "
882 "Please report this bug to the Tpetra developers.");
883 TEUCHOS_TEST_FOR_EXCEPTION(
884 pRangeMap.is_null(), std::logic_error,
885 "makeMatrix: range map is null. "
886 "Please report this bug to the Tpetra developers.");
887 TEUCHOS_TEST_FOR_EXCEPTION(
888 pRowMap.is_null(), std::logic_error,
889 "makeMatrix: row map is null. "
890 "Please report this bug to the Tpetra developers.");
894 RCP<sparse_matrix_type> A =
900 ArrayView<const GO> myRows = pRowMap->getLocalElementList();
901 const size_type myNumRows = myRows.size();
904 const GO indexBase = pRowMap->getIndexBase ();
905 for (size_type i = 0; i < myNumRows; ++i) {
906 const size_type myCurPos = myRowPtr[i];
908 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
909 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
912 for (size_type k = 0; k < curNumEntries; ++k) {
913 curColInd[k] += indexBase;
915 if (curNumEntries > 0) {
916 A->insertGlobalValues (myRows[i], curColInd, curValues);
923 myNumEntriesPerRow = null;
928 A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
936 static Teuchos::RCP<sparse_matrix_type>
937 makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
938 Teuchos::ArrayRCP<size_t>& myRowPtr,
939 Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
940 Teuchos::ArrayRCP<scalar_type>& myValues,
941 const Teuchos::RCP<const map_type>& rowMap,
942 Teuchos::RCP<const map_type>& colMap,
943 const Teuchos::RCP<const map_type>& domainMap,
944 const Teuchos::RCP<const map_type>& rangeMap,
945 const bool callFillComplete =
true)
947 using Teuchos::ArrayView;
956 RCP<sparse_matrix_type> A;
957 if (colMap.is_null ()) {
965 ArrayView<const GO> myRows = rowMap->getLocalElementList ();
966 const size_type myNumRows = myRows.size ();
969 const GO indexBase = rowMap->getIndexBase ();
970 for (size_type i = 0; i < myNumRows; ++i) {
971 const size_type myCurPos = myRowPtr[i];
972 const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
973 ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
974 ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
977 for (size_type k = 0; k < curNumEntries; ++k) {
978 curColInd[k] += indexBase;
980 if (curNumEntries > 0) {
981 A->insertGlobalValues (myRows[i], curColInd, curValues);
988 myNumEntriesPerRow = null;
993 if (callFillComplete) {
994 A->fillComplete (domainMap, rangeMap);
995 if (colMap.is_null ()) {
996 colMap = A->getColMap ();
1020 static Teuchos::RCP<const Teuchos::MatrixMarket::Banner>
1021 readBanner (std::istream& in,
1023 const bool tolerant=
false,
1025 const bool isGraph=
false)
1027 using Teuchos::MatrixMarket::Banner;
1032 typedef Teuchos::ScalarTraits<scalar_type> STS;
1034 RCP<Banner> pBanner;
1038 const bool readFailed = ! getline(in, line);
1039 TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1040 "Failed to get Matrix Market banner line from input.");
1047 pBanner = rcp (
new Banner (line, tolerant));
1048 }
catch (std::exception& e) {
1050 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1051 "Matrix Market banner line contains syntax error(s): "
1055 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() !=
"matrix",
1056 std::invalid_argument,
"The Matrix Market file does not contain "
1057 "matrix data. Its Banner (first) line says that its object type is \""
1058 << pBanner->matrixType() <<
"\", rather than the required \"matrix\".");
1062 TEUCHOS_TEST_FOR_EXCEPTION(
1063 ! STS::isComplex && pBanner->dataType() ==
"complex",
1064 std::invalid_argument,
1065 "The Matrix Market file contains complex-valued data, but you are "
1066 "trying to read it into a matrix containing entries of the real-"
1067 "valued Scalar type \""
1068 << Teuchos::TypeNameTraits<scalar_type>::name() <<
"\".");
1069 TEUCHOS_TEST_FOR_EXCEPTION(
1071 pBanner->dataType() !=
"real" &&
1072 pBanner->dataType() !=
"complex" &&
1073 pBanner->dataType() !=
"integer",
1074 std::invalid_argument,
1075 "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
1076 "Matrix Market file may not contain a \"pattern\" matrix. A "
1077 "pattern matrix is really just a graph with no weights. It "
1078 "should be stored in a CrsGraph, not a CrsMatrix.");
1080 TEUCHOS_TEST_FOR_EXCEPTION(
1082 pBanner->dataType() !=
"pattern",
1083 std::invalid_argument,
1084 "When reading Matrix Market data into a Tpetra::CrsGraph, the "
1085 "Matrix Market file must contain a \"pattern\" matrix.");
1112 static Teuchos::Tuple<global_ordinal_type, 3>
1113 readCoordDims (std::istream& in,
1115 const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1116 const Teuchos::RCP<const comm_type>& pComm,
1117 const bool tolerant =
false,
1120 using Teuchos::MatrixMarket::readCoordinateDimensions;
1121 using Teuchos::Tuple;
1126 Tuple<global_ordinal_type, 3> dims;
1132 bool success =
false;
1133 if (pComm->getRank() == 0) {
1134 TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() !=
"coordinate",
1135 std::invalid_argument,
"The Tpetra::CrsMatrix Matrix Market reader "
1136 "only accepts \"coordinate\" (sparse) matrix data.");
1140 success = readCoordinateDimensions (in, numRows, numCols,
1141 numNonzeros, lineNumber,
1147 dims[2] = numNonzeros;
1155 int the_success = success ? 1 : 0;
1156 Teuchos::broadcast (*pComm, 0, &the_success);
1157 success = (the_success == 1);
1162 Teuchos::broadcast (*pComm, 0, dims);
1170 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1171 "Error reading Matrix Market sparse matrix: failed to read "
1172 "coordinate matrix dimensions.");
1187 typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1189 typedef Teuchos::MatrixMarket::SymmetrizingGraphAdder<Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> > graph_adder_type;
1216 static Teuchos::RCP<adder_type>
1217 makeAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1218 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1219 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1220 const bool tolerant=
false,
1221 const bool debug=
false)
1223 if (pComm->getRank () == 0) {
1224 typedef Teuchos::MatrixMarket::Raw::Adder<
scalar_type,
1227 Teuchos::RCP<raw_adder_type> pRaw =
1228 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1230 return Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
1233 return Teuchos::null;
1262 static Teuchos::RCP<graph_adder_type>
1263 makeGraphAdder (
const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1264 Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1265 const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1266 const bool tolerant=
false,
1267 const bool debug=
false)
1269 if (pComm->getRank () == 0) {
1270 typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
1271 Teuchos::RCP<raw_adder_type> pRaw =
1272 Teuchos::rcp (
new raw_adder_type (dims[0], dims[1], dims[2],
1274 return Teuchos::rcp (
new graph_adder_type (pRaw, pBanner->symmType ()));
1277 return Teuchos::null;
1282 static Teuchos::RCP<sparse_graph_type>
1283 readSparseGraphHelper (std::istream& in,
1284 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1285 const Teuchos::RCP<const map_type>& rowMap,
1286 Teuchos::RCP<const map_type>& colMap,
1287 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1288 const bool tolerant,
1291 using Teuchos::MatrixMarket::Banner;
1294 using Teuchos::Tuple;
1298 const int myRank = pComm->getRank ();
1299 const int rootRank = 0;
1304 size_t lineNumber = 1;
1306 if (debug && myRank == rootRank) {
1307 cerr <<
"Matrix Market reader: readGraph:" << endl
1308 <<
"-- Reading banner line" << endl;
1317 RCP<const Banner> pBanner;
1323 int bannerIsCorrect = 1;
1324 std::ostringstream errMsg;
1326 if (myRank == rootRank) {
1329 pBanner = readBanner (in, lineNumber, tolerant, debug,
true);
1331 catch (std::exception& e) {
1332 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
1333 "threw an exception: " << e.what();
1334 bannerIsCorrect = 0;
1337 if (bannerIsCorrect) {
1342 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
1343 bannerIsCorrect = 0;
1344 errMsg <<
"The Matrix Market input file must contain a "
1345 "\"coordinate\"-format sparse graph in order to create a "
1346 "Tpetra::CrsGraph object from it, but the file's matrix "
1347 "type is \"" << pBanner->matrixType() <<
"\" instead.";
1352 if (tolerant && pBanner->matrixType() ==
"array") {
1353 bannerIsCorrect = 0;
1354 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
1355 "format sparse graph in order to create a Tpetra::CrsGraph "
1356 "object from it, but the file's matrix type is \"array\" "
1357 "instead. That probably means the file contains dense matrix "
1364 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1371 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1372 std::invalid_argument, errMsg.str ());
1374 if (debug && myRank == rootRank) {
1375 cerr <<
"-- Reading dimensions line" << endl;
1383 Tuple<global_ordinal_type, 3> dims =
1384 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1386 if (debug && myRank == rootRank) {
1387 cerr <<
"-- Making Adder for collecting graph data" << endl;
1394 RCP<graph_adder_type> pAdder =
1395 makeGraphAdder (pComm, pBanner, dims, tolerant, debug);
1397 if (debug && myRank == rootRank) {
1398 cerr <<
"-- Reading graph data" << endl;
1408 int readSuccess = 1;
1409 std::ostringstream errMsg;
1410 if (myRank == rootRank) {
1413 typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
1415 reader_type reader (pAdder);
1418 std::pair<bool, std::vector<size_t> > results =
1419 reader.read (in, lineNumber, tolerant, debug);
1420 readSuccess = results.first ? 1 : 0;
1422 catch (std::exception& e) {
1427 broadcast (*pComm, rootRank, ptr (&readSuccess));
1436 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1437 "Failed to read the Matrix Market sparse graph file: "
1441 if (debug && myRank == rootRank) {
1442 cerr <<
"-- Successfully read the Matrix Market data" << endl;
1455 if (debug && myRank == rootRank) {
1456 cerr <<
"-- Tolerant mode: rebroadcasting graph dimensions"
1458 <<
"----- Dimensions before: "
1459 << dims[0] <<
" x " << dims[1]
1463 Tuple<global_ordinal_type, 2> updatedDims;
1464 if (myRank == rootRank) {
1471 std::max (dims[0], pAdder->getAdder()->numRows());
1472 updatedDims[1] = pAdder->getAdder()->numCols();
1474 broadcast (*pComm, rootRank, updatedDims);
1475 dims[0] = updatedDims[0];
1476 dims[1] = updatedDims[1];
1477 if (debug && myRank == rootRank) {
1478 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
1492 if (myRank == rootRank) {
1499 if (dims[0] < pAdder->getAdder ()->numRows ()) {
1503 broadcast (*pComm, 0, ptr (&dimsMatch));
1504 if (dimsMatch == 0) {
1511 Tuple<global_ordinal_type, 2> addersDims;
1512 if (myRank == rootRank) {
1513 addersDims[0] = pAdder->getAdder()->numRows();
1514 addersDims[1] = pAdder->getAdder()->numCols();
1516 broadcast (*pComm, 0, addersDims);
1517 TEUCHOS_TEST_FOR_EXCEPTION(
1518 dimsMatch == 0, std::runtime_error,
1519 "The graph metadata says that the graph is " << dims[0] <<
" x "
1520 << dims[1] <<
", but the actual data says that the graph is "
1521 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
1522 "data includes more rows than reported in the metadata. This "
1523 "is not allowed when parsing in strict mode. Parse the graph in "
1524 "tolerant mode to ignore the metadata when it disagrees with the "
1530 RCP<map_type> proc0Map;
1532 if(Teuchos::is_null(rowMap)) {
1536 indexBase = rowMap->getIndexBase();
1538 if(myRank == rootRank) {
1539 proc0Map = rcp(
new map_type(dims[0],dims[0],indexBase,pComm));
1542 proc0Map = rcp(
new map_type(dims[0],0,indexBase,pComm));
1546 std::map<global_ordinal_type, size_t> numEntriesPerRow_map;
1547 if (myRank == rootRank) {
1548 const auto& entries = pAdder()->getAdder()->getEntries();
1553 for (
const auto& entry : entries) {
1555 ++numEntriesPerRow_map[gblRow];
1559 Teuchos::Array<size_t> numEntriesPerRow (proc0Map->getLocalNumElements ());
1560 for (
const auto& ent : numEntriesPerRow_map) {
1562 numEntriesPerRow[lclRow] = ent.second;
1569 std::map<global_ordinal_type, size_t> empty_map;
1570 std::swap (numEntriesPerRow_map, empty_map);
1573 RCP<sparse_graph_type> proc0Graph =
1575 constructorParams));
1576 if(myRank == rootRank) {
1577 typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
1580 const std::vector<element_type>& entries =
1581 pAdder->getAdder()->getEntries();
1584 for(
size_t curPos=0; curPos<entries.size(); curPos++) {
1585 const element_type& curEntry = entries[curPos];
1588 Teuchos::ArrayView<const global_ordinal_type> colView(&curCol,1);
1589 proc0Graph->insertGlobalIndices(curRow,colView);
1592 proc0Graph->fillComplete();
1594 RCP<sparse_graph_type> distGraph;
1595 if(Teuchos::is_null(rowMap))
1598 RCP<map_type> distMap =
1599 rcp(
new map_type(dims[0],0,pComm,GloballyDistributed));
1605 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1606 import_type importer (proc0Map, distMap);
1609 distGraph->doImport(*proc0Graph,importer,
INSERT);
1615 typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1616 import_type importer (proc0Map, rowMap);
1619 distGraph->doImport(*proc0Graph,importer,
INSERT);
1649 static Teuchos::RCP<sparse_graph_type>
1651 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1652 const bool callFillComplete=
true,
1653 const bool tolerant=
false,
1654 const bool debug=
false)
1656 using Teuchos::broadcast;
1657 using Teuchos::outArg;
1665 if (comm->getRank () == 0) {
1667 in.open (filename.c_str ());
1668 opened = in.is_open();
1674 broadcast<int, int> (*comm, 0, outArg (opened));
1675 TEUCHOS_TEST_FOR_EXCEPTION
1676 (opened == 0, std::runtime_error,
"readSparseGraphFile: "
1677 "Failed to open file \"" << filename <<
"\" on Process 0.");
1715 static Teuchos::RCP<sparse_graph_type>
1717 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1718 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1719 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1720 const bool tolerant=
false,
1721 const bool debug=
false)
1723 using Teuchos::broadcast;
1724 using Teuchos::outArg;
1732 if (pComm->getRank () == 0) {
1734 in.open (filename.c_str ());
1735 opened = in.is_open();
1741 broadcast<int, int> (*pComm, 0, outArg (opened));
1742 TEUCHOS_TEST_FOR_EXCEPTION
1743 (opened == 0, std::runtime_error,
"readSparseGraphFile: "
1744 "Failed to open file \"" << filename <<
"\" on Process 0.");
1745 if (pComm->getRank () == 0) {
1746 in.open (filename.c_str ());
1750 fillCompleteParams, tolerant, debug);
1794 static Teuchos::RCP<sparse_graph_type>
1796 const Teuchos::RCP<const map_type>& rowMap,
1797 Teuchos::RCP<const map_type>& colMap,
1798 const Teuchos::RCP<const map_type>& domainMap,
1799 const Teuchos::RCP<const map_type>& rangeMap,
1800 const bool callFillComplete=
true,
1801 const bool tolerant=
false,
1802 const bool debug=
false)
1804 using Teuchos::broadcast;
1805 using Teuchos::Comm;
1806 using Teuchos::outArg;
1809 TEUCHOS_TEST_FOR_EXCEPTION
1810 (rowMap.is_null (), std::invalid_argument,
1811 "Input rowMap must be nonnull.");
1812 RCP<const Comm<int> > comm = rowMap->getComm ();
1813 if (comm.is_null ()) {
1816 return Teuchos::null;
1825 if (comm->getRank () == 0) {
1827 in.open (filename.c_str ());
1828 opened = in.is_open();
1834 broadcast<int, int> (*comm, 0, outArg (opened));
1835 TEUCHOS_TEST_FOR_EXCEPTION
1836 (opened == 0, std::runtime_error,
"readSparseGraphFile: "
1837 "Failed to open file \"" << filename <<
"\" on Process 0.");
1839 callFillComplete, tolerant, debug);
1867 static Teuchos::RCP<sparse_graph_type>
1869 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1870 const bool callFillComplete=
true,
1871 const bool tolerant=
false,
1872 const bool debug=
false)
1874 Teuchos::RCP<const map_type> fakeRowMap;
1875 Teuchos::RCP<const map_type> fakeColMap;
1876 Teuchos::RCP<Teuchos::ParameterList> fakeCtorParams;
1878 Teuchos::RCP<sparse_graph_type> graph =
1879 readSparseGraphHelper (in, pComm,
1880 fakeRowMap, fakeColMap,
1881 fakeCtorParams, tolerant, debug);
1882 if (callFillComplete) {
1883 graph->fillComplete ();
1918 static Teuchos::RCP<sparse_graph_type>
1920 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
1921 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1922 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1923 const bool tolerant=
false,
1924 const bool debug=
false)
1926 Teuchos::RCP<const map_type> fakeRowMap;
1927 Teuchos::RCP<const map_type> fakeColMap;
1928 Teuchos::RCP<sparse_graph_type> graph =
1929 readSparseGraphHelper (in, pComm,
1930 fakeRowMap, fakeColMap,
1931 constructorParams, tolerant, debug);
1932 graph->fillComplete (fillCompleteParams);
1977 static Teuchos::RCP<sparse_graph_type>
1979 const Teuchos::RCP<const map_type>& rowMap,
1980 Teuchos::RCP<const map_type>& colMap,
1981 const Teuchos::RCP<const map_type>& domainMap,
1982 const Teuchos::RCP<const map_type>& rangeMap,
1983 const bool callFillComplete=
true,
1984 const bool tolerant=
false,
1985 const bool debug=
false)
1987 Teuchos::RCP<sparse_graph_type> graph =
1988 readSparseGraphHelper (in, rowMap->getComm (),
1989 rowMap, colMap, Teuchos::null, tolerant,
1991 if (callFillComplete) {
1992 graph->fillComplete (domainMap, rangeMap);
1997#include "MatrixMarket_TpetraNew.hpp"
2022 static Teuchos::RCP<sparse_matrix_type>
2024 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2025 const bool callFillComplete=
true,
2026 const bool tolerant=
false,
2027 const bool debug=
false)
2029 const int myRank = pComm->getRank ();
2034 in.open (filename.c_str ());
2042 return readSparse (in, pComm, callFillComplete, tolerant, debug);
2077 static Teuchos::RCP<sparse_matrix_type>
2079 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2080 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2081 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2082 const bool tolerant=
false,
2083 const bool debug=
false)
2086 if (pComm->getRank () == 0) {
2087 in.open (filename.c_str ());
2089 return readSparse (in, pComm, constructorParams,
2090 fillCompleteParams, tolerant, debug);
2131 static Teuchos::RCP<sparse_matrix_type>
2133 const Teuchos::RCP<const map_type>& rowMap,
2134 Teuchos::RCP<const map_type>& colMap,
2135 const Teuchos::RCP<const map_type>& domainMap,
2136 const Teuchos::RCP<const map_type>& rangeMap,
2137 const bool callFillComplete=
true,
2138 const bool tolerant=
false,
2139 const bool debug=
false)
2141 using Teuchos::broadcast;
2142 using Teuchos::Comm;
2143 using Teuchos::outArg;
2146 TEUCHOS_TEST_FOR_EXCEPTION(
2147 rowMap.is_null (), std::invalid_argument,
2148 "Row Map must be nonnull.");
2150 RCP<const Comm<int> > comm = rowMap->getComm ();
2151 const int myRank = comm->getRank ();
2161 in.open (filename.c_str ());
2162 opened = in.is_open();
2168 broadcast<int, int> (*comm, 0, outArg (opened));
2169 TEUCHOS_TEST_FOR_EXCEPTION(
2170 opened == 0, std::runtime_error,
2171 "readSparseFile: Failed to open file \"" << filename <<
"\" on "
2173 return readSparse (in, rowMap, colMap, domainMap, rangeMap,
2174 callFillComplete, tolerant, debug);
2202 static Teuchos::RCP<sparse_matrix_type>
2204 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2205 const bool callFillComplete=
true,
2206 const bool tolerant=
false,
2207 const bool debug=
false)
2209 using Teuchos::MatrixMarket::Banner;
2210 using Teuchos::arcp;
2211 using Teuchos::ArrayRCP;
2212 using Teuchos::broadcast;
2213 using Teuchos::null;
2216 using Teuchos::REDUCE_MAX;
2217 using Teuchos::reduceAll;
2218 using Teuchos::Tuple;
2221 typedef Teuchos::ScalarTraits<scalar_type> STS;
2223 const bool extraDebug =
false;
2224 const int myRank = pComm->getRank ();
2225 const int rootRank = 0;
2230 size_t lineNumber = 1;
2232 if (debug && myRank == rootRank) {
2233 cerr <<
"Matrix Market reader: readSparse:" << endl
2234 <<
"-- Reading banner line" << endl;
2243 RCP<const Banner> pBanner;
2249 int bannerIsCorrect = 1;
2250 std::ostringstream errMsg;
2252 if (myRank == rootRank) {
2255 pBanner = readBanner (in, lineNumber, tolerant, debug);
2257 catch (std::exception& e) {
2258 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2259 "threw an exception: " << e.what();
2260 bannerIsCorrect = 0;
2263 if (bannerIsCorrect) {
2268 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2269 bannerIsCorrect = 0;
2270 errMsg <<
"The Matrix Market input file must contain a "
2271 "\"coordinate\"-format sparse matrix in order to create a "
2272 "Tpetra::CrsMatrix object from it, but the file's matrix "
2273 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2278 if (tolerant && pBanner->matrixType() ==
"array") {
2279 bannerIsCorrect = 0;
2280 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2281 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2282 "object from it, but the file's matrix type is \"array\" "
2283 "instead. That probably means the file contains dense matrix "
2290 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2297 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2298 std::invalid_argument, errMsg.str ());
2300 if (debug && myRank == rootRank) {
2301 cerr <<
"-- Reading dimensions line" << endl;
2309 Tuple<global_ordinal_type, 3> dims =
2310 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2312 if (debug && myRank == rootRank) {
2313 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2318 RCP<adder_type> pAdder =
2319 makeAdder (pComm, pBanner, dims, tolerant, debug);
2321 if (debug && myRank == rootRank) {
2322 cerr <<
"-- Reading matrix data" << endl;
2332 int readSuccess = 1;
2333 std::ostringstream errMsg;
2334 if (myRank == rootRank) {
2337 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2339 reader_type reader (pAdder);
2342 std::pair<bool, std::vector<size_t> > results =
2343 reader.read (in, lineNumber, tolerant, debug);
2344 readSuccess = results.first ? 1 : 0;
2346 catch (std::exception& e) {
2351 broadcast (*pComm, rootRank, ptr (&readSuccess));
2360 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2361 "Failed to read the Matrix Market sparse matrix file: "
2365 if (debug && myRank == rootRank) {
2366 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2379 if (debug && myRank == rootRank) {
2380 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2382 <<
"----- Dimensions before: "
2383 << dims[0] <<
" x " << dims[1]
2387 Tuple<global_ordinal_type, 2> updatedDims;
2388 if (myRank == rootRank) {
2395 std::max (dims[0], pAdder->getAdder()->numRows());
2396 updatedDims[1] = pAdder->getAdder()->numCols();
2398 broadcast (*pComm, rootRank, updatedDims);
2399 dims[0] = updatedDims[0];
2400 dims[1] = updatedDims[1];
2401 if (debug && myRank == rootRank) {
2402 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
2416 if (myRank == rootRank) {
2423 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2427 broadcast (*pComm, 0, ptr (&dimsMatch));
2428 if (dimsMatch == 0) {
2435 Tuple<global_ordinal_type, 2> addersDims;
2436 if (myRank == rootRank) {
2437 addersDims[0] = pAdder->getAdder()->numRows();
2438 addersDims[1] = pAdder->getAdder()->numCols();
2440 broadcast (*pComm, 0, addersDims);
2441 TEUCHOS_TEST_FOR_EXCEPTION(
2442 dimsMatch == 0, std::runtime_error,
2443 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
2444 << dims[1] <<
", but the actual data says that the matrix is "
2445 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
2446 "data includes more rows than reported in the metadata. This "
2447 "is not allowed when parsing in strict mode. Parse the matrix in "
2448 "tolerant mode to ignore the metadata when it disagrees with the "
2453 if (debug && myRank == rootRank) {
2454 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
2466 ArrayRCP<size_t> numEntriesPerRow;
2467 ArrayRCP<size_t> rowPtr;
2468 ArrayRCP<global_ordinal_type> colInd;
2469 ArrayRCP<scalar_type> values;
2474 int mergeAndConvertSucceeded = 1;
2475 std::ostringstream errMsg;
2477 if (myRank == rootRank) {
2479 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
2489 const size_type numRows = dims[0];
2492 pAdder->getAdder()->merge ();
2495 const std::vector<element_type>& entries =
2496 pAdder->getAdder()->getEntries();
2499 const size_t numEntries = (size_t)entries.size();
2502 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
2503 <<
" rows and numEntries=" << numEntries
2504 <<
" entries." << endl;
2510 numEntriesPerRow = arcp<size_t> (numRows);
2511 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2512 rowPtr = arcp<size_t> (numRows+1);
2513 std::fill (rowPtr.begin(), rowPtr.end(), 0);
2514 colInd = arcp<global_ordinal_type> (numEntries);
2515 values = arcp<scalar_type> (numEntries);
2522 for (curPos = 0; curPos < numEntries; ++curPos) {
2523 const element_type& curEntry = entries[curPos];
2525 TEUCHOS_TEST_FOR_EXCEPTION(
2526 curRow < prvRow, std::logic_error,
2527 "Row indices are out of order, even though they are supposed "
2528 "to be sorted. curRow = " << curRow <<
", prvRow = "
2529 << prvRow <<
", at curPos = " << curPos <<
". Please report "
2530 "this bug to the Tpetra developers.");
2531 if (curRow > prvRow) {
2537 numEntriesPerRow[curRow]++;
2538 colInd[curPos] = curEntry.colIndex();
2539 values[curPos] = curEntry.value();
2544 rowPtr[numRows] = numEntries;
2546 catch (std::exception& e) {
2547 mergeAndConvertSucceeded = 0;
2548 errMsg <<
"Failed to merge sparse matrix entries and convert to "
2549 "CSR format: " << e.what();
2552 if (debug && mergeAndConvertSucceeded) {
2554 const size_type numRows = dims[0];
2555 const size_type maxToDisplay = 100;
2557 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
2558 << (numEntriesPerRow.size()-1) <<
"] ";
2559 if (numRows > maxToDisplay) {
2560 cerr <<
"(only showing first and last few entries) ";
2564 if (numRows > maxToDisplay) {
2565 for (size_type k = 0; k < 2; ++k) {
2566 cerr << numEntriesPerRow[k] <<
" ";
2569 for (size_type k = numRows-2; k < numRows-1; ++k) {
2570 cerr << numEntriesPerRow[k] <<
" ";
2574 for (size_type k = 0; k < numRows-1; ++k) {
2575 cerr << numEntriesPerRow[k] <<
" ";
2578 cerr << numEntriesPerRow[numRows-1];
2580 cerr <<
"]" << endl;
2582 cerr <<
"----- Proc 0: rowPtr ";
2583 if (numRows > maxToDisplay) {
2584 cerr <<
"(only showing first and last few entries) ";
2587 if (numRows > maxToDisplay) {
2588 for (size_type k = 0; k < 2; ++k) {
2589 cerr << rowPtr[k] <<
" ";
2592 for (size_type k = numRows-2; k < numRows; ++k) {
2593 cerr << rowPtr[k] <<
" ";
2597 for (size_type k = 0; k < numRows; ++k) {
2598 cerr << rowPtr[k] <<
" ";
2601 cerr << rowPtr[numRows] <<
"]" << endl;
2612 if (debug && myRank == rootRank) {
2613 cerr <<
"-- Making range, domain, and row maps" << endl;
2620 RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
2621 RCP<const map_type> pDomainMap =
2622 makeDomainMap (pRangeMap, dims[0], dims[1]);
2623 RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
2625 if (debug && myRank == rootRank) {
2626 cerr <<
"-- Distributing the matrix data" << endl;
2640 ArrayRCP<size_t> myNumEntriesPerRow;
2641 ArrayRCP<size_t> myRowPtr;
2642 ArrayRCP<global_ordinal_type> myColInd;
2643 ArrayRCP<scalar_type> myValues;
2645 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
2646 numEntriesPerRow, rowPtr, colInd, values, debug);
2648 if (debug && myRank == rootRank) {
2649 cerr <<
"-- Inserting matrix entries on each processor";
2650 if (callFillComplete) {
2651 cerr <<
" and calling fillComplete()";
2662 RCP<sparse_matrix_type> pMatrix =
2663 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
2664 pRowMap, pRangeMap, pDomainMap, callFillComplete);
2669 int localIsNull = pMatrix.is_null () ? 1 : 0;
2670 int globalIsNull = 0;
2671 reduceAll (*pComm, REDUCE_MAX, localIsNull, ptr (&globalIsNull));
2672 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
2673 "Reader::makeMatrix() returned a null pointer on at least one "
2674 "process. Please report this bug to the Tpetra developers.");
2677 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
2678 "Reader::makeMatrix() returned a null pointer. "
2679 "Please report this bug to the Tpetra developers.");
2691 if (callFillComplete) {
2692 const int numProcs = pComm->getSize ();
2694 if (extraDebug && debug) {
2696 pRangeMap->getGlobalNumElements ();
2698 pDomainMap->getGlobalNumElements ();
2699 if (myRank == rootRank) {
2700 cerr <<
"-- Matrix is "
2701 << globalNumRows <<
" x " << globalNumCols
2702 <<
" with " << pMatrix->getGlobalNumEntries()
2703 <<
" entries, and index base "
2704 << pMatrix->getIndexBase() <<
"." << endl;
2707 for (
int p = 0; p < numProcs; ++p) {
2709 cerr <<
"-- Proc " << p <<
" owns "
2710 << pMatrix->getLocalNumCols() <<
" columns, and "
2711 << pMatrix->getLocalNumEntries() <<
" entries." << endl;
2718 if (debug && myRank == rootRank) {
2719 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
2754 static Teuchos::RCP<sparse_matrix_type>
2756 const Teuchos::RCP<
const Teuchos::Comm<int> >& pComm,
2757 const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2758 const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2759 const bool tolerant=
false,
2760 const bool debug=
false)
2762 using Teuchos::MatrixMarket::Banner;
2763 using Teuchos::arcp;
2764 using Teuchos::ArrayRCP;
2765 using Teuchos::broadcast;
2766 using Teuchos::null;
2769 using Teuchos::reduceAll;
2770 using Teuchos::Tuple;
2773 typedef Teuchos::ScalarTraits<scalar_type> STS;
2775 const bool extraDebug =
false;
2776 const int myRank = pComm->getRank ();
2777 const int rootRank = 0;
2782 size_t lineNumber = 1;
2784 if (debug && myRank == rootRank) {
2785 cerr <<
"Matrix Market reader: readSparse:" << endl
2786 <<
"-- Reading banner line" << endl;
2795 RCP<const Banner> pBanner;
2801 int bannerIsCorrect = 1;
2802 std::ostringstream errMsg;
2804 if (myRank == rootRank) {
2807 pBanner = readBanner (in, lineNumber, tolerant, debug);
2809 catch (std::exception& e) {
2810 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
2811 "threw an exception: " << e.what();
2812 bannerIsCorrect = 0;
2815 if (bannerIsCorrect) {
2820 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
2821 bannerIsCorrect = 0;
2822 errMsg <<
"The Matrix Market input file must contain a "
2823 "\"coordinate\"-format sparse matrix in order to create a "
2824 "Tpetra::CrsMatrix object from it, but the file's matrix "
2825 "type is \"" << pBanner->matrixType() <<
"\" instead.";
2830 if (tolerant && pBanner->matrixType() ==
"array") {
2831 bannerIsCorrect = 0;
2832 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
2833 "format sparse matrix in order to create a Tpetra::CrsMatrix "
2834 "object from it, but the file's matrix type is \"array\" "
2835 "instead. That probably means the file contains dense matrix "
2842 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2849 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2850 std::invalid_argument, errMsg.str ());
2852 if (debug && myRank == rootRank) {
2853 cerr <<
"-- Reading dimensions line" << endl;
2861 Tuple<global_ordinal_type, 3> dims =
2862 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2864 if (debug && myRank == rootRank) {
2865 cerr <<
"-- Making Adder for collecting matrix data" << endl;
2870 RCP<adder_type> pAdder =
2871 makeAdder (pComm, pBanner, dims, tolerant, debug);
2873 if (debug && myRank == rootRank) {
2874 cerr <<
"-- Reading matrix data" << endl;
2884 int readSuccess = 1;
2885 std::ostringstream errMsg;
2886 if (myRank == rootRank) {
2889 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2891 reader_type reader (pAdder);
2894 std::pair<bool, std::vector<size_t> > results =
2895 reader.read (in, lineNumber, tolerant, debug);
2896 readSuccess = results.first ? 1 : 0;
2898 catch (std::exception& e) {
2903 broadcast (*pComm, rootRank, ptr (&readSuccess));
2912 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2913 "Failed to read the Matrix Market sparse matrix file: "
2917 if (debug && myRank == rootRank) {
2918 cerr <<
"-- Successfully read the Matrix Market data" << endl;
2931 if (debug && myRank == rootRank) {
2932 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
2934 <<
"----- Dimensions before: "
2935 << dims[0] <<
" x " << dims[1]
2939 Tuple<global_ordinal_type, 2> updatedDims;
2940 if (myRank == rootRank) {
2947 std::max (dims[0], pAdder->getAdder()->numRows());
2948 updatedDims[1] = pAdder->getAdder()->numCols();
2950 broadcast (*pComm, rootRank, updatedDims);
2951 dims[0] = updatedDims[0];
2952 dims[1] = updatedDims[1];
2953 if (debug && myRank == rootRank) {
2954 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
2968 if (myRank == rootRank) {
2975 if (dims[0] < pAdder->getAdder ()->numRows ()) {
2979 broadcast (*pComm, 0, ptr (&dimsMatch));
2980 if (dimsMatch == 0) {
2987 Tuple<global_ordinal_type, 2> addersDims;
2988 if (myRank == rootRank) {
2989 addersDims[0] = pAdder->getAdder()->numRows();
2990 addersDims[1] = pAdder->getAdder()->numCols();
2992 broadcast (*pComm, 0, addersDims);
2993 TEUCHOS_TEST_FOR_EXCEPTION(
2994 dimsMatch == 0, std::runtime_error,
2995 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
2996 << dims[1] <<
", but the actual data says that the matrix is "
2997 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
2998 "data includes more rows than reported in the metadata. This "
2999 "is not allowed when parsing in strict mode. Parse the matrix in "
3000 "tolerant mode to ignore the metadata when it disagrees with the "
3005 if (debug && myRank == rootRank) {
3006 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
3018 ArrayRCP<size_t> numEntriesPerRow;
3019 ArrayRCP<size_t> rowPtr;
3020 ArrayRCP<global_ordinal_type> colInd;
3021 ArrayRCP<scalar_type> values;
3026 int mergeAndConvertSucceeded = 1;
3027 std::ostringstream errMsg;
3029 if (myRank == rootRank) {
3031 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3041 const size_type numRows = dims[0];
3044 pAdder->getAdder()->merge ();
3047 const std::vector<element_type>& entries =
3048 pAdder->getAdder()->getEntries();
3051 const size_t numEntries = (size_t)entries.size();
3054 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
3055 <<
" rows and numEntries=" << numEntries
3056 <<
" entries." << endl;
3062 numEntriesPerRow = arcp<size_t> (numRows);
3063 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3064 rowPtr = arcp<size_t> (numRows+1);
3065 std::fill (rowPtr.begin(), rowPtr.end(), 0);
3066 colInd = arcp<global_ordinal_type> (numEntries);
3067 values = arcp<scalar_type> (numEntries);
3074 for (curPos = 0; curPos < numEntries; ++curPos) {
3075 const element_type& curEntry = entries[curPos];
3077 TEUCHOS_TEST_FOR_EXCEPTION(
3078 curRow < prvRow, std::logic_error,
3079 "Row indices are out of order, even though they are supposed "
3080 "to be sorted. curRow = " << curRow <<
", prvRow = "
3081 << prvRow <<
", at curPos = " << curPos <<
". Please report "
3082 "this bug to the Tpetra developers.");
3083 if (curRow > prvRow) {
3089 numEntriesPerRow[curRow]++;
3090 colInd[curPos] = curEntry.colIndex();
3091 values[curPos] = curEntry.value();
3096 rowPtr[numRows] = numEntries;
3098 catch (std::exception& e) {
3099 mergeAndConvertSucceeded = 0;
3100 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3101 "CSR format: " << e.what();
3104 if (debug && mergeAndConvertSucceeded) {
3106 const size_type numRows = dims[0];
3107 const size_type maxToDisplay = 100;
3109 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3110 << (numEntriesPerRow.size()-1) <<
"] ";
3111 if (numRows > maxToDisplay) {
3112 cerr <<
"(only showing first and last few entries) ";
3116 if (numRows > maxToDisplay) {
3117 for (size_type k = 0; k < 2; ++k) {
3118 cerr << numEntriesPerRow[k] <<
" ";
3121 for (size_type k = numRows-2; k < numRows-1; ++k) {
3122 cerr << numEntriesPerRow[k] <<
" ";
3126 for (size_type k = 0; k < numRows-1; ++k) {
3127 cerr << numEntriesPerRow[k] <<
" ";
3130 cerr << numEntriesPerRow[numRows-1];
3132 cerr <<
"]" << endl;
3134 cerr <<
"----- Proc 0: rowPtr ";
3135 if (numRows > maxToDisplay) {
3136 cerr <<
"(only showing first and last few entries) ";
3139 if (numRows > maxToDisplay) {
3140 for (size_type k = 0; k < 2; ++k) {
3141 cerr << rowPtr[k] <<
" ";
3144 for (size_type k = numRows-2; k < numRows; ++k) {
3145 cerr << rowPtr[k] <<
" ";
3149 for (size_type k = 0; k < numRows; ++k) {
3150 cerr << rowPtr[k] <<
" ";
3153 cerr << rowPtr[numRows] <<
"]" << endl;
3164 if (debug && myRank == rootRank) {
3165 cerr <<
"-- Making range, domain, and row maps" << endl;
3172 RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
3173 RCP<const map_type> pDomainMap =
3174 makeDomainMap (pRangeMap, dims[0], dims[1]);
3175 RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
3177 if (debug && myRank == rootRank) {
3178 cerr <<
"-- Distributing the matrix data" << endl;
3192 ArrayRCP<size_t> myNumEntriesPerRow;
3193 ArrayRCP<size_t> myRowPtr;
3194 ArrayRCP<global_ordinal_type> myColInd;
3195 ArrayRCP<scalar_type> myValues;
3197 distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
3198 numEntriesPerRow, rowPtr, colInd, values, debug);
3200 if (debug && myRank == rootRank) {
3201 cerr <<
"-- Inserting matrix entries on each process "
3202 "and calling fillComplete()" << endl;
3211 Teuchos::RCP<sparse_matrix_type> pMatrix =
3212 makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
3213 pRowMap, pRangeMap, pDomainMap, constructorParams,
3214 fillCompleteParams);
3219 int localIsNull = pMatrix.is_null () ? 1 : 0;
3220 int globalIsNull = 0;
3221 reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
3222 TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
3223 "Reader::makeMatrix() returned a null pointer on at least one "
3224 "process. Please report this bug to the Tpetra developers.");
3227 TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
3228 "Reader::makeMatrix() returned a null pointer. "
3229 "Please report this bug to the Tpetra developers.");
3238 if (extraDebug && debug) {
3239 const int numProcs = pComm->getSize ();
3241 pRangeMap->getGlobalNumElements();
3243 pDomainMap->getGlobalNumElements();
3244 if (myRank == rootRank) {
3245 cerr <<
"-- Matrix is "
3246 << globalNumRows <<
" x " << globalNumCols
3247 <<
" with " << pMatrix->getGlobalNumEntries()
3248 <<
" entries, and index base "
3249 << pMatrix->getIndexBase() <<
"." << endl;
3252 for (
int p = 0; p < numProcs; ++p) {
3254 cerr <<
"-- Proc " << p <<
" owns "
3255 << pMatrix->getLocalNumCols() <<
" columns, and "
3256 << pMatrix->getLocalNumEntries() <<
" entries." << endl;
3262 if (debug && myRank == rootRank) {
3263 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3310 static Teuchos::RCP<sparse_matrix_type>
3312 const Teuchos::RCP<const map_type>& rowMap,
3313 Teuchos::RCP<const map_type>& colMap,
3314 const Teuchos::RCP<const map_type>& domainMap,
3315 const Teuchos::RCP<const map_type>& rangeMap,
3316 const bool callFillComplete=
true,
3317 const bool tolerant=
false,
3318 const bool debug=
false)
3320 using Teuchos::MatrixMarket::Banner;
3321 using Teuchos::arcp;
3322 using Teuchos::ArrayRCP;
3323 using Teuchos::ArrayView;
3325 using Teuchos::broadcast;
3326 using Teuchos::Comm;
3327 using Teuchos::null;
3330 using Teuchos::reduceAll;
3331 using Teuchos::Tuple;
3334 typedef Teuchos::ScalarTraits<scalar_type> STS;
3336 RCP<const Comm<int> > pComm = rowMap->getComm ();
3337 const int myRank = pComm->getRank ();
3338 const int rootRank = 0;
3339 const bool extraDebug =
false;
3344 TEUCHOS_TEST_FOR_EXCEPTION(
3345 rowMap.is_null (), std::invalid_argument,
3346 "Row Map must be nonnull.");
3347 TEUCHOS_TEST_FOR_EXCEPTION(
3348 rangeMap.is_null (), std::invalid_argument,
3349 "Range Map must be nonnull.");
3350 TEUCHOS_TEST_FOR_EXCEPTION(
3351 domainMap.is_null (), std::invalid_argument,
3352 "Domain Map must be nonnull.");
3353 TEUCHOS_TEST_FOR_EXCEPTION(
3354 rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
3355 std::invalid_argument,
3356 "The specified row Map's communicator (rowMap->getComm())"
3357 "differs from the given separately supplied communicator pComm.");
3358 TEUCHOS_TEST_FOR_EXCEPTION(
3359 domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
3360 std::invalid_argument,
3361 "The specified domain Map's communicator (domainMap->getComm())"
3362 "differs from the given separately supplied communicator pComm.");
3363 TEUCHOS_TEST_FOR_EXCEPTION(
3364 rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
3365 std::invalid_argument,
3366 "The specified range Map's communicator (rangeMap->getComm())"
3367 "differs from the given separately supplied communicator pComm.");
3372 size_t lineNumber = 1;
3374 if (debug && myRank == rootRank) {
3375 cerr <<
"Matrix Market reader: readSparse:" << endl
3376 <<
"-- Reading banner line" << endl;
3385 RCP<const Banner> pBanner;
3391 int bannerIsCorrect = 1;
3392 std::ostringstream errMsg;
3394 if (myRank == rootRank) {
3397 pBanner = readBanner (in, lineNumber, tolerant, debug);
3399 catch (std::exception& e) {
3400 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
3401 "threw an exception: " << e.what();
3402 bannerIsCorrect = 0;
3405 if (bannerIsCorrect) {
3410 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
3411 bannerIsCorrect = 0;
3412 errMsg <<
"The Matrix Market input file must contain a "
3413 "\"coordinate\"-format sparse matrix in order to create a "
3414 "Tpetra::CrsMatrix object from it, but the file's matrix "
3415 "type is \"" << pBanner->matrixType() <<
"\" instead.";
3420 if (tolerant && pBanner->matrixType() ==
"array") {
3421 bannerIsCorrect = 0;
3422 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
3423 "format sparse matrix in order to create a Tpetra::CrsMatrix "
3424 "object from it, but the file's matrix type is \"array\" "
3425 "instead. That probably means the file contains dense matrix "
3432 broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
3439 TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
3440 std::invalid_argument, errMsg.str ());
3442 if (debug && myRank == rootRank) {
3443 cerr <<
"-- Reading dimensions line" << endl;
3451 Tuple<global_ordinal_type, 3> dims =
3452 readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
3454 if (debug && myRank == rootRank) {
3455 cerr <<
"-- Making Adder for collecting matrix data" << endl;
3462 RCP<adder_type> pAdder =
3463 makeAdder (pComm, pBanner, dims, tolerant, debug);
3465 if (debug && myRank == rootRank) {
3466 cerr <<
"-- Reading matrix data" << endl;
3476 int readSuccess = 1;
3477 std::ostringstream errMsg;
3478 if (myRank == rootRank) {
3481 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3483 reader_type reader (pAdder);
3486 std::pair<bool, std::vector<size_t> > results =
3487 reader.read (in, lineNumber, tolerant, debug);
3488 readSuccess = results.first ? 1 : 0;
3490 catch (std::exception& e) {
3495 broadcast (*pComm, rootRank, ptr (&readSuccess));
3504 TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
3505 "Failed to read the Matrix Market sparse matrix file: "
3509 if (debug && myRank == rootRank) {
3510 cerr <<
"-- Successfully read the Matrix Market data" << endl;
3523 if (debug && myRank == rootRank) {
3524 cerr <<
"-- Tolerant mode: rebroadcasting matrix dimensions"
3526 <<
"----- Dimensions before: "
3527 << dims[0] <<
" x " << dims[1]
3531 Tuple<global_ordinal_type, 2> updatedDims;
3532 if (myRank == rootRank) {
3539 std::max (dims[0], pAdder->getAdder()->numRows());
3540 updatedDims[1] = pAdder->getAdder()->numCols();
3542 broadcast (*pComm, rootRank, updatedDims);
3543 dims[0] = updatedDims[0];
3544 dims[1] = updatedDims[1];
3545 if (debug && myRank == rootRank) {
3546 cerr <<
"----- Dimensions after: " << dims[0] <<
" x "
3560 if (myRank == rootRank) {
3567 if (dims[0] < pAdder->getAdder ()->numRows ()) {
3571 broadcast (*pComm, 0, ptr (&dimsMatch));
3572 if (dimsMatch == 0) {
3579 Tuple<global_ordinal_type, 2> addersDims;
3580 if (myRank == rootRank) {
3581 addersDims[0] = pAdder->getAdder()->numRows();
3582 addersDims[1] = pAdder->getAdder()->numCols();
3584 broadcast (*pComm, 0, addersDims);
3585 TEUCHOS_TEST_FOR_EXCEPTION(
3586 dimsMatch == 0, std::runtime_error,
3587 "The matrix metadata says that the matrix is " << dims[0] <<
" x "
3588 << dims[1] <<
", but the actual data says that the matrix is "
3589 << addersDims[0] <<
" x " << addersDims[1] <<
". That means the "
3590 "data includes more rows than reported in the metadata. This "
3591 "is not allowed when parsing in strict mode. Parse the matrix in "
3592 "tolerant mode to ignore the metadata when it disagrees with the "
3597 if (debug && myRank == rootRank) {
3598 cerr <<
"-- Converting matrix data into CSR format on Proc 0" << endl;
3610 ArrayRCP<size_t> numEntriesPerRow;
3611 ArrayRCP<size_t> rowPtr;
3612 ArrayRCP<global_ordinal_type> colInd;
3613 ArrayRCP<scalar_type> values;
3618 int mergeAndConvertSucceeded = 1;
3619 std::ostringstream errMsg;
3621 if (myRank == rootRank) {
3623 typedef Teuchos::MatrixMarket::Raw::Element<
scalar_type,
3633 const size_type numRows = dims[0];
3636 pAdder->getAdder()->merge ();
3639 const std::vector<element_type>& entries =
3640 pAdder->getAdder()->getEntries();
3643 const size_t numEntries = (size_t)entries.size();
3646 cerr <<
"----- Proc 0: Matrix has numRows=" << numRows
3647 <<
" rows and numEntries=" << numEntries
3648 <<
" entries." << endl;
3654 numEntriesPerRow = arcp<size_t> (numRows);
3655 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3656 rowPtr = arcp<size_t> (numRows+1);
3657 std::fill (rowPtr.begin(), rowPtr.end(), 0);
3658 colInd = arcp<global_ordinal_type> (numEntries);
3659 values = arcp<scalar_type> (numEntries);
3666 for (curPos = 0; curPos < numEntries; ++curPos) {
3667 const element_type& curEntry = entries[curPos];
3669 TEUCHOS_TEST_FOR_EXCEPTION(
3670 curRow < prvRow, std::logic_error,
3671 "Row indices are out of order, even though they are supposed "
3672 "to be sorted. curRow = " << curRow <<
", prvRow = "
3673 << prvRow <<
", at curPos = " << curPos <<
". Please report "
3674 "this bug to the Tpetra developers.");
3675 if (curRow > prvRow) {
3681 numEntriesPerRow[curRow]++;
3682 colInd[curPos] = curEntry.colIndex();
3683 values[curPos] = curEntry.value();
3688 rowPtr[numRows] = numEntries;
3690 catch (std::exception& e) {
3691 mergeAndConvertSucceeded = 0;
3692 errMsg <<
"Failed to merge sparse matrix entries and convert to "
3693 "CSR format: " << e.what();
3696 if (debug && mergeAndConvertSucceeded) {
3698 const size_type numRows = dims[0];
3699 const size_type maxToDisplay = 100;
3701 cerr <<
"----- Proc 0: numEntriesPerRow[0.."
3702 << (numEntriesPerRow.size()-1) <<
"] ";
3703 if (numRows > maxToDisplay) {
3704 cerr <<
"(only showing first and last few entries) ";
3708 if (numRows > maxToDisplay) {
3709 for (size_type k = 0; k < 2; ++k) {
3710 cerr << numEntriesPerRow[k] <<
" ";
3713 for (size_type k = numRows-2; k < numRows-1; ++k) {
3714 cerr << numEntriesPerRow[k] <<
" ";
3718 for (size_type k = 0; k < numRows-1; ++k) {
3719 cerr << numEntriesPerRow[k] <<
" ";
3722 cerr << numEntriesPerRow[numRows-1];
3724 cerr <<
"]" << endl;
3726 cerr <<
"----- Proc 0: rowPtr ";
3727 if (numRows > maxToDisplay) {
3728 cerr <<
"(only showing first and last few entries) ";
3731 if (numRows > maxToDisplay) {
3732 for (size_type k = 0; k < 2; ++k) {
3733 cerr << rowPtr[k] <<
" ";
3736 for (size_type k = numRows-2; k < numRows; ++k) {
3737 cerr << rowPtr[k] <<
" ";
3741 for (size_type k = 0; k < numRows; ++k) {
3742 cerr << rowPtr[k] <<
" ";
3745 cerr << rowPtr[numRows] <<
"]" << endl;
3747 cerr <<
"----- Proc 0: colInd = [";
3748 for (
size_t k = 0; k < rowPtr[numRows]; ++k) {
3749 cerr << colInd[k] <<
" ";
3751 cerr <<
"]" << endl;
3765 if (debug && myRank == rootRank) {
3766 cerr <<
"-- Verifying Maps" << endl;
3768 TEUCHOS_TEST_FOR_EXCEPTION(
3769 as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
3770 std::invalid_argument,
3771 "The range Map has " << rangeMap->getGlobalNumElements ()
3772 <<
" entries, but the matrix has a global number of rows " << dims[0]
3774 TEUCHOS_TEST_FOR_EXCEPTION(
3775 as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
3776 std::invalid_argument,
3777 "The domain Map has " << domainMap->getGlobalNumElements ()
3778 <<
" entries, but the matrix has a global number of columns "
3782 RCP<Teuchos::FancyOStream> err = debug ?
3783 Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
3785 RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
3786 ArrayView<const global_ordinal_type> myRows =
3787 gatherRowMap->getLocalElementList ();
3788 const size_type myNumRows = myRows.size ();
3791 ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
3792 for (size_type i_ = 0; i_ < myNumRows; i_++) {
3793 gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_]-indexBase];
3799 RCP<sparse_matrix_type> A_proc0 =
3801 if (myRank == rootRank) {
3803 cerr <<
"-- Proc 0: Filling gather matrix" << endl;
3806 cerr <<
"---- Rows: " << Teuchos::toString (myRows) << endl;
3810 for (size_type i_ = 0; i_ < myNumRows; ++i_) {
3811 size_type i = myRows[i_] - indexBase;
3813 const size_type curPos = as<size_type> (rowPtr[i]);
3815 ArrayView<global_ordinal_type> curColInd =
3816 colInd.view (curPos, curNumEntries);
3817 ArrayView<scalar_type> curValues =
3818 values.view (curPos, curNumEntries);
3821 for (size_type k = 0; k < curNumEntries; ++k) {
3822 curColInd[k] += indexBase;
3825 cerr <<
"------ Columns: " << Teuchos::toString (curColInd) << endl;
3826 cerr <<
"------ Values: " << Teuchos::toString (curValues) << endl;
3829 if (curNumEntries > 0) {
3830 A_proc0->insertGlobalValues (myRows[i_], curColInd, curValues);
3836 numEntriesPerRow = null;
3842 RCP<sparse_matrix_type> A;
3844 export_type exp (gatherRowMap, rowMap);
3851 mv_type_go target_nnzPerRow(rowMap,1);
3852 mv_type_go source_nnzPerRow(gatherRowMap,1);
3853 Teuchos::ArrayRCP<GO> srcData = source_nnzPerRow.getDataNonConst(0);
3854 for (
int i=0; i<myNumRows; i++)
3855 srcData[i] = gatherNumEntriesPerRow[i];
3856 srcData = Teuchos::null;
3858 Teuchos::ArrayRCP<GO> targetData = target_nnzPerRow.getDataNonConst(0);
3859 ArrayRCP<size_t> targetData_size_t = arcp<size_t>(targetData.size());
3860 for (
int i=0; i<targetData.size(); i++)
3861 targetData_size_t[i] = targetData[i];
3863 if (colMap.is_null ()) {
3868 A->doExport (*A_proc0, exp,
INSERT);
3869 if (callFillComplete) {
3870 A->fillComplete (domainMap, rangeMap);
3882 if (callFillComplete) {
3883 const int numProcs = pComm->getSize ();
3885 if (extraDebug && debug) {
3886 const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
3887 const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
3888 if (myRank == rootRank) {
3889 cerr <<
"-- Matrix is "
3890 << globalNumRows <<
" x " << globalNumCols
3891 <<
" with " << A->getGlobalNumEntries()
3892 <<
" entries, and index base "
3893 << A->getIndexBase() <<
"." << endl;
3896 for (
int p = 0; p < numProcs; ++p) {
3898 cerr <<
"-- Proc " << p <<
" owns "
3899 << A->getLocalNumCols() <<
" columns, and "
3900 << A->getLocalNumEntries() <<
" entries." << endl;
3907 if (debug && myRank == rootRank) {
3908 cerr <<
"-- Done creating the CrsMatrix from the Matrix Market data"
3943 static Teuchos::RCP<multivector_type>
3945 const Teuchos::RCP<const comm_type>& comm,
3946 Teuchos::RCP<const map_type>& map,
3947 const bool tolerant=
false,
3948 const bool debug=
false,
3949 const bool binary=
false)
3951 using Teuchos::broadcast;
3952 using Teuchos::outArg;
3956 if (comm->getRank() == 0) {
3959 in.open (filename.c_str ());
3961 in.open (filename.c_str (), std::ios::binary);
3962 opened = in.is_open();
3968 broadcast<int, int> (*comm, 0, outArg (opened));
3969 TEUCHOS_TEST_FOR_EXCEPTION(
3970 opened == 0, std::runtime_error,
3971 "readDenseFile: Failed to open file \"" << filename <<
"\" on "
3973 return readDense (in, comm, map, tolerant, debug, binary);
4006 static Teuchos::RCP<vector_type>
4008 const Teuchos::RCP<const comm_type>& comm,
4009 Teuchos::RCP<const map_type>& map,
4010 const bool tolerant=
false,
4011 const bool debug=
false)
4013 using Teuchos::broadcast;
4014 using Teuchos::outArg;
4018 if (comm->getRank() == 0) {
4020 in.open (filename.c_str ());
4021 opened = in.is_open();
4027 broadcast<int, int> (*comm, 0, outArg (opened));
4028 TEUCHOS_TEST_FOR_EXCEPTION(
4029 opened == 0, std::runtime_error,
4030 "readVectorFile: Failed to open file \"" << filename <<
"\" on "
4032 return readVector (in, comm, map, tolerant, debug);
4104 static Teuchos::RCP<multivector_type>
4106 const Teuchos::RCP<const comm_type>& comm,
4107 Teuchos::RCP<const map_type>& map,
4108 const bool tolerant=
false,
4109 const bool debug=
false,
4110 const bool binary=
false)
4112 Teuchos::RCP<Teuchos::FancyOStream> err =
4113 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4114 return readDenseImpl<scalar_type> (in, comm, map, err, tolerant, debug, binary);
4119 static Teuchos::RCP<vector_type>
4121 const Teuchos::RCP<const comm_type>& comm,
4122 Teuchos::RCP<const map_type>& map,
4123 const bool tolerant=
false,
4124 const bool debug=
false)
4126 Teuchos::RCP<Teuchos::FancyOStream> err =
4127 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4128 return readVectorImpl<scalar_type> (in, comm, map, err, tolerant, debug);
4153 static Teuchos::RCP<const map_type>
4155 const Teuchos::RCP<const comm_type>& comm,
4156 const bool tolerant=
false,
4157 const bool debug=
false,
4158 const bool binary=
false)
4160 using Teuchos::inOutArg;
4161 using Teuchos::broadcast;
4165 if (comm->getRank () == 0) {
4167 in.open (filename.c_str (), std::ios::binary);
4169 in.open (filename.c_str ());
4174 broadcast<int, int> (*comm, 0, inOutArg (success));
4175 TEUCHOS_TEST_FOR_EXCEPTION(
4176 success == 0, std::runtime_error,
4177 "Tpetra::MatrixMarket::Reader::readMapFile: "
4178 "Failed to read file \"" << filename <<
"\" on Process 0.");
4179 return readMap (in, comm, tolerant, debug, binary);
4184 template<
class MultiVectorScalarType>
4189 readDenseImpl (std::istream& in,
4190 const Teuchos::RCP<const comm_type>& comm,
4191 Teuchos::RCP<const map_type>& map,
4192 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4193 const bool tolerant=
false,
4194 const bool debug=
false,
4195 const bool binary=
false)
4197 using Teuchos::MatrixMarket::Banner;
4198 using Teuchos::MatrixMarket::checkCommentLine;
4199 using Teuchos::ArrayRCP;
4201 using Teuchos::broadcast;
4202 using Teuchos::outArg;
4204 using Teuchos::Tuple;
4206 typedef MultiVectorScalarType ST;
4210 typedef Teuchos::ScalarTraits<ST> STS;
4211 typedef typename STS::magnitudeType MT;
4212 typedef Teuchos::ScalarTraits<MT> STM;
4217 const int myRank = comm->getRank ();
4219 if (! err.is_null ()) {
4223 *err << myRank <<
": readDenseImpl" << endl;
4225 if (! err.is_null ()) {
4260 size_t lineNumber = 1;
4263 std::ostringstream exMsg;
4264 int localBannerReadSuccess = 1;
4265 int localDimsReadSuccess = 1;
4271 *err << myRank <<
": readDenseImpl: Reading banner line (dense)" << endl;
4277 RCP<const Banner> pBanner;
4279 pBanner = readBanner (in, lineNumber, tolerant, debug);
4280 }
catch (std::exception& e) {
4282 localBannerReadSuccess = 0;
4285 if (localBannerReadSuccess) {
4286 if (pBanner->matrixType () !=
"array") {
4287 exMsg <<
"The Matrix Market file does not contain dense matrix "
4288 "data. Its banner (first) line says that its matrix type is \""
4289 << pBanner->matrixType () <<
"\", rather that the required "
4291 localBannerReadSuccess = 0;
4292 }
else if (pBanner->dataType() ==
"pattern") {
4293 exMsg <<
"The Matrix Market file's banner (first) "
4294 "line claims that the matrix's data type is \"pattern\". This does "
4295 "not make sense for a dense matrix, yet the file reports the matrix "
4296 "as dense. The only valid data types for a dense matrix are "
4297 "\"real\", \"complex\", and \"integer\".";
4298 localBannerReadSuccess = 0;
4302 dims[2] = encodeDataType (pBanner->dataType ());
4308 if (localBannerReadSuccess) {
4310 *err << myRank <<
": readDenseImpl: Reading dimensions line (dense)" << endl;
4319 bool commentLine =
true;
4321 while (commentLine) {
4324 if (in.eof () || in.fail ()) {
4325 exMsg <<
"Unable to get array dimensions line (at all) from line "
4326 << lineNumber <<
" of input stream. The input stream "
4327 <<
"claims that it is "
4328 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4329 localDimsReadSuccess = 0;
4332 if (getline (in, line)) {
4339 size_t start = 0, size = 0;
4340 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4347 std::istringstream istr (line);
4351 if (istr.eof () || istr.fail ()) {
4352 exMsg <<
"Unable to read any data from line " << lineNumber
4353 <<
" of input; the line should contain the matrix dimensions "
4354 <<
"\"<numRows> <numCols>\".";
4355 localDimsReadSuccess = 0;
4360 exMsg <<
"Failed to get number of rows from line "
4361 << lineNumber <<
" of input; the line should contains the "
4362 <<
"matrix dimensions \"<numRows> <numCols>\".";
4363 localDimsReadSuccess = 0;
4365 dims[0] = theNumRows;
4367 exMsg <<
"No more data after number of rows on line "
4368 << lineNumber <<
" of input; the line should contain the "
4369 <<
"matrix dimensions \"<numRows> <numCols>\".";
4370 localDimsReadSuccess = 0;
4375 exMsg <<
"Failed to get number of columns from line "
4376 << lineNumber <<
" of input; the line should contain "
4377 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4378 localDimsReadSuccess = 0;
4380 dims[1] = theNumCols;
4389 in.read(
reinterpret_cast<char*
>(&numRows),
sizeof(numRows));
4390 in.read(
reinterpret_cast<char*
>(&numCols),
sizeof(numCols));
4391 dims[0] = Teuchos::as<GO>(numRows);
4392 dims[1] = Teuchos::as<GO>(numCols);
4393 if ((
typeid(ST) ==
typeid(
double)) || Teuchos::ScalarTraits<ST>::isOrdinal) {
4395 }
else if (Teuchos::ScalarTraits<ST>::isComplex) {
4398 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
4399 "Unrecognized Matrix Market data type. "
4400 "We should never get here. "
4401 "Please report this bug to the Tpetra developers.");
4403 localDimsReadSuccess =
true;
4404 localDimsReadSuccess =
true;
4411 Tuple<GO, 5> bannerDimsReadResult;
4413 bannerDimsReadResult[0] = dims[0];
4414 bannerDimsReadResult[1] = dims[1];
4415 bannerDimsReadResult[2] = dims[2];
4416 bannerDimsReadResult[3] = localBannerReadSuccess;
4417 bannerDimsReadResult[4] = localDimsReadSuccess;
4421 broadcast (*comm, 0, bannerDimsReadResult);
4423 TEUCHOS_TEST_FOR_EXCEPTION(
4424 bannerDimsReadResult[3] == 0, std::runtime_error,
4425 "Failed to read banner line: " << exMsg.str ());
4426 TEUCHOS_TEST_FOR_EXCEPTION(
4427 bannerDimsReadResult[4] == 0, std::runtime_error,
4428 "Failed to read matrix dimensions line: " << exMsg.str ());
4430 dims[0] = bannerDimsReadResult[0];
4431 dims[1] = bannerDimsReadResult[1];
4432 dims[2] = bannerDimsReadResult[2];
4437 const size_t numCols =
static_cast<size_t> (dims[1]);
4442 RCP<const map_type> proc0Map;
4443 if (map.is_null ()) {
4447 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4448 const size_t localNumRows = (myRank == 0) ? numRows : 0;
4450 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4454 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4458 RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
4464 int localReadDataSuccess = 1;
4469 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4474 TEUCHOS_TEST_FOR_EXCEPTION(
4475 ! X->isConstantStride (), std::logic_error,
4476 "Can't get a 1-D view of the entries of the MultiVector X on "
4477 "Process 0, because the stride between the columns of X is not "
4478 "constant. This shouldn't happen because we just created X and "
4479 "haven't filled it in yet. Please report this bug to the Tpetra "
4486 ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4487 TEUCHOS_TEST_FOR_EXCEPTION(
4488 as<global_size_t> (X_view.size ()) < numRows * numCols,
4490 "The view of X has size " << X_view.size() <<
" which is not enough to "
4491 "accommodate the expected number of entries numRows*numCols = "
4492 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
4493 "Please report this bug to the Tpetra developers.");
4494 const size_t stride = X->getStride ();
4501 const bool isComplex = (dims[2] == 1);
4502 size_type count = 0, curRow = 0, curCol = 0;
4505 while (getline (in, line)) {
4509 size_t start = 0, size = 0;
4510 const bool commentLine =
4511 checkCommentLine (line, start, size, lineNumber, tolerant);
4512 if (! commentLine) {
4518 if (count >= X_view.size()) {
4523 TEUCHOS_TEST_FOR_EXCEPTION(
4524 count >= X_view.size(),
4526 "The Matrix Market input stream has more data in it than "
4527 "its metadata reported. Current line number is "
4528 << lineNumber <<
".");
4537 const size_t pos = line.substr (start, size).find (
':');
4538 if (pos != std::string::npos) {
4542 std::istringstream istr (line.substr (start, size));
4546 if (istr.eof() || istr.fail()) {
4553 TEUCHOS_TEST_FOR_EXCEPTION(
4554 ! tolerant && (istr.eof() || istr.fail()),
4556 "Line " << lineNumber <<
" of the Matrix Market file is "
4557 "empty, or we cannot read from it for some other reason.");
4560 ST val = STS::zero();
4563 MT real = STM::zero(), imag = STM::zero();
4580 TEUCHOS_TEST_FOR_EXCEPTION(
4581 ! tolerant && istr.eof(), std::runtime_error,
4582 "Failed to get the real part of a complex-valued matrix "
4583 "entry from line " << lineNumber <<
" of the Matrix Market "
4589 }
else if (istr.eof()) {
4590 TEUCHOS_TEST_FOR_EXCEPTION(
4591 ! tolerant && istr.eof(), std::runtime_error,
4592 "Missing imaginary part of a complex-valued matrix entry "
4593 "on line " << lineNumber <<
" of the Matrix Market file.");
4599 TEUCHOS_TEST_FOR_EXCEPTION(
4600 ! tolerant && istr.fail(), std::runtime_error,
4601 "Failed to get the imaginary part of a complex-valued "
4602 "matrix entry from line " << lineNumber <<
" of the "
4603 "Matrix Market file.");
4610 TEUCHOS_TEST_FOR_EXCEPTION(
4611 ! tolerant && istr.fail(), std::runtime_error,
4612 "Failed to get a real-valued matrix entry from line "
4613 << lineNumber <<
" of the Matrix Market file.");
4616 if (istr.fail() && tolerant) {
4622 TEUCHOS_TEST_FOR_EXCEPTION(
4623 ! tolerant && istr.fail(), std::runtime_error,
4624 "Failed to read matrix data from line " << lineNumber
4625 <<
" of the Matrix Market file.");
4628 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4630 curRow = count % numRows;
4631 curCol = count / numRows;
4632 X_view[curRow + curCol*stride] = val;
4637 TEUCHOS_TEST_FOR_EXCEPTION(
4638 ! tolerant &&
static_cast<global_size_t> (count) < numRows * numCols,
4640 "The Matrix Market metadata reports that the dense matrix is "
4641 << numRows <<
" x " << numCols <<
", and thus has "
4642 << numRows*numCols <<
" total entries, but we only found " << count
4643 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
4644 }
catch (std::exception& e) {
4646 localReadDataSuccess = 0;
4651 *err << myRank <<
": readDenseImpl: Reading matrix data (dense)"
4656 TEUCHOS_TEST_FOR_EXCEPTION(
4657 ! X->isConstantStride (), std::logic_error,
4658 "Can't get a 1-D view of the entries of the MultiVector X on "
4659 "Process 0, because the stride between the columns of X is not "
4660 "constant. This shouldn't happen because we just created X and "
4661 "haven't filled it in yet. Please report this bug to the Tpetra "
4668 auto X_view = X->getLocalViewHost (Access::OverwriteAll);
4670 TEUCHOS_TEST_FOR_EXCEPTION(
4671 as<global_size_t> (X_view.extent(0)) < numRows,
4673 "The view of X has " << X_view.extent(0) <<
" rows which is not enough to "
4674 "accommodate the expected number of entries numRows = "
4676 "Please report this bug to the Tpetra developers.");
4677 TEUCHOS_TEST_FOR_EXCEPTION(
4678 as<global_size_t> (X_view.extent(1)) < numCols,
4680 "The view of X has " << X_view.extent(1) <<
" colums which is not enough to "
4681 "accommodate the expected number of entries numRows = "
4683 "Please report this bug to the Tpetra developers.");
4685 for (
size_t curRow = 0; curRow < numRows; ++curRow) {
4686 for (
size_t curCol = 0; curCol < numCols; ++curCol) {
4687 if (Teuchos::ScalarTraits<ST>::isOrdinal){
4689 in.read(
reinterpret_cast<char*
>(&val),
sizeof(val));
4690 X_view(curRow, curCol) = val;
4693 in.read(
reinterpret_cast<char*
>(&val),
sizeof(val));
4694 X_view(curRow, curCol) = val;
4702 *err << myRank <<
": readDenseImpl: done reading data" << endl;
4706 int globalReadDataSuccess = localReadDataSuccess;
4707 broadcast (*comm, 0, outArg (globalReadDataSuccess));
4708 TEUCHOS_TEST_FOR_EXCEPTION(
4709 globalReadDataSuccess == 0, std::runtime_error,
4710 "Failed to read the multivector's data: " << exMsg.str ());
4715 if (comm->getSize () == 1 && map.is_null ()) {
4717 if (! err.is_null ()) {
4721 *err << myRank <<
": readDenseImpl: done" << endl;
4723 if (! err.is_null ()) {
4730 *err << myRank <<
": readDenseImpl: Creating target MV" << endl;
4734 RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
4737 *err << myRank <<
": readDenseImpl: Creating Export" << endl;
4743 Export<LO, GO, NT> exporter (proc0Map, map, err);
4746 *err << myRank <<
": readDenseImpl: Exporting" << endl;
4749 Y->doExport (*X, exporter,
INSERT);
4751 if (! err.is_null ()) {
4755 *err << myRank <<
": readDenseImpl: done" << endl;
4757 if (! err.is_null ()) {
4766 template<
class VectorScalarType>
4771 readVectorImpl (std::istream& in,
4772 const Teuchos::RCP<const comm_type>& comm,
4773 Teuchos::RCP<const map_type>& map,
4774 const Teuchos::RCP<Teuchos::FancyOStream>& err,
4775 const bool tolerant=
false,
4776 const bool debug=
false)
4778 using Teuchos::MatrixMarket::Banner;
4779 using Teuchos::MatrixMarket::checkCommentLine;
4781 using Teuchos::broadcast;
4782 using Teuchos::outArg;
4784 using Teuchos::Tuple;
4786 typedef VectorScalarType ST;
4790 typedef Teuchos::ScalarTraits<ST> STS;
4791 typedef typename STS::magnitudeType MT;
4792 typedef Teuchos::ScalarTraits<MT> STM;
4797 const int myRank = comm->getRank ();
4799 if (! err.is_null ()) {
4803 *err << myRank <<
": readVectorImpl" << endl;
4805 if (! err.is_null ()) {
4839 size_t lineNumber = 1;
4842 std::ostringstream exMsg;
4843 int localBannerReadSuccess = 1;
4844 int localDimsReadSuccess = 1;
4849 *err << myRank <<
": readVectorImpl: Reading banner line (dense)" << endl;
4855 RCP<const Banner> pBanner;
4857 pBanner = readBanner (in, lineNumber, tolerant, debug);
4858 }
catch (std::exception& e) {
4860 localBannerReadSuccess = 0;
4863 if (localBannerReadSuccess) {
4864 if (pBanner->matrixType () !=
"array") {
4865 exMsg <<
"The Matrix Market file does not contain dense matrix "
4866 "data. Its banner (first) line says that its matrix type is \""
4867 << pBanner->matrixType () <<
"\", rather that the required "
4869 localBannerReadSuccess = 0;
4870 }
else if (pBanner->dataType() ==
"pattern") {
4871 exMsg <<
"The Matrix Market file's banner (first) "
4872 "line claims that the matrix's data type is \"pattern\". This does "
4873 "not make sense for a dense matrix, yet the file reports the matrix "
4874 "as dense. The only valid data types for a dense matrix are "
4875 "\"real\", \"complex\", and \"integer\".";
4876 localBannerReadSuccess = 0;
4880 dims[2] = encodeDataType (pBanner->dataType ());
4886 if (localBannerReadSuccess) {
4888 *err << myRank <<
": readVectorImpl: Reading dimensions line (dense)" << endl;
4897 bool commentLine =
true;
4899 while (commentLine) {
4902 if (in.eof () || in.fail ()) {
4903 exMsg <<
"Unable to get array dimensions line (at all) from line "
4904 << lineNumber <<
" of input stream. The input stream "
4905 <<
"claims that it is "
4906 << (in.eof() ?
"at end-of-file." :
"in a failed state.");
4907 localDimsReadSuccess = 0;
4910 if (getline (in, line)) {
4917 size_t start = 0, size = 0;
4918 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4925 std::istringstream istr (line);
4929 if (istr.eof () || istr.fail ()) {
4930 exMsg <<
"Unable to read any data from line " << lineNumber
4931 <<
" of input; the line should contain the matrix dimensions "
4932 <<
"\"<numRows> <numCols>\".";
4933 localDimsReadSuccess = 0;
4938 exMsg <<
"Failed to get number of rows from line "
4939 << lineNumber <<
" of input; the line should contains the "
4940 <<
"matrix dimensions \"<numRows> <numCols>\".";
4941 localDimsReadSuccess = 0;
4943 dims[0] = theNumRows;
4945 exMsg <<
"No more data after number of rows on line "
4946 << lineNumber <<
" of input; the line should contain the "
4947 <<
"matrix dimensions \"<numRows> <numCols>\".";
4948 localDimsReadSuccess = 0;
4953 exMsg <<
"Failed to get number of columns from line "
4954 << lineNumber <<
" of input; the line should contain "
4955 <<
"the matrix dimensions \"<numRows> <numCols>\".";
4956 localDimsReadSuccess = 0;
4958 dims[1] = theNumCols;
4968 exMsg <<
"File does not contain a 1D Vector";
4969 localDimsReadSuccess = 0;
4975 Tuple<GO, 5> bannerDimsReadResult;
4977 bannerDimsReadResult[0] = dims[0];
4978 bannerDimsReadResult[1] = dims[1];
4979 bannerDimsReadResult[2] = dims[2];
4980 bannerDimsReadResult[3] = localBannerReadSuccess;
4981 bannerDimsReadResult[4] = localDimsReadSuccess;
4986 broadcast (*comm, 0, bannerDimsReadResult);
4988 TEUCHOS_TEST_FOR_EXCEPTION(
4989 bannerDimsReadResult[3] == 0, std::runtime_error,
4990 "Failed to read banner line: " << exMsg.str ());
4991 TEUCHOS_TEST_FOR_EXCEPTION(
4992 bannerDimsReadResult[4] == 0, std::runtime_error,
4993 "Failed to read matrix dimensions line: " << exMsg.str ());
4995 dims[0] = bannerDimsReadResult[0];
4996 dims[1] = bannerDimsReadResult[1];
4997 dims[2] = bannerDimsReadResult[2];
5002 const size_t numCols =
static_cast<size_t> (dims[1]);
5007 RCP<const map_type> proc0Map;
5008 if (map.is_null ()) {
5012 map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
5013 const size_t localNumRows = (myRank == 0) ? numRows : 0;
5015 proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
5019 proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
5023 RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
5029 int localReadDataSuccess = 1;
5033 *err << myRank <<
": readVectorImpl: Reading matrix data (dense)"
5038 TEUCHOS_TEST_FOR_EXCEPTION(
5039 ! X->isConstantStride (), std::logic_error,
5040 "Can't get a 1-D view of the entries of the MultiVector X on "
5041 "Process 0, because the stride between the columns of X is not "
5042 "constant. This shouldn't happen because we just created X and "
5043 "haven't filled it in yet. Please report this bug to the Tpetra "
5050 Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst ();
5051 TEUCHOS_TEST_FOR_EXCEPTION(
5052 as<global_size_t> (X_view.size ()) < numRows * numCols,
5054 "The view of X has size " << X_view <<
" which is not enough to "
5055 "accommodate the expected number of entries numRows*numCols = "
5056 << numRows <<
"*" << numCols <<
" = " << numRows*numCols <<
". "
5057 "Please report this bug to the Tpetra developers.");
5058 const size_t stride = X->getStride ();
5065 const bool isComplex = (dims[2] == 1);
5066 size_type count = 0, curRow = 0, curCol = 0;
5069 while (getline (in, line)) {
5073 size_t start = 0, size = 0;
5074 const bool commentLine =
5075 checkCommentLine (line, start, size, lineNumber, tolerant);
5076 if (! commentLine) {
5082 if (count >= X_view.size()) {
5087 TEUCHOS_TEST_FOR_EXCEPTION(
5088 count >= X_view.size(),
5090 "The Matrix Market input stream has more data in it than "
5091 "its metadata reported. Current line number is "
5092 << lineNumber <<
".");
5101 const size_t pos = line.substr (start, size).find (
':');
5102 if (pos != std::string::npos) {
5106 std::istringstream istr (line.substr (start, size));
5110 if (istr.eof() || istr.fail()) {
5117 TEUCHOS_TEST_FOR_EXCEPTION(
5118 ! tolerant && (istr.eof() || istr.fail()),
5120 "Line " << lineNumber <<
" of the Matrix Market file is "
5121 "empty, or we cannot read from it for some other reason.");
5124 ST val = STS::zero();
5127 MT real = STM::zero(), imag = STM::zero();
5144 TEUCHOS_TEST_FOR_EXCEPTION(
5145 ! tolerant && istr.eof(), std::runtime_error,
5146 "Failed to get the real part of a complex-valued matrix "
5147 "entry from line " << lineNumber <<
" of the Matrix Market "
5153 }
else if (istr.eof()) {
5154 TEUCHOS_TEST_FOR_EXCEPTION(
5155 ! tolerant && istr.eof(), std::runtime_error,
5156 "Missing imaginary part of a complex-valued matrix entry "
5157 "on line " << lineNumber <<
" of the Matrix Market file.");
5163 TEUCHOS_TEST_FOR_EXCEPTION(
5164 ! tolerant && istr.fail(), std::runtime_error,
5165 "Failed to get the imaginary part of a complex-valued "
5166 "matrix entry from line " << lineNumber <<
" of the "
5167 "Matrix Market file.");
5174 TEUCHOS_TEST_FOR_EXCEPTION(
5175 ! tolerant && istr.fail(), std::runtime_error,
5176 "Failed to get a real-valued matrix entry from line "
5177 << lineNumber <<
" of the Matrix Market file.");
5180 if (istr.fail() && tolerant) {
5186 TEUCHOS_TEST_FOR_EXCEPTION(
5187 ! tolerant && istr.fail(), std::runtime_error,
5188 "Failed to read matrix data from line " << lineNumber
5189 <<
" of the Matrix Market file.");
5192 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
5194 curRow = count % numRows;
5195 curCol = count / numRows;
5196 X_view[curRow + curCol*stride] = val;
5201 TEUCHOS_TEST_FOR_EXCEPTION(
5202 ! tolerant &&
static_cast<global_size_t> (count) < numRows * numCols,
5204 "The Matrix Market metadata reports that the dense matrix is "
5205 << numRows <<
" x " << numCols <<
", and thus has "
5206 << numRows*numCols <<
" total entries, but we only found " << count
5207 <<
" entr" << (count == 1 ?
"y" :
"ies") <<
" in the file.");
5208 }
catch (std::exception& e) {
5210 localReadDataSuccess = 0;
5215 *err << myRank <<
": readVectorImpl: done reading data" << endl;
5219 int globalReadDataSuccess = localReadDataSuccess;
5220 broadcast (*comm, 0, outArg (globalReadDataSuccess));
5221 TEUCHOS_TEST_FOR_EXCEPTION(
5222 globalReadDataSuccess == 0, std::runtime_error,
5223 "Failed to read the multivector's data: " << exMsg.str ());
5228 if (comm->getSize () == 1 && map.is_null ()) {
5230 if (! err.is_null ()) {
5234 *err << myRank <<
": readVectorImpl: done" << endl;
5236 if (! err.is_null ()) {
5243 *err << myRank <<
": readVectorImpl: Creating target MV" << endl;
5247 RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
5250 *err << myRank <<
": readVectorImpl: Creating Export" << endl;
5256 Export<LO, GO, NT> exporter (proc0Map, map, err);
5259 *err << myRank <<
": readVectorImpl: Exporting" << endl;
5262 Y->doExport (*X, exporter,
INSERT);
5264 if (! err.is_null ()) {
5268 *err << myRank <<
": readVectorImpl: done" << endl;
5270 if (! err.is_null ()) {
5299 static Teuchos::RCP<const map_type>
5301 const Teuchos::RCP<const comm_type>& comm,
5302 const bool tolerant=
false,
5303 const bool debug=
false,
5304 const bool binary=
false)
5306 Teuchos::RCP<Teuchos::FancyOStream> err =
5307 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
5308 return readMap (in, comm, err, tolerant, debug, binary);
5338 static Teuchos::RCP<const map_type>
5340 const Teuchos::RCP<const comm_type>& comm,
5341 const Teuchos::RCP<Teuchos::FancyOStream>& err,
5342 const bool tolerant=
false,
5343 const bool debug=
false,
5344 const bool binary=
false)
5346 using Teuchos::arcp;
5347 using Teuchos::Array;
5348 using Teuchos::ArrayRCP;
5350 using Teuchos::broadcast;
5351 using Teuchos::Comm;
5352 using Teuchos::CommRequest;
5353 using Teuchos::inOutArg;
5354 using Teuchos::ireceive;
5355 using Teuchos::outArg;
5357 using Teuchos::receive;
5358 using Teuchos::reduceAll;
5359 using Teuchos::REDUCE_MIN;
5360 using Teuchos::isend;
5361 using Teuchos::SerialComm;
5362 using Teuchos::toString;
5363 using Teuchos::wait;
5366 typedef ptrdiff_t int_type;
5372 const int numProcs = comm->getSize ();
5373 const int myRank = comm->getRank ();
5375 if (err.is_null ()) {
5379 std::ostringstream os;
5380 os << myRank <<
": readMap: " << endl;
5383 if (err.is_null ()) {
5389 const int sizeTag = 1339;
5391 const int dataTag = 1340;
5397 RCP<CommRequest<int> > sizeReq;
5398 RCP<CommRequest<int> > dataReq;
5403 ArrayRCP<int_type> numGidsToRecv (1);
5404 numGidsToRecv[0] = 0;
5406 sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
5409 int readSuccess = 1;
5410 std::ostringstream exMsg;
5419 RCP<const Comm<int> > proc0Comm (
new SerialComm<int> ());
5421 RCP<const map_type> dataMap;
5425 data = readDenseImpl<GO> (in, proc0Comm, dataMap, err, tolerant, debug, binary);
5427 if (data.is_null ()) {
5429 exMsg <<
"readDenseImpl() returned null." << endl;
5431 }
catch (std::exception& e) {
5433 exMsg << e.what () << endl;
5439 std::map<int, Array<GO> > pid2gids;
5444 int_type globalNumGIDs = 0;
5454 if (myRank == 0 && readSuccess == 1) {
5455 if (data->getNumVectors () == 2) {
5456 ArrayRCP<const GO> GIDs = data->getData (0);
5457 ArrayRCP<const GO> PIDs = data->getData (1);
5458 globalNumGIDs = GIDs.size ();
5462 if (globalNumGIDs > 0) {
5463 const int pid =
static_cast<int> (PIDs[0]);
5465 if (pid < 0 || pid >= numProcs) {
5467 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5468 <<
"Encountered invalid PID " << pid <<
"." << endl;
5471 const GO gid = GIDs[0];
5472 pid2gids[pid].push_back (gid);
5476 if (readSuccess == 1) {
5479 for (size_type k = 1; k < globalNumGIDs; ++k) {
5480 const int pid =
static_cast<int> (PIDs[k]);
5481 if (pid < 0 || pid >= numProcs) {
5483 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5484 <<
"Encountered invalid PID " << pid <<
"." << endl;
5487 const int_type gid = GIDs[k];
5488 pid2gids[pid].push_back (gid);
5489 if (gid < indexBase) {
5496 else if (data->getNumVectors () == 1) {
5497 if (data->getGlobalLength () % 2 !=
static_cast<GST
> (0)) {
5499 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data has the "
5500 "wrong format (for Map format 2.0). The global number of rows "
5501 "in the MultiVector must be even (divisible by 2)." << endl;
5504 ArrayRCP<const GO> theData = data->getData (0);
5505 globalNumGIDs =
static_cast<GO
> (data->getGlobalLength ()) /
5506 static_cast<GO
> (2);
5510 if (globalNumGIDs > 0) {
5511 const int pid =
static_cast<int> (theData[1]);
5512 if (pid < 0 || pid >= numProcs) {
5514 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5515 <<
"Encountered invalid PID " << pid <<
"." << endl;
5518 const GO gid = theData[0];
5519 pid2gids[pid].push_back (gid);
5525 for (int_type k = 1; k < globalNumGIDs; ++k) {
5526 const int pid =
static_cast<int> (theData[2*k + 1]);
5527 if (pid < 0 || pid >= numProcs) {
5529 exMsg <<
"Tpetra::MatrixMarket::readMap: "
5530 <<
"Encountered invalid PID " << pid <<
"." << endl;
5533 const GO gid = theData[2*k];
5534 pid2gids[pid].push_back (gid);
5535 if (gid < indexBase) {
5544 exMsg <<
"Tpetra::MatrixMarket::readMap: Input data must have "
5545 "either 1 column (for the new Map format 2.0) or 2 columns (for "
5546 "the old Map format 1.0).";
5553 int_type readResults[3];
5554 readResults[0] =
static_cast<int_type
> (indexBase);
5555 readResults[1] =
static_cast<int_type
> (globalNumGIDs);
5556 readResults[2] =
static_cast<int_type
> (readSuccess);
5557 broadcast<int, int_type> (*comm, 0, 3, readResults);
5559 indexBase =
static_cast<GO
> (readResults[0]);
5560 globalNumGIDs =
static_cast<int_type
> (readResults[1]);
5561 readSuccess =
static_cast<int> (readResults[2]);
5567 TEUCHOS_TEST_FOR_EXCEPTION(
5568 readSuccess != 1, std::runtime_error,
5569 "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
5570 "following exception message: " << exMsg.str ());
5574 for (
int p = 1; p < numProcs; ++p) {
5575 ArrayRCP<int_type> numGidsToSend (1);
5577 auto it = pid2gids.find (p);
5578 if (it == pid2gids.end ()) {
5579 numGidsToSend[0] = 0;
5581 numGidsToSend[0] = it->second.size ();
5583 sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
5584 wait<int> (*comm, outArg (sizeReq));
5589 wait<int> (*comm, outArg (sizeReq));
5594 ArrayRCP<GO> myGids;
5595 int_type myNumGids = 0;
5597 GO* myGidsRaw = NULL;
5599 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
5600 if (it != pid2gids.end ()) {
5601 myGidsRaw = it->second.getRawPtr ();
5602 myNumGids = it->second.size ();
5604 myGids = arcp<GO> (myGidsRaw, 0, myNumGids,
false);
5608 myNumGids = numGidsToRecv[0];
5609 myGids = arcp<GO> (myNumGids);
5616 if (myNumGids > 0) {
5617 dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
5621 for (
int p = 1; p < numProcs; ++p) {
5623 ArrayRCP<GO> sendGids;
5624 GO* sendGidsRaw = NULL;
5625 int_type numSendGids = 0;
5627 typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
5628 if (it != pid2gids.end ()) {
5629 numSendGids = it->second.size ();
5630 sendGidsRaw = it->second.getRawPtr ();
5631 sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids,
false);
5634 if (numSendGids > 0) {
5635 dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
5637 wait<int> (*comm, outArg (dataReq));
5639 else if (myRank == p) {
5641 wait<int> (*comm, outArg (dataReq));
5646 std::ostringstream os;
5647 os << myRank <<
": readMap: creating Map" << endl;
5650 const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
5651 RCP<const map_type> newMap;
5658 std::ostringstream errStrm;
5660 newMap = rcp (
new map_type (INVALID, myGids (), indexBase, comm));
5662 catch (std::exception& e) {
5664 errStrm <<
"Process " << comm->getRank () <<
" threw an exception: "
5665 << e.what () << std::endl;
5669 errStrm <<
"Process " << comm->getRank () <<
" threw an exception "
5670 "not a subclass of std::exception" << std::endl;
5672 Teuchos::reduceAll<int, int> (*comm, Teuchos::REDUCE_MIN,
5673 lclSuccess, Teuchos::outArg (gblSuccess));
5674 if (gblSuccess != 1) {
5677 TEUCHOS_TEST_FOR_EXCEPTION(gblSuccess != 1, std::runtime_error,
"Map constructor failed!");
5679 if (err.is_null ()) {
5683 std::ostringstream os;
5684 os << myRank <<
": readMap: done" << endl;
5687 if (err.is_null ()) {
5707 encodeDataType (
const std::string& dataType)
5709 if (dataType ==
"real" || dataType ==
"integer") {
5711 }
else if (dataType ==
"complex") {
5713 }
else if (dataType ==
"pattern") {
5719 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
5720 "Unrecognized Matrix Market data type \"" << dataType
5721 <<
"\". We should never get here. "
5722 "Please report this bug to the Tpetra developers.");
5757 static Teuchos::RCP<sparse_matrix_type>
5759 const std::string& filename_suffix,
5760 const Teuchos::RCP<const map_type>& rowMap,
5761 Teuchos::RCP<const map_type>& colMap,
5762 const Teuchos::RCP<const map_type>& domainMap,
5763 const Teuchos::RCP<const map_type>& rangeMap,
5764 const bool callFillComplete=
true,
5765 const bool tolerant=
false,
5766 const int ranksToReadAtOnce=8,
5767 const bool debug=
false)
5772 using STS =
typename Teuchos::ScalarTraits<ST>;
5774 using Teuchos::ArrayRCP;
5775 using Teuchos::arcp;
5782 TEUCHOS_TEST_FOR_EXCEPTION(
5783 rowMap.is_null (), std::invalid_argument,
5784 "Row Map must be nonnull.");
5785 Teuchos::RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
5786 TEUCHOS_TEST_FOR_EXCEPTION
5787 (comm.is_null (), std::invalid_argument,
5788 "The input row map's communicator (Teuchos::Comm object) is null.");
5789 TEUCHOS_TEST_FOR_EXCEPTION(
5790 rangeMap.is_null (), std::invalid_argument,
5791 "Range Map must be nonnull.");
5792 TEUCHOS_TEST_FOR_EXCEPTION(
5793 domainMap.is_null (), std::invalid_argument,
5794 "Domain Map must be nonnull.");
5795 TEUCHOS_TEST_FOR_EXCEPTION(
5796 domainMap->getComm().getRawPtr() != comm.getRawPtr(),
5797 std::invalid_argument,
5798 "The specified domain Map's communicator (domainMap->getComm())"
5799 "differs from row Map's communicator");
5800 TEUCHOS_TEST_FOR_EXCEPTION(
5801 rangeMap->getComm().getRawPtr() != comm.getRawPtr(),
5802 std::invalid_argument,
5803 "The specified range Map's communicator (rangeMap->getComm())"
5804 "differs from row Map's communicator");
5807 const int myRank = comm->getRank();
5808 const int numProc = comm->getSize();
5809 std::string filename = filename_prefix + std::to_string(myRank) + filename_suffix;
5812 int rank_limit = std::min(std::max(ranksToReadAtOnce,1),numProc);
5815 ArrayRCP<size_t> numEntriesPerRow;
5816 ArrayRCP<size_t> rowPtr;
5817 ArrayRCP<global_ordinal_type> colInd;
5818 ArrayRCP<scalar_type> values;
5819 std::ostringstream errMsg;
5826 bool success =
true;
5827 int bannerIsCorrect = 1, readSuccess = 1;
5828 LO numRows, numCols, numNonzeros;
5829 for(
int base_rank = 0; base_rank < numProc; base_rank += rank_limit) {
5830 int stop = std::min(base_rank+rank_limit,numProc);
5833 if(base_rank <= myRank && myRank < stop) {
5835 std::ifstream in(filename);
5836 using Teuchos::MatrixMarket::Banner;
5837 size_t lineNumber = 1;
5838 RCP<const Banner> pBanner;
5840 pBanner = readBanner (in, lineNumber, tolerant, debug);
5842 catch (std::exception& e) {
5843 errMsg <<
"Attempt to read the Matrix Market file's Banner line "
5844 "threw an exception: " << e.what();
5845 bannerIsCorrect = 0;
5847 if (bannerIsCorrect) {
5852 if (! tolerant && pBanner->matrixType() !=
"coordinate") {
5853 bannerIsCorrect = 0;
5854 errMsg <<
"The Matrix Market input file must contain a "
5855 "\"coordinate\"-format sparse matrix in order to create a "
5856 "Tpetra::CrsMatrix object from it, but the file's matrix "
5857 "type is \"" << pBanner->matrixType() <<
"\" instead.";
5862 if (tolerant && pBanner->matrixType() ==
"array") {
5863 bannerIsCorrect = 0;
5864 errMsg <<
"Matrix Market file must contain a \"coordinate\"-"
5865 "format sparse matrix in order to create a Tpetra::CrsMatrix "
5866 "object from it, but the file's matrix type is \"array\" "
5867 "instead. That probably means the file contains dense matrix "
5873 using Teuchos::MatrixMarket::readCoordinateDimensions;
5874 success = readCoordinateDimensions (in, numRows, numCols,
5875 numNonzeros, lineNumber,
5879 TEUCHOS_TEST_FOR_EXCEPTION(numRows != (LO)rowMap->getLocalNumElements(), std::invalid_argument,
5880 "# rows in file does not match rowmap.");
5881 TEUCHOS_TEST_FOR_EXCEPTION(!colMap.is_null() && numCols != (LO)colMap->getLocalNumElements(), std::invalid_argument,
5882 "# rows in file does not match colmap.");
5886 typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type,global_ordinal_type> raw_adder_type;
5887 bool tolerant_required =
true;
5888 Teuchos::RCP<raw_adder_type> pRaw =
5889 Teuchos::rcp (
new raw_adder_type (numRows,numCols,numNonzeros,tolerant_required,debug));
5890 RCP<adder_type> pAdder = Teuchos::rcp (
new adder_type (pRaw, pBanner->symmType ()));
5893 std::cerr <<
"-- Reading matrix data" << std::endl;
5898 typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
5900 reader_type reader (pAdder);
5903 std::pair<bool, std::vector<size_t> > results = reader.read (in, lineNumber, tolerant_required, debug);
5906 readSuccess = results.first ? 1 : 0;
5908 catch (std::exception& e) {
5915 typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,global_ordinal_type> element_type;
5918 pAdder->getAdder()->merge ();
5921 const std::vector<element_type>& entries = pAdder->getAdder()->getEntries();
5924 const size_t numEntries = (size_t)entries.size();
5927 std::cerr <<
"----- Proc "<<myRank<<
": Matrix has numRows=" << numRows
5928 <<
" rows and numEntries=" << numEntries
5929 <<
" entries." << std::endl;
5936 numEntriesPerRow = arcp<size_t> (numRows);
5937 std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
5938 rowPtr = arcp<size_t> (numRows+1);
5939 std::fill (rowPtr.begin(), rowPtr.end(), 0);
5940 colInd = arcp<global_ordinal_type> (numEntries);
5941 values = arcp<scalar_type> (numEntries);
5947 LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
5949 LO indexBase = rowMap->getIndexBase();
5950 for (curPos = 0; curPos < numEntries; ++curPos) {
5951 const element_type& curEntry = entries[curPos];
5953 LO l_curRow = rowMap->getLocalElement(curRow);
5956 TEUCHOS_TEST_FOR_EXCEPTION(l_curRow == INVALID,std::logic_error,
5957 "Current global row "<< curRow <<
" is invalid.");
5959 TEUCHOS_TEST_FOR_EXCEPTION(l_curRow < l_prvRow, std::logic_error,
5960 "Row indices are out of order, even though they are supposed "
5961 "to be sorted. curRow = " << l_curRow <<
", prvRow = "
5962 << l_prvRow <<
", at curPos = " << curPos <<
". Please report "
5963 "this bug to the Tpetra developers.");
5964 if (l_curRow > l_prvRow) {
5965 for (LO r = l_prvRow+1; r <= l_curRow; ++r) {
5968 l_prvRow = l_curRow;
5970 numEntriesPerRow[l_curRow]++;
5971 colInd[curPos] = curEntry.colIndex() + indexBase;
5972 values[curPos] = curEntry.value();
5978 rowPtr[numRows] = numEntries;
5989 RCP<sparse_matrix_type> A;
5990 if(colMap.is_null()) {
5992 for(
size_t i=0; i<rowMap->getLocalNumElements(); i++) {
5993 GO g_row = rowMap->getGlobalElement(i);
5994 size_t start = rowPtr[i];
5995 size_t size = rowPtr[i+1] - rowPtr[i];
5997 A->insertGlobalValues(g_row,size,&values[start],&colInd[start]);
6002 throw std::runtime_error(
"Reading with a column map is not yet implemented");
6004 RCP<const map_type> myDomainMap = domainMap.is_null() ? rowMap : domainMap;
6005 RCP<const map_type> myRangeMap = rangeMap.is_null() ? rowMap : rangeMap;
6007 A->fillComplete(myDomainMap,myRangeMap);
6011 TEUCHOS_TEST_FOR_EXCEPTION(success ==
false, std::runtime_error,
6048 template<
class SparseMatrixType>
6053 typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
6115 const std::string& matrixName,
6116 const std::string& matrixDescription,
6117 const bool debug=
false)
6119 Teuchos::RCP<const Teuchos::Comm<int> > comm = matrix.getComm ();
6120 TEUCHOS_TEST_FOR_EXCEPTION
6121 (comm.is_null (), std::invalid_argument,
6122 "The input matrix's communicator (Teuchos::Comm object) is null.");
6123 const int myRank = comm->getRank ();
6128 out.open (filename.c_str ());
6130 writeSparse (out, matrix, matrixName, matrixDescription, debug);
6139 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6140 const std::string& matrixName,
6141 const std::string& matrixDescription,
6142 const bool debug=
false)
6144 TEUCHOS_TEST_FOR_EXCEPTION
6145 (pMatrix.is_null (), std::invalid_argument,
6146 "The input matrix is null.");
6148 matrixDescription, debug);
6173 const bool debug=
false)
6181 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6182 const bool debug=
false)
6220 const std::string& matrixName,
6221 const std::string& matrixDescription,
6222 const bool debug=
false)
6224 using Teuchos::ArrayView;
6225 using Teuchos::Comm;
6226 using Teuchos::FancyOStream;
6227 using Teuchos::getFancyOStream;
6228 using Teuchos::null;
6230 using Teuchos::rcpFromRef;
6236 using STS =
typename Teuchos::ScalarTraits<ST>;
6242 Teuchos::SetScientific<ST> sci (out);
6245 RCP<const Comm<int> > comm = matrix.getComm ();
6246 TEUCHOS_TEST_FOR_EXCEPTION(
6247 comm.is_null (), std::invalid_argument,
6248 "The input matrix's communicator (Teuchos::Comm object) is null.");
6249 const int myRank = comm->getRank ();
6252 RCP<FancyOStream> err =
6253 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6255 std::ostringstream os;
6256 os << myRank <<
": writeSparse" << endl;
6259 os <<
"-- " << myRank <<
": past barrier" << endl;
6264 const bool debugPrint = debug && myRank == 0;
6266 RCP<const map_type> rowMap = matrix.getRowMap ();
6267 RCP<const map_type> colMap = matrix.getColMap ();
6268 RCP<const map_type> domainMap = matrix.getDomainMap ();
6269 RCP<const map_type> rangeMap = matrix.getRangeMap ();
6271 const global_size_t numRows = rangeMap->getGlobalNumElements ();
6272 const global_size_t numCols = domainMap->getGlobalNumElements ();
6274 if (debug && myRank == 0) {
6275 std::ostringstream os;
6276 os <<
"-- Input sparse matrix is:"
6277 <<
"---- " << numRows <<
" x " << numCols << endl
6279 << (matrix.isGloballyIndexed() ?
"Globally" :
"Locally")
6280 <<
" indexed." << endl
6281 <<
"---- Its row map has " << rowMap->getGlobalNumElements ()
6282 <<
" elements." << endl
6283 <<
"---- Its col map has " << colMap->getGlobalNumElements ()
6284 <<
" elements." << endl;
6289 const size_t localNumRows = (myRank == 0) ? numRows : 0;
6291 std::ostringstream os;
6292 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6295 RCP<const map_type> gatherRowMap =
6296 Details::computeGatherMap (rowMap, err, debug);
6304 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6305 RCP<const map_type> gatherColMap =
6306 Details::computeGatherMap (colMap, err, debug);
6310 import_type importer (rowMap, gatherRowMap);
6315 RCP<sparse_matrix_type> newMatrix =
6317 static_cast<size_t> (0)));
6319 newMatrix->doImport (matrix, importer,
INSERT);
6324 RCP<const map_type> gatherDomainMap =
6325 rcp (
new map_type (numCols, localNumCols,
6326 domainMap->getIndexBase (),
6328 RCP<const map_type> gatherRangeMap =
6329 rcp (
new map_type (numRows, localNumRows,
6330 rangeMap->getIndexBase (),
6332 newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
6336 cerr <<
"-- Output sparse matrix is:"
6337 <<
"---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
6339 << newMatrix->getDomainMap ()->getGlobalNumElements ()
6341 << newMatrix->getGlobalNumEntries () <<
" entries;" << endl
6343 << (newMatrix->isGloballyIndexed () ?
"Globally" :
"Locally")
6344 <<
" indexed." << endl
6345 <<
"---- Its row map has "
6346 << newMatrix->getRowMap ()->getGlobalNumElements ()
6347 <<
" elements, with index base "
6348 << newMatrix->getRowMap ()->getIndexBase () <<
"." << endl
6349 <<
"---- Its col map has "
6350 << newMatrix->getColMap ()->getGlobalNumElements ()
6351 <<
" elements, with index base "
6352 << newMatrix->getColMap ()->getIndexBase () <<
"." << endl
6353 <<
"---- Element count of output matrix's column Map may differ "
6354 <<
"from that of the input matrix's column Map, if some columns "
6355 <<
"of the matrix contain no entries." << endl;
6368 out <<
"%%MatrixMarket matrix coordinate "
6369 << (STS::isComplex ?
"complex" :
"real")
6370 <<
" general" << endl;
6373 if (matrixName !=
"") {
6374 printAsComment (out, matrixName);
6376 if (matrixDescription !=
"") {
6377 printAsComment (out, matrixDescription);
6387 out << newMatrix->getRangeMap ()->getGlobalNumElements () <<
" "
6388 << newMatrix->getDomainMap ()->getGlobalNumElements () <<
" "
6389 << newMatrix->getGlobalNumEntries () << endl;
6394 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6395 const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
6403 if (newMatrix->isGloballyIndexed()) {
6406 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6407 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6408 for (GO globalRowIndex = minAllGlobalIndex;
6409 globalRowIndex <= maxAllGlobalIndex;
6411 typename sparse_matrix_type::global_inds_host_view_type ind;
6412 typename sparse_matrix_type::values_host_view_type val;
6413 newMatrix->getGlobalRowView (globalRowIndex, ind, val);
6414 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6415 const GO globalColIndex = ind(ii);
6418 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6419 << (globalColIndex + 1 - colIndexBase) <<
" ";
6420 if (STS::isComplex) {
6421 out << STS::real (val(ii)) <<
" " << STS::imag (val(ii));
6430 using OTG = Teuchos::OrdinalTraits<GO>;
6431 for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
6432 localRowIndex <= gatherRowMap->getMaxLocalIndex();
6435 const GO globalRowIndex =
6436 gatherRowMap->getGlobalElement (localRowIndex);
6437 TEUCHOS_TEST_FOR_EXCEPTION(
6438 globalRowIndex == OTG::invalid(), std::logic_error,
6439 "Failed to convert the supposed local row index "
6440 << localRowIndex <<
" into a global row index. "
6441 "Please report this bug to the Tpetra developers.");
6442 typename sparse_matrix_type::local_inds_host_view_type ind;
6443 typename sparse_matrix_type::values_host_view_type val;
6444 newMatrix->getLocalRowView (localRowIndex, ind, val);
6445 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6447 const GO globalColIndex =
6448 newMatrix->getColMap()->getGlobalElement (ind(ii));
6449 TEUCHOS_TEST_FOR_EXCEPTION(
6450 globalColIndex == OTG::invalid(), std::logic_error,
6451 "On local row " << localRowIndex <<
" of the sparse matrix: "
6452 "Failed to convert the supposed local column index "
6453 << ind(ii) <<
" into a global column index. Please report "
6454 "this bug to the Tpetra developers.");
6457 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6458 << (globalColIndex + 1 - colIndexBase) <<
" ";
6459 if (STS::isComplex) {
6460 out << STS::real (val(ii)) <<
" " << STS::imag (val(ii));
6474 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6475 const std::string& matrixName,
6476 const std::string& matrixDescription,
6477 const bool debug=
false)
6479 TEUCHOS_TEST_FOR_EXCEPTION
6480 (pMatrix.is_null (), std::invalid_argument,
6481 "The input matrix is null.");
6482 writeSparse (out, *pMatrix, matrixName, matrixDescription, debug);
6518 const std::string& graphName,
6519 const std::string& graphDescription,
6520 const bool debug=
false)
6522 using Teuchos::ArrayView;
6523 using Teuchos::Comm;
6524 using Teuchos::FancyOStream;
6525 using Teuchos::getFancyOStream;
6526 using Teuchos::null;
6528 using Teuchos::rcpFromRef;
6539 if (rowMap.is_null ()) {
6542 auto comm = rowMap->getComm ();
6543 if (comm.is_null ()) {
6546 const int myRank = comm->getRank ();
6549 RCP<FancyOStream> err =
6550 debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6552 std::ostringstream os;
6553 os << myRank <<
": writeSparseGraph" << endl;
6556 os <<
"-- " << myRank <<
": past barrier" << endl;
6561 const bool debugPrint = debug && myRank == 0;
6568 const global_size_t numRows = rangeMap->getGlobalNumElements ();
6569 const global_size_t numCols = domainMap->getGlobalNumElements ();
6571 if (debug && myRank == 0) {
6572 std::ostringstream os;
6573 os <<
"-- Input sparse graph is:"
6574 <<
"---- " << numRows <<
" x " << numCols <<
" with "
6578 <<
" indexed." << endl
6579 <<
"---- Its row Map has " << rowMap->getGlobalNumElements ()
6580 <<
" elements." << endl
6581 <<
"---- Its col Map has " << colMap->getGlobalNumElements ()
6582 <<
" elements." << endl;
6587 const size_t localNumRows = (myRank == 0) ? numRows : 0;
6589 std::ostringstream os;
6590 os <<
"-- " << myRank <<
": making gatherRowMap" << endl;
6593 auto gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
6601 const size_t localNumCols = (myRank == 0) ? numCols : 0;
6602 auto gatherColMap = Details::computeGatherMap (colMap, err, debug);
6611 static_cast<size_t> (0));
6618 RCP<const map_type> gatherDomainMap =
6619 rcp (
new map_type (numCols, localNumCols,
6620 domainMap->getIndexBase (),
6622 RCP<const map_type> gatherRangeMap =
6623 rcp (
new map_type (numRows, localNumRows,
6624 rangeMap->getIndexBase (),
6626 newGraph.
fillComplete (gatherDomainMap, gatherRangeMap);
6630 cerr <<
"-- Output sparse graph is:"
6631 <<
"---- " << newGraph.
getRangeMap ()->getGlobalNumElements ()
6638 <<
" indexed." << endl
6639 <<
"---- Its row map has "
6640 << newGraph.
getRowMap ()->getGlobalNumElements ()
6641 <<
" elements, with index base "
6642 << newGraph.
getRowMap ()->getIndexBase () <<
"." << endl
6643 <<
"---- Its col map has "
6644 << newGraph.
getColMap ()->getGlobalNumElements ()
6645 <<
" elements, with index base "
6646 << newGraph.
getColMap ()->getIndexBase () <<
"." << endl
6647 <<
"---- Element count of output graph's column Map may differ "
6648 <<
"from that of the input matrix's column Map, if some columns "
6649 <<
"of the matrix contain no entries." << endl;
6663 out <<
"%%MatrixMarket matrix coordinate pattern general" << endl;
6666 if (graphName !=
"") {
6667 printAsComment (out, graphName);
6669 if (graphDescription !=
"") {
6670 printAsComment (out, graphDescription);
6681 out << newGraph.
getRangeMap ()->getGlobalNumElements () <<
" "
6682 << newGraph.
getDomainMap ()->getGlobalNumElements () <<
" "
6688 const GO rowIndexBase = gatherRowMap->getIndexBase ();
6689 const GO colIndexBase = newGraph.
getColMap()->getIndexBase ();
6700 const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6701 const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6702 for (GO globalRowIndex = minAllGlobalIndex;
6703 globalRowIndex <= maxAllGlobalIndex;
6705 typename crs_graph_type::global_inds_host_view_type ind;
6707 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6708 const GO globalColIndex = ind(ii);
6711 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6712 << (globalColIndex + 1 - colIndexBase) <<
" ";
6718 typedef Teuchos::OrdinalTraits<GO> OTG;
6719 for (LO localRowIndex = gatherRowMap->getMinLocalIndex ();
6720 localRowIndex <= gatherRowMap->getMaxLocalIndex ();
6723 const GO globalRowIndex =
6724 gatherRowMap->getGlobalElement (localRowIndex);
6725 TEUCHOS_TEST_FOR_EXCEPTION
6726 (globalRowIndex == OTG::invalid (), std::logic_error,
"Failed "
6727 "to convert the supposed local row index " << localRowIndex <<
6728 " into a global row index. Please report this bug to the "
6729 "Tpetra developers.");
6730 typename crs_graph_type::local_inds_host_view_type ind;
6732 for (
size_t ii = 0; ii < ind.extent(0); ii++) {
6734 const GO globalColIndex =
6735 newGraph.
getColMap ()->getGlobalElement (ind(ii));
6736 TEUCHOS_TEST_FOR_EXCEPTION(
6737 globalColIndex == OTG::invalid(), std::logic_error,
6738 "On local row " << localRowIndex <<
" of the sparse graph: "
6739 "Failed to convert the supposed local column index "
6740 << ind(ii) <<
" into a global column index. Please report "
6741 "this bug to the Tpetra developers.");
6744 out << (globalRowIndex + 1 - rowIndexBase) <<
" "
6745 << (globalColIndex + 1 - colIndexBase) <<
" ";
6761 const bool debug=
false)
6803 const std::string& graphName,
6804 const std::string& graphDescription,
6805 const bool debug=
false)
6808 if (comm.is_null ()) {
6816 const int myRank = comm->getRank ();
6821 out.open (filename.c_str ());
6836 const bool debug=
false)
6851 const Teuchos::RCP<const crs_graph_type>& pGraph,
6852 const std::string& graphName,
6853 const std::string& graphDescription,
6854 const bool debug=
false)
6870 const Teuchos::RCP<const crs_graph_type>& pGraph,
6871 const bool debug=
false)
6901 const bool debug=
false)
6909 const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6910 const bool debug=
false)
6946 const std::string& matrixName,
6947 const std::string& matrixDescription,
6948 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6949 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6951 const int myRank = X.
getMap ().is_null () ? 0 :
6952 (X.
getMap ()->getComm ().is_null () ? 0 :
6953 X.
getMap ()->getComm ()->getRank ());
6957 out.open (filename.c_str());
6960 writeDense (out, X, matrixName, matrixDescription, err, dbg);
6973 const Teuchos::RCP<const multivector_type>& X,
6974 const std::string& matrixName,
6975 const std::string& matrixDescription,
6976 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6977 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6979 TEUCHOS_TEST_FOR_EXCEPTION(
6980 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
6981 "writeDenseFile: The input MultiVector X is null.");
6982 writeDenseFile (filename, *X, matrixName, matrixDescription, err, dbg);
6993 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6994 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7006 const Teuchos::RCP<const multivector_type>& X,
7007 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7008 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7010 TEUCHOS_TEST_FOR_EXCEPTION(
7011 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7012 "writeDenseFile: The input MultiVector X is null.");
7050 const std::string& matrixName,
7051 const std::string& matrixDescription,
7052 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7053 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7055 using Teuchos::Comm;
7056 using Teuchos::outArg;
7057 using Teuchos::REDUCE_MAX;
7058 using Teuchos::reduceAll;
7062 RCP<const Comm<int> > comm = X.
getMap ().is_null () ?
7063 Teuchos::null : X.
getMap ()->getComm ();
7064 const int myRank = comm.is_null () ? 0 : comm->getRank ();
7069 const bool debug = ! dbg.is_null ();
7072 std::ostringstream os;
7073 os << myRank <<
": writeDense" << endl;
7078 writeDenseHeader (out, X, matrixName, matrixDescription, err, dbg);
7083 for (
size_t j = 0; j < numVecs; ++j) {
7084 writeDenseColumn (out, * (X.
getVector (j)), err, dbg);
7089 std::ostringstream os;
7090 os << myRank <<
": writeDense: Done" << endl;
7124 writeDenseHeader (std::ostream& out,
7126 const std::string& matrixName,
7127 const std::string& matrixDescription,
7128 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7129 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7131 using Teuchos::Comm;
7132 using Teuchos::outArg;
7134 using Teuchos::REDUCE_MAX;
7135 using Teuchos::reduceAll;
7137 typedef Teuchos::ScalarTraits<scalar_type> STS;
7138 const char prefix[] =
"Tpetra::MatrixMarket::writeDenseHeader: ";
7140 RCP<const Comm<int> > comm = X.getMap ().is_null () ?
7141 Teuchos::null : X.getMap ()->getComm ();
7142 const int myRank = comm.is_null () ? 0 : comm->getRank ();
7149 const bool debug = ! dbg.is_null ();
7153 std::ostringstream os;
7154 os << myRank <<
": writeDenseHeader" << endl;
7173 std::ostringstream hdr;
7175 std::string dataType;
7176 if (STS::isComplex) {
7177 dataType =
"complex";
7178 }
else if (STS::isOrdinal) {
7179 dataType =
"integer";
7183 hdr <<
"%%MatrixMarket matrix array " << dataType <<
" general"
7188 if (matrixName !=
"") {
7189 printAsComment (hdr, matrixName);
7191 if (matrixDescription !=
"") {
7192 printAsComment (hdr, matrixDescription);
7195 hdr << X.getGlobalLength () <<
" " << X.getNumVectors () << endl;
7199 }
catch (std::exception& e) {
7200 if (! err.is_null ()) {
7201 *err << prefix <<
"While writing the Matrix Market header, "
7202 "Process 0 threw an exception: " << e.what () << endl;
7211 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
7212 TEUCHOS_TEST_FOR_EXCEPTION(
7213 gblErr == 1, std::runtime_error, prefix <<
"Some error occurred "
7214 "which prevented this method from completing.");
7218 *dbg << myRank <<
": writeDenseHeader: Done" << endl;
7241 writeDenseColumn (std::ostream& out,
7243 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7244 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7246 using Teuchos::arcp;
7247 using Teuchos::Array;
7248 using Teuchos::ArrayRCP;
7249 using Teuchos::ArrayView;
7250 using Teuchos::Comm;
7251 using Teuchos::CommRequest;
7252 using Teuchos::ireceive;
7253 using Teuchos::isend;
7254 using Teuchos::outArg;
7255 using Teuchos::REDUCE_MAX;
7256 using Teuchos::reduceAll;
7258 using Teuchos::TypeNameTraits;
7259 using Teuchos::wait;
7261 typedef Teuchos::ScalarTraits<scalar_type> STS;
7263 const Comm<int>& comm = * (X.getMap ()->getComm ());
7264 const int myRank = comm.getRank ();
7265 const int numProcs = comm.getSize ();
7272 const bool debug = ! dbg.is_null ();
7276 std::ostringstream os;
7277 os << myRank <<
": writeDenseColumn" << endl;
7286 Teuchos::SetScientific<scalar_type> sci (out);
7288 const size_t myNumRows = X.getLocalLength ();
7289 const size_t numCols = X.getNumVectors ();
7292 const int sizeTag = 1337;
7293 const int dataTag = 1338;
7354 RCP<CommRequest<int> > sendReqSize, sendReqData;
7360 Array<ArrayRCP<size_t> > recvSizeBufs (3);
7361 Array<ArrayRCP<scalar_type> > recvDataBufs (3);
7362 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7363 Array<RCP<CommRequest<int> > > recvDataReqs (3);
7366 ArrayRCP<size_t> sendDataSize (1);
7367 sendDataSize[0] = myNumRows;
7371 std::ostringstream os;
7372 os << myRank <<
": Post receive-size receives from "
7373 "Procs 1 and 2: tag = " << sizeTag << endl;
7377 recvSizeBufs[0].resize (1);
7381 (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7382 recvSizeBufs[1].resize (1);
7383 (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7384 recvSizeBufs[2].resize (1);
7385 (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7388 ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
7392 ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
7395 else if (myRank == 1 || myRank == 2) {
7397 std::ostringstream os;
7398 os << myRank <<
": Post send-size send: size = "
7399 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7406 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7410 std::ostringstream os;
7411 os << myRank <<
": Not posting my send-size send yet" << endl;
7420 std::ostringstream os;
7421 os << myRank <<
": Pack my entries" << endl;
7424 ArrayRCP<scalar_type> sendDataBuf;
7426 sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
7427 X.get1dCopy (sendDataBuf (), myNumRows);
7429 catch (std::exception& e) {
7431 if (! err.is_null ()) {
7432 std::ostringstream os;
7433 os <<
"Process " << myRank <<
": Attempt to pack my MultiVector "
7434 "entries threw an exception: " << e.what () << endl;
7439 std::ostringstream os;
7440 os << myRank <<
": Done packing my entries" << endl;
7449 *dbg << myRank <<
": Post send-data send: tag = " << dataTag
7452 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7460 std::ostringstream os;
7461 os << myRank <<
": Write my entries" << endl;
7467 const size_t printNumRows = myNumRows;
7468 ArrayView<const scalar_type> printData = sendDataBuf ();
7469 const size_t printStride = printNumRows;
7470 if (
static_cast<size_t> (printData.size ()) < printStride * numCols) {
7472 if (! err.is_null ()) {
7473 std::ostringstream os;
7474 os <<
"Process " << myRank <<
": My MultiVector data's size "
7475 << printData.size () <<
" does not match my local dimensions "
7476 << printStride <<
" x " << numCols <<
"." << endl;
7484 for (
size_t col = 0; col < numCols; ++col) {
7485 for (
size_t row = 0; row < printNumRows; ++row) {
7486 if (STS::isComplex) {
7487 out << STS::real (printData[row + col * printStride]) <<
" "
7488 << STS::imag (printData[row + col * printStride]) << endl;
7490 out << printData[row + col * printStride] << endl;
7499 const int recvRank = 1;
7500 const int circBufInd = recvRank % 3;
7502 std::ostringstream os;
7503 os << myRank <<
": Wait on receive-size receive from Process "
7504 << recvRank << endl;
7508 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7512 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7513 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7515 if (! err.is_null ()) {
7516 std::ostringstream os;
7517 os << myRank <<
": Result of receive-size receive from Process "
7518 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7519 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7520 "This should never happen, and suggests that the receive never "
7521 "got posted. Please report this bug to the Tpetra developers."
7536 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7538 std::ostringstream os;
7539 os << myRank <<
": Post receive-data receive from Process "
7540 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7541 << recvDataBufs[circBufInd].size () << endl;
7544 if (! recvSizeReqs[circBufInd].is_null ()) {
7546 if (! err.is_null ()) {
7547 std::ostringstream os;
7548 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7549 "null, before posting the receive-data receive from Process "
7550 << recvRank <<
". This should never happen. Please report "
7551 "this bug to the Tpetra developers." << endl;
7555 recvDataReqs[circBufInd] =
7556 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7557 recvRank, dataTag, comm);
7560 else if (myRank == 1) {
7563 std::ostringstream os;
7564 os << myRank <<
": Wait on my send-size send" << endl;
7567 wait<int> (comm, outArg (sendReqSize));
7573 for (
int p = 1; p < numProcs; ++p) {
7575 if (p + 2 < numProcs) {
7577 const int recvRank = p + 2;
7578 const int circBufInd = recvRank % 3;
7580 std::ostringstream os;
7581 os << myRank <<
": Post receive-size receive from Process "
7582 << recvRank <<
": tag = " << sizeTag << endl;
7585 if (! recvSizeReqs[circBufInd].is_null ()) {
7587 if (! err.is_null ()) {
7588 std::ostringstream os;
7589 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
7590 <<
"null, for the receive-size receive from Process "
7591 << recvRank <<
"! This may mean that this process never "
7592 <<
"finished waiting for the receive from Process "
7593 << (recvRank - 3) <<
"." << endl;
7597 recvSizeReqs[circBufInd] =
7598 ireceive<int, size_t> (recvSizeBufs[circBufInd],
7599 recvRank, sizeTag, comm);
7602 if (p + 1 < numProcs) {
7603 const int recvRank = p + 1;
7604 const int circBufInd = recvRank % 3;
7608 std::ostringstream os;
7609 os << myRank <<
": Wait on receive-size receive from Process "
7610 << recvRank << endl;
7613 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7617 size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7618 if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7620 if (! err.is_null ()) {
7621 std::ostringstream os;
7622 os << myRank <<
": Result of receive-size receive from Process "
7623 << recvRank <<
" is Teuchos::OrdinalTraits<size_t>::invalid() "
7624 <<
"= " << Teuchos::OrdinalTraits<size_t>::invalid () <<
". "
7625 "This should never happen, and suggests that the receive never "
7626 "got posted. Please report this bug to the Tpetra developers."
7640 recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7642 std::ostringstream os;
7643 os << myRank <<
": Post receive-data receive from Process "
7644 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
7645 << recvDataBufs[circBufInd].size () << endl;
7648 if (! recvDataReqs[circBufInd].is_null ()) {
7650 if (! err.is_null ()) {
7651 std::ostringstream os;
7652 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
7653 <<
"null, for the receive-data receive from Process "
7654 << recvRank <<
"! This may mean that this process never "
7655 <<
"finished waiting for the receive from Process "
7656 << (recvRank - 3) <<
"." << endl;
7660 recvDataReqs[circBufInd] =
7661 ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7662 recvRank, dataTag, comm);
7666 const int recvRank = p;
7667 const int circBufInd = recvRank % 3;
7669 std::ostringstream os;
7670 os << myRank <<
": Wait on receive-data receive from Process "
7671 << recvRank << endl;
7674 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
7681 std::ostringstream os;
7682 os << myRank <<
": Write entries from Process " << recvRank
7684 *dbg << os.str () << endl;
7686 size_t printNumRows = (recvSizeBufs[circBufInd])[0];
7687 if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7689 if (! err.is_null ()) {
7690 std::ostringstream os;
7691 os << myRank <<
": Result of receive-size receive from Process "
7692 << recvRank <<
" was Teuchos::OrdinalTraits<size_t>::"
7693 "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
7694 <<
". This should never happen, and suggests that its "
7695 "receive-size receive was never posted. "
7696 "Please report this bug to the Tpetra developers." << endl;
7703 if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
7705 if (! err.is_null ()) {
7706 std::ostringstream os;
7707 os << myRank <<
": Result of receive-size receive from Proc "
7708 << recvRank <<
" was " << printNumRows <<
" > 0, but "
7709 "recvDataBufs[" << circBufInd <<
"] is null. This should "
7710 "never happen. Please report this bug to the Tpetra "
7711 "developers." << endl;
7721 ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
7722 const size_t printStride = printNumRows;
7726 for (
size_t col = 0; col < numCols; ++col) {
7727 for (
size_t row = 0; row < printNumRows; ++row) {
7728 if (STS::isComplex) {
7729 out << STS::real (printData[row + col * printStride]) <<
" "
7730 << STS::imag (printData[row + col * printStride]) << endl;
7732 out << printData[row + col * printStride] << endl;
7737 else if (myRank == p) {
7740 std::ostringstream os;
7741 os << myRank <<
": Wait on my send-data send" << endl;
7744 wait<int> (comm, outArg (sendReqData));
7746 else if (myRank == p + 1) {
7749 std::ostringstream os;
7750 os << myRank <<
": Post send-data send: tag = " << dataTag
7754 sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7757 std::ostringstream os;
7758 os << myRank <<
": Wait on my send-size send" << endl;
7761 wait<int> (comm, outArg (sendReqSize));
7763 else if (myRank == p + 2) {
7766 std::ostringstream os;
7767 os << myRank <<
": Post send-size send: size = "
7768 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
7771 sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7776 reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
7777 TEUCHOS_TEST_FOR_EXCEPTION(
7778 gblErr == 1, std::runtime_error,
"Tpetra::MatrixMarket::writeDense "
7779 "experienced some kind of error and was unable to complete.");
7783 *dbg << myRank <<
": writeDenseColumn: Done" << endl;
7797 const Teuchos::RCP<const multivector_type>& X,
7798 const std::string& matrixName,
7799 const std::string& matrixDescription,
7800 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7801 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7803 TEUCHOS_TEST_FOR_EXCEPTION(
7804 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7805 "writeDense: The input MultiVector X is null.");
7806 writeDense (out, *X, matrixName, matrixDescription, err, dbg);
7817 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7818 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7830 const Teuchos::RCP<const multivector_type>& X,
7831 const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7832 const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7834 TEUCHOS_TEST_FOR_EXCEPTION(
7835 X.is_null (), std::invalid_argument,
"Tpetra::MatrixMarket::"
7836 "writeDense: The input MultiVector X is null.");
7862 Teuchos::RCP<Teuchos::FancyOStream> err =
7863 Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
7878 const Teuchos::RCP<Teuchos::FancyOStream>& err,
7879 const bool debug=
false)
7881 using Teuchos::Array;
7882 using Teuchos::ArrayRCP;
7883 using Teuchos::ArrayView;
7884 using Teuchos::Comm;
7885 using Teuchos::CommRequest;
7886 using Teuchos::ireceive;
7887 using Teuchos::isend;
7889 using Teuchos::TypeNameTraits;
7890 using Teuchos::wait;
7893 typedef int pid_type;
7908 typedef ptrdiff_t int_type;
7909 TEUCHOS_TEST_FOR_EXCEPTION(
7910 sizeof (GO) >
sizeof (int_type), std::logic_error,
7911 "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
7912 <<
" is too big for ptrdiff_t. sizeof(GO) = " <<
sizeof (GO)
7913 <<
" > sizeof(ptrdiff_t) = " <<
sizeof (ptrdiff_t) <<
".");
7914 TEUCHOS_TEST_FOR_EXCEPTION(
7915 sizeof (pid_type) >
sizeof (int_type), std::logic_error,
7916 "The (MPI) process rank type pid_type=" <<
7917 TypeNameTraits<pid_type>::name () <<
" is too big for ptrdiff_t. "
7918 "sizeof(pid_type) = " <<
sizeof (pid_type) <<
" > sizeof(ptrdiff_t)"
7919 " = " <<
sizeof (ptrdiff_t) <<
".");
7921 const Comm<int>& comm = * (map.
getComm ());
7922 const int myRank = comm.getRank ();
7923 const int numProcs = comm.getSize ();
7925 if (! err.is_null ()) {
7929 std::ostringstream os;
7930 os << myRank <<
": writeMap" << endl;
7933 if (! err.is_null ()) {
7940 const int sizeTag = 1337;
7941 const int dataTag = 1338;
8000 RCP<CommRequest<int> > sendReqSize, sendReqData;
8006 Array<ArrayRCP<int_type> > recvSizeBufs (3);
8007 Array<ArrayRCP<int_type> > recvDataBufs (3);
8008 Array<RCP<CommRequest<int> > > recvSizeReqs (3);
8009 Array<RCP<CommRequest<int> > > recvDataReqs (3);
8012 ArrayRCP<int_type> sendDataSize (1);
8013 sendDataSize[0] = myNumRows;
8017 std::ostringstream os;
8018 os << myRank <<
": Post receive-size receives from "
8019 "Procs 1 and 2: tag = " << sizeTag << endl;
8023 recvSizeBufs[0].resize (1);
8024 (recvSizeBufs[0])[0] = -1;
8025 recvSizeBufs[1].resize (1);
8026 (recvSizeBufs[1])[0] = -1;
8027 recvSizeBufs[2].resize (1);
8028 (recvSizeBufs[2])[0] = -1;
8031 ireceive<int, int_type> (recvSizeBufs[1], 1, sizeTag, comm);
8035 ireceive<int, int_type> (recvSizeBufs[2], 2, sizeTag, comm);
8038 else if (myRank == 1 || myRank == 2) {
8040 std::ostringstream os;
8041 os << myRank <<
": Post send-size send: size = "
8042 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
8049 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
8053 std::ostringstream os;
8054 os << myRank <<
": Not posting my send-size send yet" << endl;
8065 std::ostringstream os;
8066 os << myRank <<
": Pack my GIDs and PIDs" << endl;
8070 ArrayRCP<int_type> sendDataBuf (myNumRows * 2);
8073 const int_type myMinGblIdx =
8075 for (
size_t k = 0; k < myNumRows; ++k) {
8076 const int_type gid = myMinGblIdx +
static_cast<int_type
> (k);
8077 const int_type pid =
static_cast<int_type
> (myRank);
8078 sendDataBuf[2*k] = gid;
8079 sendDataBuf[2*k+1] = pid;
8084 for (
size_t k = 0; k < myNumRows; ++k) {
8085 const int_type gid =
static_cast<int_type
> (myGblInds[k]);
8086 const int_type pid =
static_cast<int_type
> (myRank);
8087 sendDataBuf[2*k] = gid;
8088 sendDataBuf[2*k+1] = pid;
8093 std::ostringstream os;
8094 os << myRank <<
": Done packing my GIDs and PIDs" << endl;
8101 *err << myRank <<
": Post send-data send: tag = " << dataTag
8104 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
8109 *err << myRank <<
": Write MatrixMarket header" << endl;
8114 std::ostringstream hdr;
8118 hdr <<
"%%MatrixMarket matrix array integer general" << endl
8119 <<
"% Format: Version 2.0" << endl
8121 <<
"% This file encodes a Tpetra::Map." << endl
8122 <<
"% It is stored as a dense vector, with twice as many " << endl
8123 <<
"% entries as the global number of GIDs (global indices)." << endl
8124 <<
"% (GID, PID) pairs are stored contiguously, where the PID " << endl
8125 <<
"% is the rank of the process owning that GID." << endl
8130 std::ostringstream os;
8131 os << myRank <<
": Write my GIDs and PIDs" << endl;
8137 const int_type printNumRows = myNumRows;
8138 ArrayView<const int_type> printData = sendDataBuf ();
8139 for (int_type k = 0; k < printNumRows; ++k) {
8140 const int_type gid = printData[2*k];
8141 const int_type pid = printData[2*k+1];
8142 out << gid << endl << pid << endl;
8148 const int recvRank = 1;
8149 const int circBufInd = recvRank % 3;
8151 std::ostringstream os;
8152 os << myRank <<
": Wait on receive-size receive from Process "
8153 << recvRank << endl;
8157 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
8161 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
8162 if (debug && recvNumRows == -1) {
8163 std::ostringstream os;
8164 os << myRank <<
": Result of receive-size receive from Process "
8165 << recvRank <<
" is -1. This should never happen, and "
8166 "suggests that the receive never got posted. Please report "
8167 "this bug to the Tpetra developers." << endl;
8172 recvDataBufs[circBufInd].resize (recvNumRows * 2);
8174 std::ostringstream os;
8175 os << myRank <<
": Post receive-data receive from Process "
8176 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
8177 << recvDataBufs[circBufInd].size () << endl;
8180 if (! recvSizeReqs[circBufInd].is_null ()) {
8181 std::ostringstream os;
8182 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
8183 "null, before posting the receive-data receive from Process "
8184 << recvRank <<
". This should never happen. Please report "
8185 "this bug to the Tpetra developers." << endl;
8188 recvDataReqs[circBufInd] =
8189 ireceive<int, int_type> (recvDataBufs[circBufInd],
8190 recvRank, dataTag, comm);
8193 else if (myRank == 1) {
8196 std::ostringstream os;
8197 os << myRank <<
": Wait on my send-size send" << endl;
8200 wait<int> (comm, outArg (sendReqSize));
8206 for (
int p = 1; p < numProcs; ++p) {
8208 if (p + 2 < numProcs) {
8210 const int recvRank = p + 2;
8211 const int circBufInd = recvRank % 3;
8213 std::ostringstream os;
8214 os << myRank <<
": Post receive-size receive from Process "
8215 << recvRank <<
": tag = " << sizeTag << endl;
8218 if (! recvSizeReqs[circBufInd].is_null ()) {
8219 std::ostringstream os;
8220 os << myRank <<
": recvSizeReqs[" << circBufInd <<
"] is not "
8221 <<
"null, for the receive-size receive from Process "
8222 << recvRank <<
"! This may mean that this process never "
8223 <<
"finished waiting for the receive from Process "
8224 << (recvRank - 3) <<
"." << endl;
8227 recvSizeReqs[circBufInd] =
8228 ireceive<int, int_type> (recvSizeBufs[circBufInd],
8229 recvRank, sizeTag, comm);
8232 if (p + 1 < numProcs) {
8233 const int recvRank = p + 1;
8234 const int circBufInd = recvRank % 3;
8238 std::ostringstream os;
8239 os << myRank <<
": Wait on receive-size receive from Process "
8240 << recvRank << endl;
8243 wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
8247 const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
8248 if (debug && recvNumRows == -1) {
8249 std::ostringstream os;
8250 os << myRank <<
": Result of receive-size receive from Process "
8251 << recvRank <<
" is -1. This should never happen, and "
8252 "suggests that the receive never got posted. Please report "
8253 "this bug to the Tpetra developers." << endl;
8258 recvDataBufs[circBufInd].resize (recvNumRows * 2);
8260 std::ostringstream os;
8261 os << myRank <<
": Post receive-data receive from Process "
8262 << recvRank <<
": tag = " << dataTag <<
", buffer size = "
8263 << recvDataBufs[circBufInd].size () << endl;
8266 if (! recvDataReqs[circBufInd].is_null ()) {
8267 std::ostringstream os;
8268 os << myRank <<
": recvDataReqs[" << circBufInd <<
"] is not "
8269 <<
"null, for the receive-data receive from Process "
8270 << recvRank <<
"! This may mean that this process never "
8271 <<
"finished waiting for the receive from Process "
8272 << (recvRank - 3) <<
"." << endl;
8275 recvDataReqs[circBufInd] =
8276 ireceive<int, int_type> (recvDataBufs[circBufInd],
8277 recvRank, dataTag, comm);
8281 const int recvRank = p;
8282 const int circBufInd = recvRank % 3;
8284 std::ostringstream os;
8285 os << myRank <<
": Wait on receive-data receive from Process "
8286 << recvRank << endl;
8289 wait<int> (comm, outArg (recvDataReqs[circBufInd]));
8296 std::ostringstream os;
8297 os << myRank <<
": Write GIDs and PIDs from Process "
8298 << recvRank << endl;
8299 *err << os.str () << endl;
8301 const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
8302 if (debug && printNumRows == -1) {
8303 std::ostringstream os;
8304 os << myRank <<
": Result of receive-size receive from Process "
8305 << recvRank <<
" was -1. This should never happen, and "
8306 "suggests that its receive-size receive was never posted. "
8307 "Please report this bug to the Tpetra developers." << endl;
8310 if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
8311 std::ostringstream os;
8312 os << myRank <<
": Result of receive-size receive from Proc "
8313 << recvRank <<
" was " << printNumRows <<
" > 0, but "
8314 "recvDataBufs[" << circBufInd <<
"] is null. This should "
8315 "never happen. Please report this bug to the Tpetra "
8316 "developers." << endl;
8319 ArrayView<const int_type> printData = (recvDataBufs[circBufInd]) ();
8320 for (int_type k = 0; k < printNumRows; ++k) {
8321 const int_type gid = printData[2*k];
8322 const int_type pid = printData[2*k+1];
8323 out << gid << endl << pid << endl;
8326 else if (myRank == p) {
8329 std::ostringstream os;
8330 os << myRank <<
": Wait on my send-data send" << endl;
8333 wait<int> (comm, outArg (sendReqData));
8335 else if (myRank == p + 1) {
8338 std::ostringstream os;
8339 os << myRank <<
": Post send-data send: tag = " << dataTag
8343 sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
8346 std::ostringstream os;
8347 os << myRank <<
": Wait on my send-size send" << endl;
8350 wait<int> (comm, outArg (sendReqSize));
8352 else if (myRank == p + 2) {
8355 std::ostringstream os;
8356 os << myRank <<
": Post send-size send: size = "
8357 << sendDataSize[0] <<
", tag = " << sizeTag << endl;
8360 sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
8364 if (! err.is_null ()) {
8368 *err << myRank <<
": writeMap: Done" << endl;
8370 if (! err.is_null ()) {
8380 const int myRank = map.
getComm ()->getRank ();
8383 out.open (filename.c_str());
8416 printAsComment (std::ostream& out,
const std::string& str)
8419 std::istringstream inpstream (str);
8422 while (getline (inpstream, line)) {
8423 if (! line.empty()) {
8426 if (line[0] ==
'%') {
8427 out << line << endl;
8430 out <<
"%% " << line << endl;
8458 Teuchos::ParameterList pl;
8484 Teuchos::ParameterList pl;
8527 const Teuchos::ParameterList& params)
8530 std::string tmpFile =
"__TMP__" + fileName;
8531 const int myRank = A.
getDomainMap()->getComm()->getRank();
8532 bool precisionChanged=
false;
8542 if (std::ifstream(tmpFile))
8543 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
8544 "writeOperator: temporary file " << tmpFile <<
" already exists");
8545 out.open(tmpFile.c_str());
8546 if (params.isParameter(
"precision")) {
8547 oldPrecision = out.precision(params.get<
int>(
"precision"));
8548 precisionChanged=
true;
8552 const std::string header = writeOperatorImpl(out, A, params);
8555 if (precisionChanged)
8556 out.precision(oldPrecision);
8558 out.open(fileName.c_str(), std::ios::binary);
8559 bool printMatrixMarketHeader =
true;
8560 if (params.isParameter(
"print MatrixMarket header"))
8561 printMatrixMarketHeader = params.get<
bool>(
"print MatrixMarket header");
8562 if (printMatrixMarketHeader && myRank == 0) {
8567 std::ifstream src(tmpFile, std::ios_base::binary);
8571 remove(tmpFile.c_str());
8616 const Teuchos::ParameterList& params)
8618 const int myRank = A.
getDomainMap ()->getComm ()->getRank ();
8635 std::ostringstream tmpOut;
8637 if (params.isParameter (
"precision") && params.isType<
int> (
"precision")) {
8638 (void) tmpOut.precision (params.get<
int> (
"precision"));
8642 const std::string header = writeOperatorImpl (tmpOut, A, params);
8645 bool printMatrixMarketHeader =
true;
8646 if (params.isParameter (
"print MatrixMarket header") &&
8647 params.isType<
bool> (
"print MatrixMarket header")) {
8648 printMatrixMarketHeader = params.get<
bool> (
"print MatrixMarket header");
8650 if (printMatrixMarketHeader && myRank == 0) {
8666 out << tmpOut.str ();
8680 writeOperatorImpl (std::ostream& os,
8681 const operator_type& A,
8682 const Teuchos::ParameterList& params)
8686 using Teuchos::ArrayRCP;
8687 using Teuchos::Array;
8692 typedef Teuchos::OrdinalTraits<LO> TLOT;
8693 typedef Teuchos::OrdinalTraits<GO> TGOT;
8697 const map_type& domainMap = *(A.getDomainMap());
8698 RCP<const map_type> rangeMap = A.getRangeMap();
8699 RCP<const Teuchos::Comm<int> > comm = rangeMap->getComm();
8700 const int myRank = comm->getRank();
8701 const size_t numProcs = comm->getSize();
8704 if (params.isParameter(
"probing size"))
8705 numMVs = params.get<
int>(
"probing size");
8708 GO minColGid = domainMap.getMinAllGlobalIndex();
8709 GO maxColGid = domainMap.getMaxAllGlobalIndex();
8714 GO numGlobElts = maxColGid - minColGid + TGOT::one();
8715 GO numChunks = numGlobElts / numMVs;
8716 GO rem = numGlobElts % numMVs;
8717 GO indexBase = rangeMap->getIndexBase();
8719 int offsetToUseInPrinting = 1 - indexBase;
8720 if (params.isParameter(
"zero-based indexing")) {
8721 if (params.get<
bool>(
"zero-based indexing") ==
true)
8722 offsetToUseInPrinting = -indexBase;
8726 size_t numLocalRangeEntries = rangeMap->getLocalNumElements();
8729 RCP<const map_type> allGidsMap = rcp(
new map_type(TGOT::invalid(), numLocalRangeEntries,
8732 mv_type_go allGids(allGidsMap,1);
8733 Teuchos::ArrayRCP<GO> allGidsData = allGids.getDataNonConst(0);
8735 for (
size_t i=0; i<numLocalRangeEntries; i++)
8736 allGidsData[i] = rangeMap->getGlobalElement(i);
8737 allGidsData = Teuchos::null;
8740 GO numTargetMapEntries=TGOT::zero();
8741 Teuchos::Array<GO> importGidList;
8743 numTargetMapEntries = rangeMap->getGlobalNumElements();
8744 importGidList.reserve(numTargetMapEntries);
8745 for (GO j=0; j<numTargetMapEntries; ++j) importGidList.push_back(j + indexBase);
8747 importGidList.reserve(numTargetMapEntries);
8749 RCP<map_type> importGidMap = rcp(
new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
8752 import_type gidImporter(allGidsMap, importGidMap);
8753 mv_type_go importedGids(importGidMap, 1);
8754 importedGids.doImport(allGids, gidImporter,
INSERT);
8757 ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
8758 RCP<const map_type> importMap = rcp(
new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm) );
8761 import_type importer(rangeMap, importMap);
8763 RCP<mv_type> colsOnPid0 = rcp(
new mv_type(importMap,numMVs));
8765 RCP<mv_type> ei = rcp(
new mv_type(A.getDomainMap(),numMVs));
8766 RCP<mv_type> colsA = rcp(
new mv_type(A.getRangeMap(),numMVs));
8768 Array<GO> globalColsArray, localColsArray;
8769 globalColsArray.reserve(numMVs);
8770 localColsArray.reserve(numMVs);
8772 ArrayRCP<ArrayRCP<Scalar> > eiData(numMVs);
8773 for (
size_t i=0; i<numMVs; ++i)
8774 eiData[i] = ei->getDataNonConst(i);
8779 for (GO k=0; k<numChunks; ++k) {
8780 for (
size_t j=0; j<numMVs; ++j ) {
8782 GO curGlobalCol = minColGid + k*numMVs + j;
8783 globalColsArray.push_back(curGlobalCol);
8785 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8786 if (curLocalCol != TLOT::invalid()) {
8787 eiData[j][curLocalCol] = TGOT::one();
8788 localColsArray.push_back(curLocalCol);
8794 A.apply(*ei,*colsA);
8796 colsOnPid0->doImport(*colsA,importer,
INSERT);
8799 globalNnz += writeColumns(os,*colsOnPid0, numMVs, importedGidsData(),
8800 globalColsArray, offsetToUseInPrinting);
8803 for (
size_t j=0; j<numMVs; ++j ) {
8804 for (
int i=0; i<localColsArray.size(); ++i)
8805 eiData[j][localColsArray[i]] = TGOT::zero();
8807 globalColsArray.clear();
8808 localColsArray.clear();
8816 for (
int j=0; j<rem; ++j ) {
8817 GO curGlobalCol = maxColGid - rem + j + TGOT::one();
8818 globalColsArray.push_back(curGlobalCol);
8820 LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
8821 if (curLocalCol != TLOT::invalid()) {
8822 eiData[j][curLocalCol] = TGOT::one();
8823 localColsArray.push_back(curLocalCol);
8829 A.apply(*ei,*colsA);
8831 colsOnPid0->doImport(*colsA,importer,
INSERT);
8833 globalNnz += writeColumns(os,*colsOnPid0, rem, importedGidsData(),
8834 globalColsArray, offsetToUseInPrinting);
8837 for (
int j=0; j<rem; ++j ) {
8838 for (
int i=0; i<localColsArray.size(); ++i)
8839 eiData[j][localColsArray[i]] = TGOT::zero();
8841 globalColsArray.clear();
8842 localColsArray.clear();
8851 std::ostringstream oss;
8853 oss <<
"%%MatrixMarket matrix coordinate ";
8854 if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
8859 oss <<
" general" << std::endl;
8860 oss <<
"% Tpetra::Operator" << std::endl;
8861 std::time_t now = std::time(NULL);
8862 oss <<
"% time stamp: " << ctime(&now);
8863 oss <<
"% written from " << numProcs <<
" processes" << std::endl;
8864 size_t numRows = rangeMap->getGlobalNumElements();
8865 size_t numCols = domainMap.getGlobalNumElements();
8866 oss << numRows <<
" " << numCols <<
" " << globalNnz << std::endl;
8873 writeColumns(std::ostream& os, mv_type
const &colsA,
size_t const &numCols,
8874 Teuchos::ArrayView<const global_ordinal_type>
const &rowGids,
8875 Teuchos::Array<global_ordinal_type>
const &colsArray,
8880 typedef Teuchos::ScalarTraits<Scalar> STS;
8883 const Scalar zero = STS::zero();
8884 const size_t numRows = colsA.getGlobalLength();
8885 for (
size_t j=0; j<numCols; ++j) {
8886 Teuchos::ArrayRCP<const Scalar>
const curCol = colsA.getData(j);
8887 const GO J = colsArray[j];
8888 for (
size_t i=0; i<numRows; ++i) {
8889 const Scalar val = curCol[i];
8891 os << rowGids[i]+indexBase <<
" " << J+indexBase <<
" " << val << std::endl;
8914 const std::string& filename_suffix,
8916 const std::string& matrixName,
8917 const std::string& matrixDescription,
8918 const int ranksToWriteAtOnce=8,
8919 const bool debug=
false) {
8924 using STS =
typename Teuchos::ScalarTraits<ST>;
8928 Teuchos::RCP<const Teuchos::Comm<int> > comm = matrix.getComm ();
8929 TEUCHOS_TEST_FOR_EXCEPTION
8930 (comm.is_null (), std::invalid_argument,
8931 "The input matrix's communicator (Teuchos::Comm object) is null.");
8932 TEUCHOS_TEST_FOR_EXCEPTION
8933 (matrix.isGloballyIndexed() || !matrix.isFillComplete(), std::invalid_argument,
8934 "The input matrix must not be GloballyIndexed and must be fillComplete.");
8937 const int myRank = comm->getRank();
8938 const int numProc = comm->getSize();
8939 std::string filename = filename_prefix + std::to_string(myRank) + filename_suffix;
8940 RCP<const map_type> rowMap = matrix.getRowMap();
8941 RCP<const map_type> colMap = matrix.getColMap();
8942 size_t local_nnz = matrix.getLocalNumEntries();
8943 size_t local_num_rows = rowMap->getLocalNumElements();
8944 size_t local_num_cols = colMap->getLocalNumElements();
8945 const GO rowIndexBase = rowMap->getIndexBase();
8946 const GO colIndexBase = colMap->getIndexBase();
8949 int rank_limit = std::min(std::max(ranksToWriteAtOnce,1),numProc);
8952 for(
int base_rank = 0; base_rank < numProc; base_rank += rank_limit) {
8953 int stop = std::min(base_rank+rank_limit,numProc);
8955 if(base_rank <= myRank && myRank < stop) {
8957 std::ofstream out(filename);
8960 out <<
"%%MatrixMarket matrix coordinate "
8961 << (STS::isComplex ?
"complex" :
"real")
8962 <<
" general" << std::endl;
8965 if (matrixName !=
"") {
8966 printAsComment (out, matrixName);
8968 if (matrixDescription !=
"") {
8969 printAsComment (out, matrixDescription);
8975 out << local_num_rows <<
" " << local_num_cols <<
" " << local_nnz <<std::endl;
8982 Teuchos::SetScientific<ST> sci (out);
8984 for(
size_t l_row = 0; l_row < local_num_rows; l_row++) {
8985 GO g_row = rowMap->getGlobalElement(l_row);
8987 typename sparse_matrix_type::local_inds_host_view_type indices;
8988 typename sparse_matrix_type::values_host_view_type values;
8989 matrix.getLocalRowView(l_row, indices, values);
8990 for (
size_t ii = 0; ii < indices.extent(0); ii++) {
8991 const GO g_col = colMap->getGlobalElement(indices(ii));
8994 out << (g_row + 1 - rowIndexBase) <<
" "
8995 << (g_col + 1 - colIndexBase) <<
" ";
8996 if (STS::isComplex) {
8997 out << STS::real(values(ii)) <<
" " << STS::imag(values(ii));
From a distributed map build a map with all GIDs on the root node.
Declaration of a function that prints strings from each process.
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
Teuchos::RCP< const map_type > getColMap() const override
Returns the Map that describes the column distribution in this graph.
global_size_t getGlobalNumEntries() const override
Returns the global number of entries in the graph.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const override
Returns the communicator.
void fillComplete(const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Tell the graph that you are done changing its structure.
Teuchos::RCP< const map_type > getDomainMap() const override
Returns the Map associated with the domain of this graph.
void getGlobalRowView(const global_ordinal_type gblRow, global_inds_host_view_type &gblColInds) const override
Get a const view of the given global row's global column indices.
Teuchos::RCP< const map_type > getRangeMap() const override
Returns the Map associated with the domain of this graph.
Teuchos::RCP< const map_type > getRowMap() const override
Returns the Map that describes the row distribution in this graph.
void getLocalRowView(const LocalOrdinal lclRow, local_inds_host_view_type &lclColInds) const override
Get a const view of the given local row's local column indices.
bool isGloballyIndexed() const override
Whether the graph's column indices are stored as global indices.
void doImport(const SrcDistObject &source, const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, const CombineMode CM, const bool restrictedMode=false)
Import data into this object using an Import object ("forward mode").
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
A parallel distribution of indices over processes.
Teuchos::ArrayView< const global_ordinal_type > getLocalElementList() const
Return a NONOWNING view of the global indices owned by this process.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
global_ordinal_type getMinGlobalIndex() const
The minimum global index owned by the calling process.
global_size_t getGlobalNumElements() const
The number of elements in this Map.
bool isContiguous() const
True if this Map is distributed contiguously, else false.
size_t getLocalNumElements() const
The number of elements belonging to the calling process.
Matrix Market file reader for CrsMatrix and MultiVector.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
The MultiVector specialization associated with SparseMatrixType.
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.
static Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read Map (as a MultiVector) from the given input stream, with optional debugging output stream.
SparseMatrixType::global_ordinal_type global_ordinal_type
static Teuchos::RCP< const map_type > readMap(std::istream &in, const Teuchos::RCP< const comm_type > &comm, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read Map (as a MultiVector) from the given input stream.
static Teuchos::RCP< vector_type > readVector(std::istream &in, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream.
static Teuchos::RCP< sparse_matrix_type > readSparsePerRank(const std::string &filename_prefix, const std::string &filename_suffix, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const int ranksToReadAtOnce=8, const bool debug=false)
Read a Tpetra::CrsMatrix from a file per rank setup.
Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > vector_type
The Vector specialization associated with SparseMatrixType.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream, with provided Maps.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
SparseMatrixType::node_type node_type
The fourth template parameter of CrsMatrix and MultiVector.
static Teuchos::RCP< multivector_type > readDenseFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read dense matrix (as a MultiVector) from the given Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream.
static Teuchos::RCP< vector_type > readVectorFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read a Vector from the given Matrix Market file.
static Teuchos::RCP< const map_type > readMapFile(const std::string &filename, const Teuchos::RCP< const comm_type > &comm, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read Map (as a MultiVector) from the given Matrix Market file.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream, with provided Maps.
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file, with provided Maps.
SparseMatrixType::local_ordinal_type local_ordinal_type
SparseMatrixType sparse_matrix_type
This class' template parameter; a specialization of CrsMatrix.
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market file.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market input stream.
static Teuchos::RCP< multivector_type > readDense(std::istream &in, const Teuchos::RCP< const comm_type > &comm, Teuchos::RCP< const map_type > &map, const bool tolerant=false, const bool debug=false, const bool binary=false)
Read dense matrix (as a MultiVector) from the given Matrix Market input stream.
static Teuchos::RCP< sparse_graph_type > readSparseGraph(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market input stream.
static 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.
SparseMatrixType::scalar_type scalar_type
static Teuchos::RCP< sparse_graph_type > readSparseGraphFile(const std::string &filename, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse graph from the given Matrix Market file, with provided Maps.
CrsGraph< local_ordinal_type, global_ordinal_type, node_type > sparse_graph_type
The CrsGraph specialization associated with SparseMatrixType.
Matrix Market file writer for CrsMatrix and MultiVector.
static void writeMap(std::ostream &out, const map_type &map, const bool debug=false)
Print the Map to the given output stream.
static void writeDense(std::ostream &out, const Teuchos::RCP< const multivector_type > &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and or description.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Only for backwards compatibility; prefer the overload above.
CrsGraph< local_ordinal_type, global_ordinal_type, node_type > crs_graph_type
Specialization of Tpetra::CrsGraph that matches SparseMatrixType.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Only for backwards compatibility; prefer the overload above.
static void writeSparseGraphFile(const std::string &filename, const crs_graph_type &graph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), with no comments.
static void writeDense(std::ostream &out, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
static void writeSparsePerRank(const std::string &filename_prefix, const std::string &filename_suffix, const sparse_matrix_type &matrix, const std::string &matrixName, const std::string &matrixDescription, const int ranksToWriteAtOnce=8, const bool debug=false)
Write a Tpetra::CrsMatrix to a file per rank.
static void writeSparseGraphFile(const std::string &filename, const Teuchos::RCP< const crs_graph_type > &pGraph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), with no comments,...
SparseMatrixType sparse_matrix_type
Template parameter of this class; specialization of CrsMatrix.
static void writeSparse(std::ostream &out, const sparse_matrix_type &matrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
SparseMatrixType::scalar_type scalar_type
Type of the entries of the sparse matrix.
static void writeSparseGraphFile(const std::string &filename, const Teuchos::RCP< const crs_graph_type > &pGraph, const std::string &graphName, const std::string &graphDescription, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename), taking the graph by T...
static void 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.
static void writeDense(std::ostream &out, const Teuchos::RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
SparseMatrixType::global_ordinal_type global_ordinal_type
Type of indices as read from the Matrix Market file.
SparseMatrixType::node_type node_type
The Kokkos Node type; fourth template parameter of Tpetra::CrsMatrix.
static void writeDense(std::ostream &out, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static void writeSparse(std::ostream &out, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Only for backwards compatibility; prefer the overload above.
static void writeOperator(std::ostream &out, const operator_type &A)
Write a Tpetra::Operator to an output stream.
static void writeOperator(const std::string &fileName, operator_type const &A)
Write a Tpetra::Operator to a 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.
static void writeDenseFile(const std::string &filename, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
static void writeDenseFile(const std::string &filename, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static void writeSparseGraph(std::ostream &out, const crs_graph_type &graph, const bool debug=false)
Print the sparse graph in Matrix Market format to the given output stream, with no comments.
static void writeOperator(const std::string &fileName, const operator_type &A, const Teuchos::ParameterList ¶ms)
Write a Tpetra::Operator to a file, with options.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
Specialization of Tpetra::MultiVector that matches SparseMatrixType.
static void writeSparse(std::ostream &out, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Only for backwards compatibility; prefer the overload above.
static void writeMapFile(const std::string &filename, const map_type &map)
Write the Map to the given file.
static void writeDenseFile(const std::string &filename, const Teuchos::RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Specialization of Tpetra::Map that matches SparseMatrixType.
SparseMatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the sparse matrix.
static void writeSparseFile(const std::string &filename, const sparse_matrix_type &matrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
static void writeSparseGraph(std::ostream &out, const crs_graph_type &graph, const std::string &graphName, const std::string &graphDescription, const bool debug=false)
Print the sparse graph in Matrix Market format to the given output stream.
static void writeSparseGraphFile(const std::string &filename, const crs_graph_type &graph, const std::string &graphName, const std::string &graphDescription, const bool debug=false)
Print the sparse graph in Matrix Market format to the given file (by filename).
static void writeMap(std::ostream &out, const map_type &map, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool debug=false)
Print the Map to the given output stream out.
static void writeSparseFile(const std::string &filename, const sparse_matrix_type &matrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
static void writeSparse(std::ostream &out, const sparse_matrix_type &matrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
One or more distributed dense vectors.
Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getVector(const size_t j) const
Return a Vector which is a const view of column j.
size_t getNumVectors() const
Number of columns in the multivector.
Abstract interface for operators (e.g., matrices and preconditioners).
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....
A distributed dense vector.
Matrix Market file readers and writers for sparse and dense matrices (as CrsMatrix resp....
void gathervPrint(std::ostream &out, const std::string &s, const Teuchos::Comm< int > &comm)
On Process 0 in the given communicator, print strings from each process in that communicator,...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
size_t global_size_t
Global size_t object.
@ INSERT
Insert new values that don't currently exist.