48 #ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
54 #ifdef HAVE_XPETRA_EPETRA
56 #include "Epetra_MpiComm.h"
60 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
61 #include <EpetraExt_MatrixMatrix.h>
62 #include <EpetraExt_RowMatrixOut.h>
63 #include <EpetraExt_MultiVectorOut.h>
64 #include <EpetraExt_CrsMatrixIn.h>
65 #include <EpetraExt_MultiVectorIn.h>
66 #include <EpetraExt_BlockMapIn.h>
69 #include <EpetraExt_BlockMapOut.h>
72 #ifdef HAVE_XPETRA_TPETRA
73 #include <MatrixMarket_Tpetra.hpp>
74 #include <Tpetra_RowMatrixTransposer.hpp>
75 #include <TpetraExt_MatrixMatrix.hpp>
76 #include <Xpetra_TpetraMultiVector.hpp>
77 #include <Xpetra_TpetraCrsGraph.hpp>
78 #include <Xpetra_TpetraCrsMatrix.hpp>
79 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
80 #include "Tpetra_Util.hpp"
83 #ifdef HAVE_XPETRA_EPETRA
90 #include "Xpetra_CrsMatrixWrap.hpp"
93 #include "Xpetra_Map.hpp"
94 #include "Xpetra_StridedMap.hpp"
95 #include "Xpetra_StridedMapFactory.hpp"
96 #include "Xpetra_MapExtractor.hpp"
99 #include <Teuchos_MatrixMarket_Raw_Writer.hpp>
100 #include <Teuchos_MatrixMarket_Raw_Reader.hpp>
105 #ifdef HAVE_XPETRA_EPETRA
107 template <
class SC,
class LO,
class GO,
class NO>
108 RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO>>
111 "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
112 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
117 inline RCP<Xpetra::CrsMatrixWrap<double, int, int, Xpetra::EpetraNode>> Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double, int, int, Xpetra::EpetraNode>(RCP<Epetra_CrsMatrix>& epAB) {
130 template <
class SC,
class LO,
class GO,
class NO>
131 RCP<Xpetra::MultiVector<SC, LO, GO, NO>>
134 "Convert_Epetra_MultiVector_ToXpetra_MultiVector cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
135 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
140 inline RCP<Xpetra::MultiVector<double, int, int, Xpetra::EpetraNode>> Convert_Epetra_MultiVector_ToXpetra_MultiVector<double, int, int, Xpetra::EpetraNode>(RCP<Epetra_MultiVector>& epX) {
146 RCP<Xpetra::MultiVector<SC, LO, GO, NO>> tmp = Xpetra::toXpetra<GO, NO>(epX);
156 template <
class Scalar,
157 class LocalOrdinal = int,
158 class GlobalOrdinal = LocalOrdinal,
159 class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
162 #undef XPETRA_IO_SHORT
166 #ifdef HAVE_XPETRA_EPETRA
183 if (xeMap == Teuchos::null)
184 throw Exceptions::BadCast(
"Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
185 return xeMap->getEpetra_Map();
190 #ifdef HAVE_XPETRA_TPETRA
211 if (tmp_TMap == Teuchos::null)
212 throw Exceptions::BadCast(
"Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
213 return tmp_TMap->getTpetra_Map();
221 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tmp_Map = rcpFromRef(M);
222 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
224 if (tmp_EMap != Teuchos::null) {
225 int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
230 #endif // HAVE_XPETRA_EPETRAEXT
232 #ifdef HAVE_XPETRA_TPETRA
233 const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>& tmp_TMap =
235 if (tmp_TMap != Teuchos::null) {
236 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> TMap = tmp_TMap->
getTpetra_Map();
237 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeMapFile(fileName, *TMap);
240 #endif // HAVE_XPETRA_TPETRA
248 std::string mapfile =
"map_" + fileName;
251 RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_Vec = Teuchos::rcpFromRef(vec);
252 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
254 if (tmp_EVec != Teuchos::null) {
255 int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
260 #endif // HAVE_XPETRA_EPETRA
262 #ifdef HAVE_XPETRA_TPETRA
263 const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TVec =
265 if (tmp_TVec != Teuchos::null) {
266 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TVec = tmp_TVec->
getTpetra_MultiVector();
267 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
270 #endif // HAVE_XPETRA_TPETRA
272 throw Exceptions::BadCast(
"Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
288 RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_CrsMtx = crsOp.
getCrsMatrix();
289 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
291 if (tmp_ECrsMtx != Teuchos::null) {
293 int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
300 #ifdef HAVE_XPETRA_TPETRA
301 const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TCrsMtx =
303 if (tmp_TCrsMtx != Teuchos::null) {
304 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A = tmp_TCrsMtx->
getTpetra_CrsMatrix();
305 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeSparseFile(fileName, A);
308 const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_BlockCrs =
310 if (tmp_BlockCrs != Teuchos::null) {
311 std::ofstream outstream(fileName, std::ofstream::out);
312 Teuchos::FancyOStream ofs(Teuchos::rcpFromRef(outstream));
313 tmp_BlockCrs->getTpetra_BlockCrsMatrix()->describe(ofs, Teuchos::VERB_EXTREME);
317 #endif // HAVE_XPETRA_TPETRA
319 throw Exceptions::BadCast(
"Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
326 RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_CrsMtx = crsOp.
getCrsMatrix();
328 ArrayRCP<const size_t> rowptr_RCP;
329 ArrayRCP<LocalOrdinal> rowptr2_RCP;
330 ArrayRCP<const LocalOrdinal> colind_RCP;
331 ArrayRCP<const Scalar> vals_RCP;
332 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
334 ArrayView<const size_t> rowptr = rowptr_RCP();
335 ArrayView<const LocalOrdinal> colind = colind_RCP();
336 ArrayView<const Scalar> vals = vals_RCP();
338 rowptr2_RCP.resize(rowptr.size());
339 ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
340 for (
size_t j = 0; j < rowptr.size(); j++)
341 rowptr2[j] = rowptr[j];
343 Teuchos::MatrixMarket::Raw::Writer<Scalar, LocalOrdinal> writer;
344 writer.writeFile(fileName +
"." + std::to_string(Op.
getRowMap()->getComm()->getSize()) +
"." + std::to_string(Op.
getRowMap()->getComm()->getRank()),
345 rowptr2, colind, vals,
346 rowptr.size() - 1, Op.
getColMap()->getLocalNumElements());
367 for (
size_t row = 0; row < Op.
Rows(); ++row) {
368 for (
size_t col = 0; col < Op.
Cols(); ++col) {
369 RCP<const Matrix> m = Op.
getMatrix(row, col);
370 if (m != Teuchos::null) {
372 "Sub block matrix (" << row <<
"," << col <<
") is not of type CrsMatrixWrap.");
373 XpIO::Write(fileName +
toString(row) +
toString(col) +
".m", *m, writeAllMaps);
382 for (
size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
383 RCP<const Map> map = rangeMapExtractor->getMap(row);
384 XpIO::Write(
"subRangeMap_" + fileName + XpIO::toString<size_t>(row) +
".m", *map);
386 XpIO::Write(
"fullRangeMap_" + fileName +
".m", *(rangeMapExtractor->getFullMap()));
388 for (
size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
389 RCP<const Map> map = domainMapExtractor->getMap(col);
390 XpIO::Write(
"subDomainMap_" + fileName + XpIO::toString<size_t>(col) +
".m", *map);
392 XpIO::Write(
"fullDomainMap_" + fileName +
".m", *(domainMapExtractor->getFullMap()));
396 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
Read(
const std::string& fileName,
Xpetra::UnderlyingLib lib,
const RCP<
const Teuchos::Comm<int>>& comm,
bool binary =
false) {
397 if (binary ==
false) {
400 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
401 Epetra_CrsMatrix* eA;
403 int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
407 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
409 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A =
410 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
416 #ifdef HAVE_XPETRA_TPETRA
417 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
419 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
421 bool callFillComplete =
true;
423 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
441 std::ifstream ifs(fileName.c_str(), std::ios::binary);
444 ifs.read(reinterpret_cast<char*>(&m),
sizeof(m));
445 ifs.read(reinterpret_cast<char*>(&n),
sizeof(n));
446 ifs.read(reinterpret_cast<char*>(&nnz),
sizeof(nnz));
448 int myRank = comm->getRank();
453 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A;
456 Teuchos::Array<GlobalOrdinal> inds;
457 Teuchos::Array<Scalar> vals;
459 Teuchos::ArrayRCP<size_t> numEntriesPerRow(m, (
size_t)(0));
460 for (
int i = 0; i < m; i++) {
462 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
463 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
464 numEntriesPerRow[row] = rownnz;
465 for (
int j = 0; j < rownnz; j++) {
467 ifs.read(reinterpret_cast<char*>(&index),
sizeof(index));
469 for (
int j = 0; j < rownnz; j++) {
471 ifs.read(reinterpret_cast<char*>(&value),
sizeof(value));
478 ifs.seekg(0, ifs.beg);
480 ifs.read(reinterpret_cast<char*>(&junk),
sizeof(junk));
481 ifs.read(reinterpret_cast<char*>(&junk),
sizeof(junk));
482 ifs.read(reinterpret_cast<char*>(&junk),
sizeof(junk));
483 for (
int i = 0; i < m; i++) {
485 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
486 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
489 for (
int j = 0; j < rownnz; j++) {
491 ifs.read(reinterpret_cast<char*>(&index),
sizeof(index));
492 inds[j] = Teuchos::as<GlobalOrdinal>(index);
494 for (
int j = 0; j < rownnz; j++) {
496 ifs.read(reinterpret_cast<char*>(&value),
sizeof(value));
497 vals[j] = Teuchos::as<Scalar>(value);
499 A->insertGlobalValues(row, inds, vals);
503 Teuchos::ArrayRCP<size_t> numEntriesPerRow(0, (
size_t)(0));
507 A->fillComplete(domainMap, rangeMap);
512 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
520 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
521 Read(
const std::string& filename,
526 const bool callFillComplete =
true,
527 const bool binary =
false,
528 const bool tolerant =
false,
529 const bool debug =
false) {
532 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
533 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
536 if (binary ==
false) {
538 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
539 Epetra_CrsMatrix* eA;
545 if (colMap.is_null()) {
546 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
550 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
556 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
557 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A =
558 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
565 #ifdef HAVE_XPETRA_TPETRA
566 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
567 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
568 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
570 const RCP<const map_type> tpetraRowMap =
Map2TpetraMap(*rowMap);
571 RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null :
Map2TpetraMap(*colMap));
572 const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap :
Map2TpetraMap(*rangeMap));
573 const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap :
Map2TpetraMap(*domainMap));
575 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
576 callFillComplete, tolerant, debug);
593 auto tempA =
Read(filename, lib, rowMap->getComm(), binary);
598 if (callFillComplete)
599 A->fillComplete(domainMap, rangeMap);
604 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
615 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
ReadLocal(
const std::string& filename,
620 const bool callFillComplete =
true,
621 const bool binary =
false,
622 const bool tolerant =
false,
623 const bool debug =
false) {
624 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(),
Exceptions::RuntimeError,
"Utils::ReadLocal() : rowMap cannot be null");
625 TEUCHOS_TEST_FOR_EXCEPTION(colMap.is_null(),
Exceptions::RuntimeError,
"Utils::ReadLocal() : colMap cannot be null");
631 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
632 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
634 std::string rankFilename = filename +
"." + std::to_string(rowMap->getComm()->getSize()) +
"." + std::to_string(rowMap->getComm()->getRank());
635 RCP<matrix_type> A = rcp(
new crs_wrap_type(rowMap, colMap, 0));
637 if (binary ==
false) {
638 RCP<Teuchos::ParameterList> params = rcp(
new Teuchos::ParameterList());
639 params->set(
"Parse tolerantly", tolerant);
640 params->set(
"Debug mode", debug);
642 LocalOrdinal numRows = rowMap->getLocalNumElements();
643 LocalOrdinal numCols = colMap->getLocalNumElements();
645 ArrayRCP<LocalOrdinal> rowptr2_RCP;
646 ArrayRCP<LocalOrdinal> colind2_RCP;
647 ArrayRCP<Scalar> vals2_RCP;
649 Teuchos::MatrixMarket::Raw::Reader<Scalar, LocalOrdinal> reader;
650 reader.readFile(rowptr2_RCP, colind2_RCP, vals2_RCP,
654 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
656 ArrayRCP<size_t> rowptr_RCP;
657 ArrayRCP<LocalOrdinal> colind_RCP;
658 ArrayRCP<Scalar> vals_RCP;
659 ACrs->allocateAllValues(colind2_RCP.size(), rowptr_RCP, colind_RCP, vals_RCP);
661 rowptr_RCP.assign(rowptr2_RCP.begin(), rowptr2_RCP.end());
662 colind_RCP = colind2_RCP;
663 vals_RCP = vals2_RCP;
665 ACrs->setAllValues(rowptr_RCP, colind_RCP, vals_RCP);
668 std::ifstream ifs = std::ifstream(rankFilename.c_str(), std::ios::binary);
672 ifs.read(reinterpret_cast<char*>(&m),
sizeof(m));
673 ifs.read(reinterpret_cast<char*>(&n),
sizeof(n));
674 ifs.read(reinterpret_cast<char*>(&nnz),
sizeof(nnz));
676 TEUCHOS_ASSERT_EQUALITY(Teuchos::as<int>(rowMap->getLocalNumElements()), m);
678 Teuchos::ArrayRCP<size_t> rowptrRCP;
679 Teuchos::ArrayRCP<LocalOrdinal> indicesRCP;
680 Teuchos::ArrayRCP<Scalar> valuesRCP;
682 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
684 ACrs->allocateAllValues(nnz, rowptrRCP, indicesRCP, valuesRCP);
686 Teuchos::ArrayView<size_t> rowptr = rowptrRCP();
687 Teuchos::ArrayView<LocalOrdinal> indices = indicesRCP();
688 Teuchos::ArrayView<Scalar> values = valuesRCP();
693 for (
int i = 0; i < m; i++) {
695 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
696 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
698 rowptr[row + 1] += rownnz;
699 ifs.seekg(
sizeof(
int) * rownnz +
sizeof(
double) * rownnz, ifs.cur);
701 for (
int i = 0; i < m; i++)
702 rowptr[i + 1] += rowptr[i];
703 TEUCHOS_ASSERT(Teuchos::as<int>(rowptr[m]) == nnz);
706 ifs.seekg(
sizeof(
int) * 3, ifs.beg);
709 for (
int i = 0; i < m; i++) {
711 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
712 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
713 size_t ptr = rowptr[row];
714 for (
int j = 0; j < rownnz; j++) {
716 ifs.read(reinterpret_cast<char*>(&index),
sizeof(index));
717 indices[ptr] = Teuchos::as<LocalOrdinal>(index);
719 sorted = sorted & (indices[ptr - 1] < indices[ptr]);
723 for (
int j = 0; j < rownnz; j++) {
725 ifs.read(reinterpret_cast<char*>(&value),
sizeof(value));
726 values[ptr] = Teuchos::as<Scalar>(value);
729 rowptr[row] += rownnz;
731 for (
int i = m; i > 0; i--)
732 rowptr[i] = rowptr[i - 1];
735 #ifdef HAVE_XPETRA_TPETRA
737 for (LocalOrdinal lclRow = 0; lclRow < m; lclRow++) {
738 size_t rowBegin = rowptr[lclRow];
739 size_t rowEnd = rowptr[lclRow + 1];
740 Tpetra::sort2(&indices[rowBegin], &indices[rowEnd], &values[rowBegin]);
744 TEUCHOS_ASSERT(sorted);
747 ACrs->setAllValues(rowptrRCP, indicesRCP, valuesRCP);
750 if (callFillComplete)
751 A->fillComplete(domainMap, rangeMap);
757 const RCP<const Map>& map,
758 const bool binary =
false) {
765 #ifdef HAVE_XPETRA_TPETRA
766 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
767 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
768 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
769 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
771 RCP<const map_type> temp =
toTpetra(map);
772 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp,
false,
false, binary);
782 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
785 static RCP<const Map>
ReadMap(
const std::string& fileName,
787 const RCP<
const Teuchos::Comm<int>>& comm,
788 const bool binary =
false) {
792 #ifdef HAVE_XPETRA_TPETRA
793 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
794 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
796 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tMap = reader_type::readMapFile(fileName, comm,
false,
false, binary);
808 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
827 size_t numBlocks = 2;
829 std::vector<RCP<const Map>> rangeMapVec;
830 for (
size_t row = 0; row < numBlocks; ++row) {
831 RCP<const Map> map = XpIO::ReadMap(
"subRangeMap_" + fileName + XpIO::toString<size_t>(row) +
".m", lib, comm);
832 rangeMapVec.push_back(map);
834 RCP<const Map> fullRangeMap = XpIO::ReadMap(
"fullRangeMap_" + fileName +
".m", lib, comm);
836 std::vector<RCP<const Map>> domainMapVec;
837 for (
size_t col = 0; col < numBlocks; ++col) {
838 RCP<const Map> map = XpIO::ReadMap(
"subDomainMap_" + fileName + XpIO::toString<size_t>(col) +
".m", lib, comm);
839 domainMapVec.push_back(map);
841 RCP<const Map> fullDomainMap = XpIO::ReadMap(
"fullDomainMap_" + fileName +
".m", lib, comm);
857 bool bRangeUseThyraStyleNumbering =
false;
864 RCP<const MapExtractor> rangeMapExtractor =
865 Teuchos::rcp(
new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
868 bool bDomainUseThyraStyleNumbering =
false;
875 RCP<const MapExtractor> domainMapExtractor =
876 Teuchos::rcp(
new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
878 RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 33));
881 for (
size_t row = 0; row < numBlocks; ++row) {
882 for (
size_t col = 0; col < numBlocks; ++col) {
883 RCP<const Map> rowSubMap = XpIO::ReadMap(
"rowmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
884 RCP<const Map> colSubMap = XpIO::ReadMap(
"colmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
885 RCP<const Map> domSubMap = XpIO::ReadMap(
"domainmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
886 RCP<const Map> ranSubMap = XpIO::ReadMap(
"rangemap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
887 RCP<Matrix> mat = XpIO::Read(fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
888 bOp->setMatrix(row, col, mat);
900 std::ostringstream buf;
906 #ifdef HAVE_XPETRA_EPETRA
916 template <
class Scalar>
923 #ifdef HAVE_XPETRA_EPETRA
928 if (xeMap == Teuchos::null)
929 throw Exceptions::BadCast(
"IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
930 return xeMap->getEpetra_Map();
935 #ifdef HAVE_XPETRA_TPETRA
940 if (tmp_TMap == Teuchos::null)
941 throw Exceptions::BadCast(
"IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
942 return tmp_TMap->getTpetra_Map();
950 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tmp_Map = rcpFromRef(M);
951 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
953 if (tmp_EMap != Teuchos::null) {
954 int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
959 #endif // HAVE_XPETRA_EPETRA
961 #ifdef HAVE_XPETRA_TPETRA
962 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
963 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
966 const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>& tmp_TMap =
968 if (tmp_TMap != Teuchos::null) {
969 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> TMap = tmp_TMap->
getTpetra_Map();
970 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeMapFile(fileName, *TMap);
974 #endif // HAVE_XPETRA_TPETRA
980 std::string mapfile =
"map_" + fileName;
983 RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_Vec = Teuchos::rcpFromRef(vec);
984 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
986 if (tmp_EVec != Teuchos::null) {
987 int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
992 #endif // HAVE_XPETRA_EPETRAEXT
994 #ifdef HAVE_XPETRA_TPETRA
995 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
996 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
999 const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TVec =
1001 if (tmp_TVec != Teuchos::null) {
1002 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TVec = tmp_TVec->
getTpetra_MultiVector();
1003 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
1007 #endif // HAVE_XPETRA_TPETRA
1009 throw Exceptions::BadCast(
"Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
1025 RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_CrsMtx = crsOp.
getCrsMatrix();
1026 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1028 if (tmp_ECrsMtx != Teuchos::null) {
1030 int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
1035 #endif // endif HAVE_XPETRA_EPETRA
1037 #ifdef HAVE_XPETRA_TPETRA
1038 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1039 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1042 const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TCrsMtx =
1044 if (tmp_TCrsMtx != Teuchos::null) {
1045 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A = tmp_TCrsMtx->
getTpetra_CrsMatrix();
1046 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeSparseFile(fileName, A);
1049 const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_BlockCrs =
1051 if (tmp_BlockCrs != Teuchos::null) {
1052 std::ofstream outstream(fileName, std::ofstream::out);
1053 Teuchos::FancyOStream ofs(Teuchos::rcpFromRef(outstream));
1054 tmp_BlockCrs->getTpetra_BlockCrsMatrix()->describe(ofs, Teuchos::VERB_EXTREME);
1059 #endif // HAVE_XPETRA_TPETRA
1061 throw Exceptions::BadCast(
"Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
1074 RCP<const Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>> tmp_Graph = rcpFromRef(graph);
1076 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1078 if (tmp_ECrsGraph != Teuchos::null) {
1081 #endif // endif HAVE_XPETRA_EPETRA
1083 #ifdef HAVE_XPETRA_TPETRA
1084 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1085 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1088 RCP<const Xpetra::TpetraCrsGraph<LocalOrdinal, GlobalOrdinal, Node>> tmp_TCrsGraph =
1090 if (tmp_TCrsGraph != Teuchos::null) {
1091 RCP<const Tpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>> G = tmp_TCrsGraph->
getTpetra_CrsGraph();
1092 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeSparseGraphFile(fileName, G);
1096 #endif // HAVE_XPETRA_TPETRA
1098 throw Exceptions::BadCast(
"Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
1105 RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_CrsMtx = crsOp.
getCrsMatrix();
1107 ArrayRCP<const size_t> rowptr_RCP;
1108 ArrayRCP<LocalOrdinal> rowptr2_RCP;
1109 ArrayRCP<const LocalOrdinal> colind_RCP;
1110 ArrayRCP<const Scalar> vals_RCP;
1111 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
1113 ArrayView<const size_t> rowptr = rowptr_RCP();
1114 ArrayView<const LocalOrdinal> colind = colind_RCP();
1115 ArrayView<const Scalar> vals = vals_RCP();
1117 rowptr2_RCP.resize(rowptr.size());
1118 ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
1119 for (
size_t j = 0; j < Teuchos::as<size_t>(rowptr.size()); j++)
1120 rowptr2[j] = rowptr[j];
1122 Teuchos::MatrixMarket::Raw::Writer<Scalar, LocalOrdinal> writer;
1123 writer.writeFile(fileName +
"." + std::to_string(Op.
getRowMap()->getComm()->getSize()) +
"." + std::to_string(Op.
getRowMap()->getComm()->getRank()),
1124 rowptr2, colind, vals,
1125 rowptr.size() - 1, Op.
getColMap()->getLocalNumElements());
1149 for (
size_t row = 0; row < Op.
Rows(); ++row) {
1150 for (
size_t col = 0; col < Op.
Cols(); ++col) {
1151 RCP<const Matrix> m = Op.
getMatrix(row, col);
1152 if (m != Teuchos::null) {
1154 "Sub block matrix (" << row <<
"," << col <<
") is not of type CrsMatrixWrap.");
1155 XpIO::Write(fileName +
toString(row) +
toString(col) +
".m", *m, writeAllMaps);
1164 for (
size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
1165 RCP<const Map> map = rangeMapExtractor->getMap(row);
1166 XpIO::Write(
"subRangeMap_" + fileName + XpIO::toString<size_t>(row) +
".m", *map);
1168 XpIO::Write(
"fullRangeMap_" + fileName +
".m", *(rangeMapExtractor->getFullMap()));
1170 for (
size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
1171 RCP<const Map> map = domainMapExtractor->getMap(col);
1172 XpIO::Write(
"subDomainMap_" + fileName + XpIO::toString<size_t>(col) +
".m", *map);
1174 XpIO::Write(
"fullDomainMap_" + fileName +
".m", *(domainMapExtractor->getFullMap()));
1178 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
Read(
const std::string& fileName,
Xpetra::UnderlyingLib lib,
const RCP<
const Teuchos::Comm<int>>& comm,
bool binary =
false) {
1179 if (binary ==
false) {
1182 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1183 Epetra_CrsMatrix* eA;
1185 int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
1189 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1191 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A =
1192 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1198 #ifdef HAVE_XPETRA_TPETRA
1199 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1200 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1203 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1205 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1207 bool callFillComplete =
true;
1209 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
1228 std::ifstream ifs(fileName.c_str(), std::ios::binary);
1231 ifs.read(reinterpret_cast<char*>(&m),
sizeof(m));
1232 ifs.read(reinterpret_cast<char*>(&n),
sizeof(n));
1233 ifs.read(reinterpret_cast<char*>(&nnz),
sizeof(nnz));
1235 int myRank = comm->getRank();
1241 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A;
1244 Teuchos::Array<GlobalOrdinal> inds;
1245 Teuchos::Array<Scalar> vals;
1247 Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
1248 for (
int i = 0; i < m; i++) {
1250 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
1251 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
1252 numEntriesPerRow[i] = rownnz;
1253 for (
int j = 0; j < rownnz; j++) {
1255 ifs.read(reinterpret_cast<char*>(&index),
sizeof(index));
1257 for (
int j = 0; j < rownnz; j++) {
1259 ifs.read(reinterpret_cast<char*>(&value),
sizeof(value));
1266 ifs.seekg(0, ifs.beg);
1268 ifs.read(reinterpret_cast<char*>(&m),
sizeof(junk));
1269 ifs.read(reinterpret_cast<char*>(&n),
sizeof(junk));
1270 ifs.read(reinterpret_cast<char*>(&nnz),
sizeof(junk));
1271 for (
int i = 0; i < m; i++) {
1273 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
1274 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
1275 inds.resize(rownnz);
1276 vals.resize(rownnz);
1277 for (
int j = 0; j < rownnz; j++) {
1279 ifs.read(reinterpret_cast<char*>(&index),
sizeof(index));
1280 inds[j] = Teuchos::as<GlobalOrdinal>(index);
1282 for (
int j = 0; j < rownnz; j++) {
1284 ifs.read(reinterpret_cast<char*>(&value),
sizeof(value));
1285 vals[j] = Teuchos::as<Scalar>(value);
1287 A->insertGlobalValues(row, inds, vals);
1291 Teuchos::ArrayRCP<size_t> numEntriesPerRow(0);
1295 A->fillComplete(domainMap, rangeMap);
1300 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1308 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
Read(
const std::string& filename,
1313 const bool callFillComplete =
true,
1314 const bool binary =
false,
1315 const bool tolerant =
false,
1316 const bool debug =
false) {
1319 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
1320 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
1323 if (binary ==
false) {
1325 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1326 Epetra_CrsMatrix* eA;
1332 if (colMap.is_null()) {
1333 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
1337 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
1343 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1344 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A =
1345 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1352 #ifdef HAVE_XPETRA_TPETRA
1353 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1354 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1357 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1358 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1359 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1361 const RCP<const map_type> tpetraRowMap =
Map2TpetraMap(*rowMap);
1362 RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null :
Map2TpetraMap(*colMap));
1363 const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap :
Map2TpetraMap(*rangeMap));
1364 const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap :
Map2TpetraMap(*domainMap));
1366 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
1367 callFillComplete, tolerant, debug);
1385 auto tempA =
Read(filename, lib, rowMap->getComm(), binary);
1390 if (callFillComplete)
1391 A->fillComplete(domainMap, rangeMap);
1396 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1407 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
ReadLocal(
const std::string& filename,
1412 const bool callFillComplete =
true,
1413 const bool binary =
false,
1414 const bool tolerant =
false,
1415 const bool debug =
false) {
1416 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(),
Exceptions::RuntimeError,
"Utils::ReadLocal() : rowMap cannot be null");
1417 TEUCHOS_TEST_FOR_EXCEPTION(colMap.is_null(),
Exceptions::RuntimeError,
"Utils::ReadLocal() : colMap cannot be null");
1423 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
1424 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
1426 std::string rankFilename = filename +
"." + std::to_string(rowMap->getComm()->getSize()) +
"." + std::to_string(rowMap->getComm()->getRank());
1427 RCP<matrix_type> A = rcp(
new crs_wrap_type(rowMap, colMap, 0));
1429 if (binary ==
false) {
1430 RCP<Teuchos::ParameterList> params = rcp(
new Teuchos::ParameterList());
1431 params->set(
"Parse tolerantly", tolerant);
1432 params->set(
"Debug mode", debug);
1437 ArrayRCP<LocalOrdinal> rowptr2_RCP;
1438 ArrayRCP<LocalOrdinal> colind2_RCP;
1439 ArrayRCP<Scalar> vals2_RCP;
1441 Teuchos::MatrixMarket::Raw::Reader<Scalar, LocalOrdinal> reader;
1442 reader.readFile(rowptr2_RCP, colind2_RCP, vals2_RCP,
1446 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
1448 ArrayRCP<size_t> rowptr_RCP;
1449 ArrayRCP<LocalOrdinal> colind_RCP;
1450 ArrayRCP<Scalar> vals_RCP;
1451 ACrs->allocateAllValues(colind2_RCP.size(), rowptr_RCP, colind_RCP, vals_RCP);
1453 rowptr_RCP.assign(rowptr2_RCP.begin(), rowptr2_RCP.end());
1454 colind_RCP = colind2_RCP;
1455 vals_RCP = vals2_RCP;
1457 ACrs->setAllValues(rowptr_RCP, colind_RCP, vals_RCP);
1460 std::ifstream ifs = std::ifstream(rankFilename.c_str(), std::ios::binary);
1464 ifs.read(reinterpret_cast<char*>(&m),
sizeof(m));
1465 ifs.read(reinterpret_cast<char*>(&n),
sizeof(n));
1466 ifs.read(reinterpret_cast<char*>(&nnz),
sizeof(nnz));
1468 TEUCHOS_ASSERT_EQUALITY(Teuchos::as<int>(rowMap->getLocalNumElements()), m);
1470 Teuchos::ArrayRCP<size_t> rowptrRCP;
1471 Teuchos::ArrayRCP<LocalOrdinal> indicesRCP;
1472 Teuchos::ArrayRCP<Scalar> valuesRCP;
1474 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
1476 ACrs->allocateAllValues(nnz, rowptrRCP, indicesRCP, valuesRCP);
1478 Teuchos::ArrayView<size_t> rowptr = rowptrRCP();
1479 Teuchos::ArrayView<LocalOrdinal> indices = indicesRCP();
1480 Teuchos::ArrayView<Scalar> values = valuesRCP();
1485 for (
int i = 0; i < m; i++) {
1487 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
1488 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
1490 rowptr[row + 1] += rownnz;
1491 ifs.seekg(
sizeof(
int) * rownnz +
sizeof(
double) * rownnz, ifs.cur);
1493 for (
int i = 0; i < m; i++)
1494 rowptr[i + 1] += rowptr[i];
1495 TEUCHOS_ASSERT(Teuchos::as<int>(rowptr[m]) == nnz);
1498 ifs.seekg(
sizeof(
int) * 3, ifs.beg);
1501 for (
int i = 0; i < m; i++) {
1503 ifs.read(reinterpret_cast<char*>(&row),
sizeof(row));
1504 ifs.read(reinterpret_cast<char*>(&rownnz),
sizeof(rownnz));
1505 size_t ptr = rowptr[row];
1506 for (
int j = 0; j < rownnz; j++) {
1508 ifs.read(reinterpret_cast<char*>(&index),
sizeof(index));
1509 indices[ptr] = Teuchos::as<LocalOrdinal>(index);
1511 sorted = sorted & (indices[ptr - 1] < indices[ptr]);
1515 for (
int j = 0; j < rownnz; j++) {
1517 ifs.read(reinterpret_cast<char*>(&value),
sizeof(value));
1518 values[ptr] = Teuchos::as<Scalar>(value);
1521 rowptr[row] += rownnz;
1523 for (
int i = m; i > 0; i--)
1524 rowptr[i] = rowptr[i - 1];
1527 #ifdef HAVE_XPETRA_TPETRA
1530 size_t rowBegin = rowptr[lclRow];
1531 size_t rowEnd = rowptr[lclRow + 1];
1532 Tpetra::sort2(&indices[rowBegin], &indices[rowEnd], &values[rowBegin]);
1536 TEUCHOS_ASSERT(sorted);
1539 ACrs->setAllValues(rowptrRCP, indicesRCP, valuesRCP);
1542 if (callFillComplete)
1543 A->fillComplete(domainMap, rangeMap);
1548 static RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
ReadMultiVector(
const std::string& fileName,
1550 const bool binary =
false) {
1556 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1557 TEUCHOS_ASSERT(!binary);
1558 Epetra_MultiVector* MV;
1559 int rv = EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(),
toEpetra(map), MV);
1561 RCP<Epetra_MultiVector> MVrcp = rcp(MV);
1562 return Convert_Epetra_MultiVector_ToXpetra_MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(MVrcp);
1567 #ifdef HAVE_XPETRA_TPETRA
1568 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1569 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1572 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1573 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1574 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1575 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
1577 RCP<const map_type> temp =
toTpetra(map);
1578 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp,
false,
false, binary);
1579 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> rmv =
Xpetra::toXpetra(TMV);
1589 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1592 static RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
ReadMap(
const std::string& fileName,
1594 const RCP<
const Teuchos::Comm<int>>& comm,
1595 const bool binary =
false) {
1599 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1600 TEUCHOS_ASSERT(!binary);
1602 int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(
Xpetra::toEpetra(comm)), eMap);
1606 RCP<Epetra_Map> eMap1 = rcp(
new Epetra_Map(*eMap));
1607 return Xpetra::toXpetra<int, Node>(*eMap1);
1612 #ifdef HAVE_XPETRA_TPETRA
1613 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1614 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1617 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1618 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1620 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tMap = reader_type::readMapFile(fileName, comm,
false,
false, binary);
1633 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1655 size_t numBlocks = 2;
1657 std::vector<RCP<const Map>> rangeMapVec;
1658 for (
size_t row = 0; row < numBlocks; ++row) {
1659 RCP<const Map> map = XpIO::ReadMap(
"subRangeMap_" + fileName + XpIO::toString<size_t>(row) +
".m", lib, comm);
1660 rangeMapVec.push_back(map);
1662 RCP<const Map> fullRangeMap = XpIO::ReadMap(
"fullRangeMap_" + fileName +
".m", lib, comm);
1664 std::vector<RCP<const Map>> domainMapVec;
1665 for (
size_t col = 0; col < numBlocks; ++col) {
1666 RCP<const Map> map = XpIO::ReadMap(
"subDomainMap_" + fileName + XpIO::toString<size_t>(col) +
".m", lib, comm);
1667 domainMapVec.push_back(map);
1669 RCP<const Map> fullDomainMap = XpIO::ReadMap(
"fullDomainMap_" + fileName +
".m", lib, comm);
1685 bool bRangeUseThyraStyleNumbering =
false;
1692 RCP<const MapExtractor> rangeMapExtractor =
1693 rcp(
new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
1696 bool bDomainUseThyraStyleNumbering =
false;
1702 RCP<const MapExtractor> domainMapExtractor =
1703 rcp(
new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
1705 RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 33));
1708 for (
size_t row = 0; row < numBlocks; ++row) {
1709 for (
size_t col = 0; col < numBlocks; ++col) {
1710 RCP<const Map> rowSubMap = XpIO::ReadMap(
"rowmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
1711 RCP<const Map> colSubMap = XpIO::ReadMap(
"colmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
1712 RCP<const Map> domSubMap = XpIO::ReadMap(
"domainmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
1713 RCP<const Map> ranSubMap = XpIO::ReadMap(
"rangemap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", lib, comm);
1714 RCP<Matrix> mat = XpIO::Read(fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) +
".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1715 bOp->setMatrix(row, col, mat);
1719 bOp->fillComplete();
1727 std::ostringstream buf;
1732 #endif // HAVE_XPETRA_EPETRA
1736 #define XPETRA_IO_SHORT
RCP< CrsMatrix > getCrsMatrix() const
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Read(const std::string &filename, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rowMap, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> colMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> domainMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rangeMap=Teuchos::null, const bool callFillComplete=true, const bool binary=false, const bool tolerant=false, const bool debug=false)
Read matrix from file in Matrix Market or binary format.
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadLocal(const std::string &filename, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rowMap, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> colMap, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> domainMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rangeMap=Teuchos::null, const bool callFillComplete=true, const bool binary=false, const bool tolerant=false, const bool debug=false)
Read matrix from local files in Matrix Market or binary format.
static void Write(const std::string &fileName, const Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, const bool &writeAllMaps=false)
Save CrsGraph to file in Matrix Market format.
static RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadBlockedCrsMatrix(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int >> &comm)
Read block matrix from one file per block in Matrix Market format.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const =0
Returns the Map that describes the column distribution in this graph.
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map, const bool binary=false)
static void WriteLocal(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save local parts of matrix to files in Matrix Market format.
Exception throws to report errors in the internal logical of the program.
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsGraph() const
Get the underlying Tpetra graph.
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
Exception indicating invalid cast attempted.
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadLocal(const std::string &filename, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rowMap, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> colMap, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> domainMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rangeMap=Teuchos::null, const bool callFillComplete=true, const bool binary=false, const bool tolerant=false, const bool debug=false)
Read matrix from local files in Matrix Market or binary format.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual size_t Cols() const
number of column blocks
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
static RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadBlockedCrsMatrix(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int >> &comm)
Read block matrix from one file per block in Matrix Market format.
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save block matrix to one file per block in Matrix Market format.
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
virtual const Teuchos::RCP< const map_type > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
Returns the Map associated with the domain of this graph.
static RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int >> &comm, const bool binary=false)
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save block matrix to one file per block in Matrix Market format.
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Read(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int >> &comm, bool binary=false)
Read matrix from file in Matrix Market or binary format.
RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getTpetra_Map() const
Get the underlying Tpetra map.
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
static void Write(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save matrix to file in Matrix Market format.
static void WriteLocal(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save local parts of matrix to files in Matrix Market format.
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadMultiVector(const std::string &fileName, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> &map, const bool binary=false)
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > Convert_Epetra_MultiVector_ToXpetra_MultiVector(RCP< Epetra_MultiVector > &epX)
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
Get the underlying Tpetra matrix.
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
Concrete implementation of Xpetra::Matrix.
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
virtual size_t Rows() const
number of row blocks
virtual const Teuchos::RCP< const map_type > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_MultiVector() const
Get the underlying Tpetra multivector.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
Returns the Map associated with the domain of this graph.
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
static void Write(const std::string &fileName, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save vector to file in Matrix Market format.
static void Write(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save matrix to file in Matrix Market format.
static void Write(const std::string &fileName, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save vector to file in Matrix Market format.
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Read(const std::string &filename, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rowMap, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> colMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> domainMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> rangeMap=Teuchos::null, const bool callFillComplete=true, const bool binary=false, const bool tolerant=false, const bool debug=false)
Read matrix from file in Matrix Market or binary format.
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)
Constructor specifying the number of non-zeros for all rows.
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Read(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int >> &comm, bool binary=false)
Read matrix from file in Matrix Market or binary format.
Xpetra-specific matrix class.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
Returns the Map that describes the row distribution in this graph.
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int >> &comm, const bool binary=false)
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.