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;
67 Write(mapfile, *(vec.
getMap()));
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;
85 Write(mapfile, *(vec.
getMap()));
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;
104 Write(mapfile, *(vec.
getMap()));
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>
122 Write(
"rowmap_" + fileName, *(Op.
getRowMap()));
126 Write(
"rangemap_" + fileName, *(Op.
getRangeMap()));
128 Write(
"colmap_" + fileName, *(Op.
getColMap()));
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>
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));
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));
288 inds[j] = Teuchos::as<GlobalOrdinal>(index);
290 for (
int j = 0; j < rownnz; j++) {
292 ifs.read(
reinterpret_cast<char*
>(&value),
sizeof(value));
293 vals[j] = Teuchos::as<Scalar>(value);
295 A->insertGlobalValues(row, inds, vals);
303 A->fillComplete(domainMap, rangeMap);
311template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
318 const bool callFillComplete,
328 if (binary ==
false) {
335 RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
336 const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
337 const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
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);
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));
471 indices[ptr] = Teuchos::as<LocalOrdinal>(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));
480 values[ptr] = Teuchos::as<Scalar>(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);
581 RCP<const Map> fullRangeMap = ReadMap(
"fullRangeMap_" + fileName +
".m", lib, comm);
583 std::vector<RCP<const Map>> domainMapVec;
584 for (
size_t col = 0; col < numBlocks; ++col) {
586 domainMapVec.push_back(map);
588 RCP<const Map> fullDomainMap = ReadMap(
"fullDomainMap_" + fileName +
".m", lib, comm);
604 bool bRangeUseThyraStyleNumbering =
false;
615 bool bDomainUseThyraStyleNumbering =
false;
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>
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 RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map, const bool binary=false)
Read a MultiVector from file in Matrix Matrix or binary format.
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< 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 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 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 void WriteGOMV(const std::string &fileName, const Xpetra::MultiVector< GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save multivector with Scalar=GlobalOrdinal 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< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReadMultiVectorLO(const std::string &fileName, const RCP< const Map > &map, const bool binary=false)
Read a MultiVector with Scalar=LocalOrdinal from file in Matrix Matrix or binary format.
static void WriteLOMV(const std::string &fileName, const Xpetra::MultiVector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save multivector with Scalar=LocalOrdinal to file in Matrix Market format.
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
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.
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
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)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
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)
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.