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_TpetraCrsMatrix.hpp>
78 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
79 #endif
80 
81 #ifdef HAVE_XPETRA_EPETRA
82 #include <Xpetra_EpetraMap.hpp>
83 #endif
84 
85 #include "Xpetra_Matrix.hpp"
86 #include "Xpetra_MatrixMatrix.hpp"
87 #include "Xpetra_CrsMatrixWrap.hpp"
88 #include "Xpetra_BlockedCrsMatrix.hpp"
89 
90 #include "Xpetra_Map.hpp"
91 #include "Xpetra_StridedMap.hpp"
92 #include "Xpetra_StridedMapFactory.hpp"
93 #include "Xpetra_MapExtractor.hpp"
94 #include "Xpetra_MatrixFactory.hpp"
95 
96 #include <Teuchos_MatrixMarket_Raw_Writer.hpp>
97 #include <string>
98 
99 
100 namespace Xpetra {
101 
102 
103 #ifdef HAVE_XPETRA_EPETRA
104  // This non-member templated function exists so that the matrix-matrix multiply will compile if Epetra, Tpetra, and ML are enabled.
105  template<class SC,class LO,class GO,class NO>
106  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >
107  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap (RCP<Epetra_CrsMatrix> &/* epAB */) {
108  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
109  "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
110  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
111  }
112 
113  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
114  template<>
115  inline RCP<Xpetra::CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_CrsMatrix> &epAB) {
116  typedef double SC;
117  typedef int LO;
118  typedef int GO;
119  typedef Xpetra::EpetraNode NO;
120 
121  RCP<Xpetra::EpetraCrsMatrixT<GO,NO> > tmpC1 = rcp(new Xpetra::EpetraCrsMatrixT<GO,NO>(epAB));
122  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC,LO,GO,NO> >(tmpC1);
123  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > tmpC3 = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(tmpC2));
124 
125  return tmpC3;
126  }
127 
128 
129  template<class SC,class LO,class GO,class NO>
130  RCP<Xpetra::MultiVector<SC,LO,GO,NO> >
131  Convert_Epetra_MultiVector_ToXpetra_MultiVector (RCP<Epetra_MultiVector> &epX) {
132  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
133  "Convert_Epetra_MultiVector_ToXpetra_MultiVector cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
134  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
135  }
136 
137  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
138  template<>
139  inline RCP<Xpetra::MultiVector<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_MultiVector_ToXpetra_MultiVector<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_MultiVector> &epX) {
140  typedef double SC;
141  typedef int LO;
142  typedef int GO;
143  typedef Xpetra::EpetraNode NO;
144 
145  RCP<Xpetra::MultiVector<SC,LO,GO,NO >> tmp = Xpetra::toXpetra<GO,NO>(epX);
146  return tmp;
147  }
148 
149 #endif
150 
155  template <class Scalar,
156  class LocalOrdinal = int,
157  class GlobalOrdinal = LocalOrdinal,
158  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
159  class IO {
160 
161  private:
162 #undef XPETRA_IO_SHORT
163 #include "Xpetra_UseShortNames.hpp"
164 
165  public:
166 
167 #ifdef HAVE_XPETRA_EPETRA
169  // @{
170  /*static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec);
171  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec);
172 
173  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec);
174  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec);
175 
176  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op);
177  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op);
178 
179  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op);
180  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op);*/
181 
182  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
183  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
184  if (xeMap == Teuchos::null)
185  throw Exceptions::BadCast("Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
186  return xeMap->getEpetra_Map();
187  }
188  // @}
189 #endif
190 
191 #ifdef HAVE_XPETRA_TPETRA
193  // @{
194  /*static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec);
195  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec);
196  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec);
197 
198  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec);
199  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec);
200 
201  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op);
202  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op);
203 
204  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op);
205  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op);
206 
207  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op);
208  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op);*/
209 
210 
211  static const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
212  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
213  if (tmp_TMap == Teuchos::null)
214  throw Exceptions::BadCast("Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
215  return tmp_TMap->getTpetra_Map();
216  }
217 #endif
218 
219 
221 
222 
223  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
224  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > tmp_Map = rcpFromRef(M);
225 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
226  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
227  if (tmp_EMap != Teuchos::null) {
228  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
229  if (rv != 0)
230  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
231  return;
232  }
233 #endif // HAVE_XPETRA_EPETRAEXT
234 
235 #ifdef HAVE_XPETRA_TPETRA
236  const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> > &tmp_TMap =
237  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
238  if (tmp_TMap != Teuchos::null) {
239  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
240  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
241  return;
242  }
243 #endif // HAVE_XPETRA_TPETRA
244 
245  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
246 
247  } //Write
248 
250  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
251  std::string mapfile = "map_" + fileName;
252  Write(mapfile, *(vec.getMap()));
253 
254  RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_Vec = Teuchos::rcpFromRef(vec);
255 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
256  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
257  if (tmp_EVec != Teuchos::null) {
258  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
259  if (rv != 0)
260  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
261  return;
262  }
263 #endif // HAVE_XPETRA_EPETRA
264 
265 #ifdef HAVE_XPETRA_TPETRA
266  const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &tmp_TVec =
267  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
268  if (tmp_TVec != Teuchos::null) {
269  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
270  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
271  return;
272  }
273 #endif // HAVE_XPETRA_TPETRA
274 
275  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
276 
277  } //Write
278 
279 
281  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
282 
283  Write("rowmap_" + fileName, *(Op.getRowMap()));
284  if ( !Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
285  Write("domainmap_" + fileName, *(Op.getDomainMap()));
286  if ( !Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
287  Write("rangemap_" + fileName, *(Op.getRangeMap()));
288  if ( !Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps )
289  Write("colmap_" + fileName, *(Op.getColMap()));
290 
291  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
292  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
293  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
294 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
295  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
296  if (tmp_ECrsMtx != Teuchos::null) {
297  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
298  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
299  if (rv != 0)
300  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
301  return;
302  }
303 #endif
304 
305 #ifdef HAVE_XPETRA_TPETRA
306  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& tmp_TCrsMtx =
307  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
308  if (tmp_TCrsMtx != Teuchos::null) {
309  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
310  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
311  return;
312  }
313 #endif // HAVE_XPETRA_TPETRA
314 
315  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
316  } //Write
317 
318 
320  static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
321  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
322  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
323  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
324 
325  ArrayRCP<const size_t> rowptr_RCP;
326  ArrayRCP<LocalOrdinal> rowptr2_RCP;
327  ArrayRCP<const LocalOrdinal> colind_RCP;
328  ArrayRCP<const Scalar> vals_RCP;
329  tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
330 
331  ArrayView<const size_t> rowptr = rowptr_RCP();
332  ArrayView<const LocalOrdinal> colind = colind_RCP();
333  ArrayView<const Scalar> vals = vals_RCP();
334 
335  rowptr2_RCP.resize(rowptr.size());
336  ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
337  for (size_t j = 0; j<rowptr.size(); j++)
338  rowptr2[j] = rowptr[j];
339 
340  Teuchos::MatrixMarket::Raw::Writer<Scalar,LocalOrdinal> writer;
341  writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
342  rowptr2,colind,vals,
343  rowptr.size()-1,Op.getColMap()->getNodeNumElements());
344  } //WriteLocal
345 
346 
348  static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
349  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
350  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMapExtractor;
351  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
352  typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
353  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
354  typedef Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpIO;
355 
356  // write all matrices with their maps
357  for (size_t r = 0; r < Op.Rows(); ++r) {
358  for (size_t c = 0; c < Op.Cols(); ++c) {
359  RCP<const XpMat > m = Op.getMatrix(r,c);
360  if(m != Teuchos::null) { // skip empty blocks
361  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const XpCrsMatWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
362  "Sub block matrix (" << r << "," << c << ") is not of type CrsMatrixWrap.");
363  XpIO::Write(fileName + toString(r) + toString(c) + ".m", *m, writeAllMaps);
364  }
365  }
366  }
367 
368  // write map information of map extractors
369  RCP<const XpMapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
370  RCP<const XpMapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
371 
372  for(size_t r = 0; r < rangeMapExtractor->NumMaps(); ++r) {
373  RCP<const XpMap> map = rangeMapExtractor->getMap(r);
374  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", *map);
375  }
376  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
377 
378  for(size_t c = 0; c < domainMapExtractor->NumMaps(); ++c) {
379  RCP<const XpMap> map = domainMapExtractor->getMap(c);
380  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", *map);
381  }
382  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
383  } //WriteBlockCrsMatrix
384 
386  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) {
387  if (binary == false) {
388  // Matrix Market file format (ASCII)
389  if (lib == Xpetra::UseEpetra) {
390 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
391  Epetra_CrsMatrix *eA;
392  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
393  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
394  if (rv != 0)
395  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
396 
397  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
398 
399  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A =
400  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
401  return A;
402 #else
403  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
404 #endif
405  } else if (lib == Xpetra::UseTpetra) {
406 #ifdef HAVE_XPETRA_TPETRA
407  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
408 
409  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
410 
411  bool callFillComplete = true;
412 
413  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
414 
415  if (tA.is_null())
416  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
417 
418  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
419  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
420  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
421 
422  return A;
423 #else
424  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
425 #endif
426  } else {
427  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
428  }
429  } else {
430  // Custom file format (binary)
431  std::ifstream ifs(fileName.c_str(), std::ios::binary);
432  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
433  int m, n, nnz;
434  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
435  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
436  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
437 
438  int myRank = comm->getRank();
439 
440  GO indexBase = 0;
441  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
442  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
443  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
444 
445  if (myRank == 0) {
446  Teuchos::Array<GlobalOrdinal> inds;
447  Teuchos::Array<Scalar> vals;
448  // Scan matrix to determine the exact nnz per row.
449  Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
450  for (int i = 0; i < m; i++) {
451  int row, rownnz;
452  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
453  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
454  numEntriesPerRow[i] = rownnz;
455  for (int j = 0; j < rownnz; j++) {
456  int index;
457  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
458  }
459  for (int j = 0; j < rownnz; j++) {
460  double value;
461  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
462  }
463  }
464 
465  A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
466 
467  // Now that nnz per row are known, reread and store the matrix.
468  ifs.seekg(0, ifs.beg); //rewind to beginning of file
469  int junk; //skip header info
470  ifs.read(reinterpret_cast<char*>(&m), sizeof(junk));
471  ifs.read(reinterpret_cast<char*>(&n), sizeof(junk));
472  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(junk));
473  for (int i = 0; i < m; i++) {
474  int row, rownnz;
475  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
476  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
477  inds.resize(rownnz);
478  vals.resize(rownnz);
479  for (int j = 0; j < rownnz; j++) {
480  int index;
481  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
482  inds[j] = Teuchos::as<GlobalOrdinal>(index);
483  }
484  for (int j = 0; j < rownnz; j++) {
485  double value;
486  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
487  vals[j] = Teuchos::as<Scalar>(value);
488  }
489  A->insertGlobalValues(row, inds, vals);
490  }
491  } //if (myRank == 0)
492 
493  A->fillComplete(domainMap, rangeMap);
494 
495  return A;
496  } //if (binary == false) ... else
497 
498  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
499 
500  } //Read()
501 
502 
507  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
508  Read(const std::string& filename,
509  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rowMap,
510  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
511  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
512  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
513  const bool callFillComplete = true,
514  const bool binary = false,
515  const bool tolerant = false,
516  const bool debug = false) {
517  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
518 
519  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
520  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
521 
522  const Xpetra::UnderlyingLib lib = rowMap->lib();
523  if (binary == false) {
524  if (lib == Xpetra::UseEpetra) {
525 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
526  Epetra_CrsMatrix *eA;
527  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
528  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
529  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
530  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
531  int rv;
532  if (colMap.is_null()) {
533  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
534 
535  } else {
536  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
537  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
538  }
539 
540  if (rv != 0)
541  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
542 
543  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
544  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A =
545  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
546 
547  return A;
548 #else
549  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
550 #endif
551  } else if (lib == Xpetra::UseTpetra) {
552 #ifdef HAVE_XPETRA_TPETRA
553  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
554  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
555  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
556 
557  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
558  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
559  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
560  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
561 
562  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
563  callFillComplete, tolerant, debug);
564  if (tA.is_null())
565  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
566 
567  RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tA));
568  RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
569  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tmpA2));
570 
571  return A;
572 #else
573  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
574 #endif
575  } else {
576  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
577  }
578  } else {
579  // Custom file format (binary)
580  std::ifstream ifs(filename.c_str(), std::ios::binary);
581  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
582  int m, n, nnz;
583  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
584  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
585  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
586 
587  //2020-June-05 JHU : for Tpetra, this will probably fail because Tpetra now requires staticly-sized matrix graphs.
588  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap, 1);
589 
590  //2019-06-07 JHU I don't see why this should matter.
591  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
592 
593  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
594  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
595 
596  Teuchos::Array<GlobalOrdinal> inds;
597  Teuchos::Array<Scalar> vals;
598  for (int i = 0; i < m; i++) {
599  int row, rownnz;
600  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
601  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
602  inds.resize(rownnz);
603  vals.resize(rownnz);
604  for (int j = 0; j < rownnz; j++) {
605  int index;
606  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
607  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
608  }
609  for (int j = 0; j < rownnz; j++) {
610  double value;
611  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
612  vals[j] = Teuchos::as<SC>(value);
613  }
614  //This implies that row is not a global index.
615  A->insertGlobalValues(rowElements[row], inds, vals);
616  }
617  A->fillComplete(domainMap, rangeMap);
618  return A;
619  }
620 
621  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
622  }
624 
625 
626  static RCP<MultiVector> ReadMultiVector (const std::string& fileName, const RCP<const Map>& map) {
627  Xpetra::UnderlyingLib lib = map->lib();
628 
629  if (lib == Xpetra::UseEpetra) {
630  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
631 
632  } else if (lib == Xpetra::UseTpetra) {
633 #ifdef HAVE_XPETRA_TPETRA
634  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
635  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
636  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
637  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
638 
639  RCP<const map_type> temp = toTpetra(map);
640  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
641  RCP<MultiVector> rmv = Xpetra::toXpetra(TMV);
642  return rmv;
643 #else
644  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
645 #endif
646  } else {
647  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
648  }
649 
650  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
651  }
652 
653  static RCP<const Map> ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
654  if (lib == Xpetra::UseEpetra) {
655  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
656  } else if (lib == Xpetra::UseTpetra) {
657 #ifdef HAVE_XPETRA_TPETRA
658  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
659  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
660 
661  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
662  if (tMap.is_null())
663  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
664 
665  return Xpetra::toXpetra(tMap);
666 #else
667  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
668 #endif
669  } else {
670  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
671  }
672 
673  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
674 
675  }
676 
678  static RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ReadBlockedCrsMatrix (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
679  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
680  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMapExtractor;
681  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
682  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
683  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
684  typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
685  typedef Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpIO;
686 
687  size_t numBlocks = 2; // TODO user parameter?
688 
689  std::vector<RCP<const XpMap> > rgMapVec;
690  for(size_t r = 0; r < numBlocks; ++r) {
691  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
692  rgMapVec.push_back(map);
693  }
694  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
695 
696  std::vector<RCP<const XpMap> > doMapVec;
697  for(size_t c = 0; c < numBlocks; ++c) {
698  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
699  doMapVec.push_back(map);
700  }
701  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
702 
703  /*std::vector<RCP<const XpMap> > testRgMapVec;
704  for(size_t r = 0; r < numBlocks; ++r) {
705  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
706  testRgMapVec.push_back(map);
707  }
708  std::vector<RCP<const XpMap> > testDoMapVec;
709  for(size_t c = 0; c < numBlocks; ++c) {
710  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
711  testDoMapVec.push_back(map);
712  }*/
713 
714  // create map extractors
715 
716  // range map extractor
717  bool bRangeUseThyraStyleNumbering = false;
718  /*GlobalOrdinal gMinGids = 0;
719  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
720  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
721  }
722  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
723  */
724  RCP<const XpMapExtractor> rangeMapExtractor =
725  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
726 
727 
728  // domain map extractor
729  bool bDomainUseThyraStyleNumbering = false;
730  /*gMinGids = 0;
731  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
732  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
733  }
734  if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
735  */
736  RCP<const XpMapExtractor> domainMapExtractor =
737  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
738 
739  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
740 
741  // write all matrices with their maps
742  for (size_t r = 0; r < numBlocks; ++r) {
743  for (size_t c = 0; c < numBlocks; ++c) {
744  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
745  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
746  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
747  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
748  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
749  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
750  bOp->setMatrix(r, c, mat);
751  }
752  }
753 
754  bOp->fillComplete();
755 
756  return bOp;
757  } //ReadBlockedCrsMatrix
758 
759 
761  template<class T>
762  static std::string toString(const T& what) {
763  std::ostringstream buf;
764  buf << what;
765  return buf.str();
766  }
767  };
768 
769 
770 #ifdef HAVE_XPETRA_EPETRA
780  template <class Scalar>
781  class IO<Scalar,int,int,EpetraNode> {
782  public:
783  typedef int LocalOrdinal;
784  typedef int GlobalOrdinal;
785  typedef EpetraNode Node;
786 
787 #ifdef HAVE_XPETRA_EPETRA
789  // @{
790  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
791  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
792  if (xeMap == Teuchos::null)
793  throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
794  return xeMap->getEpetra_Map();
795  }
796  // @}
797 #endif
798 
799 #ifdef HAVE_XPETRA_TPETRA
801  // @{
802  static const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
803  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
804  if (tmp_TMap == Teuchos::null)
805  throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
806  return tmp_TMap->getTpetra_Map();
807  }
808 #endif
809 
810 
812 
813 
814  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
815  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > tmp_Map = rcpFromRef(M);
816 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
817  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
818  if (tmp_EMap != Teuchos::null) {
819  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
820  if (rv != 0)
821  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
822  return;
823  }
824 #endif // HAVE_XPETRA_EPETRA
825 
826 #ifdef HAVE_XPETRA_TPETRA
827 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
828  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
829  // do nothing
830 # else
831  const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> > &tmp_TMap =
832  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
833  if (tmp_TMap != Teuchos::null) {
834  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
835  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
836  return;
837  }
838 # endif
839 #endif // HAVE_XPETRA_TPETRA
840  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
841  }
842 
844  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
845  std::string mapfile = "map_" + fileName;
846  Write(mapfile, *(vec.getMap()));
847 
848  RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_Vec = Teuchos::rcpFromRef(vec);
849 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
850  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
851  if (tmp_EVec != Teuchos::null) {
852  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
853  if (rv != 0)
854  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
855  return;
856  }
857 #endif // HAVE_XPETRA_EPETRAEXT
858 
859 #ifdef HAVE_XPETRA_TPETRA
860 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
861  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
862  // do nothin
863 # else
864  const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &tmp_TVec =
865  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
866  if (tmp_TVec != Teuchos::null) {
867  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
868  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
869  return;
870  }
871 # endif
872 #endif // HAVE_XPETRA_TPETRA
873 
874  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
875 
876  } //Write
877 
878 
879 
881  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
882 
883  Write("rowmap_" + fileName, *(Op.getRowMap()));
884  if ( !Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
885  Write("domainmap_" + fileName, *(Op.getDomainMap()));
886  if ( !Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
887  Write("rangemap_" + fileName, *(Op.getRangeMap()));
888  if ( !Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps )
889  Write("colmap_" + fileName, *(Op.getColMap()));
890 
891  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
892  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
893  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
894 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
895  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
896  if (tmp_ECrsMtx != Teuchos::null) {
897  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
898  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
899  if (rv != 0)
900  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
901  return;
902  }
903 #endif // endif HAVE_XPETRA_EPETRA
904 
905 #ifdef HAVE_XPETRA_TPETRA
906 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
907  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
908  // do nothin
909 # else
910  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& tmp_TCrsMtx =
911  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
912  if (tmp_TCrsMtx != Teuchos::null) {
913  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
914  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
915  return;
916  }
917 # endif
918 #endif // HAVE_XPETRA_TPETRA
919 
920  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
921  } //Write
922 
923 
925  static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
926  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp =
927  dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
928  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
929 
930  ArrayRCP<const size_t> rowptr_RCP;
931  ArrayRCP<LocalOrdinal> rowptr2_RCP;
932  ArrayRCP<const LocalOrdinal> colind_RCP;
933  ArrayRCP<const Scalar> vals_RCP;
934  tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
935 
936  ArrayView<const size_t> rowptr = rowptr_RCP();
937  ArrayView<const LocalOrdinal> colind = colind_RCP();
938  ArrayView<const Scalar> vals = vals_RCP();
939 
940  rowptr2_RCP.resize(rowptr.size());
941  ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
942  for (size_t j = 0; j<rowptr.size(); j++)
943  rowptr2[j] = rowptr[j];
944 
945  Teuchos::MatrixMarket::Raw::Writer<Scalar,LocalOrdinal> writer;
946  writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
947  rowptr2,colind,vals,
948  rowptr.size()-1,Op.getColMap()->getNodeNumElements());
949  } //WriteLocal
950 
952  static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
953  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
954  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMapExtractor;
955  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
956  typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
957  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
958  typedef Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpIO;
959 
960  // write all matrices with their maps
961  for (size_t r = 0; r < Op.Rows(); ++r) {
962  for (size_t c = 0; c < Op.Cols(); ++c) {
963  RCP<const XpMat > m = Op.getMatrix(r,c);
964  if(m != Teuchos::null) { // skip empty blocks
965  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const XpCrsMatWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
966  "Sub block matrix (" << r << "," << c << ") is not of type CrsMatrixWrap.");
967  XpIO::Write(fileName + toString(r) + toString(c) + ".m", *m, writeAllMaps);
968  }
969  }
970  }
971 
972  // write map information of map extractors
973  RCP<const XpMapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
974  RCP<const XpMapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
975 
976  for(size_t r = 0; r < rangeMapExtractor->NumMaps(); ++r) {
977  RCP<const XpMap> map = rangeMapExtractor->getMap(r);
978  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", *map);
979  }
980  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
981 
982  for(size_t c = 0; c < domainMapExtractor->NumMaps(); ++c) {
983  RCP<const XpMap> map = domainMapExtractor->getMap(c);
984  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", *map);
985  }
986  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
987  } //WriteBlockCrsMatrix
988 
990  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) {
991  if (binary == false) {
992  // Matrix Market file format (ASCII)
993  if (lib == Xpetra::UseEpetra) {
994 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
995  Epetra_CrsMatrix *eA;
996  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
997  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
998  if (rv != 0)
999  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1000 
1001  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1002 
1003  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A =
1004  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1005  return A;
1006 #else
1007  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1008 #endif
1009  } else if (lib == Xpetra::UseTpetra) {
1010 #ifdef HAVE_XPETRA_TPETRA
1011 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1012  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1013  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
1014 # else
1015  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1016 
1017  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1018 
1019  bool callFillComplete = true;
1020 
1021  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
1022 
1023  if (tA.is_null())
1024  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1025 
1026  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
1027  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
1028  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
1029 
1030  return A;
1031 # endif
1032 #else
1033  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1034 #endif
1035  } else {
1036  throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1037  }
1038  } else {
1039  // Custom file format (binary)
1040  std::ifstream ifs(fileName.c_str(), std::ios::binary);
1041  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
1042  int m, n, nnz;
1043  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1044  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1045  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1046 
1047  int myRank = comm->getRank();
1048 
1049  GlobalOrdinal indexBase = 0;
1050  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
1051  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
1052 
1053  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
1054 
1055  if (myRank == 0) {
1056  Teuchos::Array<GlobalOrdinal> inds;
1057  Teuchos::Array<Scalar> vals;
1058  // Scan matrix to determine the exact nnz per row.
1059  Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
1060  for (int i = 0; i < m; i++) {
1061  int row, rownnz;
1062  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1063  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1064  numEntriesPerRow[i] = rownnz;
1065  for (int j = 0; j < rownnz; j++) {
1066  int index;
1067  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1068  }
1069  for (int j = 0; j < rownnz; j++) {
1070  double value;
1071  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1072  }
1073  }
1074 
1075  A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, numEntriesPerRow);
1076 
1077  // Now that nnz per row are known, reread and store the matrix.
1078  ifs.seekg(0, ifs.beg); //rewind to beginning of file
1079  int junk; //skip header info
1080  ifs.read(reinterpret_cast<char*>(&m), sizeof(junk));
1081  ifs.read(reinterpret_cast<char*>(&n), sizeof(junk));
1082  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(junk));
1083  for (int i = 0; i < m; i++) {
1084  int row, rownnz;
1085  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1086  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1087  inds.resize(rownnz);
1088  vals.resize(rownnz);
1089  for (int j = 0; j < rownnz; j++) {
1090  int index;
1091  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1092  inds[j] = Teuchos::as<GlobalOrdinal>(index);
1093  }
1094  for (int j = 0; j < rownnz; j++) {
1095  double value;
1096  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1097  vals[j] = Teuchos::as<Scalar>(value);
1098  }
1099  A->insertGlobalValues(row, inds, vals);
1100  }
1101  } //if (myRank == 0)
1102 
1103  A->fillComplete(domainMap, rangeMap);
1104 
1105  return A;
1106  }
1107 
1108  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1109 
1110  } //Read()
1111 
1112 
1117  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Read(const std::string& filename,
1118  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rowMap,
1119  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
1120  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
1121  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
1122  const bool callFillComplete = true,
1123  const bool binary = false,
1124  const bool tolerant = false,
1125  const bool debug = false) {
1126  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
1127 
1128  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
1129  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
1130 
1131  const Xpetra::UnderlyingLib lib = rowMap->lib();
1132  if (binary == false) {
1133  if (lib == Xpetra::UseEpetra) {
1134 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1135  Epetra_CrsMatrix *eA;
1136  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
1137  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
1138  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
1139  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
1140  int rv;
1141  if (colMap.is_null()) {
1142  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
1143 
1144  } else {
1145  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
1146  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
1147  }
1148 
1149  if (rv != 0)
1150  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1151 
1152  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1153  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A =
1154  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1155 
1156  return A;
1157 #else
1158  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1159 #endif
1160  } else if (lib == Xpetra::UseTpetra) {
1161 #ifdef HAVE_XPETRA_TPETRA
1162 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1163  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1164  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1165 # else
1166  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1167  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1168  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1169 
1170  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
1171  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
1172  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
1173  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
1174 
1175  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
1176  callFillComplete, tolerant, debug);
1177  if (tA.is_null())
1178  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1179 
1180  RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tA));
1181  RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
1182  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tmpA2));
1183 
1184  return A;
1185 # endif
1186 #else
1187  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1188 #endif
1189  } else {
1190  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1191  }
1192  } else {
1193  // Custom file format (binary)
1194  std::ifstream ifs(filename.c_str(), std::ios::binary);
1195  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
1196  int m, n, nnz;
1197  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1198  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1199  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1200 
1201  //2020-June-05 JHU : for Tpetra, this will probably fail because Tpetra now requires staticly-sized matrix graphs.
1202  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap, 1);
1203 
1204  //2019-06-07 JHU I don't see why this should matter.
1205  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1206 
1207  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
1208  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
1209 
1210  Teuchos::Array<GlobalOrdinal> inds;
1211  Teuchos::Array<Scalar> vals;
1212  for (int i = 0; i < m; i++) {
1213  int row, rownnz;
1214  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1215  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1216  inds.resize(rownnz);
1217  vals.resize(rownnz);
1218  for (int j = 0; j < rownnz; j++) {
1219  int index;
1220  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1221  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
1222  }
1223  for (int j = 0; j < rownnz; j++) {
1224  double value;
1225  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1226  vals[j] = Teuchos::as<Scalar>(value);
1227  }
1228  //This implies that row is not a global index.
1229  A->insertGlobalValues(rowElements[row], inds, vals);
1230  }
1231  A->fillComplete(domainMap, rangeMap);
1232  return A;
1233  }
1234 
1235  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1236  }
1238 
1239 
1240  static RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ReadMultiVector (const std::string& fileName, const RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& map) {
1241  Xpetra::UnderlyingLib lib = map->lib();
1242 
1243  if (lib == Xpetra::UseEpetra) {
1244  // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
1245  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1246 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1247  Epetra_MultiVector * MV;
1248  EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
1249  RCP<Epetra_MultiVector> MVrcp = rcp(MV);
1250  return Convert_Epetra_MultiVector_ToXpetra_MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(MVrcp);
1251 #else
1252  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1253 #endif
1254  } else if (lib == Xpetra::UseTpetra) {
1255 #ifdef HAVE_XPETRA_TPETRA
1256 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1257  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1258  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1259 # else
1260  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1261  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1262  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1263  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
1264 
1265  RCP<const map_type> temp = toTpetra(map);
1266  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
1267  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rmv = Xpetra::toXpetra(TMV);
1268  return rmv;
1269 # endif
1270 #else
1271  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1272 #endif
1273  } else {
1274  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1275  }
1276 
1277  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1278 
1279  }
1280 
1281 
1282  static RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1283  if (lib == Xpetra::UseEpetra) {
1284  // do we need another specialization for <double,int,int> ??
1285  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1286 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1287  Epetra_Map *eMap;
1288  int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
1289  if (rv != 0)
1290  throw Exceptions::RuntimeError("Error reading map from file " + fileName + " with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
1291 
1292  RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
1293  return Xpetra::toXpetra<int,Node>(*eMap1);
1294 #else
1295  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1296 #endif
1297  } else if (lib == Xpetra::UseTpetra) {
1298 #ifdef HAVE_XPETRA_TPETRA
1299 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1300  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1301  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1302 # else
1303  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1304  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1305 
1306  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
1307  if (tMap.is_null())
1308  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1309 
1310  return Xpetra::toXpetra(tMap);
1311 # endif
1312 #else
1313  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1314 #endif
1315  } else {
1316  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1317  }
1318 
1319  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1320 
1321  }
1322 
1324  static RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ReadBlockedCrsMatrix (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1325  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
1326  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMapExtractor;
1327  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
1328  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
1329  typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
1330  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
1331  typedef Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpIO;
1332 
1333 
1334  size_t numBlocks = 2; // TODO user parameter?
1335 
1336  std::vector<RCP<const XpMap> > rgMapVec;
1337  for(size_t r = 0; r < numBlocks; ++r) {
1338  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
1339  rgMapVec.push_back(map);
1340  }
1341  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
1342 
1343  std::vector<RCP<const XpMap> > doMapVec;
1344  for(size_t c = 0; c < numBlocks; ++c) {
1345  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
1346  doMapVec.push_back(map);
1347  }
1348  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
1349 
1350  /*std::vector<RCP<const XpMap> > testRgMapVec;
1351  for(size_t r = 0; r < numBlocks; ++r) {
1352  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
1353  testRgMapVec.push_back(map);
1354  }
1355  std::vector<RCP<const XpMap> > testDoMapVec;
1356  for(size_t c = 0; c < numBlocks; ++c) {
1357  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
1358  testDoMapVec.push_back(map);
1359  }*/
1360 
1361  // create map extractors
1362 
1363  // range map extractor
1364  bool bRangeUseThyraStyleNumbering = false;
1365  /*
1366  GlobalOrdinal gMinGids = 0;
1367  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
1368  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
1369  }
1370  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;*/
1371  RCP<const XpMapExtractor> rangeMapExtractor =
1372  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
1373 
1374  // domain map extractor
1375  bool bDomainUseThyraStyleNumbering = false;
1376  /*gMinGids = 0;
1377  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
1378  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
1379  }
1380  if ( gMinGids==0 && testDoMapVec.size() > 1) bDomainUseThyraStyleNumbering = true;*/
1381  RCP<const XpMapExtractor> domainMapExtractor =
1382  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
1383 
1384  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
1385 
1386  // write all matrices with their maps
1387  for (size_t r = 0; r < numBlocks; ++r) {
1388  for (size_t c = 0; c < numBlocks; ++c) {
1389  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1390  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1391  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1392  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1393  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1394  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
1395  bOp->setMatrix(r, c, mat);
1396  }
1397  }
1398 
1399  bOp->fillComplete();
1400 
1401  return bOp;
1402  } //ReadBlockedCrsMatrix
1403 
1405  template<class T>
1406  static std::string toString(const T& what) {
1407  std::ostringstream buf;
1408  buf << what;
1409  return buf.str();
1410  }
1411  };
1412 #endif // HAVE_XPETRA_EPETRA
1413 
1414 
1415 } // end namespace Xpetra
1416 
1417 #define XPETRA_IO_SHORT
1418 
1419 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
Xpetra_BlockedCrsMatrix.hpp
Xpetra_ConfigDefs.hpp
Xpetra_EpetraMap.hpp
Xpetra_EpetraMultiVector.hpp
Xpetra_EpetraUtils.hpp
Xpetra_MatrixFactory.hpp
Xpetra_MatrixMatrix.hpp
Xpetra_Matrix.hpp
LO
LocalOrdinal LO
Definition: Xpetra_UseShortNamesOrdinal.hpp:138
GO
GlobalOrdinal GO
Definition: Xpetra_UseShortNamesOrdinal.hpp:139
NO
Node NO
Definition: Xpetra_UseShortNamesOrdinal.hpp:140
SC
Scalar SC
Definition: Xpetra_UseShortNamesScalar.hpp:178
Xpetra_UseShortNames.hpp
Epetra_CrsMatrix
Epetra_Map
Epetra_MultiVector
Teuchos::ArrayRCP::resize
void resize(const size_type n, const T &val=T())
Teuchos::ArrayView::size
size_type size() const
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
Teuchos::MatrixMarket::Raw::Writer::writeFile
void writeFile(const std::string &filename, const ArrayView< const OrdinalType > &rowptr, const ArrayView< const OrdinalType > &colind, const ArrayView< const ScalarType > &values, const OrdinalType numRows, const OrdinalType numCols)
Teuchos::RCP::is_null
bool is_null() const
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1296
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1085
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1082
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1299
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition: Xpetra_BlockedCrsMatrix.hpp:1332
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_decl.hpp:85
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap_def.hpp:457
Xpetra::DistObject::getMap
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
Xpetra::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: Xpetra_Exceptions.hpp:87
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:102
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:1117
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 matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:952
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:802
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:925
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:1282
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:881
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:790
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:1406
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:814
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 matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:1324
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:844
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:990
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:1240
Xpetra::IO
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
Definition: Xpetra_IO.hpp:159
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:653
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:508
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:211
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:386
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 matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:678
Xpetra::IO::toString
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:762
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:250
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:182
Xpetra::IO::Write
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:223
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:281
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:320
Xpetra::IO::WriteBlockedCrsMatrix
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:348
Xpetra::IO::ReadMultiVector
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
Definition: Xpetra_IO.hpp:626
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::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory.hpp:234
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix.hpp:97
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
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::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....
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....
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_UNREACHABLE_RETURN
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
Teuchos::toString
std::string toString(const T &t)
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra
Xpetra namespace
Definition: Xpetra_BlockedCrsMatrix.hpp:88
Xpetra::toEpetra
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
Definition: Xpetra_EpetraCrsGraph.cpp:57
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition: Xpetra_EpetraCrsGraph.cpp:65
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:107
Xpetra::toTpetra
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
Definition: Xpetra_TpetraCrsGraph_decl.hpp:378
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:131