Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_IO_decl.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#ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
11#define PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
12
13#include <fstream>
14#include "Xpetra_ConfigDefs.hpp"
15
16#ifdef HAVE_XPETRA_EPETRA
17#ifdef HAVE_MPI
18#include "Epetra_MpiComm.h"
19#endif
20#endif
21
22#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
23#include <EpetraExt_MatrixMatrix.h>
24#include <EpetraExt_RowMatrixOut.h>
25#include <EpetraExt_MultiVectorOut.h>
26#include <EpetraExt_CrsMatrixIn.h>
27#include <EpetraExt_MultiVectorIn.h>
28#include <EpetraExt_BlockMapIn.h>
31#include <EpetraExt_BlockMapOut.h>
32#endif
33
34#ifdef HAVE_XPETRA_TPETRA
35#include <MatrixMarket_Tpetra.hpp>
36#include <Tpetra_RowMatrixTransposer.hpp>
37#include <TpetraExt_MatrixMatrix.hpp>
38#include <Xpetra_TpetraMultiVector.hpp>
39#include <Xpetra_TpetraCrsGraph.hpp>
40#include <Xpetra_TpetraCrsMatrix.hpp>
41#include <Xpetra_TpetraBlockCrsMatrix.hpp>
42#include "Tpetra_Util.hpp"
43#endif
44
45#ifdef HAVE_XPETRA_EPETRA
46#include <Xpetra_EpetraMap.hpp>
47#endif
48
49#include "Xpetra_Matrix.hpp"
50#include "Xpetra_MatrixMatrix.hpp"
51#include "Xpetra_Helpers.hpp"
52#include "Xpetra_CrsGraph.hpp"
53#include "Xpetra_CrsMatrixWrap.hpp"
54#include "Xpetra_BlockedCrsMatrix.hpp"
55
56#include "Xpetra_Map.hpp"
57#include "Xpetra_StridedMap.hpp"
58#include "Xpetra_StridedMapFactory.hpp"
59#include "Xpetra_MapExtractor.hpp"
60#include "Xpetra_MatrixFactory.hpp"
61
63#include <Teuchos_MatrixMarket_Raw_Writer.hpp>
64#include <Teuchos_MatrixMarket_Raw_Reader.hpp>
65#include <string>
66
67namespace Xpetra {
68
69#ifdef HAVE_XPETRA_EPETRA
70// This non-member templated function exists so that the matrix-matrix multiply will compile if Epetra, Tpetra, and ML are enabled.
71template <class SC, class LO, class GO, class NO>
72RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO>>
75 "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
76 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
77}
78
79// specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
80template <>
81inline RCP<Xpetra::CrsMatrixWrap<double, int, int, Xpetra::EpetraNode>> Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double, int, int, Xpetra::EpetraNode>(RCP<Epetra_CrsMatrix>& epAB) {
82 typedef double SC;
83 typedef int LO;
84 typedef int GO;
85 typedef Xpetra::EpetraNode NO;
86
87 RCP<Xpetra::EpetraCrsMatrixT<GO, NO>> tmpC1 = rcp(new Xpetra::EpetraCrsMatrixT<GO, NO>(epAB));
88 RCP<Xpetra::CrsMatrix<SC, LO, GO, NO>> tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC, LO, GO, NO>>(tmpC1);
89 RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO>> tmpC3 = rcp(new Xpetra::CrsMatrixWrap<SC, LO, GO, NO>(tmpC2));
90
91 return tmpC3;
92}
93
94template <class SC, class LO, class GO, class NO>
95RCP<Xpetra::MultiVector<SC, LO, GO, NO>>
98 "Convert_Epetra_MultiVector_ToXpetra_MultiVector cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
99 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
100}
101
102// specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
103template <>
104inline RCP<Xpetra::MultiVector<double, int, int, Xpetra::EpetraNode>> Convert_Epetra_MultiVector_ToXpetra_MultiVector<double, int, int, Xpetra::EpetraNode>(RCP<Epetra_MultiVector>& epX) {
105 typedef double SC;
106 typedef int LO;
107 typedef int GO;
108 typedef Xpetra::EpetraNode NO;
109
110 RCP<Xpetra::MultiVector<SC, LO, GO, NO>> tmp = Xpetra::toXpetra<GO, NO>(epX);
111 return tmp;
112}
113
114#endif
115
120template <class Scalar,
121 class LocalOrdinal = int,
122 class GlobalOrdinal = LocalOrdinal,
123 class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
124class IO {
125 private:
126#undef XPETRA_IO_SHORT
128
129 public:
130#ifdef HAVE_XPETRA_EPETRA
132 // @{
133 /*static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec);
134 static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec);
135
136 static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec);
137 static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec);
138
139 static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op);
140 static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op);
141
142 static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op);
143 static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op);*/
144
146 // @}
147#endif
148
149#ifdef HAVE_XPETRA_TPETRA
151 // @{
152 /*static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec);
153 static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec);
154 static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec);
155
156 static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec);
157 static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec);
158
159 static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op);
160 static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op);
161
162 static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op);
163 static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op);
164
165 static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op);
166 static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op);*/
167
169#endif
170
172
173
174 static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& M);
175
177 static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec);
178
179 // CAG:
180 // The class is templated on the usual SC-LO-GO-NO.
181 // Instead of instantiating the entire class using Scalar=LO or GO and then dealing with the headaches that integer-valued CrsMatrix creates, we use these two methods.
182
184 static void WriteLOMV(const std::string& fileName, const Xpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& vec);
185
187 static void WriteGOMV(const std::string& fileName, const Xpetra::MultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& vec);
188
190 static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const bool& writeAllMaps = false);
191
193 static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op);
194
209 static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const bool& writeAllMaps = false);
210
212 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);
213
219 Read(const std::string& filename,
221 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = Teuchos::null,
222 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domainMap = Teuchos::null,
223 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rangeMap = Teuchos::null,
224 const bool callFillComplete = true,
225 const bool binary = false,
226 const bool tolerant = false,
227 const bool debug = false);
228
240 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domainMap = Teuchos::null,
241 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rangeMap = Teuchos::null,
242 const bool callFillComplete = true,
243 const bool binary = false,
244 const bool tolerant = false,
245 const bool debug = false);
247
249 static RCP<MultiVector> ReadMultiVector(const std::string& fileName,
250 const RCP<const Map>& map,
251 const bool binary = false);
252
255 const RCP<const Map>& map,
256 const bool binary = false);
257
258 static RCP<const Map> ReadMap(const std::string& fileName,
260 const RCP<const Teuchos::Comm<int>>& comm,
261 const bool binary = false);
262
277
279 template <class T>
280 static std::string toString(const T& what);
281};
282
283#ifdef HAVE_XPETRA_EPETRA
293template <class Scalar>
294class IO<Scalar, int, int, EpetraNode> {
295 public:
296 typedef int LocalOrdinal;
297 typedef int GlobalOrdinal;
299
300#ifdef HAVE_XPETRA_EPETRA
302 // @{
304 RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>> xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(Teuchos::rcpFromRef(map));
305 if (xeMap == Teuchos::null)
306 throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
307 return xeMap->getEpetra_Map();
308 }
309 // @}
310#endif
311
312#ifdef HAVE_XPETRA_TPETRA
314 // @{
316 const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(rcpFromRef(map));
317 if (tmp_TMap == Teuchos::null)
318 throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
319 return tmp_TMap->getTpetra_Map();
320 }
321#endif
322
324
325
326 static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& M) {
328#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
329 const RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(tmp_Map);
330 if (tmp_EMap != Teuchos::null) {
331 int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
332 if (rv != 0)
333 throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
334 return;
335 }
336#endif // HAVE_XPETRA_EPETRA
337
338#ifdef HAVE_XPETRA_TPETRA
339#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
340 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
341 // do nothing
342#else
344 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(tmp_Map);
345 if (tmp_TMap != Teuchos::null) {
346 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> TMap = tmp_TMap->getTpetra_Map();
347 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeMapFile(fileName, *TMap);
348 return;
349 }
350#endif
351#endif // HAVE_XPETRA_TPETRA
352 throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
353 }
354
355 static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) {
356 std::string mapfile = "map_" + fileName;
357 Write(mapfile, *(vec.getMap()));
358
360#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
361 const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>>& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>>(tmp_Vec);
362 if (tmp_EVec != Teuchos::null) {
363 int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
364 if (rv != 0)
365 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
366 return;
367 }
368#endif // HAVE_XPETRA_EPETRAEXT
369
370#ifdef HAVE_XPETRA_TPETRA
371#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
372 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
373 // do nothin
374#else
376 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
377 if (tmp_TVec != Teuchos::null) {
378 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TVec = tmp_TVec->getTpetra_MultiVector();
379 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
380 return;
381 }
382#endif
383#endif // HAVE_XPETRA_TPETRA
384
385 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
386 }
387
388 static void WriteLOMV(const std::string& fileName, const Xpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& vec) {
389 std::string mapfile = "map_" + fileName;
390 Write(mapfile, *(vec.getMap()));
391
393#ifdef HAVE_XPETRA_TPETRA
395 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
396 if (tmp_TVec != Teuchos::null) {
398 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
399 return;
400 } else
401#endif // HAVE_XPETRA_TPETRA
402 {
403 throw Exceptions::RuntimeError("Xpetra cannot write MV<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
404 }
405
406 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
407 }
408
409 static void WriteGOMV(const std::string& fileName, const Xpetra::MultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& vec) {
410 std::string mapfile = "map_" + fileName;
411 Write(mapfile, *(vec.getMap()));
412
414#ifdef HAVE_XPETRA_TPETRA
416 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>(tmp_Vec);
417 if (tmp_TVec != Teuchos::null) {
419 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>::writeDenseFile(fileName, TVec);
420 return;
421 } else
422#endif // HAVE_XPETRA_TPETRA
423 {
424 throw Exceptions::RuntimeError("Xpetra cannot write MV<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node> when the underlying library is Epetra.");
425 }
426
427 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
428 }
429
430 static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const bool& writeAllMaps = false) {
431 Write("rowmap_" + fileName, *(Op.getRowMap()));
432 if (!Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps)
433 Write("domainmap_" + fileName, *(Op.getDomainMap()));
434 if (!Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps)
435 Write("rangemap_" + fileName, *(Op.getRangeMap()));
436 if (!Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps)
437 Write("colmap_" + fileName, *(Op.getColMap()));
438
442#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
443 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(tmp_CrsMtx);
444 if (tmp_ECrsMtx != Teuchos::null) {
445 RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
446 int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
447 if (rv != 0)
448 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
449 return;
450 }
451#endif // endif HAVE_XPETRA_EPETRA
452
453#ifdef HAVE_XPETRA_TPETRA
454#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
455 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
456 // do nothin
457#else
459 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_CrsMtx);
460 if (tmp_TCrsMtx != Teuchos::null) {
462 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeSparseFile(fileName, A);
463 return;
464 }
466 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_CrsMtx);
467 if (tmp_BlockCrs != Teuchos::null) {
468 std::ofstream outstream(fileName, std::ofstream::out);
469 Teuchos::FancyOStream ofs(Teuchos::rcpFromRef(outstream));
470 tmp_BlockCrs->getTpetra_BlockCrsMatrix()->describe(ofs, Teuchos::VERB_EXTREME);
471 return;
472 }
473
474#endif
475#endif // HAVE_XPETRA_TPETRA
476
477 throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
478 }
479
480 static void Write(const std::string& fileName, const Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>& graph, const bool& writeAllMaps = false) {
481 Write("rowmap_" + fileName, *(graph.getRowMap()));
482 if (!graph.getDomainMap()->isSameAs(*(graph.getRowMap())) || writeAllMaps)
483 Write("domainmap_" + fileName, *(graph.getDomainMap()));
484 if (!graph.getRangeMap()->isSameAs(*(graph.getRowMap())) || writeAllMaps)
485 Write("rangemap_" + fileName, *(graph.getRangeMap()));
486 if (!graph.getColMap()->isSameAs(*(graph.getDomainMap())) || writeAllMaps)
487 Write("colmap_" + fileName, *(graph.getColMap()));
488
490
491#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
492 const RCP<const Xpetra::EpetraCrsGraphT<GlobalOrdinal, Node>>& tmp_ECrsGraph = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsGraphT<GlobalOrdinal, Node>>(tmp_Graph);
493 if (tmp_ECrsGraph != Teuchos::null) {
494 throw Exceptions::BadCast("Writing not implemented for EpetraCrsGraphT");
495 }
496#endif // endif HAVE_XPETRA_EPETRA
497
498#ifdef HAVE_XPETRA_TPETRA
499#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
500 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
501 // do nothin
502#else
504 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsGraph<LocalOrdinal, GlobalOrdinal, Node>>(tmp_Graph);
505 if (tmp_TCrsGraph != Teuchos::null) {
506 RCP<const Tpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>> G = tmp_TCrsGraph->getTpetra_CrsGraph();
507 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::writeSparseGraphFile(fileName, G);
508 return;
509 }
510#endif
511#endif // HAVE_XPETRA_TPETRA
512
513 throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
514 }
515
516 static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
520
521 ArrayRCP<const size_t> rowptr_RCP;
522 ArrayRCP<LocalOrdinal> rowptr2_RCP;
524 ArrayRCP<const Scalar> vals_RCP;
525 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
526
527 ArrayView<const size_t> rowptr = rowptr_RCP();
528 ArrayView<const LocalOrdinal> colind = colind_RCP();
529 ArrayView<const Scalar> vals = vals_RCP();
530
531 rowptr2_RCP.resize(rowptr.size());
532 ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
533 for (size_t j = 0; j < Teuchos::as<size_t>(rowptr.size()); j++)
534 rowptr2[j] = rowptr[j];
535
537 writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
538 rowptr2, colind, vals,
539 rowptr.size() - 1, Op.getColMap()->getLocalNumElements());
540 }
541
542 static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const bool& writeAllMaps = false) {
543 // write all matrices with their maps
544 for (size_t row = 0; row < Op.Rows(); ++row) {
545 for (size_t col = 0; col < Op.Cols(); ++col) {
546 auto m = Op.getMatrix(row, col);
547 if (m != Teuchos::null) { // skip empty blocks
548 const bool cond = Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(m) == Teuchos::null;
551 "Sub block matrix (" << row << "," << col << ") is not of type CrsMatrixWrap.");
552 Write(fileName + toString(row) + toString(col) + ".m", *m, writeAllMaps);
553 }
554 }
555 }
556
557 // write map information of map extractors
558 auto rangeMapExtractor = Op.getRangeMapExtractor();
559 auto domainMapExtractor = Op.getDomainMapExtractor();
560
561 for (size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
562 auto map = rangeMapExtractor->getMap(row);
563 Write("subRangeMap_" + fileName + toString(row) + ".m", *map);
564 }
565 Write("fullRangeMap_" + fileName + ".m", *(rangeMapExtractor->getFullMap()));
566
567 for (size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
568 auto map = domainMapExtractor->getMap(col);
569 Write("subDomainMap_" + fileName + toString(col) + ".m", *map);
570 }
571 Write("fullDomainMap_" + fileName + ".m", *(domainMapExtractor->getFullMap()));
572 }
573
574 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) {
575 if (!binary) {
576 // Matrix Market file format (ASCII)
577 if (lib == Xpetra::UseEpetra) {
578#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
580 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
581 int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
582 if (rv != 0)
583 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
584
585 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
586
588 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
589 return A;
590#else
591 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
592#endif
593 } else if (lib == Xpetra::UseTpetra) {
594#ifdef HAVE_XPETRA_TPETRA
595#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
596 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
597 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
598#else
599 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
600
601 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
602
603 bool callFillComplete = true;
604
605 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
606
607 if (tA.is_null())
608 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
609
611 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmpA1);
613
614 return A;
615#endif
616#else
617 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
618#endif
619 } else {
620 throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
621 }
622 } else {
623 // Custom file format (binary)
624 std::ifstream ifs(fileName.c_str(), std::ios::binary);
625 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
626 int m, n, nnz;
627 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
628 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
629 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
630
631 int myRank = comm->getRank();
632
633 GlobalOrdinal indexBase = 0;
634 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
635 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
636
638
639 if (myRank == 0) {
642 // Scan matrix to determine the exact nnz per row.
643 Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
644 for (int i = 0; i < m; i++) {
645 int row, rownnz;
646 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
647 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
648 numEntriesPerRow[i] = rownnz;
649 for (int j = 0; j < rownnz; j++) {
650 int index;
651 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
652 }
653 for (int j = 0; j < rownnz; j++) {
654 double value;
655 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
656 }
657 }
658
660
661 // Now that nnz per row are known, reread and store the matrix.
662 ifs.seekg(0, ifs.beg); // rewind to beginning of file
663 int junk; // skip header info
664 ifs.read(reinterpret_cast<char*>(&m), sizeof(junk));
665 ifs.read(reinterpret_cast<char*>(&n), sizeof(junk));
666 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(junk));
667 for (int i = 0; i < m; i++) {
668 int row, rownnz;
669 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
670 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
671 inds.resize(rownnz);
672 vals.resize(rownnz);
673 for (int j = 0; j < rownnz; j++) {
674 int index;
675 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
676 inds[j] = Teuchos::as<GlobalOrdinal>(index);
677 }
678 for (int j = 0; j < rownnz; j++) {
679 double value;
680 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
681 vals[j] = Teuchos::as<Scalar>(value);
682 }
683 A->insertGlobalValues(row, inds, vals);
684 }
685 } // if (myRank == 0)
686 else {
687 Teuchos::ArrayRCP<size_t> numEntriesPerRow(0);
689 }
690
691 A->fillComplete(domainMap, rangeMap);
692
693 return A;
694 }
695
696 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
697 }
698
701 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = Teuchos::null,
702 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domainMap = Teuchos::null,
703 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rangeMap = Teuchos::null,
704 const bool callFillComplete = true,
705 const bool binary = false,
706 const bool tolerant = false,
707 const bool debug = false) {
708 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
709
710 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
711 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
712
713 const Xpetra::UnderlyingLib lib = rowMap->lib();
714 if (binary == false) {
715 if (lib == Xpetra::UseEpetra) {
716#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
718 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
720 const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(*domainMap));
721 const Epetra_Map& epetraRangeMap = (rangeMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(*rangeMap));
722 int rv;
723 if (colMap.is_null()) {
724 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
725
726 } else {
727 const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
728 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
729 }
730
731 if (rv != 0)
732 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
733
734 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
736 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
737
738 return A;
739#else
740 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
741#endif
742 } else if (lib == Xpetra::UseTpetra) {
743#ifdef HAVE_XPETRA_TPETRA
744#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
745 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
746 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
747#else
748 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
749 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
750 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
751
752 const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
753 RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
754 const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
755 const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
756
757 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
758 callFillComplete, tolerant, debug);
759 if (tA.is_null())
760 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
761
763 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmpA1);
765
766 return A;
767#endif
768#else
769 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
770#endif
771 } else {
772 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
773 }
774 } else {
775 // Read in on rank 0.
776 auto tempA = Read(filename, lib, rowMap->getComm(), binary);
777
779 auto importer = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(tempA->getRowMap(), rowMap);
780 A->doImport(*tempA, *importer, Xpetra::INSERT);
781 if (callFillComplete)
782 A->fillComplete(domainMap, rangeMap);
783
784 return A;
785 }
786
787 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
788 }
789
793 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domainMap = Teuchos::null,
794 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rangeMap = Teuchos::null,
795 const bool callFillComplete = true,
796 const bool binary = false,
797 const bool tolerant = false,
798 const bool debug = false) {
799 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::ReadLocal() : rowMap cannot be null");
800 TEUCHOS_TEST_FOR_EXCEPTION(colMap.is_null(), Exceptions::RuntimeError, "Utils::ReadLocal() : colMap cannot be null");
801
805
806 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domain = (domainMap.is_null() ? rowMap : domainMap);
807 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> range = (rangeMap.is_null() ? rowMap : rangeMap);
808
809 std::string rankFilename = filename + "." + std::to_string(rowMap->getComm()->getSize()) + "." + std::to_string(rowMap->getComm()->getRank());
810 RCP<matrix_type> A = rcp(new crs_wrap_type(rowMap, colMap, 0));
811
812 if (binary == false) {
814 params->set("Parse tolerantly", tolerant);
815 params->set("Debug mode", debug);
816
817 LocalOrdinal numRows = rowMap->getLocalNumElements();
818 LocalOrdinal numCols = colMap->getLocalNumElements();
819
820 ArrayRCP<LocalOrdinal> rowptr2_RCP;
821 ArrayRCP<LocalOrdinal> colind2_RCP;
822 ArrayRCP<Scalar> vals2_RCP;
823
825 reader.readFile(rowptr2_RCP, colind2_RCP, vals2_RCP,
826 numRows, numCols,
827 rankFilename);
828
829 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
830
831 ArrayRCP<size_t> rowptr_RCP;
832 ArrayRCP<LocalOrdinal> colind_RCP;
833 ArrayRCP<Scalar> vals_RCP;
834 ACrs->allocateAllValues(colind2_RCP.size(), rowptr_RCP, colind_RCP, vals_RCP);
835
836 rowptr_RCP.assign(rowptr2_RCP.begin(), rowptr2_RCP.end());
837 colind_RCP = colind2_RCP;
838 vals_RCP = vals2_RCP;
839
840 ACrs->setAllValues(rowptr_RCP, colind_RCP, vals_RCP);
841 } else {
842 // Custom file format (binary)
843 std::ifstream ifs = std::ifstream(rankFilename.c_str(), std::ios::binary);
844 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
845
846 int m, n, nnz;
847 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
848 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
849 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
850
851 TEUCHOS_ASSERT_EQUALITY(Teuchos::as<int>(rowMap->getLocalNumElements()), m);
852
856
857 RCP<crs_type> ACrs = Teuchos::rcp_dynamic_cast<crs_wrap_type>(A)->getCrsMatrix();
858
859 ACrs->allocateAllValues(nnz, rowptrRCP, indicesRCP, valuesRCP);
860
861 Teuchos::ArrayView<size_t> rowptr = rowptrRCP();
862 Teuchos::ArrayView<LocalOrdinal> indices = indicesRCP();
863 Teuchos::ArrayView<Scalar> values = valuesRCP();
864
865 bool sorted = true;
866
867 // Read in rowptr
868 for (int i = 0; i < m; i++) {
869 int row, rownnz;
870 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
871 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
872
873 rowptr[row + 1] += rownnz;
874 ifs.seekg(sizeof(int) * rownnz + sizeof(double) * rownnz, ifs.cur);
875 }
876 for (int i = 0; i < m; i++)
877 rowptr[i + 1] += rowptr[i];
878 TEUCHOS_ASSERT(Teuchos::as<int>(rowptr[m]) == nnz);
879
880 // reset to where the data starts
881 ifs.seekg(sizeof(int) * 3, ifs.beg);
882
883 // read in entries
884 for (int i = 0; i < m; i++) {
885 int row, rownnz;
886 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
887 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
888 size_t ptr = rowptr[row];
889 for (int j = 0; j < rownnz; j++) {
890 int index;
891 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
892 indices[ptr] = Teuchos::as<LocalOrdinal>(index);
893 if (j > 0)
894 sorted = sorted & (indices[ptr - 1] < indices[ptr]);
895 ++ptr;
896 }
897 ptr = rowptr[row];
898 for (int j = 0; j < rownnz; j++) {
899 double value;
900 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
901 values[ptr] = Teuchos::as<Scalar>(value);
902 ++ptr;
903 }
904 rowptr[row] += rownnz;
905 }
906 for (int i = m; i > 0; i--)
907 rowptr[i] = rowptr[i - 1];
908 rowptr[0] = 0;
909
910#ifdef HAVE_XPETRA_TPETRA
911 if (!sorted) {
912 for (LocalOrdinal lclRow = 0; lclRow < m; lclRow++) {
913 size_t rowBegin = rowptr[lclRow];
914 size_t rowEnd = rowptr[lclRow + 1];
915 Tpetra::sort2(&indices[rowBegin], &indices[rowEnd], &values[rowBegin]);
916 }
917 }
918#else
919 TEUCHOS_ASSERT(sorted);
920#endif
921
922 ACrs->setAllValues(rowptrRCP, indicesRCP, valuesRCP);
923 }
924
925 if (callFillComplete)
926 A->fillComplete(domainMap, rangeMap);
927 return A;
928 }
929
932 const bool binary = false) {
933 Xpetra::UnderlyingLib lib = map->lib();
934
935 if (lib == Xpetra::UseEpetra) {
936 // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
937 // TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
938#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
939 TEUCHOS_ASSERT(!binary);
941 int rv = EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
942 if (rv != 0) throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToMultiVector failed");
943 RCP<Epetra_MultiVector> MVrcp = rcp(MV);
944 return Convert_Epetra_MultiVector_ToXpetra_MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(MVrcp);
945#else
946 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
947#endif
948 } else if (lib == Xpetra::UseTpetra) {
949#ifdef HAVE_XPETRA_TPETRA
950#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
951 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
952 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
953#else
954 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
955 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
956 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
957 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
958
959 RCP<const map_type> temp = toTpetra(map);
960 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp, false, false, binary);
962 return rmv;
963#endif
964#else
965 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
966#endif
967 } else {
968 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
969 }
970
971 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
972 }
973
976 const bool binary = false) {
977 Xpetra::UnderlyingLib lib = map->lib();
978
979 if (lib == Xpetra::UseTpetra) {
980#ifdef HAVE_XPETRA_TPETRA
981 typedef Tpetra::CrsMatrix<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
982 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
983 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
984 typedef Tpetra::MultiVector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
985
986 RCP<const map_type> temp = toTpetra(map);
987 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName, map->getComm(), temp, false, false, binary);
989 return rmv;
990#else
991 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
992#endif
993 } else {
994 throw Exceptions::RuntimeError("Utils::ReadMultiVectorLO : only implemented for Tpetra");
995 }
996
997 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
998 }
999
1002 const RCP<const Teuchos::Comm<int>>& comm,
1003 const bool binary = false) {
1004 if (lib == Xpetra::UseEpetra) {
1005 // do we need another specialization for <double,int,int> ??
1006 // TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1007#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1008 TEUCHOS_ASSERT(!binary);
1009 Epetra_Map* eMap;
1010 int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
1011 if (rv != 0)
1012 throw Exceptions::RuntimeError("Error reading map from file " + fileName + " with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
1013
1014 RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
1015 return Xpetra::toXpetra<int, Node>(*eMap1);
1016#else
1017 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1018#endif
1019 } else if (lib == Xpetra::UseTpetra) {
1020#ifdef HAVE_XPETRA_TPETRA
1021#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1022 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1023 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1024#else
1025 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1026 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1027
1028 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tMap = reader_type::readMapFile(fileName, comm, false, false, binary);
1029 if (tMap.is_null())
1030 throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1031
1032 return Xpetra::toXpetra(tMap);
1033#endif
1034#else
1035 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1036#endif
1037 } else {
1038 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1039 }
1040
1041 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1042 }
1043
1045 size_t numBlocks = 2; // TODO user parameter?
1046
1047 std::vector<RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>> rangeMapVec;
1048 for (size_t row = 0; row < numBlocks; ++row) {
1049 auto map = ReadMap("subRangeMap_" + fileName + toString(row) + ".m", lib, comm);
1050 rangeMapVec.push_back(map);
1051 }
1052 auto fullRangeMap = ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
1053
1054 std::vector<RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>> domainMapVec;
1055 for (size_t col = 0; col < numBlocks; ++col) {
1056 auto map = ReadMap("subDomainMap_" + fileName + toString(col) + ".m", lib, comm);
1057 domainMapVec.push_back(map);
1058 }
1059 auto fullDomainMap = ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
1060
1061 /*std::vector<RCP<const XpMap> > testRgMapVec;
1062 for(size_t r = 0; r < numBlocks; ++r) {
1063 RCP<const XpMap> map = ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
1064 testRgMapVec.push_back(map);
1065 }
1066 std::vector<RCP<const XpMap> > testDoMapVec;
1067 for(size_t c = 0; c < numBlocks; ++c) {
1068 RCP<const XpMap> map = ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
1069 testDoMapVec.push_back(map);
1070 }*/
1071
1072 // create map extractors
1073
1074 // range map extractor
1075 bool bRangeUseThyraStyleNumbering = false;
1076 /*
1077 GlobalOrdinal gMinGids = 0;
1078 for(size_t v = 0; v < testRgMapVec.size(); ++v) {
1079 gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
1080 }
1081 if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;*/
1083 rcp(new Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node>(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
1084
1085 // domain map extractor
1086 bool bDomainUseThyraStyleNumbering = false;
1087 /*gMinGids = 0;
1088 for(size_t v = 0; v < testDoMapVec.size(); ++v) {
1089 gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
1090 }
1091 if ( gMinGids==0 && testDoMapVec.size() > 1) bDomainUseThyraStyleNumbering = true;*/
1093 rcp(new Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node>(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
1094
1095 auto bOp = Teuchos::rcp(new Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(rangeMapExtractor, domainMapExtractor, 33));
1096
1097 // Read all matrices with their maps and create the BlockedCrsMatrix
1098 for (size_t row = 0; row < numBlocks; ++row) {
1099 for (size_t col = 0; col < numBlocks; ++col) {
1100 auto rowSubMap = ReadMap("rowmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
1101 auto colSubMap = ReadMap("colmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
1102 auto domSubMap = ReadMap("domainmap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
1103 auto ranSubMap = ReadMap("rangemap_" + fileName + toString(row) + toString(col) + ".m", lib, comm);
1104 auto mat = Read(fileName + toString(row) + toString(col) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1105 bOp->setMatrix(row, col, mat);
1106 }
1107 }
1108
1109 bOp->fillComplete();
1110
1111 return bOp;
1112 }
1113
1115 template <class T>
1116 static std::string toString(const T& what) {
1117 std::ostringstream buf;
1118 buf << what;
1119 return buf.str();
1120 }
1121};
1122
1123#endif // HAVE_XPETRA_EPETRA
1124
1125} // end namespace Xpetra
1126
1127#define XPETRA_IO_SHORT
1128
1129#endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
LocalOrdinal LO
GlobalOrdinal GO
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
void resize(size_type new_size, const value_type &x=value_type())
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.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
Returns the Map associated with the domain of this graph.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
Returns the Map that describes the row distribution in this graph.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
Returns the Map associated with the domain of this graph.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const =0
Returns the Map that describes the column distribution in this graph.
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, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > rowMap, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > colMap=Teuchos::null, 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)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadMultiVector(const std::string &fileName, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map, const bool binary=false)
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)
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
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.
static void WriteLocal(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
static void Write(const std::string &fileName, const Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, const bool &writeAllMaps=false)
static void Write(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
static RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadBlockedCrsMatrix(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
static void WriteGOMV(const std::string &fileName, const Xpetra::MultiVector< GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &vec)
static RCP< Xpetra::MultiVector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReadMultiVectorLO(const std::string &fileName, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map, const bool binary=false)
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
static void WriteLOMV(const std::string &fileName, const Xpetra::MultiVector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &vec)
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)
static RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm, const bool binary=false)
static void Write(const std::string &fileName, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
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 const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
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)
std::string toString(const T &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > Convert_Epetra_MultiVector_ToXpetra_MultiVector(RCP< Epetra_MultiVector > &epX)
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
RCP< Xpetra::MultiVector< double, int, int, Xpetra::EpetraNode > > Convert_Epetra_MultiVector_ToXpetra_MultiVector< double, int, int, Xpetra::EpetraNode >(RCP< Epetra_MultiVector > &epX)
RCP< Xpetra::CrsMatrixWrap< double, int, int, Xpetra::EpetraNode > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap< double, int, int, Xpetra::EpetraNode >(RCP< Epetra_CrsMatrix > &epAB)
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &)