Xpetra  Version of the Day
Xpetra_IO.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 
48 #ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
50 
51 #include <fstream>
52 #include "Xpetra_ConfigDefs.hpp"
53 
54 #ifdef HAVE_XPETRA_EPETRA
55 # ifdef HAVE_MPI
56 # include "Epetra_MpiComm.h"
57 # endif
58 #endif
59 
60 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
61 #include <EpetraExt_MatrixMatrix.h>
62 #include <EpetraExt_RowMatrixOut.h>
63 #include <EpetraExt_MultiVectorOut.h>
64 #include <EpetraExt_CrsMatrixIn.h>
65 #include <EpetraExt_MultiVectorIn.h>
66 #include <EpetraExt_BlockMapIn.h>
67 #include <Xpetra_EpetraUtils.hpp>
68 #include <Xpetra_EpetraMultiVector.hpp>
69 #include <EpetraExt_BlockMapOut.h>
70 #endif
71 
72 #ifdef HAVE_XPETRA_TPETRA
73 #include <MatrixMarket_Tpetra.hpp>
74 #include <Tpetra_RowMatrixTransposer.hpp>
75 #include <TpetraExt_MatrixMatrix.hpp>
76 #include <Xpetra_TpetraMultiVector.hpp>
77 #include <Xpetra_TpetraCrsGraph.hpp>
78 #include <Xpetra_TpetraCrsMatrix.hpp>
79 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
80 #endif
81 
82 #ifdef HAVE_XPETRA_EPETRA
83 #include <Xpetra_EpetraMap.hpp>
84 #endif
85 
86 #include "Xpetra_Matrix.hpp"
87 #include "Xpetra_MatrixMatrix.hpp"
88 #include "Xpetra_CrsGraph.hpp"
89 #include "Xpetra_CrsMatrixWrap.hpp"
90 #include "Xpetra_BlockedCrsMatrix.hpp"
91 
92 #include "Xpetra_Map.hpp"
93 #include "Xpetra_StridedMap.hpp"
94 #include "Xpetra_StridedMapFactory.hpp"
95 #include "Xpetra_MapExtractor.hpp"
96 #include "Xpetra_MatrixFactory.hpp"
97 
98 #include <Teuchos_MatrixMarket_Raw_Writer.hpp>
99 #include <string>
100 
101 
102 namespace Xpetra {
103 
104 
105 #ifdef HAVE_XPETRA_EPETRA
106  // This non-member templated function exists so that the matrix-matrix multiply will compile if Epetra, Tpetra, and ML are enabled.
107  template<class SC,class LO,class GO,class NO>
108  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >
109  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap (RCP<Epetra_CrsMatrix> &/* epAB */) {
110  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
111  "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
112  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
113  }
114 
115  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
116  template<>
117  inline RCP<Xpetra::CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_CrsMatrix> &epAB) {
118  typedef double SC;
119  typedef int LO;
120  typedef int GO;
121  typedef Xpetra::EpetraNode NO;
122 
123  RCP<Xpetra::EpetraCrsMatrixT<GO,NO> > tmpC1 = rcp(new Xpetra::EpetraCrsMatrixT<GO,NO>(epAB));
124  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC,LO,GO,NO> >(tmpC1);
125  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > tmpC3 = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(tmpC2));
126 
127  return tmpC3;
128  }
129 
130 
131  template<class SC,class LO,class GO,class NO>
132  RCP<Xpetra::MultiVector<SC,LO,GO,NO> >
133  Convert_Epetra_MultiVector_ToXpetra_MultiVector (RCP<Epetra_MultiVector> &epX) {
134  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
135  "Convert_Epetra_MultiVector_ToXpetra_MultiVector cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
136  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
137  }
138 
139  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
140  template<>
141  inline RCP<Xpetra::MultiVector<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_MultiVector_ToXpetra_MultiVector<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_MultiVector> &epX) {
142  typedef double SC;
143  typedef int LO;
144  typedef int GO;
145  typedef Xpetra::EpetraNode NO;
146 
147  RCP<Xpetra::MultiVector<SC,LO,GO,NO >> tmp = Xpetra::toXpetra<GO,NO>(epX);
148  return tmp;
149  }
150 
151 #endif
152 
157  template <class Scalar,
158  class LocalOrdinal = int,
159  class GlobalOrdinal = LocalOrdinal,
160  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
161  class IO {
162 
163  private:
164 #undef XPETRA_IO_SHORT
165 #include "Xpetra_UseShortNames.hpp"
166 
167  public:
168 
169 #ifdef HAVE_XPETRA_EPETRA
170  // @{
172  /*static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec);
173  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec);
174 
175  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec);
176  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec);
177 
178  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op);
179  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op);
180 
181  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op);
182  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op);*/
183 
184  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
185  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
186  if (xeMap == Teuchos::null)
187  throw Exceptions::BadCast("Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
188  return xeMap->getEpetra_Map();
189  }
190  // @}
191 #endif
192 
193 #ifdef HAVE_XPETRA_TPETRA
194  // @{
196  /*static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec);
197  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec);
198  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec);
199 
200  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec);
201  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec);
202 
203  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op);
204  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op);
205 
206  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op);
207  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op);
208 
209  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op);
210  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op);*/
211 
212 
213  static const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
214  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
215  if (tmp_TMap == Teuchos::null)
216  throw Exceptions::BadCast("Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
217  return tmp_TMap->getTpetra_Map();
218  }
219 #endif
220 
221 
223 
224 
225  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
226  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > tmp_Map = rcpFromRef(M);
227 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
228  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
229  if (tmp_EMap != Teuchos::null) {
230  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
231  if (rv != 0)
232  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
233  return;
234  }
235 #endif // HAVE_XPETRA_EPETRAEXT
236 
237 #ifdef HAVE_XPETRA_TPETRA
238  const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> > &tmp_TMap =
239  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
240  if (tmp_TMap != Teuchos::null) {
241  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
242  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
243  return;
244  }
245 #endif // HAVE_XPETRA_TPETRA
246 
247  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
248 
249  } //Write
250 
252  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
253  std::string mapfile = "map_" + fileName;
254  Write(mapfile, *(vec.getMap()));
255 
256  RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_Vec = Teuchos::rcpFromRef(vec);
257 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
258  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
259  if (tmp_EVec != Teuchos::null) {
260  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
261  if (rv != 0)
262  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
263  return;
264  }
265 #endif // HAVE_XPETRA_EPETRA
266 
267 #ifdef HAVE_XPETRA_TPETRA
268  const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &tmp_TVec =
269  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
270  if (tmp_TVec != Teuchos::null) {
271  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
272  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
273  return;
274  }
275 #endif // HAVE_XPETRA_TPETRA
276 
277  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
278 
279  } //Write
280 
281 
283  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
284 
285  Write("rowmap_" + fileName, *(Op.getRowMap()));
286  if ( !Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
287  Write("domainmap_" + fileName, *(Op.getDomainMap()));
288  if ( !Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
289  Write("rangemap_" + fileName, *(Op.getRangeMap()));
290  if ( !Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps )
291  Write("colmap_" + fileName, *(Op.getColMap()));
292 
293  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
294  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
295  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
296 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
297  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
298  if (tmp_ECrsMtx != Teuchos::null) {
299  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
300  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
301  if (rv != 0)
302  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
303  return;
304  }
305 #endif
306 
307 #ifdef HAVE_XPETRA_TPETRA
308  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& tmp_TCrsMtx =
309  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
310  if (tmp_TCrsMtx != Teuchos::null) {
311  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
312  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
313  return;
314  }
315 #endif // HAVE_XPETRA_TPETRA
316 
317  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
318  } //Write
319 
320 
322  static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
323  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
324  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
325  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
326 
327  ArrayRCP<const size_t> rowptr_RCP;
328  ArrayRCP<LocalOrdinal> rowptr2_RCP;
329  ArrayRCP<const LocalOrdinal> colind_RCP;
330  ArrayRCP<const Scalar> vals_RCP;
331  tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
332 
333  ArrayView<const size_t> rowptr = rowptr_RCP();
334  ArrayView<const LocalOrdinal> colind = colind_RCP();
335  ArrayView<const Scalar> vals = vals_RCP();
336 
337  rowptr2_RCP.resize(rowptr.size());
338  ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
339  for (size_t j = 0; j<rowptr.size(); j++)
340  rowptr2[j] = rowptr[j];
341 
342  Teuchos::MatrixMarket::Raw::Writer<Scalar,LocalOrdinal> writer;
343  writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
344  rowptr2,colind,vals,
345  rowptr.size()-1,Op.getColMap()->getNodeNumElements());
346  } //WriteLocal
347 
348 
363  static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
364  using XpIO = Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
365 
366  // Write all matrix blocks with their maps
367  for (size_t row = 0; row < Op.Rows(); ++row) {
368  for (size_t col = 0; col < Op.Cols(); ++col) {
369  RCP<const Matrix > m = Op.getMatrix(row,col);
370  if(m != Teuchos::null) { // skip empty blocks
371  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
372  "Sub block matrix (" << row << "," << col << ") is not of type CrsMatrixWrap.");
373  XpIO::Write(fileName + toString(row) + toString(col) + ".m", *m, writeAllMaps);
374  }
375  }
376  }
377 
378  // write map information of map extractors
379  RCP<const MapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
380  RCP<const MapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
381 
382  for(size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
383  RCP<const Map> map = rangeMapExtractor->getMap(row);
384  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", *map);
385  }
386  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
387 
388  for(size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
389  RCP<const Map> map = domainMapExtractor->getMap(col);
390  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", *map);
391  }
392  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
393  } // WriteBlockedCrsMatrix
394 
396  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) {
397  if (binary == false) {
398  // Matrix Market file format (ASCII)
399  if (lib == Xpetra::UseEpetra) {
400 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
401  Epetra_CrsMatrix *eA;
402  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
403  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
404  if (rv != 0)
405  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
406 
407  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
408 
409  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A =
410  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
411  return A;
412 #else
413  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
414 #endif
415  } else if (lib == Xpetra::UseTpetra) {
416 #ifdef HAVE_XPETRA_TPETRA
417  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
418 
419  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
420 
421  bool callFillComplete = true;
422 
423  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
424 
425  if (tA.is_null())
426  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
427 
428  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
429  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
430  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
431 
432  return A;
433 #else
434  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
435 #endif
436  } else {
437  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
438  }
439  } else {
440  // Custom file format (binary)
441  std::ifstream ifs(fileName.c_str(), std::ios::binary);
442  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
443  int m, n, nnz;
444  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
445  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
446  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
447 
448  int myRank = comm->getRank();
449 
450  GO indexBase = 0;
451  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
452  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
453  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
454 
455  if (myRank == 0) {
456  Teuchos::Array<GlobalOrdinal> inds;
457  Teuchos::Array<Scalar> vals;
458  // Scan matrix to determine the exact nnz per row.
459  Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
460  for (int i = 0; i < m; i++) {
461  int row, rownnz;
462  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
463  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
464  numEntriesPerRow[i] = rownnz;
465  for (int j = 0; j < rownnz; j++) {
466  int index;
467  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
468  }
469  for (int j = 0; j < rownnz; j++) {
470  double value;
471  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
472  }
473  }
474 
475  A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
476 
477  // Now that nnz per row are known, reread and store the matrix.
478  ifs.seekg(0, ifs.beg); //rewind to beginning of file
479  int junk; //skip header info
480  ifs.read(reinterpret_cast<char*>(&m), sizeof(junk));
481  ifs.read(reinterpret_cast<char*>(&n), sizeof(junk));
482  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(junk));
483  for (int i = 0; i < m; i++) {
484  int row, rownnz;
485  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
486  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
487  inds.resize(rownnz);
488  vals.resize(rownnz);
489  for (int j = 0; j < rownnz; j++) {
490  int index;
491  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
492  inds[j] = Teuchos::as<GlobalOrdinal>(index);
493  }
494  for (int j = 0; j < rownnz; j++) {
495  double value;
496  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
497  vals[j] = Teuchos::as<Scalar>(value);
498  }
499  A->insertGlobalValues(row, inds, vals);
500  }
501  } //if (myRank == 0)
502 
503  A->fillComplete(domainMap, rangeMap);
504 
505  return A;
506  } //if (binary == false) ... else
507 
508  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
509 
510  } //Read()
511 
512 
517  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
518  Read(const std::string& filename,
519  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rowMap,
520  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
521  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
522  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
523  const bool callFillComplete = true,
524  const bool binary = false,
525  const bool tolerant = false,
526  const bool debug = false) {
527  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
528 
529  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
530  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
531 
532  const Xpetra::UnderlyingLib lib = rowMap->lib();
533  if (binary == false) {
534  if (lib == Xpetra::UseEpetra) {
535 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
536  Epetra_CrsMatrix *eA;
537  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
538  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
539  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
540  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
541  int rv;
542  if (colMap.is_null()) {
543  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
544 
545  } else {
546  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
547  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
548  }
549 
550  if (rv != 0)
551  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
552 
553  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
554  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A =
555  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
556 
557  return A;
558 #else
559  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
560 #endif
561  } else if (lib == Xpetra::UseTpetra) {
562 #ifdef HAVE_XPETRA_TPETRA
563  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
564  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
565  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
566 
567  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
568  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
569  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
570  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
571 
572  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
573  callFillComplete, tolerant, debug);
574  if (tA.is_null())
575  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
576 
577  RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tA));
578  RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
579  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tmpA2));
580 
581  return A;
582 #else
583  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
584 #endif
585  } else {
586  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
587  }
588  } else {
589  // Custom file format (binary)
590  std::ifstream ifs(filename.c_str(), std::ios::binary);
591  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
592  int m, n, nnz;
593  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
594  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
595  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
596 
597  //2020-June-05 JHU : for Tpetra, this will probably fail because Tpetra now requires staticly-sized matrix graphs.
598  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap, 1);
599 
600  //2019-06-07 JHU I don't see why this should matter.
601  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
602 
603  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
604  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
605 
606  Teuchos::Array<GlobalOrdinal> inds;
607  Teuchos::Array<Scalar> vals;
608  for (int i = 0; i < m; i++) {
609  int row, rownnz;
610  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
611  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
612  inds.resize(rownnz);
613  vals.resize(rownnz);
614  for (int j = 0; j < rownnz; j++) {
615  int index;
616  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
617  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
618  }
619  for (int j = 0; j < rownnz; j++) {
620  double value;
621  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
622  vals[j] = Teuchos::as<SC>(value);
623  }
624  //This implies that row is not a global index.
625  A->insertGlobalValues(rowElements[row], inds, vals);
626  }
627  A->fillComplete(domainMap, rangeMap);
628  return A;
629  }
630 
631  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
632  }
634 
635 
636  static RCP<MultiVector> ReadMultiVector (const std::string& fileName, const RCP<const Map>& map) {
637  Xpetra::UnderlyingLib lib = map->lib();
638 
639  if (lib == Xpetra::UseEpetra) {
640  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
641 
642  } else if (lib == Xpetra::UseTpetra) {
643 #ifdef HAVE_XPETRA_TPETRA
644  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
645  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
646  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
647  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
648 
649  RCP<const map_type> temp = toTpetra(map);
650  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
651  RCP<MultiVector> rmv = Xpetra::toXpetra(TMV);
652  return rmv;
653 #else
654  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
655 #endif
656  } else {
657  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
658  }
659 
660  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
661  }
662 
663  static RCP<const Map> ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
664  if (lib == Xpetra::UseEpetra) {
665  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
666  } else if (lib == Xpetra::UseTpetra) {
667 #ifdef HAVE_XPETRA_TPETRA
668  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
669  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
670 
671  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
672  if (tMap.is_null())
673  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
674 
675  return Xpetra::toXpetra(tMap);
676 #else
677  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
678 #endif
679  } else {
680  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
681  }
682 
683  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
684 
685  }
686 
700  static RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ReadBlockedCrsMatrix (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
701  using XpIO = Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
702 
703  size_t numBlocks = 2; // TODO user parameter?
704 
705  std::vector<RCP<const Map> > rangeMapVec;
706  for(size_t row = 0; row < numBlocks; ++row) {
707  RCP<const Map> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", lib, comm);
708  rangeMapVec.push_back(map);
709  }
710  RCP<const Map> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
711 
712  std::vector<RCP<const Map> > domainMapVec;
713  for(size_t col = 0; col < numBlocks; ++col) {
714  RCP<const Map> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", lib, comm);
715  domainMapVec.push_back(map);
716  }
717  RCP<const Map> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
718 
719  /*std::vector<RCP<const XpMap> > testRgMapVec;
720  for(size_t r = 0; r < numBlocks; ++r) {
721  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
722  testRgMapVec.push_back(map);
723  }
724  std::vector<RCP<const XpMap> > testDoMapVec;
725  for(size_t c = 0; c < numBlocks; ++c) {
726  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
727  testDoMapVec.push_back(map);
728  }*/
729 
730  // create map extractors
731 
732  // range map extractor
733  bool bRangeUseThyraStyleNumbering = false;
734  /*GlobalOrdinal gMinGids = 0;
735  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
736  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
737  }
738  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
739  */
740  RCP<const MapExtractor> rangeMapExtractor =
741  Teuchos::rcp(new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
742 
743 
744  // domain map extractor
745  bool bDomainUseThyraStyleNumbering = false;
746  /*gMinGids = 0;
747  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
748  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
749  }
750  if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
751  */
752  RCP<const MapExtractor> domainMapExtractor =
753  Teuchos::rcp(new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
754 
755  RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor,33));
756 
757  // Read all sub-matrices with their maps and place into blocked operator
758  for (size_t row = 0; row < numBlocks; ++row) {
759  for (size_t col = 0; col < numBlocks; ++col) {
760  RCP<const Map> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
761  RCP<const Map> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
762  RCP<const Map> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
763  RCP<const Map> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
764  RCP<Matrix> mat = XpIO::Read(fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
765  bOp->setMatrix(row, col, mat);
766  }
767  }
768 
769  bOp->fillComplete();
770 
771  return bOp;
772  } // ReadBlockedCrsMatrix
773 
774 
776  template<class T>
777  static std::string toString(const T& what) {
778  std::ostringstream buf;
779  buf << what;
780  return buf.str();
781  }
782  };
783 
784 
785 #ifdef HAVE_XPETRA_EPETRA
786 
795  template <class Scalar>
796  class IO<Scalar,int,int,EpetraNode> {
797  public:
798  typedef int LocalOrdinal;
799  typedef int GlobalOrdinal;
800  typedef EpetraNode Node;
801 
802 #ifdef HAVE_XPETRA_EPETRA
803  // @{
805  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
806  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
807  if (xeMap == Teuchos::null)
808  throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
809  return xeMap->getEpetra_Map();
810  }
811  // @}
812 #endif
813 
814 #ifdef HAVE_XPETRA_TPETRA
815  // @{
817  static const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
818  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
819  if (tmp_TMap == Teuchos::null)
820  throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
821  return tmp_TMap->getTpetra_Map();
822  }
823 #endif
824 
825 
827 
828 
829  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
830  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > tmp_Map = rcpFromRef(M);
831 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
832  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
833  if (tmp_EMap != Teuchos::null) {
834  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
835  if (rv != 0)
836  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
837  return;
838  }
839 #endif // HAVE_XPETRA_EPETRA
840 
841 #ifdef HAVE_XPETRA_TPETRA
842 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
843  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
844  // do nothing
845 # else
846  const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> > &tmp_TMap =
847  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
848  if (tmp_TMap != Teuchos::null) {
849  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
850  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
851  return;
852  }
853 # endif
854 #endif // HAVE_XPETRA_TPETRA
855  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
856  }
857 
859  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
860  std::string mapfile = "map_" + fileName;
861  Write(mapfile, *(vec.getMap()));
862 
863  RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_Vec = Teuchos::rcpFromRef(vec);
864 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
865  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
866  if (tmp_EVec != Teuchos::null) {
867  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
868  if (rv != 0)
869  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
870  return;
871  }
872 #endif // HAVE_XPETRA_EPETRAEXT
873 
874 #ifdef HAVE_XPETRA_TPETRA
875 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
876  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
877  // do nothin
878 # else
879  const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &tmp_TVec =
880  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
881  if (tmp_TVec != Teuchos::null) {
882  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
883  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
884  return;
885  }
886 # endif
887 #endif // HAVE_XPETRA_TPETRA
888 
889  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
890 
891  } //Write
892 
893 
894 
896  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
897 
898  Write("rowmap_" + fileName, *(Op.getRowMap()));
899  if ( !Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
900  Write("domainmap_" + fileName, *(Op.getDomainMap()));
901  if ( !Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
902  Write("rangemap_" + fileName, *(Op.getRangeMap()));
903  if ( !Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps )
904  Write("colmap_" + fileName, *(Op.getColMap()));
905 
906  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
907  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
908  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
909 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
910  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
911  if (tmp_ECrsMtx != Teuchos::null) {
912  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
913  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
914  if (rv != 0)
915  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
916  return;
917  }
918 #endif // endif HAVE_XPETRA_EPETRA
919 
920 #ifdef HAVE_XPETRA_TPETRA
921 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
922  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
923  // do nothin
924 # else
925  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& tmp_TCrsMtx =
926  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
927  if (tmp_TCrsMtx != Teuchos::null) {
928  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
929  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
930  return;
931  }
932 # endif
933 #endif // HAVE_XPETRA_TPETRA
934 
935  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
936  } //Write
937 
938 
940  static void Write(const std::string& fileName, const Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> & graph, const bool &writeAllMaps = false) {
941 
942  Write("rowmap_" + fileName, *(graph.getRowMap()));
943  if ( !graph.getDomainMap()->isSameAs(*(graph.getRowMap())) || writeAllMaps )
944  Write("domainmap_" + fileName, *(graph.getDomainMap()));
945  if ( !graph.getRangeMap()->isSameAs(*(graph.getRowMap())) || writeAllMaps )
946  Write("rangemap_" + fileName, *(graph.getRangeMap()));
947  if ( !graph.getColMap()->isSameAs(*(graph.getDomainMap())) || writeAllMaps )
948  Write("colmap_" + fileName, *(graph.getColMap()));
949 
950  RCP<const Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > tmp_Graph = rcpFromRef(graph);
951 
952 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
953  const RCP<const Xpetra::EpetraCrsGraphT<GlobalOrdinal,Node> >& tmp_ECrsGraph = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsGraphT<GlobalOrdinal,Node> >(tmp_Graph);
954  if (tmp_ECrsGraph != Teuchos::null) {
955  throw Exceptions::BadCast("Writing not implemented for EpetraCrsGraphT");
956  }
957 #endif // endif HAVE_XPETRA_EPETRA
958 
959 #ifdef HAVE_XPETRA_TPETRA
960 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
961  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
962  // do nothin
963 # else
964  RCP<const Xpetra::TpetraCrsGraph<LocalOrdinal, GlobalOrdinal, Node> > tmp_TCrsGraph =
965  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsGraph<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Graph);
966  if (tmp_TCrsGraph != Teuchos::null) {
967  RCP<const Tpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > G = tmp_TCrsGraph->getTpetra_CrsGraph();
968  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseGraphFile(fileName, G);
969  return;
970  }
971 # endif
972 #endif // HAVE_XPETRA_TPETRA
973 
974  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
975  } //Write
976 
977 
979  static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
980  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
981  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
982  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
983 
984  ArrayRCP<const size_t> rowptr_RCP;
985  ArrayRCP<LocalOrdinal> rowptr2_RCP;
986  ArrayRCP<const LocalOrdinal> colind_RCP;
987  ArrayRCP<const Scalar> vals_RCP;
988  tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
989 
990  ArrayView<const size_t> rowptr = rowptr_RCP();
991  ArrayView<const LocalOrdinal> colind = colind_RCP();
992  ArrayView<const Scalar> vals = vals_RCP();
993 
994  rowptr2_RCP.resize(rowptr.size());
995  ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
996  for (size_t j = 0; j<Teuchos::as<size_t>(rowptr.size()); j++)
997  rowptr2[j] = rowptr[j];
998 
999  Teuchos::MatrixMarket::Raw::Writer<Scalar,LocalOrdinal> writer;
1000  writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
1001  rowptr2,colind,vals,
1002  rowptr.size()-1,Op.getColMap()->getNodeNumElements());
1003  } //WriteLocal
1004 
1021  static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
1022 #include "Xpetra_UseShortNames.hpp"
1023  using XpIO = Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
1024 
1025  // write all matrices with their maps
1026  for (size_t row = 0; row < Op.Rows(); ++row) {
1027  for (size_t col = 0; col < Op.Cols(); ++col) {
1028  RCP<const Matrix> m = Op.getMatrix(row, col);
1029  if(m != Teuchos::null) { // skip empty blocks
1030  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
1031  "Sub block matrix (" << row << "," << col << ") is not of type CrsMatrixWrap.");
1032  XpIO::Write(fileName + toString(row) + toString(col) + ".m", *m, writeAllMaps);
1033  }
1034  }
1035  }
1036 
1037  // write map information of map extractors
1038  RCP<const MapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
1039  RCP<const MapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
1040 
1041  for(size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
1042  RCP<const Map> map = rangeMapExtractor->getMap(row);
1043  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", *map);
1044  }
1045  XpIO::Write("fullRangeMap_" + fileName +".m", *(rangeMapExtractor->getFullMap()));
1046 
1047  for(size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
1048  RCP<const Map> map = domainMapExtractor->getMap(col);
1049  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", *map);
1050  }
1051  XpIO::Write("fullDomainMap_" + fileName+ ".m", *(domainMapExtractor->getFullMap()));
1052  } // WriteBlockedCrsMatrix
1053 
1055  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) {
1056  if (binary == false) {
1057  // Matrix Market file format (ASCII)
1058  if (lib == Xpetra::UseEpetra) {
1059 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1060  Epetra_CrsMatrix *eA;
1061  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
1062  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
1063  if (rv != 0)
1064  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1065 
1066  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1067 
1068  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A =
1069  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1070  return A;
1071 #else
1072  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1073 #endif
1074  } else if (lib == Xpetra::UseTpetra) {
1075 #ifdef HAVE_XPETRA_TPETRA
1076 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1077  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1078  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
1079 # else
1080  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1081 
1082  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1083 
1084  bool callFillComplete = true;
1085 
1086  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
1087 
1088  if (tA.is_null())
1089  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1090 
1091  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
1092  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
1093  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
1094 
1095  return A;
1096 # endif
1097 #else
1098  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1099 #endif
1100  } else {
1101  throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1102  }
1103  } else {
1104  // Custom file format (binary)
1105  std::ifstream ifs(fileName.c_str(), std::ios::binary);
1106  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
1107  int m, n, nnz;
1108  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1109  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1110  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1111 
1112  int myRank = comm->getRank();
1113 
1114  GlobalOrdinal indexBase = 0;
1115  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
1116  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
1117 
1118  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
1119 
1120  if (myRank == 0) {
1121  Teuchos::Array<GlobalOrdinal> inds;
1122  Teuchos::Array<Scalar> vals;
1123  // Scan matrix to determine the exact nnz per row.
1124  Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
1125  for (int i = 0; i < m; i++) {
1126  int row, rownnz;
1127  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1128  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1129  numEntriesPerRow[i] = rownnz;
1130  for (int j = 0; j < rownnz; j++) {
1131  int index;
1132  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1133  }
1134  for (int j = 0; j < rownnz; j++) {
1135  double value;
1136  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1137  }
1138  }
1139 
1140  A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
1141 
1142  // Now that nnz per row are known, reread and store the matrix.
1143  ifs.seekg(0, ifs.beg); //rewind to beginning of file
1144  int junk; //skip header info
1145  ifs.read(reinterpret_cast<char*>(&m), sizeof(junk));
1146  ifs.read(reinterpret_cast<char*>(&n), sizeof(junk));
1147  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(junk));
1148  for (int i = 0; i < m; i++) {
1149  int row, rownnz;
1150  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1151  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1152  inds.resize(rownnz);
1153  vals.resize(rownnz);
1154  for (int j = 0; j < rownnz; j++) {
1155  int index;
1156  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1157  inds[j] = Teuchos::as<GlobalOrdinal>(index);
1158  }
1159  for (int j = 0; j < rownnz; j++) {
1160  double value;
1161  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1162  vals[j] = Teuchos::as<Scalar>(value);
1163  }
1164  A->insertGlobalValues(row, inds, vals);
1165  }
1166  } //if (myRank == 0)
1167 
1168  A->fillComplete(domainMap, rangeMap);
1169 
1170  return A;
1171  }
1172 
1173  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1174 
1175  } //Read()
1176 
1177 
1182  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Read(const std::string& filename,
1183  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rowMap,
1184  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
1185  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
1186  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
1187  const bool callFillComplete = true,
1188  const bool binary = false,
1189  const bool tolerant = false,
1190  const bool debug = false) {
1191  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
1192 
1193  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
1194  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
1195 
1196  const Xpetra::UnderlyingLib lib = rowMap->lib();
1197  if (binary == false) {
1198  if (lib == Xpetra::UseEpetra) {
1199 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1200  Epetra_CrsMatrix *eA;
1201  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
1202  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
1203  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
1204  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
1205  int rv;
1206  if (colMap.is_null()) {
1207  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
1208 
1209  } else {
1210  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
1211  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
1212  }
1213 
1214  if (rv != 0)
1215  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1216 
1217  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1218  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A =
1219  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1220 
1221  return A;
1222 #else
1223  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1224 #endif
1225  } else if (lib == Xpetra::UseTpetra) {
1226 #ifdef HAVE_XPETRA_TPETRA
1227 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1228  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1229  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1230 # else
1231  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1232  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1233  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1234 
1235  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
1236  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
1237  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
1238  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
1239 
1240  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
1241  callFillComplete, tolerant, debug);
1242  if (tA.is_null())
1243  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1244 
1245  RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tA));
1246  RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
1247  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tmpA2));
1248 
1249  return A;
1250 # endif
1251 #else
1252  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1253 #endif
1254  } else {
1255  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1256  }
1257  } else {
1258  // Custom file format (binary)
1259  std::ifstream ifs(filename.c_str(), std::ios::binary);
1260  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
1261  int m, n, nnz;
1262  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1263  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1264  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1265 
1266  //2020-June-05 JHU : for Tpetra, this will probably fail because Tpetra now requires staticly-sized matrix graphs.
1267  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap, 1);
1268 
1269  //2019-06-07 JHU I don't see why this should matter.
1270  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1271 
1272  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
1273  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
1274 
1275  Teuchos::Array<GlobalOrdinal> inds;
1276  Teuchos::Array<Scalar> vals;
1277  for (int i = 0; i < m; i++) {
1278  int row, rownnz;
1279  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1280  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1281  inds.resize(rownnz);
1282  vals.resize(rownnz);
1283  for (int j = 0; j < rownnz; j++) {
1284  int index;
1285  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1286  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
1287  }
1288  for (int j = 0; j < rownnz; j++) {
1289  double value;
1290  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1291  vals[j] = Teuchos::as<Scalar>(value);
1292  }
1293  //This implies that row is not a global index.
1294  A->insertGlobalValues(rowElements[row], inds, vals);
1295  }
1296  A->fillComplete(domainMap, rangeMap);
1297  return A;
1298  }
1299 
1300  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1301  }
1303 
1304 
1305  static RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ReadMultiVector (const std::string& fileName, const RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& map) {
1306  Xpetra::UnderlyingLib lib = map->lib();
1307 
1308  if (lib == Xpetra::UseEpetra) {
1309  // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
1310  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1311 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1312  Epetra_MultiVector * MV;
1313  int rv = EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
1314  if(rv != 0) throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToMultiVector failed");
1315  RCP<Epetra_MultiVector> MVrcp = rcp(MV);
1316  return Convert_Epetra_MultiVector_ToXpetra_MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(MVrcp);
1317 #else
1318  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1319 #endif
1320  } else if (lib == Xpetra::UseTpetra) {
1321 #ifdef HAVE_XPETRA_TPETRA
1322 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1323  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1324  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1325 # else
1326  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1327  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1328  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1329  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
1330 
1331  RCP<const map_type> temp = toTpetra(map);
1332  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
1333  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rmv = Xpetra::toXpetra(TMV);
1334  return rmv;
1335 # endif
1336 #else
1337  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1338 #endif
1339  } else {
1340  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1341  }
1342 
1343  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1344 
1345  }
1346 
1347 
1348  static RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1349  if (lib == Xpetra::UseEpetra) {
1350  // do we need another specialization for <double,int,int> ??
1351  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1352 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1353  Epetra_Map *eMap;
1354  int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
1355  if (rv != 0)
1356  throw Exceptions::RuntimeError("Error reading map from file " + fileName + " with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
1357 
1358  RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
1359  return Xpetra::toXpetra<int,Node>(*eMap1);
1360 #else
1361  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1362 #endif
1363  } else if (lib == Xpetra::UseTpetra) {
1364 #ifdef HAVE_XPETRA_TPETRA
1365 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1366  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1367  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1368 # else
1369  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1370  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1371 
1372  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
1373  if (tMap.is_null())
1374  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1375 
1376  return Xpetra::toXpetra(tMap);
1377 # endif
1378 #else
1379  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1380 #endif
1381  } else {
1382  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1383  }
1384 
1385  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1386 
1387  }
1388 
1404  static RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ReadBlockedCrsMatrix (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1405 #include "Xpetra_UseShortNames.hpp"
1406  using XpIO = Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
1407 
1408  size_t numBlocks = 2; // TODO user parameter?
1409 
1410  std::vector<RCP<const Map> > rangeMapVec;
1411  for(size_t row = 0; row < numBlocks; ++row) {
1412  RCP<const Map> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", lib, comm);
1413  rangeMapVec.push_back(map);
1414  }
1415  RCP<const Map> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
1416 
1417  std::vector<RCP<const Map> > domainMapVec;
1418  for(size_t col = 0; col < numBlocks; ++col) {
1419  RCP<const Map> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", lib, comm);
1420  domainMapVec.push_back(map);
1421  }
1422  RCP<const Map> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
1423 
1424  /*std::vector<RCP<const XpMap> > testRgMapVec;
1425  for(size_t r = 0; r < numBlocks; ++r) {
1426  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
1427  testRgMapVec.push_back(map);
1428  }
1429  std::vector<RCP<const XpMap> > testDoMapVec;
1430  for(size_t c = 0; c < numBlocks; ++c) {
1431  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
1432  testDoMapVec.push_back(map);
1433  }*/
1434 
1435  // create map extractors
1436 
1437  // range map extractor
1438  bool bRangeUseThyraStyleNumbering = false;
1439  /*
1440  GlobalOrdinal gMinGids = 0;
1441  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
1442  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
1443  }
1444  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;*/
1445  RCP<const MapExtractor> rangeMapExtractor =
1446  rcp(new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
1447 
1448  // domain map extractor
1449  bool bDomainUseThyraStyleNumbering = false;
1450  /*gMinGids = 0;
1451  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
1452  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
1453  }
1454  if ( gMinGids==0 && testDoMapVec.size() > 1) bDomainUseThyraStyleNumbering = true;*/
1455  RCP<const MapExtractor> domainMapExtractor =
1456  rcp(new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
1457 
1458  RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 33));
1459 
1460  // Read all matrices with their maps and create the BlockedCrsMatrix
1461  for (size_t row = 0; row < numBlocks; ++row) {
1462  for (size_t col = 0; col < numBlocks; ++col) {
1463  RCP<const Map> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1464  RCP<const Map> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1465  RCP<const Map> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1466  RCP<const Map> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1467  RCP<Matrix> mat = XpIO::Read(fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1468  bOp->setMatrix(row, col, mat);
1469  }
1470  }
1471 
1472  bOp->fillComplete();
1473 
1474  return bOp;
1475  } // ReadBlockedCrsMatrix
1476 
1478  template<class T>
1479  static std::string toString(const T& what) {
1480  std::ostringstream buf;
1481  buf << what;
1482  return buf.str();
1483  }
1484  };
1485 #endif // HAVE_XPETRA_EPETRA
1486 
1487 
1488 } // end namespace Xpetra
1489 
1490 #define XPETRA_IO_SHORT
1491 
1492 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
Xpetra::toString
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
Definition: Xpetra_Utils.cpp:51
Xpetra::IO< Scalar, int, int, EpetraNode >::Write
static void Write(const std::string &fileName, const Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, const bool &writeAllMaps=false)
Save CrsGraph to file in Matrix Market format.
Definition: Xpetra_IO.hpp:940
Xpetra::IO::Read
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)
Read matrix from file in Matrix Market or binary format.
Definition: Xpetra_IO.hpp:518
Xpetra::TpetraMap::getTpetra_Map
RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getTpetra_Map() const
Get the underlying Tpetra map.
Definition: Xpetra_TpetraMap_def.hpp:242
Xpetra::Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &)
Definition: Xpetra_IO.hpp:109
Xpetra::IO< Scalar, int, int, EpetraNode >::Read
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)
Read matrix from file in Matrix Market or binary format.
Definition: Xpetra_IO.hpp:1182
Xpetra::Matrix::getRowMap
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Xpetra_Matrix.hpp:320
GO
GlobalOrdinal GO
Definition: Xpetra_UseShortNamesOrdinal.hpp:139
Xpetra_UseShortNames.hpp
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:540
Xpetra::IO::Write
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:225
Xpetra
Xpetra namespace
Definition: Xpetra_BlockedCrsMatrix.hpp:89
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory.hpp:234
Xpetra_EpetraUtils.hpp
Xpetra::IO< Scalar, int, int, EpetraNode >::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.hpp:979
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:101
LO
LocalOrdinal LO
Definition: Xpetra_UseShortNamesOrdinal.hpp:138
Xpetra::toEpetra
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
Definition: Xpetra_EpetraCrsGraph.cpp:57
Xpetra::IO< Scalar, int, int, EpetraNode >::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.hpp:1404
Xpetra::IO< Scalar, int, int, EpetraNode >::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.hpp:817
Xpetra_Matrix.hpp
Xpetra_MatrixMatrix.hpp
Xpetra_CrsGraph.hpp
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1090
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.hpp:213
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.hpp:330
Xpetra::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: Xpetra_Exceptions.hpp:86
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:78
Xpetra::IO::ReadMap
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
Definition: Xpetra_IO.hpp:663
Xpetra::CrsGraph::getDomainMap
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
Returns the Map associated with the domain of this graph.
Xpetra::IO< Scalar, int, int, EpetraNode >::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.hpp:1055
Xpetra::IO::ReadMultiVector
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
Definition: Xpetra_IO.hpp:636
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::IO< Scalar, int, int, EpetraNode >::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.hpp:1021
Xpetra::TpetraMultiVector::getTpetra_MultiVector
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_MultiVector() const
Get the underlying Tpetra multivector.
Definition: Xpetra_TpetraMultiVector_def.hpp:450
Xpetra::CrsGraph::getRangeMap
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
Returns the Map associated with the domain of this graph.
Xpetra::Convert_Epetra_MultiVector_ToXpetra_MultiVector
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > Convert_Epetra_MultiVector_ToXpetra_MultiVector(RCP< Epetra_MultiVector > &epX)
Definition: Xpetra_IO.hpp:133
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.hpp:363
Xpetra::IO< Scalar, int, int, EpetraNode >::toString
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:1479
Xpetra_EpetraMap.hpp
Xpetra::IO< Scalar, int, int, EpetraNode >::ReadMultiVector
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReadMultiVector(const std::string &fileName, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map)
Definition: Xpetra_IO.hpp:1305
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
Xpetra::toTpetra
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
Definition: Xpetra_TpetraCrsGraph_decl.hpp:425
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition: Xpetra_EpetraCrsGraph.cpp:65
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrix
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
Definition: Xpetra_EpetraCrsMatrix.hpp:239
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.hpp:283
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.hpp:322
Xpetra_EpetraMultiVector.hpp
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.hpp:396
Xpetra_ConfigDefs.hpp
SC
Scalar SC
Definition: Xpetra_UseShortNamesScalar.hpp:182
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.hpp:184
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_fwd.hpp:51
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition: Xpetra_BlockedCrsMatrix.hpp:1340
Xpetra::Operator::getRangeMap
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
NO
Node NO
Definition: Xpetra_UseShortNamesOrdinal.hpp:140
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1093
Xpetra::CrsGraph::getRowMap
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
Returns the Map that describes the row distribution in this graph.
Xpetra::IO< Scalar, int, int, EpetraNode >::ReadMap
static RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
Definition: Xpetra_IO.hpp:1348
Xpetra::IO< Scalar, int, int, EpetraNode >::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.hpp:805
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1304
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1307
Xpetra::TpetraCrsGraph::getTpetra_CrsGraph
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsGraph() const
Get the underlying Tpetra graph.
Definition: Xpetra_TpetraCrsGraph_def.hpp:413
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap_def.hpp:482
Xpetra::IO
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
Definition: Xpetra_IO.hpp:161
Xpetra::IO< Scalar, int, int, EpetraNode >::Write
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:829
Xpetra::IO::Write
static void Write(const std::string &fileName, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save vector to file in Matrix Market format.
Definition: Xpetra_IO.hpp:252
Xpetra::IO< Scalar, int, int, EpetraNode >::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.hpp:896
Xpetra::IO< Scalar, int, int, EpetraNode >::Write
static void Write(const std::string &fileName, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save vector to file in Matrix Market format.
Definition: Xpetra_IO.hpp:859
Xpetra::Operator::getDomainMap
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_fwd.hpp:51
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.hpp:700
Xpetra::IO::toString
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:777
Xpetra::CrsGraph::getColMap
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const =0
Returns the Map that describes the column distribution in this graph.