Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
63 namespace Xpetra {
64 
65 #ifdef HAVE_XPETRA_EPETRA
66 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
67 const 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
76 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77 const 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 
85 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86 void 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 
111 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112 void 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 
140 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141 void 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 
162 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
163 void 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 
184 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
185 void 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 
230 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
231 void 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 
257 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
258 void 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 
288 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
289 Teuchos::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 
408 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
409 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
410 IO<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 
496 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
497 Teuchos::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 
637 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
638 RCP<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 
667 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
668 RCP<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 
694 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
695 RCP<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 
721 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
722 RCP<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 
793 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
794 template <class T>
795 std::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::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::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap_def.hpp:433
Xpetra::toString
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
Definition: Xpetra_Utils.cpp:15
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::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix_def.hpp:873
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_EpetraUtils.hpp
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::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:65
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::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix_def.hpp:1067
Xpetra_CrsGraph.hpp
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::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: Xpetra_Exceptions.hpp:50
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::DistObject< Scalar, LocalOrdinal, GlobalOrdinal, Node >::getMap
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
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:29
Xpetra::toTpetra
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
Definition: Xpetra_TpetraCrsGraph_decl.hpp:372
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrix
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
Definition: Xpetra_EpetraCrsMatrix.hpp:207
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix_def.hpp:1061
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...
Xpetra::TpetraMap::getTpetra_Map
RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getTpetra_Map() const
Get the underlying Tpetra map.
Definition: Xpetra_TpetraMap_def.hpp:238
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix_def.hpp:867
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::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition: Xpetra_BlockedCrsMatrix_def.hpp:1103
Xpetra_EpetraMap.hpp
Xpetra_IO_decl.hpp
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory_def.hpp:22
Xpetra::TpetraCrsMatrix::getTpetra_CrsMatrix
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
Get the underlying Tpetra matrix.
Definition: Xpetra_TpetraCrsMatrix_def.hpp:611
Xpetra_EpetraMultiVector.hpp
Xpetra_ConfigDefs.hpp
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_decl.hpp:46
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::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::TpetraMultiVector::getTpetra_MultiVector
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_MultiVector() const
Get the underlying Tpetra multivector.
Definition: Xpetra_TpetraMultiVector_def.hpp:519
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::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::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::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::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::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::toEpetra
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
Definition: Xpetra_EpetraCrsGraph.cpp:21
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_decl.hpp:60
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::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::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition: Xpetra_EpetraCrsGraph.cpp:29