Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_IO_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Xpetra: A linear algebra interface package
4//
5// Copyright 2012 NTESS and the Xpetra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#include "Xpetra_IO_decl.hpp"
11#include <fstream>
12#include "Xpetra_ConfigDefs.hpp"
13
14#ifdef HAVE_XPETRA_EPETRA
15#ifdef HAVE_MPI
16#include "Epetra_MpiComm.h"
17#endif
18#endif
19
20#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
21#include <EpetraExt_MatrixMatrix.h>
22#include <EpetraExt_RowMatrixOut.h>
23#include <EpetraExt_MultiVectorOut.h>
24#include <EpetraExt_CrsMatrixIn.h>
25#include <EpetraExt_MultiVectorIn.h>
26#include <EpetraExt_BlockMapIn.h>
27#include <Xpetra_EpetraUtils.hpp>
28#include <Xpetra_EpetraMultiVector.hpp>
29#include <EpetraExt_BlockMapOut.h>
30#endif
31
32#ifdef HAVE_XPETRA_TPETRA
33#include <MatrixMarket_Tpetra.hpp>
34#include <Tpetra_RowMatrixTransposer.hpp>
35#include <TpetraExt_MatrixMatrix.hpp>
36#include <Xpetra_TpetraMultiVector.hpp>
37#include <Xpetra_TpetraCrsGraph.hpp>
38#include <Xpetra_TpetraCrsMatrix.hpp>
39#include <Xpetra_TpetraBlockCrsMatrix.hpp>
40#include "Tpetra_Util.hpp"
41#endif
42
43#ifdef HAVE_XPETRA_EPETRA
44#include <Xpetra_EpetraMap.hpp>
45#endif
46
47#include "Xpetra_Matrix.hpp"
48#include "Xpetra_MatrixMatrix.hpp"
49#include "Xpetra_CrsGraph.hpp"
50#include "Xpetra_CrsMatrixWrap.hpp"
51#include "Xpetra_BlockedCrsMatrix.hpp"
52
53#include "Xpetra_Map.hpp"
54#include "Xpetra_StridedMap.hpp"
55#include "Xpetra_StridedMapFactory.hpp"
56#include "Xpetra_MapExtractor.hpp"
57#include "Xpetra_MatrixFactory.hpp"
58
59#include <Teuchos_MatrixMarket_Raw_Writer.hpp>
60#include <Teuchos_MatrixMarket_Raw_Reader.hpp>
61#include <string>
62
63namespace Xpetra {
64
65#ifdef HAVE_XPETRA_EPETRA
66template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
67const Epetra_Map& IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& map) {
68 RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>> xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(Teuchos::rcpFromRef(map));
69 if (xeMap == Teuchos::null)
70 throw Exceptions::BadCast("Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
71 return xeMap->getEpetra_Map();
72}
73#endif
74
75#ifdef HAVE_XPETRA_TPETRA
76template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77const RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2TpetraMap(const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& map) {
78 const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(rcpFromRef(map));
79 if (tmp_TMap == Teuchos::null)
80 throw Exceptions::BadCast("Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
81 return tmp_TMap->getTpetra_Map();
82}
83#endif
84
85template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& M) {
87 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tmp_Map = rcpFromRef(M);
88#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
89 const RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(tmp_Map);
90 if (tmp_EMap != Teuchos::null) {
91 int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
92 if (rv != 0)
93 throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
94 return;
95 }
96#endif // HAVE_XPETRA_EPETRAEXT
97
98#ifdef HAVE_XPETRA_TPETRA
99 const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>& tmp_TMap =
100 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(tmp_Map);
101 if (tmp_TMap != Teuchos::null) {
102 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> TMap = tmp_TMap->getTpetra_Map();
103 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeMapFile(fileName, *TMap);
104 return;
105 }
106#endif // HAVE_XPETRA_TPETRA
107
108 throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
109}
110
111template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) {
113 std::string mapfile = "map_" + fileName;
114 Write(mapfile, *(vec.getMap()));
115
116 RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_Vec = Teuchos::rcpFromRef(vec);
117#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
118 const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>>& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>>(tmp_Vec);
119 if (tmp_EVec != Teuchos::null) {
120 int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
121 if (rv != 0)
122 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
123 return;
124 }
125#endif // HAVE_XPETRA_EPETRA
126
127#ifdef HAVE_XPETRA_TPETRA
128 const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TVec =
129 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
130 if (tmp_TVec != Teuchos::null) {
131 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TVec = tmp_TVec->getTpetra_MultiVector();
132 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
133 return;
134 }
135#endif // HAVE_XPETRA_TPETRA
136
137 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
138}
139
140template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::WriteLOMV(const std::string& fileName, const Xpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& vec) {
142 std::string mapfile = "map_" + fileName;
143 Write(mapfile, *(vec.getMap()));
144
145 RCP<const Xpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> tmp_Vec = Teuchos::rcpFromRef(vec);
146#ifdef HAVE_XPETRA_TPETRA
147 const RCP<const Xpetra::TpetraMultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TVec =
148 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
149 if (tmp_TVec != Teuchos::null) {
150 RCP<const Tpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> TVec = tmp_TVec->getTpetra_MultiVector();
151 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
152 return;
153 } else
154#endif // HAVE_XPETRA_TPETRA
155 {
156 throw Exceptions::RuntimeError("Xpetra cannot write MV<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
157 }
158
159 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
160}
161
162template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
163void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::WriteGOMV(const std::string& fileName, const Xpetra::MultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& vec) {
164 std::string mapfile = "map_" + fileName;
165 Write(mapfile, *(vec.getMap()));
166
167 RCP<const Xpetra::MultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> tmp_Vec = Teuchos::rcpFromRef(vec);
168#ifdef HAVE_XPETRA_TPETRA
169 const RCP<const Xpetra::TpetraMultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TVec =
170 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
171 if (tmp_TVec != Teuchos::null) {
172 RCP<const Tpetra::MultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> TVec = tmp_TVec->getTpetra_MultiVector();
173 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
174 return;
175 } else
176#endif // HAVE_XPETRA_TPETRA
177 {
178 throw Exceptions::RuntimeError("Xpetra cannot write MV<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
179 }
180
181 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
182}
183
184template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
185void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const bool& writeAllMaps) {
186 Write("rowmap_" + fileName, *(Op.getRowMap()));
187 if (!Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps)
188 Write("domainmap_" + fileName, *(Op.getDomainMap()));
189 if (!Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps)
190 Write("rangemap_" + fileName, *(Op.getRangeMap()));
191 if (!Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps)
192 Write("colmap_" + fileName, *(Op.getColMap()));
193
194 const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
195 dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
196 RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_CrsMtx = crsOp.getCrsMatrix();
197#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
198 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(tmp_CrsMtx);
199 if (tmp_ECrsMtx != Teuchos::null) {
200 RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
201 int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
202 if (rv != 0)
203 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
204 return;
205 }
206#endif
207
208#ifdef HAVE_XPETRA_TPETRA
209 const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_TCrsMtx =
210 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_CrsMtx);
211 if (tmp_TCrsMtx != Teuchos::null) {
212 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A = tmp_TCrsMtx->getTpetra_CrsMatrix();
213 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeSparseFile(fileName, A);
214 return;
215 }
216 const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_BlockCrs =
217 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_CrsMtx);
218 if (tmp_BlockCrs != Teuchos::null) {
219 std::ofstream outstream(fileName, std::ofstream::out);
220 Teuchos::FancyOStream ofs(Teuchos::rcpFromRef(outstream));
221 tmp_BlockCrs->getTpetra_BlockCrsMatrix()->describe(ofs, Teuchos::VERB_EXTREME);
222 return;
223 }
224
225#endif // HAVE_XPETRA_TPETRA
226
227 throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
228}
229
230template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
231void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
232 const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
233 dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
234 RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_CrsMtx = crsOp.getCrsMatrix();
235
236 ArrayRCP<const size_t> rowptr_RCP;
237 ArrayRCP<LocalOrdinal> rowptr2_RCP;
238 ArrayRCP<const LocalOrdinal> colind_RCP;
239 ArrayRCP<const Scalar> vals_RCP;
240 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
241
242 ArrayView<const size_t> rowptr = rowptr_RCP();
243 ArrayView<const LocalOrdinal> colind = colind_RCP();
244 ArrayView<const Scalar> vals = vals_RCP();
245
246 rowptr2_RCP.resize(rowptr.size());
247 ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
248 for (LocalOrdinal j = 0; j < rowptr.size(); j++)
249 rowptr2[j] = rowptr[j];
250
251 Teuchos::MatrixMarket::Raw::Writer<Scalar, LocalOrdinal> writer;
252 writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
253 rowptr2, colind, vals,
254 rowptr.size() - 1, Op.getColMap()->getLocalNumElements());
255}
256
257template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
258void IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const bool& writeAllMaps) {
259 // Write all matrix blocks with their maps
260 for (size_t row = 0; row < Op.Rows(); ++row) {
261 for (size_t col = 0; col < Op.Cols(); ++col) {
262 RCP<const Matrix> m = Op.getMatrix(row, col);
263 if (m != Teuchos::null) { // skip empty blocks
264 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
265 "Sub block matrix (" << row << "," << col << ") is not of type CrsMatrixWrap.");
266 Write(fileName + toString(row) + toString(col) + ".m", *m, writeAllMaps);
267 }
268 }
269 }
270
271 // write map information of map extractors
272 RCP<const MapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
273 RCP<const MapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
274
275 for (size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
276 RCP<const Map> map = rangeMapExtractor->getMap(row);
277 Write("subRangeMap_" + fileName + toString(row) + ".m", *map);
278 }
279 Write("fullRangeMap_" + fileName + ".m", *(rangeMapExtractor->getFullMap()));
280
281 for (size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
282 RCP<const Map> map = domainMapExtractor->getMap(col);
283 Write("subDomainMap_" + fileName + toString(col) + ".m", *map);
284 }
285 Write("fullDomainMap_" + fileName + ".m", *(domainMapExtractor->getFullMap()));
286}
287
288template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
289Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Read(const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int>>& comm, bool binary) {
290 if (binary == false) {
291 // Matrix Market file format (ASCII)
292 if (lib == Xpetra::UseEpetra) {
293#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
294 Epetra_CrsMatrix* eA;
295 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
296 int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
297 if (rv != 0)
298 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
299
300 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
301
302 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A =
303 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
304 return A;
305#else
306 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
307#endif
308 } else if (lib == Xpetra::UseTpetra) {
309#ifdef HAVE_XPETRA_TPETRA
310 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
311
312 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
313
314 bool callFillComplete = true;
315
316 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
317
318 if (tA.is_null())
319 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
320
321 RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
322 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmpA1);
323 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
324
325 return A;
326#else
327 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
328#endif
329 } else {
330 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
331 }
332 } else {
333 // Custom file format (binary)
334 std::ifstream ifs(fileName.c_str(), std::ios::binary);
335 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
336 int m, n, nnz;
337 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
338 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
339 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
340
341 int myRank = comm->getRank();
342
343 GlobalOrdinal indexBase = 0;
344 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
345 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
346 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A;
347
348 if (myRank == 0) {
349 Teuchos::Array<GlobalOrdinal> inds;
350 Teuchos::Array<Scalar> vals;
351 // Scan matrix to determine the exact nnz per row.
352 Teuchos::ArrayRCP<size_t> numEntriesPerRow(m, (size_t)(0));
353 for (int i = 0; i < m; i++) {
354 int row, rownnz;
355 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
356 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
357 numEntriesPerRow[row] = rownnz;
358 for (int j = 0; j < rownnz; j++) {
359 int index;
360 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
361 }
362 for (int j = 0; j < rownnz; j++) {
363 double value;
364 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
365 }
366 }
367
368 A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
369
370 // Now that nnz per row are known, reread and store the matrix.
371 ifs.seekg(0, ifs.beg); // rewind to beginning of file
372 int junk; // skip header info
373 ifs.read(reinterpret_cast<char*>(&junk), sizeof(junk));
374 ifs.read(reinterpret_cast<char*>(&junk), sizeof(junk));
375 ifs.read(reinterpret_cast<char*>(&junk), sizeof(junk));
376 for (int i = 0; i < m; i++) {
377 int row, rownnz;
378 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
379 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
380 inds.resize(rownnz);
381 vals.resize(rownnz);
382 for (int j = 0; j < rownnz; j++) {
383 int index;
384 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
385 inds[j] = Teuchos::as<GlobalOrdinal>(index);
386 }
387 for (int j = 0; j < rownnz; j++) {
388 double value;
389 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
390 vals[j] = Teuchos::as<Scalar>(value);
391 }
392 A->insertGlobalValues(row, inds, vals);
393 }
394 } // if (myRank == 0)
395 else {
396 Teuchos::ArrayRCP<size_t> numEntriesPerRow(0, (size_t)(0));
397 A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
398 }
399
400 A->fillComplete(domainMap, rangeMap);
401
402 return A;
403 } // if (binary == false) ... else
404
405 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
406}
407
408template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
409Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
410IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Read(const std::string& filename,
411 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap,
412 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap,
413 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domainMap,
414 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rangeMap,
415 const bool callFillComplete,
416 const bool binary,
417 const bool tolerant,
418 const bool debug) {
419 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
420
421 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
422 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
423
424 const Xpetra::UnderlyingLib lib = rowMap->lib();
425 if (binary == false) {
426 if (lib == Xpetra::UseEpetra) {
427#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
428 Epetra_CrsMatrix* eA;
429 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
430 const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(*rowMap);
431 const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(*domainMap));
432 const Epetra_Map& epetraRangeMap = (rangeMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(*rangeMap));
433 int rv;
434 if (colMap.is_null()) {
435 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
436
437 } else {
438 const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
439 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
440 }
441
442 if (rv != 0)
443 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
444
445 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
446 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A =
447 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
448
449 return A;
450#else
451 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
452#endif
453 } else if (lib == Xpetra::UseTpetra) {
454#ifdef HAVE_XPETRA_TPETRA
455 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
456 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
457 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
458
459 const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
460 RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
461 const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
462 const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
463
464 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
465 callFillComplete, tolerant, debug);
466 if (tA.is_null())
467 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
468
469 RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
470 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmpA1);
471 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
472
473 return A;
474#else
475 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
476#endif
477 } else {
478 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
479 }
480 } else {
481 // Read in on rank 0.
482 auto tempA = Read(filename, lib, rowMap->getComm(), binary);
483
484 auto A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, 0);
485 auto importer = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(tempA->getRowMap(), rowMap);
486 A->doImport(*tempA, *importer, Xpetra::INSERT);
487 if (callFillComplete)
488 A->fillComplete(domainMap, rangeMap);
489
490 return A;
491 }
492
493 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
494}
495
496template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
497Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReadLocal(const std::string& filename,
498 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap,
499 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap,
500 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domainMap,
501 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rangeMap,
502 const bool callFillComplete,
503 const bool binary,
504 const bool tolerant,
505 const bool debug) {
506 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::ReadLocal() : rowMap cannot be null");
507 TEUCHOS_TEST_FOR_EXCEPTION(colMap.is_null(), Exceptions::RuntimeError, "Utils::ReadLocal() : colMap cannot be null");
508
509 using matrix_type = Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
510 using crs_wrap_type = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
511 using crs_type = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
512
513 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
514 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
515
516 std::string rankFilename = filename + "." + std::to_string(rowMap->getComm()->getSize()) + "." + std::to_string(rowMap->getComm()->getRank());
517 RCP<matrix_type> A = rcp(new crs_wrap_type(rowMap, colMap, 0));
518
519 if (binary == false) {
520 RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
521 params->set("Parse tolerantly", tolerant);
522 params->set("Debug mode", debug);
523
524 LocalOrdinal numRows = rowMap->getLocalNumElements();
525 LocalOrdinal numCols = colMap->getLocalNumElements();
526
527 ArrayRCP<LocalOrdinal> rowptr2_RCP;
528 ArrayRCP<LocalOrdinal> colind2_RCP;
529 ArrayRCP<Scalar> vals2_RCP;
530
531 Teuchos::MatrixMarket::Raw::Reader<Scalar, LocalOrdinal> reader;
532 reader.readFile(rowptr2_RCP, colind2_RCP, vals2_RCP,
533 numRows, numCols,
534 rankFilename);
535
536 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
537
538 ArrayRCP<size_t> rowptr_RCP;
539 ArrayRCP<LocalOrdinal> colind_RCP;
540 ArrayRCP<Scalar> vals_RCP;
541 ACrs->allocateAllValues(colind2_RCP.size(), rowptr_RCP, colind_RCP, vals_RCP);
542
543 rowptr_RCP.assign(rowptr2_RCP.begin(), rowptr2_RCP.end());
544 colind_RCP = colind2_RCP;
545 vals_RCP = vals2_RCP;
546
547 ACrs->setAllValues(rowptr_RCP, colind_RCP, vals_RCP);
548 } else {
549 // Custom file format (binary)
550 std::ifstream ifs = std::ifstream(rankFilename.c_str(), std::ios::binary);
551 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
552
553 int m, n, nnz;
554 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
555 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
556 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
557
558 TEUCHOS_ASSERT_EQUALITY(Teuchos::as<int>(rowMap->getLocalNumElements()), m);
559
560 Teuchos::ArrayRCP<size_t> rowptrRCP;
561 Teuchos::ArrayRCP<LocalOrdinal> indicesRCP;
562 Teuchos::ArrayRCP<Scalar> valuesRCP;
563
564 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
565
566 ACrs->allocateAllValues(nnz, rowptrRCP, indicesRCP, valuesRCP);
567
568 Teuchos::ArrayView<size_t> rowptr = rowptrRCP();
569 Teuchos::ArrayView<LocalOrdinal> indices = indicesRCP();
570 Teuchos::ArrayView<Scalar> values = valuesRCP();
571
572 bool sorted = true;
573
574 // Read in rowptr
575 for (int i = 0; i < m; i++) {
576 int row, rownnz;
577 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
578 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
579
580 rowptr[row + 1] += rownnz;
581 ifs.seekg(sizeof(int) * rownnz + sizeof(double) * rownnz, ifs.cur);
582 }
583 for (int i = 0; i < m; i++)
584 rowptr[i + 1] += rowptr[i];
585 TEUCHOS_ASSERT(Teuchos::as<int>(rowptr[m]) == nnz);
586
587 // reset to where the data starts
588 ifs.seekg(sizeof(int) * 3, ifs.beg);
589
590 // read in entries
591 for (int i = 0; i < m; i++) {
592 int row, rownnz;
593 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
594 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
595 size_t ptr = rowptr[row];
596 for (int j = 0; j < rownnz; j++) {
597 int index;
598 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
599 indices[ptr] = Teuchos::as<LocalOrdinal>(index);
600 if (j > 0)
601 sorted = sorted & (indices[ptr - 1] < indices[ptr]);
602 ++ptr;
603 }
604 ptr = rowptr[row];
605 for (int j = 0; j < rownnz; j++) {
606 double value;
607 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
608 values[ptr] = Teuchos::as<Scalar>(value);
609 ++ptr;
610 }
611 rowptr[row] += rownnz;
612 }
613 for (int i = m; i > 0; i--)
614 rowptr[i] = rowptr[i - 1];
615 rowptr[0] = 0;
616
617#ifdef HAVE_XPETRA_TPETRA
618 if (!sorted) {
619 for (LocalOrdinal lclRow = 0; lclRow < m; lclRow++) {
620 size_t rowBegin = rowptr[lclRow];
621 size_t rowEnd = rowptr[lclRow + 1];
622 Tpetra::sort2(&indices[rowBegin], &indices[rowEnd], &values[rowBegin]);
623 }
624 }
625#else
626 TEUCHOS_ASSERT(sorted);
627#endif
628
629 ACrs->setAllValues(rowptrRCP, indicesRCP, valuesRCP);
630 }
631
632 if (callFillComplete)
633 A->fillComplete(domainMap, rangeMap);
634 return A;
635}
636
637template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
638RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReadMultiVector(const std::string& fileName,
639 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& map,
640 const bool binary) {
641 Xpetra::UnderlyingLib lib = map->lib();
642
643 if (lib == Xpetra::UseEpetra) {
644 TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
645
646 } else if (lib == Xpetra::UseTpetra) {
647#ifdef HAVE_XPETRA_TPETRA
648 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
649 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
650 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
651 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
652
653 RCP<const map_type> temp = toTpetra(map);
654 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp, false, false, binary);
655 RCP<MultiVector> rmv = Xpetra::toXpetra(TMV);
656 return rmv;
657#else
658 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
659#endif
660 } else {
661 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
662 }
663
664 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
665}
666
667template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
668RCP<Xpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReadMultiVectorLO(const std::string& fileName,
669 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& map,
670 const bool binary) {
671 Xpetra::UnderlyingLib lib = map->lib();
672
673 if (lib == Xpetra::UseTpetra) {
674#ifdef HAVE_XPETRA_TPETRA
675 typedef Tpetra::CrsMatrix<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
676 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
677 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
678 typedef Tpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
679
680 RCP<const map_type> temp = toTpetra(map);
681 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp, false, false, binary);
682 RCP<Xpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> rmv = Xpetra::toXpetra(TMV);
683 return rmv;
684#else
685 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
686#endif
687 } else {
688 throw Exceptions::RuntimeError("Utils::ReadMultiVectorLO : only implemented for Tpetra");
689 }
690
691 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
692}
693
694template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
695RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReadMap(const std::string& fileName,
696 Xpetra::UnderlyingLib lib,
697 const RCP<const Teuchos::Comm<int>>& comm,
698 const bool binary) {
699 if (lib == Xpetra::UseEpetra) {
700 TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
701 } else if (lib == Xpetra::UseTpetra) {
702#ifdef HAVE_XPETRA_TPETRA
703 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
704 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
705
706 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tMap = reader_type::readMapFile(fileName, comm, false, false, binary);
707 if (tMap.is_null())
708 throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
709
710 return Xpetra::toXpetra(tMap);
711#else
712 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
713#endif
714 } else {
715 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
716 }
717
718 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
719}
720
721template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
722RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReadBlockedCrsMatrix(const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int>>& comm) {
723 size_t numBlocks = 2; // TODO user parameter?
724
725 std::vector<RCP<const Map>> rangeMapVec;
726 for (size_t row = 0; row < numBlocks; ++row) {
727 RCP<const Map> map = ReadMap("subRangeMap_" + fileName + toString(row) + ".m", lib, comm);
728 rangeMapVec.push_back(map);
729 }
730 RCP<const Map> fullRangeMap = ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
731
732 std::vector<RCP<const Map>> domainMapVec;
733 for (size_t col = 0; col < numBlocks; ++col) {
734 RCP<const Map> map = ReadMap("subDomainMap_" + fileName + toString(col) + ".m", lib, comm);
735 domainMapVec.push_back(map);
736 }
737 RCP<const Map> fullDomainMap = ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
738
739 /*std::vector<RCP<const XpMap> > testRgMapVec;
740 for(size_t r = 0; r < numBlocks; ++r) {
741 RCP<const XpMap> map = ReadMap("rangemap_" + fileName + toString<size_t>(r) + "0.m", lib, comm);
742 testRgMapVec.push_back(map);
743 }
744 std::vector<RCP<const XpMap> > testDoMapVec;
745 for(size_t c = 0; c < numBlocks; ++c) {
746 RCP<const XpMap> map = ReadMap("domainmap_" + fileName + "0" + toString<size_t>(c) + ".m", lib, comm);
747 testDoMapVec.push_back(map);
748 }*/
749
750 // create map extractors
751
752 // range map extractor
753 bool bRangeUseThyraStyleNumbering = false;
754 /*GlobalOrdinal gMinGids = 0;
755 for(size_t v = 0; v < testRgMapVec.size(); ++v) {
756 gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
757 }
758 if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
759 */
760 RCP<const MapExtractor> rangeMapExtractor =
761 Teuchos::rcp(new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
762
763 // domain map extractor
764 bool bDomainUseThyraStyleNumbering = false;
765 /*gMinGids = 0;
766 for(size_t v = 0; v < testDoMapVec.size(); ++v) {
767 gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
768 }
769 if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
770 */
771 RCP<const MapExtractor> domainMapExtractor =
772 Teuchos::rcp(new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
773
774 RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 33));
775
776 // Read all sub-matrices with their maps and place into blocked operator
777 for (size_t row = 0; row < numBlocks; ++row) {
778 for (size_t col = 0; col < numBlocks; ++col) {
779 RCP<const Map> rowSubMap = ReadMap("rowmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
780 RCP<const Map> colSubMap = ReadMap("colmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
781 RCP<const Map> domSubMap = ReadMap("domainmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
782 RCP<const Map> ranSubMap = ReadMap("rangemap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
783 RCP<Matrix> mat = Read(fileName + toString(row) + toString(col) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
784 bOp->setMatrix(row, col, mat);
785 }
786 }
787
788 bOp->fillComplete();
789
790 return bOp;
791} // ReadBlockedCrsMatrix
792
793template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
794template <class T>
795std::string IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toString(const T& what) {
796 std::ostringstream buf;
797 buf << what;
798 return buf.str();
799}
800} // namespace Xpetra
Xpetra_ConfigDefs.hpp
Xpetra_CrsGraph.hpp
Xpetra_EpetraMap.hpp
Xpetra_EpetraMultiVector.hpp
Xpetra_EpetraUtils.hpp
Xpetra_IO_decl.hpp
Epetra_CrsMatrix
Epetra_Map
Teuchos::ArrayRCP::size
size_type size() const
Teuchos::ArrayRCP::resize
void resize(const size_type n, const T &val=T())
Teuchos::ArrayRCP::assign
void assign(size_type n, const T &val)
Teuchos::ArrayRCP::begin
iterator begin() const
Teuchos::ArrayRCP::end
iterator end() const
Teuchos::ArrayView::size
size_type size() const
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
Teuchos::MatrixMarket::Raw::Reader::readFile
bool readFile(ArrayRCP< Ordinal > &rowptr, ArrayRCP< Ordinal > &colind, ArrayRCP< Scalar > &values, Ordinal &numRows, Ordinal &numCols, const std::string &filename)
Teuchos::MatrixMarket::Raw::Writer::writeFile
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)
Teuchos::RCP::is_null
bool is_null() const
Teuchos::basic_FancyOStream
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition Xpetra_BlockedCrsMatrix_def.hpp:1067
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition Xpetra_BlockedCrsMatrix_def.hpp:1103
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition Xpetra_BlockedCrsMatrix_def.hpp:1061
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition Xpetra_BlockedCrsMatrix_def.hpp:867
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition Xpetra_BlockedCrsMatrix_def.hpp:873
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition Xpetra_CrsMatrixWrap_decl.hpp:50
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition Xpetra_CrsMatrixWrap_def.hpp:437
Xpetra::DistObject::getMap
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
Xpetra::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition Xpetra_Exceptions.hpp:50
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition Xpetra_Exceptions.hpp:65
Xpetra::IO::Write
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.
Definition Xpetra_IO_def.hpp:185
Xpetra::IO::Read
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.
Definition Xpetra_IO_def.hpp:289
Xpetra::IO::ReadMap
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm, const bool binary=false)
Definition Xpetra_IO_def.hpp:695
Xpetra::IO::ReadMultiVector
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.
Definition Xpetra_IO_def.hpp:638
Xpetra::IO::Map2TpetraMap
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.
Definition Xpetra_IO_def.hpp:77
Xpetra::IO::Write
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition Xpetra_IO_def.hpp:86
Xpetra::IO::ReadBlockedCrsMatrix
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.
Definition Xpetra_IO_def.hpp:722
Xpetra::IO::ReadLocal
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.
Definition Xpetra_IO_def.hpp:497
Xpetra::IO::WriteBlockedCrsMatrix
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.
Definition Xpetra_IO_def.hpp:258
Xpetra::IO::WriteGOMV
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.
Definition Xpetra_IO_def.hpp:163
Xpetra::IO::Map2EpetraMap
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.
Definition Xpetra_IO_def.hpp:67
Xpetra::IO::WriteLocal
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.
Definition Xpetra_IO_def.hpp:231
Xpetra::IO::ReadMultiVectorLO
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.
Definition Xpetra_IO_def.hpp:668
Xpetra::IO::WriteLOMV
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.
Definition Xpetra_IO_def.hpp:141
Xpetra::IO::toString
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition Xpetra_IO_def.hpp:795
Xpetra::ImportFactory::Build
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.
Definition Xpetra_ImportFactory.hpp:38
Xpetra::MapFactory::Build
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.
Definition Xpetra_MapFactory_def.hpp:31
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition Xpetra_MatrixFactory_def.hpp:22
Xpetra::Matrix
Xpetra-specific matrix class.
Definition Xpetra_Matrix_decl.hpp:60
Xpetra::Matrix::getRowMap
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition Xpetra_Matrix_def.hpp:134
Xpetra::Matrix::getColMap
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
Definition Xpetra_Matrix_def.hpp:143
Xpetra::Operator::getDomainMap
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....
Xpetra::Operator::getRangeMap
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....
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_ASSERT_EQUALITY
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
TEUCHOS_UNREACHABLE_RETURN
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
Teuchos::toString
std::string toString(const T &t)
Teuchos::VERB_EXTREME
VERB_EXTREME
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra::toEpetra
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
Definition Xpetra_EpetraCrsGraph.cpp:21
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition Xpetra_EpetraCrsGraph.cpp:29
Xpetra::toTpetra
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
Definition Xpetra_TpetraCrsGraph_decl.hpp:376
Xpetra::toString
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
Definition Xpetra_Utils.cpp:15