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
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>
21#include "Tpetra_Util.hpp"
22
23#include "Xpetra_Matrix.hpp"
24#include "Xpetra_MatrixMatrix.hpp"
25#include "Xpetra_CrsGraph.hpp"
26#include "Xpetra_CrsMatrixWrap.hpp"
27#include "Xpetra_BlockedCrsMatrix.hpp"
28
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"
34
35#include <Teuchos_MatrixMarket_Raw_Writer.hpp>
36#include <Teuchos_MatrixMarket_Raw_Reader.hpp>
37#include <string>
38
39namespace Xpetra {
40
41template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
43 const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(rcpFromRef(map));
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();
47}
48
49template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
52
54 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(tmp_Map);
55 if (tmp_TMap != Teuchos::null) {
56 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> TMap = tmp_TMap->getTpetra_Map();
58 return;
59 }
60
61 throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
62}
63
64template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
66 std::string mapfile = "map_" + fileName;
67 Write(mapfile, *(vec.getMap()));
68
70
72 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
73 if (tmp_TVec != Teuchos::null) {
76 return;
77 }
78
79 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
80}
81
82template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84 std::string mapfile = "map_" + fileName;
85 Write(mapfile, *(vec.getMap()));
86
89 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
90 if (tmp_TVec != Teuchos::null) {
93 return;
94 } else {
95 throw Exceptions::RuntimeError("Xpetra cannot write MV<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
96 }
97
98 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
99}
100
101template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
103 std::string mapfile = "map_" + fileName;
104 Write(mapfile, *(vec.getMap()));
105
108 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
109 if (tmp_TVec != Teuchos::null) {
112 return;
113 } else {
114 throw Exceptions::RuntimeError("Xpetra cannot write MV<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
115 }
116
117 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
118}
119
120template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
122 Write("rowmap_" + fileName, *(Op.getRowMap()));
123 if (!Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps)
124 Write("domainmap_" + fileName, *(Op.getDomainMap()));
125 if (!Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps)
126 Write("rangemap_" + fileName, *(Op.getRangeMap()));
127 if (!Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps)
128 Write("colmap_" + fileName, *(Op.getColMap()));
129
133
135 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_CrsMtx);
136 if (tmp_TCrsMtx != Teuchos::null) {
139 return;
140 }
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);
145 Teuchos::FancyOStream ofs(Teuchos::rcpFromRef(outstream));
146 tmp_BlockCrs->getTpetra_BlockCrsMatrix()->describe(ofs, Teuchos::VERB_EXTREME);
147 return;
148 }
149
150 throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
151}
152
153template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
158
159 ArrayRCP<const size_t> rowptr_RCP;
160 ArrayRCP<LocalOrdinal> rowptr2_RCP;
162 ArrayRCP<const Scalar> vals_RCP;
163 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
164
165 ArrayView<const size_t> rowptr = rowptr_RCP();
166 ArrayView<const LocalOrdinal> colind = colind_RCP();
167 ArrayView<const Scalar> vals = vals_RCP();
168
169 rowptr2_RCP.resize(rowptr.size());
170 ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
171 for (LocalOrdinal j = 0; j < rowptr.size(); j++)
172 rowptr2[j] = rowptr[j];
173
175 writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
176 rowptr2, colind, vals,
177 rowptr.size() - 1, Op.getColMap()->getLocalNumElements());
178}
179
180template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
182 // Write all matrix blocks with their maps
183 for (size_t row = 0; row < Op.Rows(); ++row) {
184 for (size_t col = 0; col < Op.Cols(); ++col) {
185 RCP<const Matrix> m = Op.getMatrix(row, col);
186 if (m != Teuchos::null) { // skip empty blocks
187 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
188 "Sub block matrix (" << row << "," << col << ") is not of type CrsMatrixWrap.");
189 Write(fileName + toString(row) + toString(col) + ".m", *m, writeAllMaps);
190 }
191 }
192 }
193
194 // write map information of map extractors
195 RCP<const MapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
196 RCP<const MapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
197
198 for (size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
199 RCP<const Map> map = rangeMapExtractor->getMap(row);
200 Write("subRangeMap_" + fileName + toString(row) + ".m", *map);
201 }
202 Write("fullRangeMap_" + fileName + ".m", *(rangeMapExtractor->getFullMap()));
203
204 for (size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
205 RCP<const Map> map = domainMapExtractor->getMap(col);
206 Write("subDomainMap_" + fileName + toString(col) + ".m", *map);
207 }
208 Write("fullDomainMap_" + fileName + ".m", *(domainMapExtractor->getFullMap()));
209}
210
211template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
213 if (binary == false) {
214 // Matrix Market file format (ASCII)
215 if (lib == Xpetra::UseTpetra) {
217
219
220 bool callFillComplete = true;
221
222 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
223
224 if (tA.is_null())
225 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
226
228 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmpA1);
230
231 return A;
232 } else {
233 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
234 }
235 } else {
236 // Custom file format (binary)
237 std::ifstream ifs(fileName.c_str(), std::ios::binary);
238 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
239 int m, n, nnz;
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));
243
244 int myRank = comm->getRank();
245
246 GlobalOrdinal indexBase = 0;
247 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
248 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
250
251 if (myRank == 0) {
254 // Scan matrix to determine the exact nnz per row.
255 Teuchos::ArrayRCP<size_t> numEntriesPerRow(m, (size_t)(0));
256 for (int i = 0; i < m; i++) {
257 int row, rownnz;
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++) {
262 int index;
263 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
264 }
265 for (int j = 0; j < rownnz; j++) {
266 double value;
267 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
268 }
269 }
270
272
273 // Now that nnz per row are known, reread and store the matrix.
274 ifs.seekg(0, ifs.beg); // rewind to beginning of file
275 int junk; // skip header info
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++) {
280 int row, rownnz;
281 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
282 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
283 inds.resize(rownnz);
284 vals.resize(rownnz);
285 for (int j = 0; j < rownnz; j++) {
286 int index;
287 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
288 inds[j] = Teuchos::as<GlobalOrdinal>(index);
289 }
290 for (int j = 0; j < rownnz; j++) {
291 double value;
292 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
293 vals[j] = Teuchos::as<Scalar>(value);
294 }
295 A->insertGlobalValues(row, inds, vals);
296 }
297 } // if (myRank == 0)
298 else {
299 Teuchos::ArrayRCP<size_t> numEntriesPerRow(0, (size_t)(0));
301 }
302
303 A->fillComplete(domainMap, rangeMap);
304
305 return A;
306 } // if (binary == false) ... else
307
308 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
309}
310
311template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
318 const bool callFillComplete,
319 const bool binary,
320 const bool tolerant,
321 const bool debug) {
322 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
323
324 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
325 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
326
327 const Xpetra::UnderlyingLib lib = rowMap->lib();
328 if (binary == false) {
329 if (lib == Xpetra::UseTpetra) {
333
334 const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
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));
338
339 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
340 callFillComplete, tolerant, debug);
341 if (tA.is_null())
342 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
343
345 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmpA1);
347
348 return A;
349 } else {
350 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
351 }
352 } else {
353 // Read in on rank 0.
354 auto tempA = Read(filename, lib, rowMap->getComm(), binary);
355
357 auto importer = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(tempA->getRowMap(), rowMap);
358 A->doImport(*tempA, *importer, Xpetra::INSERT);
359 if (callFillComplete)
360 A->fillComplete(domainMap, rangeMap);
361
362 return A;
363 }
364
365 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
366}
367
368template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
374 const bool callFillComplete,
375 const bool binary,
376 const bool tolerant,
377 const bool debug) {
378 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::ReadLocal() : rowMap cannot be null");
379 TEUCHOS_TEST_FOR_EXCEPTION(colMap.is_null(), Exceptions::RuntimeError, "Utils::ReadLocal() : colMap cannot be null");
380
384
385 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
386 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
387
388 std::string rankFilename = filename + "." + std::to_string(rowMap->getComm()->getSize()) + "." + std::to_string(rowMap->getComm()->getRank());
389 RCP<matrix_type> A = rcp(new crs_wrap_type(rowMap, colMap, 0));
390
391 if (binary == false) {
393 params->set("Parse tolerantly", tolerant);
394 params->set("Debug mode", debug);
395
396 LocalOrdinal numRows = rowMap->getLocalNumElements();
397 LocalOrdinal numCols = colMap->getLocalNumElements();
398
399 ArrayRCP<LocalOrdinal> rowptr2_RCP;
400 ArrayRCP<LocalOrdinal> colind2_RCP;
401 ArrayRCP<Scalar> vals2_RCP;
402
404 reader.readFile(rowptr2_RCP, colind2_RCP, vals2_RCP,
405 numRows, numCols,
406 rankFilename);
407
408 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
409
410 ArrayRCP<size_t> rowptr_RCP;
411 ArrayRCP<LocalOrdinal> colind_RCP;
412 ArrayRCP<Scalar> vals_RCP;
413 ACrs->allocateAllValues(colind2_RCP.size(), rowptr_RCP, colind_RCP, vals_RCP);
414
415 rowptr_RCP.assign(rowptr2_RCP.begin(), rowptr2_RCP.end());
416 colind_RCP = colind2_RCP;
417 vals_RCP = vals2_RCP;
418
419 ACrs->setAllValues(rowptr_RCP, colind_RCP, vals_RCP);
420 } else {
421 // Custom file format (binary)
422 std::ifstream ifs = std::ifstream(rankFilename.c_str(), std::ios::binary);
423 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
424
425 int m, n, nnz;
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));
429
430 TEUCHOS_ASSERT_EQUALITY(Teuchos::as<int>(rowMap->getLocalNumElements()), m);
431
435
436 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
437
438 ACrs->allocateAllValues(nnz, rowptrRCP, indicesRCP, valuesRCP);
439
440 Teuchos::ArrayView<size_t> rowptr = rowptrRCP();
441 Teuchos::ArrayView<LocalOrdinal> indices = indicesRCP();
442 Teuchos::ArrayView<Scalar> values = valuesRCP();
443
444 bool sorted = true;
445
446 // Read in rowptr
447 for (int i = 0; i < m; i++) {
448 int row, rownnz;
449 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
450 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
451
452 rowptr[row + 1] += rownnz;
453 ifs.seekg(sizeof(int) * rownnz + sizeof(double) * rownnz, ifs.cur);
454 }
455 for (int i = 0; i < m; i++)
456 rowptr[i + 1] += rowptr[i];
457 TEUCHOS_ASSERT(Teuchos::as<int>(rowptr[m]) == nnz);
458
459 // reset to where the data starts
460 ifs.seekg(sizeof(int) * 3, ifs.beg);
461
462 // read in entries
463 for (int i = 0; i < m; i++) {
464 int row, rownnz;
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++) {
469 int index;
470 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
471 indices[ptr] = Teuchos::as<LocalOrdinal>(index);
472 if (j > 0)
473 sorted = sorted & (indices[ptr - 1] < indices[ptr]);
474 ++ptr;
475 }
476 ptr = rowptr[row];
477 for (int j = 0; j < rownnz; j++) {
478 double value;
479 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
480 values[ptr] = Teuchos::as<Scalar>(value);
481 ++ptr;
482 }
483 rowptr[row] += rownnz;
484 }
485 for (int i = m; i > 0; i--)
486 rowptr[i] = rowptr[i - 1];
487 rowptr[0] = 0;
488
489 if (!sorted) {
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]);
494 }
495 }
496
497 ACrs->setAllValues(rowptrRCP, indicesRCP, valuesRCP);
498 }
499
500 if (callFillComplete)
501 A->fillComplete(domainMap, rangeMap);
502 return A;
503}
504
505template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
508 const bool binary) {
509 Xpetra::UnderlyingLib lib = map->lib();
510
511 if (lib == Xpetra::UseTpetra) {
516
517 RCP<const map_type> temp = toTpetra(map);
518 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp, false, false, binary);
520 return rmv;
521 } else {
522 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
523 }
524
525 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
526}
527
528template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
531 const bool binary) {
532 Xpetra::UnderlyingLib lib = map->lib();
533
534 if (lib == Xpetra::UseTpetra) {
539
540 RCP<const map_type> temp = toTpetra(map);
541 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp, false, false, binary);
543 return rmv;
544 } else {
545 throw Exceptions::RuntimeError("Utils::ReadMultiVectorLO : only implemented for Tpetra");
546 }
547
548 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
549}
550
551template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
554 const RCP<const Teuchos::Comm<int>>& comm,
555 const bool binary) {
556 if (lib == Xpetra::UseTpetra) {
559
560 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tMap = reader_type::readMapFile(fileName, comm, false, false, binary);
561 if (tMap.is_null())
562 throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
563
564 return Xpetra::toXpetra(tMap);
565 } else {
566 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
567 }
568
569 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
570}
571
572template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
574 size_t numBlocks = 2; // TODO user parameter?
575
576 std::vector<RCP<const Map>> rangeMapVec;
577 for (size_t row = 0; row < numBlocks; ++row) {
578 RCP<const Map> map = ReadMap("subRangeMap_" + fileName + toString(row) + ".m", lib, comm);
579 rangeMapVec.push_back(map);
580 }
581 RCP<const Map> fullRangeMap = ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
582
583 std::vector<RCP<const Map>> domainMapVec;
584 for (size_t col = 0; col < numBlocks; ++col) {
585 RCP<const Map> map = ReadMap("subDomainMap_" + fileName + toString(col) + ".m", lib, comm);
586 domainMapVec.push_back(map);
587 }
588 RCP<const Map> fullDomainMap = ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
589
590 /*std::vector<RCP<const XpMap> > testRgMapVec;
591 for(size_t r = 0; r < numBlocks; ++r) {
592 RCP<const XpMap> map = ReadMap("rangemap_" + fileName + toString<size_t>(r) + "0.m", lib, comm);
593 testRgMapVec.push_back(map);
594 }
595 std::vector<RCP<const XpMap> > testDoMapVec;
596 for(size_t c = 0; c < numBlocks; ++c) {
597 RCP<const XpMap> map = ReadMap("domainmap_" + fileName + "0" + toString<size_t>(c) + ".m", lib, comm);
598 testDoMapVec.push_back(map);
599 }*/
600
601 // create map extractors
602
603 // range map extractor
604 bool bRangeUseThyraStyleNumbering = false;
605 /*GlobalOrdinal gMinGids = 0;
606 for(size_t v = 0; v < testRgMapVec.size(); ++v) {
607 gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
608 }
609 if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
610 */
611 RCP<const MapExtractor> rangeMapExtractor =
612 Teuchos::rcp(new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
613
614 // domain map extractor
615 bool bDomainUseThyraStyleNumbering = false;
616 /*gMinGids = 0;
617 for(size_t v = 0; v < testDoMapVec.size(); ++v) {
618 gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
619 }
620 if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
621 */
622 RCP<const MapExtractor> domainMapExtractor =
623 Teuchos::rcp(new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
624
625 RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 33));
626
627 // Read all sub-matrices with their maps and place into blocked operator
628 for (size_t row = 0; row < numBlocks; ++row) {
629 for (size_t col = 0; col < numBlocks; ++col) {
630 RCP<const Map> rowSubMap = ReadMap("rowmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
631 RCP<const Map> colSubMap = ReadMap("colmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
632 RCP<const Map> domSubMap = ReadMap("domainmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
633 RCP<const Map> ranSubMap = ReadMap("rangemap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
634 RCP<Matrix> mat = Read(fileName + toString(row) + toString(col) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
635 bOp->setMatrix(row, col, mat);
636 }
637 }
638
639 bOp->fillComplete();
640
641 return bOp;
642} // ReadBlockedCrsMatrix
643
644template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
645template <class T>
647 std::ostringstream buf;
648 buf << what;
649 return buf.str();
650}
651} // namespace Xpetra
void resize(size_type new_size, const value_type &x=value_type())
size_type size() const
void resize(const size_type n, const T &val=T())
void assign(size_type n, const T &val)
iterator begin() const
iterator end() const
size_type size() const
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)
bool is_null() const
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.