15#include <Tpetra_RowMatrixTransposer.hpp>
16#include <TpetraExt_MatrixMatrix.hpp>
17#include <Xpetra_TpetraMultiVector.hpp>
18#include <Xpetra_TpetraCrsGraph.hpp>
19#include <Xpetra_TpetraCrsMatrix.hpp>
20#include <Xpetra_TpetraBlockCrsMatrix.hpp>
23#include "Xpetra_Matrix.hpp"
24#include "Xpetra_MatrixMatrix.hpp"
26#include "Xpetra_CrsMatrixWrap.hpp"
27#include "Xpetra_BlockedCrsMatrix.hpp"
29#include "Xpetra_Map.hpp"
30#include "Xpetra_StridedMap.hpp"
31#include "Xpetra_StridedMapFactory.hpp"
32#include "Xpetra_MapExtractor.hpp"
33#include "Xpetra_MatrixFactory.hpp"
35#include <Teuchos_MatrixMarket_Raw_Writer.hpp>
36#include <Teuchos_MatrixMarket_Raw_Reader.hpp>
41template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
44 if (tmp_TMap == Teuchos::null)
45 throw Exceptions::BadCast(
"Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
46 return tmp_TMap->getTpetra_Map();
49template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
54 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(tmp_Map);
55 if (tmp_TMap != Teuchos::null) {
64template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 std::string mapfile =
"map_" + fileName;
72 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
73 if (tmp_TVec != Teuchos::null) {
79 throw Exceptions::BadCast(
"Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
82template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 std::string mapfile =
"map_" + fileName;
89 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
90 if (tmp_TVec != Teuchos::null) {
95 throw Exceptions::RuntimeError(
"Xpetra cannot write MV<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
98 throw Exceptions::BadCast(
"Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
101template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 std::string mapfile =
"map_" + fileName;
108 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
109 if (tmp_TVec != Teuchos::null) {
114 throw Exceptions::RuntimeError(
"Xpetra cannot write MV<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
117 throw Exceptions::BadCast(
"Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
120template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
135 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_CrsMtx);
136 if (tmp_TCrsMtx != Teuchos::null) {
142 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_CrsMtx);
143 if (tmp_BlockCrs != Teuchos::null) {
144 std::ofstream outstream(fileName, std::ofstream::out);
150 throw Exceptions::BadCast(
"Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
153template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
163 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
171 for (LocalOrdinal j = 0; j < rowptr.
size(); j++)
172 rowptr2[j] = rowptr[j];
175 writer.
writeFile(fileName +
"." + std::to_string(Op.
getRowMap()->getComm()->getSize()) +
"." + std::to_string(Op.
getRowMap()->getComm()->getRank()),
176 rowptr2, colind, vals,
180template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
181void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::WriteBlockedCrsMatrix(
const std::string& fileName,
const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op,
const bool& writeAllMaps) {
183 for (
size_t row = 0; row < Op.
Rows(); ++row) {
184 for (
size_t col = 0; col < Op.
Cols(); ++col) {
186 if (m != Teuchos::null) {
188 "Sub block matrix (" << row <<
"," << col <<
") is not of type CrsMatrixWrap.");
198 for (
size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
200 Write(
"subRangeMap_" + fileName +
toString(row) +
".m", *map);
202 Write(
"fullRangeMap_" + fileName +
".m", *(rangeMapExtractor->getFullMap()));
204 for (
size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
206 Write(
"subDomainMap_" + fileName +
toString(col) +
".m", *map);
208 Write(
"fullDomainMap_" + fileName +
".m", *(domainMapExtractor->getFullMap()));
211template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
213 if (binary ==
false) {
220 bool callFillComplete =
true;
237 std::ifstream ifs(fileName.c_str(), std::ios::binary);
240 ifs.read(
reinterpret_cast<char*
>(&m),
sizeof(m));
241 ifs.read(
reinterpret_cast<char*
>(&n),
sizeof(n));
242 ifs.read(
reinterpret_cast<char*
>(&nnz),
sizeof(nnz));
244 int myRank = comm->getRank();
246 GlobalOrdinal indexBase = 0;
256 for (
int i = 0; i < m; i++) {
258 ifs.read(
reinterpret_cast<char*
>(&row),
sizeof(row));
259 ifs.read(
reinterpret_cast<char*
>(&rownnz),
sizeof(rownnz));
260 numEntriesPerRow[row] = rownnz;
261 for (
int j = 0; j < rownnz; j++) {
263 ifs.read(
reinterpret_cast<char*
>(&index),
sizeof(index));
265 for (
int j = 0; j < rownnz; j++) {
267 ifs.read(
reinterpret_cast<char*
>(&value),
sizeof(value));
271 A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
274 ifs.seekg(0, ifs.beg);
276 ifs.read(
reinterpret_cast<char*
>(&junk),
sizeof(junk));
277 ifs.read(
reinterpret_cast<char*
>(&junk),
sizeof(junk));
278 ifs.read(
reinterpret_cast<char*
>(&junk),
sizeof(junk));
279 for (
int i = 0; i < m; i++) {
281 ifs.read(
reinterpret_cast<char*
>(&row),
sizeof(row));
282 ifs.read(
reinterpret_cast<char*
>(&rownnz),
sizeof(rownnz));
285 for (
int j = 0; j < rownnz; j++) {
287 ifs.read(
reinterpret_cast<char*
>(&index),
sizeof(index));
290 for (
int j = 0; j < rownnz; j++) {
292 ifs.read(
reinterpret_cast<char*
>(&value),
sizeof(value));
295 A->insertGlobalValues(row, inds, vals);
300 A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
303 A->fillComplete(domainMap, rangeMap);
311template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
313IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Read(
const std::string& filename,
318 const bool callFillComplete,
328 if (binary ==
false) {
339 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
340 callFillComplete, tolerant, debug);
354 auto tempA =
Read(filename, lib, rowMap->getComm(), binary);
356 auto A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, 0);
359 if (callFillComplete)
360 A->fillComplete(domainMap, rangeMap);
368template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
374 const bool callFillComplete,
388 std::string rankFilename = filename +
"." + std::to_string(rowMap->getComm()->getSize()) +
"." + std::to_string(rowMap->getComm()->getRank());
391 if (binary ==
false) {
393 params->set(
"Parse tolerantly", tolerant);
394 params->set(
"Debug mode", debug);
396 LocalOrdinal numRows = rowMap->getLocalNumElements();
397 LocalOrdinal numCols = colMap->getLocalNumElements();
404 reader.
readFile(rowptr2_RCP, colind2_RCP, vals2_RCP,
408 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
413 ACrs->allocateAllValues(colind2_RCP.
size(), rowptr_RCP, colind_RCP, vals_RCP);
416 colind_RCP = colind2_RCP;
417 vals_RCP = vals2_RCP;
419 ACrs->setAllValues(rowptr_RCP, colind_RCP, vals_RCP);
422 std::ifstream ifs = std::ifstream(rankFilename.c_str(), std::ios::binary);
426 ifs.read(
reinterpret_cast<char*
>(&m),
sizeof(m));
427 ifs.read(
reinterpret_cast<char*
>(&n),
sizeof(n));
428 ifs.read(
reinterpret_cast<char*
>(&nnz),
sizeof(nnz));
436 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
438 ACrs->allocateAllValues(nnz, rowptrRCP, indicesRCP, valuesRCP);
447 for (
int i = 0; i < m; i++) {
449 ifs.read(
reinterpret_cast<char*
>(&row),
sizeof(row));
450 ifs.read(
reinterpret_cast<char*
>(&rownnz),
sizeof(rownnz));
452 rowptr[row + 1] += rownnz;
453 ifs.seekg(
sizeof(
int) * rownnz +
sizeof(
double) * rownnz, ifs.cur);
455 for (
int i = 0; i < m; i++)
456 rowptr[i + 1] += rowptr[i];
460 ifs.seekg(
sizeof(
int) * 3, ifs.beg);
463 for (
int i = 0; i < m; i++) {
465 ifs.read(
reinterpret_cast<char*
>(&row),
sizeof(row));
466 ifs.read(
reinterpret_cast<char*
>(&rownnz),
sizeof(rownnz));
467 size_t ptr = rowptr[row];
468 for (
int j = 0; j < rownnz; j++) {
470 ifs.read(
reinterpret_cast<char*
>(&index),
sizeof(index));
473 sorted = sorted & (indices[ptr - 1] < indices[ptr]);
477 for (
int j = 0; j < rownnz; j++) {
479 ifs.read(
reinterpret_cast<char*
>(&value),
sizeof(value));
483 rowptr[row] += rownnz;
485 for (
int i = m; i > 0; i--)
486 rowptr[i] = rowptr[i - 1];
490 for (LocalOrdinal lclRow = 0; lclRow < m; lclRow++) {
491 size_t rowBegin = rowptr[lclRow];
492 size_t rowEnd = rowptr[lclRow + 1];
493 Tpetra::sort2(&indices[rowBegin], &indices[rowEnd], &values[rowBegin]);
497 ACrs->setAllValues(rowptrRCP, indicesRCP, valuesRCP);
500 if (callFillComplete)
501 A->fillComplete(domainMap, rangeMap);
505template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
518 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp,
false,
false, binary);
528template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
541 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp,
false,
false, binary);
551template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
572template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
574 size_t numBlocks = 2;
576 std::vector<RCP<const Map>> rangeMapVec;
577 for (
size_t row = 0; row < numBlocks; ++row) {
579 rangeMapVec.push_back(map);
583 std::vector<RCP<const Map>> domainMapVec;
584 for (
size_t col = 0; col < numBlocks; ++col) {
586 domainMapVec.push_back(map);
604 bool bRangeUseThyraStyleNumbering =
false;
612 Teuchos::rcp(
new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
615 bool bDomainUseThyraStyleNumbering =
false;
623 Teuchos::rcp(
new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
628 for (
size_t row = 0; row < numBlocks; ++row) {
629 for (
size_t col = 0; col < numBlocks; ++col) {
635 bOp->setMatrix(row, col, mat);
644template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
646std::string IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toString(
const T& what) {
647 std::ostringstream buf;
void resize(size_type new_size, const value_type &x=value_type())
void resize(const size_type n, const T &val=T())
void assign(size_type n, const T &val)
bool readFile(ArrayRCP< Ordinal > &rowptr, ArrayRCP< Ordinal > &colind, ArrayRCP< Scalar > &values, Ordinal &numRows, Ordinal &numCols, const std::string &filename)
void writeFile(const std::string &filename, const ArrayView< const OrdinalType > &rowptr, const ArrayView< const OrdinalType > &colind, const ArrayView< const ScalarType > &values, const OrdinalType numRows, const OrdinalType numCols)
virtual size_t Cols() const
number of column blocks
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual size_t Rows() const
number of row blocks
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Concrete implementation of Xpetra::Matrix.
RCP< CrsMatrix > getCrsMatrix() const
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
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.
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm, const bool binary=false)
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.
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
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.
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.
Xpetra-specific matrix class.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
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....
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....
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
TypeTo as(const TypeFrom &t)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
basic_FancyOStream< char > FancyOStream
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2, const bool stableSort=false)
RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toXpetra(RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > graph)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.