All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
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"
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> >
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 
122  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC,LO,GO,NO> >(tmpC1);
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> >
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,
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
168  // @{
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 
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
192  // @{
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 
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) {
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
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 
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
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 
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
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) {
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 
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) {
351  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
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 
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 
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;
444 
445  //2019-06-07 JHU I don't see why this should matter.
446  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
447 
448  if (myRank == 0) {
451  for (int i = 0; i < m; i++) {
452  int row, rownnz;
453  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
454  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
455  inds.resize(rownnz);
456  vals.resize(rownnz);
457  for (int j = 0; j < rownnz; j++) {
458  int index;
459  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
460  inds[j] = Teuchos::as<GlobalOrdinal>(index);
461  }
462  for (int j = 0; j < rownnz; j++) {
463  double value;
464  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
465  vals[j] = Teuchos::as<SC>(value);
466  }
467  A->insertGlobalValues(row, inds, vals);
468  }
469  }
470 
471  A->fillComplete(domainMap, rangeMap);
472 
473  return A;
474  }
475 
476  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
477 
478  } //Read()
479 
480 
486  Read(const std::string& filename,
488  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
489  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
490  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
491  const bool callFillComplete = true,
492  const bool binary = false,
493  const bool tolerant = false,
494  const bool debug = false) {
495  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
496 
497  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
498  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
499 
500  const Xpetra::UnderlyingLib lib = rowMap->lib();
501  if (binary == false) {
502  if (lib == Xpetra::UseEpetra) {
503 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
504  Epetra_CrsMatrix *eA;
505  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
507  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
508  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
509  int rv;
510  if (colMap.is_null()) {
511  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
512 
513  } else {
514  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
515  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
516  }
517 
518  if (rv != 0)
519  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
520 
521  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
523  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
524 
525  return A;
526 #else
527  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
528 #endif
529  } else if (lib == Xpetra::UseTpetra) {
530 #ifdef HAVE_XPETRA_TPETRA
531  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
532  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
533  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
534 
535  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
536  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
537  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
538  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
539 
540  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
541  callFillComplete, tolerant, debug);
542  if (tA.is_null())
543  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
544 
548 
549  return A;
550 #else
551  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
552 #endif
553  } else {
554  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
555  }
556  } else {
557  // Custom file format (binary)
558  std::ifstream ifs(filename.c_str(), std::ios::binary);
559  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
560  int m, n, nnz;
561  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
562  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
563  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
564 
566 
567  //2019-06-07 JHU I don't see why this should matter.
568  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
569 
570  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
571  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
572 
575  for (int i = 0; i < m; i++) {
576  int row, rownnz;
577  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
578  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
579  inds.resize(rownnz);
580  vals.resize(rownnz);
581  for (int j = 0; j < rownnz; j++) {
582  int index;
583  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
584  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
585  }
586  for (int j = 0; j < rownnz; j++) {
587  double value;
588  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
589  vals[j] = Teuchos::as<SC>(value);
590  }
591  A->insertGlobalValues(rowElements[row], inds, vals);
592  }
593  A->fillComplete(domainMap, rangeMap);
594  return A;
595  }
596 
597  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
598  }
600 
601 
602  static RCP<MultiVector> ReadMultiVector (const std::string& fileName, const RCP<const Map>& map) {
603  Xpetra::UnderlyingLib lib = map->lib();
604 
605  if (lib == Xpetra::UseEpetra) {
606  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
607 
608  } else if (lib == Xpetra::UseTpetra) {
609 #ifdef HAVE_XPETRA_TPETRA
610  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
611  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
612  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
613  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
614 
615  RCP<const map_type> temp = toTpetra(map);
616  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
618  return rmv;
619 #else
620  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
621 #endif
622  } else {
623  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
624  }
625 
626  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
627  }
628 
629  static RCP<const Map> ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
630  if (lib == Xpetra::UseEpetra) {
631  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
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 
637  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
638  if (tMap.is_null())
639  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
640 
641  return Xpetra::toXpetra(tMap);
642 #else
643  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
644 #endif
645  } else {
646  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
647  }
648 
649  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
650 
651  }
652 
657  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
658  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
662 
663  size_t numBlocks = 2; // TODO user parameter?
664 
665  std::vector<RCP<const XpMap> > rgMapVec;
666  for(size_t r = 0; r < numBlocks; ++r) {
667  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
668  rgMapVec.push_back(map);
669  }
670  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
671 
672  std::vector<RCP<const XpMap> > doMapVec;
673  for(size_t c = 0; c < numBlocks; ++c) {
674  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
675  doMapVec.push_back(map);
676  }
677  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
678 
679  /*std::vector<RCP<const XpMap> > testRgMapVec;
680  for(size_t r = 0; r < numBlocks; ++r) {
681  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
682  testRgMapVec.push_back(map);
683  }
684  std::vector<RCP<const XpMap> > testDoMapVec;
685  for(size_t c = 0; c < numBlocks; ++c) {
686  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
687  testDoMapVec.push_back(map);
688  }*/
689 
690  // create map extractors
691 
692  // range map extractor
693  bool bRangeUseThyraStyleNumbering = false;
694  /*GlobalOrdinal gMinGids = 0;
695  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
696  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
697  }
698  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
699  */
700  RCP<const XpMapExtractor> rangeMapExtractor =
701  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
702 
703 
704  // domain map extractor
705  bool bDomainUseThyraStyleNumbering = false;
706  /*gMinGids = 0;
707  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
708  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
709  }
710  if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
711  */
712  RCP<const XpMapExtractor> domainMapExtractor =
713  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
714 
715  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
716 
717  // write all matrices with their maps
718  for (size_t r = 0; r < numBlocks; ++r) {
719  for (size_t c = 0; c < numBlocks; ++c) {
720  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
721  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
722  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
723  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
724  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
725  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
726  bOp->setMatrix(r, c, mat);
727  }
728  }
729 
730  bOp->fillComplete();
731 
732  return bOp;
733  } //ReadBlockedCrsMatrix
734 
735 
737  template<class T>
738  static std::string toString(const T& what) {
739  std::ostringstream buf;
740  buf << what;
741  return buf.str();
742  }
743  };
744 
745 
746 #ifdef HAVE_XPETRA_EPETRA
747 
756  template <class Scalar>
757  class IO<Scalar,int,int,EpetraNode> {
758  public:
759  typedef int LocalOrdinal;
760  typedef int GlobalOrdinal;
761  typedef EpetraNode Node;
762 
763 #ifdef HAVE_XPETRA_EPETRA
764  // @{
767  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
768  if (xeMap == Teuchos::null)
769  throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
770  return xeMap->getEpetra_Map();
771  }
772  // @}
773 #endif
774 
775 #ifdef HAVE_XPETRA_TPETRA
776  // @{
779  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
780  if (tmp_TMap == Teuchos::null)
781  throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
782  return tmp_TMap->getTpetra_Map();
783  }
784 #endif
785 
786 
788 
789 
790  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
792 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
793  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
794  if (tmp_EMap != Teuchos::null) {
795  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
796  if (rv != 0)
797  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
798  return;
799  }
800 #endif // HAVE_XPETRA_EPETRA
801 
802 #ifdef HAVE_XPETRA_TPETRA
803 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
804  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
805  // do nothing
806 # else
808  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
809  if (tmp_TMap != Teuchos::null) {
810  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
811  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
812  return;
813  }
814 # endif
815 #endif // HAVE_XPETRA_TPETRA
816  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
817  }
818 
820  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
821  std::string mapfile = "map_" + fileName;
822  Write(mapfile, *(vec.getMap()));
823 
825 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
826  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
827  if (tmp_EVec != Teuchos::null) {
828  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
829  if (rv != 0)
830  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
831  return;
832  }
833 #endif // HAVE_XPETRA_EPETRAEXT
834 
835 #ifdef HAVE_XPETRA_TPETRA
836 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
837  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
838  // do nothin
839 # else
841  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
842  if (tmp_TVec != Teuchos::null) {
843  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
844  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
845  return;
846  }
847 # endif
848 #endif // HAVE_XPETRA_TPETRA
849 
850  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
851 
852  } //Write
853 
854 
855 
857  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
858 
859  Write("rowmap_" + fileName, *(Op.getRowMap()));
860  if ( !Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
861  Write("domainmap_" + fileName, *(Op.getDomainMap()));
862  if ( !Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
863  Write("rangemap_" + fileName, *(Op.getRangeMap()));
864  if ( !Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps )
865  Write("colmap_" + fileName, *(Op.getColMap()));
866 
870 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
871  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
872  if (tmp_ECrsMtx != Teuchos::null) {
873  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
874  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
875  if (rv != 0)
876  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
877  return;
878  }
879 #endif // endif HAVE_XPETRA_EPETRA
880 
881 #ifdef HAVE_XPETRA_TPETRA
882 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
883  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
884  // do nothin
885 # else
887  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
888  if (tmp_TCrsMtx != Teuchos::null) {
889  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
890  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
891  return;
892  }
893 # endif
894 #endif // HAVE_XPETRA_TPETRA
895 
896  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
897  } //Write
898 
899 
901  static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
905 
906  ArrayRCP<const size_t> rowptr_RCP;
907  ArrayRCP<LocalOrdinal> rowptr2_RCP;
908  ArrayRCP<const LocalOrdinal> colind_RCP;
909  ArrayRCP<const Scalar> vals_RCP;
910  tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
911 
912  ArrayView<const size_t> rowptr = rowptr_RCP();
913  ArrayView<const LocalOrdinal> colind = colind_RCP();
914  ArrayView<const Scalar> vals = vals_RCP();
915 
916  rowptr2_RCP.resize(rowptr.size());
917  ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
918  for (size_t j = 0; j<rowptr.size(); j++)
919  rowptr2[j] = rowptr[j];
920 
922  writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
923  rowptr2,colind,vals,
924  rowptr.size()-1,Op.getColMap()->getNodeNumElements());
925  } //WriteLocal
926 
928  static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
931  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
935 
936  // write all matrices with their maps
937  for (size_t r = 0; r < Op.Rows(); ++r) {
938  for (size_t c = 0; c < Op.Cols(); ++c) {
939  RCP<const XpMat > m = Op.getMatrix(r,c);
940  if(m != Teuchos::null) { // skip empty blocks
941  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const XpCrsMatWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
942  "Sub block matrix (" << r << "," << c << ") is not of type CrsMatrixWrap.");
943  XpIO::Write(fileName + toString(r) + toString(c) + ".m", *m, writeAllMaps);
944  }
945  }
946  }
947 
948  // write map information of map extractors
949  RCP<const XpMapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
950  RCP<const XpMapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
951 
952  for(size_t r = 0; r < rangeMapExtractor->NumMaps(); ++r) {
953  RCP<const XpMap> map = rangeMapExtractor->getMap(r);
954  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", *map);
955  }
956  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
957 
958  for(size_t c = 0; c < domainMapExtractor->NumMaps(); ++c) {
959  RCP<const XpMap> map = domainMapExtractor->getMap(c);
960  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", *map);
961  }
962  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
963  } //WriteBlockCrsMatrix
964 
966  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) {
967  if (binary == false) {
968  // Matrix Market file format (ASCII)
969  if (lib == Xpetra::UseEpetra) {
970 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
971  Epetra_CrsMatrix *eA;
972  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
973  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
974  if (rv != 0)
975  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
976 
977  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
978 
980  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
981  return A;
982 #else
983  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
984 #endif
985  } else if (lib == Xpetra::UseTpetra) {
986 #ifdef HAVE_XPETRA_TPETRA
987 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
988  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
989  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
990 # else
991  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
992 
993  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
994 
995  bool callFillComplete = true;
996 
997  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
998 
999  if (tA.is_null())
1000  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1001 
1005 
1006  return A;
1007 # endif
1008 #else
1009  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1010 #endif
1011  } else {
1012  throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1013  }
1014  } else {
1015  // Custom file format (binary)
1016  std::ifstream ifs(fileName.c_str(), std::ios::binary);
1017  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
1018  int m, n, nnz;
1019  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1020  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1021  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1022 
1023  int myRank = comm->getRank();
1024 
1025  GlobalOrdinal indexBase = 0;
1026  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
1027  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
1029 
1030  //2019-06-07 JHU I don't see why this should matter.
1031  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1032 
1033  if (myRank == 0) {
1036  for (int i = 0; i < m; i++) {
1037  int row, rownnz;
1038  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1039  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1040  inds.resize(rownnz);
1041  vals.resize(rownnz);
1042  for (int j = 0; j < rownnz; j++) {
1043  int index;
1044  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1045  inds[j] = Teuchos::as<GlobalOrdinal>(index);
1046  }
1047  for (int j = 0; j < rownnz; j++) {
1048  double value;
1049  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1050  vals[j] = Teuchos::as<Scalar>(value);
1051  }
1052  A->insertGlobalValues(row, inds, vals);
1053  }
1054  }
1055 
1056  A->fillComplete(domainMap, rangeMap);
1057 
1058  return A;
1059  }
1060 
1061  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1062 
1063  } //Read()
1064 
1065 
1072  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
1073  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
1074  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
1075  const bool callFillComplete = true,
1076  const bool binary = false,
1077  const bool tolerant = false,
1078  const bool debug = false) {
1079  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
1080 
1081  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
1082  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
1083 
1084  const Xpetra::UnderlyingLib lib = rowMap->lib();
1085  if (binary == false) {
1086  if (lib == Xpetra::UseEpetra) {
1087 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1088  Epetra_CrsMatrix *eA;
1089  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
1091  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
1092  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
1093  int rv;
1094  if (colMap.is_null()) {
1095  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
1096 
1097  } else {
1098  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
1099  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
1100  }
1101 
1102  if (rv != 0)
1103  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1104 
1105  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1107  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1108 
1109  return A;
1110 #else
1111  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1112 #endif
1113  } else if (lib == Xpetra::UseTpetra) {
1114 #ifdef HAVE_XPETRA_TPETRA
1115 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1116  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1117  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1118 # else
1119  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1120  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1121  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1122 
1123  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
1124  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
1125  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
1126  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
1127 
1128  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
1129  callFillComplete, tolerant, debug);
1130  if (tA.is_null())
1131  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1132 
1136 
1137  return A;
1138 # endif
1139 #else
1140  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1141 #endif
1142  } else {
1143  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1144  }
1145  } else {
1146  // Custom file format (binary)
1147  std::ifstream ifs(filename.c_str(), std::ios::binary);
1148  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
1149  int m, n, nnz;
1150  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1151  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1152  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1153 
1155 
1156  //2019-06-07 JHU I don't see why this should matter.
1157  //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1158 
1159  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
1160  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
1161 
1164  for (int i = 0; i < m; i++) {
1165  int row, rownnz;
1166  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1167  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1168  inds.resize(rownnz);
1169  vals.resize(rownnz);
1170  for (int j = 0; j < rownnz; j++) {
1171  int index;
1172  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1173  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
1174  }
1175  for (int j = 0; j < rownnz; j++) {
1176  double value;
1177  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1178  vals[j] = Teuchos::as<Scalar>(value);
1179  }
1180  A->insertGlobalValues(rowElements[row], inds, vals);
1181  }
1182  A->fillComplete(domainMap, rangeMap);
1183  return A;
1184  }
1185 
1186  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1187  }
1189 
1190 
1192  Xpetra::UnderlyingLib lib = map->lib();
1193 
1194  if (lib == Xpetra::UseEpetra) {
1195  // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
1196  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1197 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1198  Epetra_MultiVector * MV;
1199  EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
1200  RCP<Epetra_MultiVector> MVrcp = rcp(MV);
1201  return Convert_Epetra_MultiVector_ToXpetra_MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(MVrcp);
1202 #else
1203  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1204 #endif
1205  } else if (lib == Xpetra::UseTpetra) {
1206 #ifdef HAVE_XPETRA_TPETRA
1207 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1208  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1209  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1210 # else
1211  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1212  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1213  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1214  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
1215 
1216  RCP<const map_type> temp = toTpetra(map);
1217  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
1219  return rmv;
1220 # endif
1221 #else
1222  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1223 #endif
1224  } else {
1225  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1226  }
1227 
1228  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1229 
1230  }
1231 
1232 
1233  static RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1234  if (lib == Xpetra::UseEpetra) {
1235  // do we need another specialization for <double,int,int> ??
1236  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1237 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1238  Epetra_Map *eMap;
1239  int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
1240  if (rv != 0)
1241  throw Exceptions::RuntimeError("Error reading map from file " + fileName + " with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
1242 
1243  RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
1244  return Xpetra::toXpetra<int,Node>(*eMap1);
1245 #else
1246  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1247 #endif
1248  } else if (lib == Xpetra::UseTpetra) {
1249 #ifdef HAVE_XPETRA_TPETRA
1250 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1251  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1252  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1253 # else
1254  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1255  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1256 
1257  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
1258  if (tMap.is_null())
1259  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1260 
1261  return Xpetra::toXpetra(tMap);
1262 # endif
1263 #else
1264  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1265 #endif
1266  } else {
1267  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1268  }
1269 
1270  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1271 
1272  }
1273 
1278  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
1279  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
1283 
1284 
1285  size_t numBlocks = 2; // TODO user parameter?
1286 
1287  std::vector<RCP<const XpMap> > rgMapVec;
1288  for(size_t r = 0; r < numBlocks; ++r) {
1289  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
1290  rgMapVec.push_back(map);
1291  }
1292  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
1293 
1294  std::vector<RCP<const XpMap> > doMapVec;
1295  for(size_t c = 0; c < numBlocks; ++c) {
1296  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
1297  doMapVec.push_back(map);
1298  }
1299  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
1300 
1301  /*std::vector<RCP<const XpMap> > testRgMapVec;
1302  for(size_t r = 0; r < numBlocks; ++r) {
1303  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
1304  testRgMapVec.push_back(map);
1305  }
1306  std::vector<RCP<const XpMap> > testDoMapVec;
1307  for(size_t c = 0; c < numBlocks; ++c) {
1308  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
1309  testDoMapVec.push_back(map);
1310  }*/
1311 
1312  // create map extractors
1313 
1314  // range map extractor
1315  bool bRangeUseThyraStyleNumbering = false;
1316  /*
1317  GlobalOrdinal gMinGids = 0;
1318  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
1319  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
1320  }
1321  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;*/
1322  RCP<const XpMapExtractor> rangeMapExtractor =
1323  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
1324 
1325  // domain map extractor
1326  bool bDomainUseThyraStyleNumbering = false;
1327  /*gMinGids = 0;
1328  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
1329  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
1330  }
1331  if ( gMinGids==0 && testDoMapVec.size() > 1) bDomainUseThyraStyleNumbering = true;*/
1332  RCP<const XpMapExtractor> domainMapExtractor =
1333  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
1334 
1335  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
1336 
1337  // write all matrices with their maps
1338  for (size_t r = 0; r < numBlocks; ++r) {
1339  for (size_t c = 0; c < numBlocks; ++c) {
1340  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1341  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1342  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1343  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1344  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1345  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
1346  bOp->setMatrix(r, c, mat);
1347  }
1348  }
1349 
1350  bOp->fillComplete();
1351 
1352  return bOp;
1353  } //ReadBlockedCrsMatrix
1354 
1356  template<class T>
1357  static std::string toString(const T& what) {
1358  std::ostringstream buf;
1359  buf << what;
1360  return buf.str();
1361  }
1362  };
1363 #endif // HAVE_XPETRA_EPETRA
1364 
1365 
1366 } // end namespace Xpetra
1367 
1368 #define XPETRA_IO_SHORT
1369 
1370 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
RCP< CrsMatrix > getCrsMatrix() const
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:486
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &)
Definition: Xpetra_IO.hpp:107
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:1070
GlobalOrdinal GO
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:223
size_type size() const
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:901
Exception throws to report errors in the internal logical of the program.
LocalOrdinal LO
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
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:1275
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:778
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
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
Exception indicating invalid cast attempted.
void resize(const size_type n, const T &val=T())
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
Definition: Xpetra_IO.hpp:629
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
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:966
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
Definition: Xpetra_IO.hpp:602
virtual size_t Cols() const
number of column blocks
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
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:928
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > Convert_Epetra_MultiVector_ToXpetra_MultiVector(RCP< Epetra_MultiVector > &epX)
Definition: Xpetra_IO.hpp:131
void resize(size_type new_size, const value_type &x=value_type())
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
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:1357
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:1191
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
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)
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
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
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
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
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
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
Concrete implementation of Xpetra::Matrix.
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
virtual size_t Rows() const
number of row blocks
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
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:1233
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:766
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
Definition: Xpetra_IO.hpp:159
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:790
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
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:857
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:820
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra-specific matrix class.
std::string toString(const T &t)
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:654
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
bool is_null() const
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:738