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 "Xpetra_ConfigDefs.hpp"
52 
53 #ifdef HAVE_XPETRA_EPETRA
54 # ifdef HAVE_MPI
55 # include "Epetra_MpiComm.h"
56 # endif
57 #endif
58 
59 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
60 #include <EpetraExt_MatrixMatrix.h>
61 #include <EpetraExt_RowMatrixOut.h>
62 #include <EpetraExt_MultiVectorOut.h>
63 #include <EpetraExt_CrsMatrixIn.h>
64 #include <EpetraExt_MultiVectorIn.h>
65 #include <EpetraExt_BlockMapIn.h>
66 #include <Xpetra_EpetraUtils.hpp>
68 #include <EpetraExt_BlockMapOut.h>
69 #endif
70 
71 #ifdef HAVE_XPETRA_TPETRA
72 #include <MatrixMarket_Tpetra.hpp>
73 #include <Tpetra_RowMatrixTransposer.hpp>
74 #include <TpetraExt_MatrixMatrix.hpp>
78 #endif
79 
80 #ifdef HAVE_XPETRA_EPETRA
81 #include <Xpetra_EpetraMap.hpp>
82 #endif
83 
84 #include "Xpetra_Matrix.hpp"
85 #include "Xpetra_MatrixMatrix.hpp"
86 #include "Xpetra_CrsMatrixWrap.hpp"
88 
89 #include "Xpetra_Map.hpp"
90 #include "Xpetra_StridedMap.hpp"
92 #include "Xpetra_MapExtractor.hpp"
93 #include "Xpetra_MatrixFactory.hpp"
94 
95 
96 
97 namespace Xpetra {
98 
99 
100 #ifdef HAVE_XPETRA_EPETRA
101  // This non-member templated function exists so that the matrix-matrix multiply will compile if Epetra, Tpetra, and ML are enabled.
102  template<class SC,class LO,class GO,class NO>
103  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >
104  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap (RCP<Epetra_CrsMatrix> &epAB) {
105  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
106  "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
107  return Teuchos::null;
108  }
109 
110  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
111  template<>
112  inline RCP<Xpetra::CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_CrsMatrix> &epAB) {
113  typedef double SC;
114  typedef int LO;
115  typedef int GO;
116  typedef Xpetra::EpetraNode NO;
117 
118  RCP<Xpetra::EpetraCrsMatrixT<GO,NO> > tmpC1 = rcp(new Xpetra::EpetraCrsMatrixT<GO,NO>(epAB));
119  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC,LO,GO,NO> >(tmpC1);
120  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > tmpC3 = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(tmpC2));
121 
122  return tmpC3;
123  }
124 #endif
125 
130  template <class Scalar,
131  class LocalOrdinal = int,
132  class GlobalOrdinal = LocalOrdinal,
134  class IO {
135 
136  private:
137 #undef XPETRA_IO_SHORT
138 #include "Xpetra_UseShortNames.hpp"
139 
140  public:
141 
142 #ifdef HAVE_XPETRA_EPETRA
143  // @{
145  /*static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec);
146  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec);
147 
148  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec);
149  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec);
150 
151  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op);
152  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op);
153 
154  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op);
155  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op);*/
156 
157  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
158  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
159  if (xeMap == Teuchos::null)
160  throw Exceptions::BadCast("Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
161  return xeMap->getEpetra_Map();
162  }
163  // @}
164 #endif
165 
166 #ifdef HAVE_XPETRA_TPETRA
167  // @{
169  /*static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec);
170  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec);
171  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec);
172 
173  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec);
174  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec);
175 
176  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op);
177  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op);
178 
179  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op);
180  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op);
181 
182  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op);
183  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op);*/
184 
185 
186  static const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
187  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
188  if (tmp_TMap == Teuchos::null)
189  throw Exceptions::BadCast("Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
190  return tmp_TMap->getTpetra_Map();
191  }
192 #endif
193 
194 
196 
197 
198  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
199  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > tmp_Map = rcpFromRef(M);
200 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
201  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
202  if (tmp_EMap != Teuchos::null) {
203  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
204  if (rv != 0)
205  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
206  return;
207  }
208 #endif // HAVE_XPETRA_EPETRAEXT
209 
210 #ifdef HAVE_XPETRA_TPETRA
211  const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> > &tmp_TMap =
212  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
213  if (tmp_TMap != Teuchos::null) {
214  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
215  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
216  return;
217  }
218 #endif // HAVE_XPETRA_TPETRA
219 
220  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
221 
222  } //Write
223 
225  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
226  std::string mapfile = "map_" + fileName;
227  Write(mapfile, *(vec.getMap()));
228 
229  RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_Vec = Teuchos::rcpFromRef(vec);
230 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
231  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
232  if (tmp_EVec != Teuchos::null) {
233  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
234  if (rv != 0)
235  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
236  return;
237  }
238 #endif // HAVE_XPETRA_EPETRA
239 
240 #ifdef HAVE_XPETRA_TPETRA
241  const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &tmp_TVec =
242  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
243  if (tmp_TVec != Teuchos::null) {
244  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
245  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
246  return;
247  }
248 #endif // HAVE_XPETRA_TPETRA
249 
250  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
251 
252  } //Write
253 
254 
256  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
257 
258  Write("rowmap_" + fileName, *(Op.getRowMap()));
259  Write("colmap_" + fileName, *(Op.getColMap()));
260  Write("domainmap_" + fileName, *(Op.getDomainMap()));
261  Write("rangemap_" + fileName, *(Op.getRangeMap()));
262 
265  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
266 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
267  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
268  if (tmp_ECrsMtx != Teuchos::null) {
269  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
270  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
271  if (rv != 0)
272  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
273  return;
274  }
275 #endif
276 
277 #ifdef HAVE_XPETRA_TPETRA
278  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& tmp_TCrsMtx =
279  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
280  if (tmp_TCrsMtx != Teuchos::null) {
281  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
282  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
283  return;
284  }
285 #endif // HAVE_XPETRA_TPETRA
286 
287  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
288  } //Write
289 
294  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
298 
299  // write all matrices with their maps
300  for (size_t r = 0; r < Op.Rows(); ++r) {
301  for (size_t c = 0; c < Op.Cols(); ++c) {
302  RCP<const XpMat > m = Op.getMatrix(r,c);
303  if(m != Teuchos::null) { // skip empty blocks
304  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const XpCrsMatWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
305  "Sub block matrix (" << r << "," << c << ") is not of type CrsMatrixWrap.");
306  XpIO::Write(fileName + toString(r) + toString(c) + ".m", *m);
307  }
308  }
309  }
310 
311  // write map information of map extractors
312  RCP<const XpMapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
313  RCP<const XpMapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
314 
315  for(size_t r = 0; r < rangeMapExtractor->NumMaps(); ++r) {
316  RCP<const XpMap> map = rangeMapExtractor->getMap(r);
317  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", *map);
318  }
319  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
320 
321  for(size_t c = 0; c < domainMapExtractor->NumMaps(); ++c) {
322  RCP<const XpMap> map = domainMapExtractor->getMap(c);
323  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", *map);
324  }
325  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
326  } //WriteBlockCrsMatrix
327 
329  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) {
330  if (binary == false) {
331  // Matrix Market file format (ASCII)
332  if (lib == Xpetra::UseEpetra) {
333 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
334  Epetra_CrsMatrix *eA;
335  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
336  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
337  if (rv != 0)
338  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
339 
340  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
341 
342  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A =
343  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
344  return A;
345 #else
346  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
347 #endif
348  } else if (lib == Xpetra::UseTpetra) {
349 #ifdef HAVE_XPETRA_TPETRA
350  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
351 
352  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
353 
354  //RCP<Node> node = Xpetra::DefaultPlatform::getDefaultPlatform().getNode();
355  Teuchos::ParameterList pl = Teuchos::ParameterList();
356  RCP<Node> node = rcp(new Node(pl));
357  bool callFillComplete = true;
358 
359  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, node, callFillComplete);
360 
361  if (tA.is_null())
362  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
363 
364  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
365  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
366  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
367 
368  return A;
369 #else
370  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
371 #endif
372  } else {
373  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
374  }
375  } else {
376  // Custom file format (binary)
377  std::ifstream ifs(fileName.c_str(), std::ios::binary);
378  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
379  int m, n, nnz;
380  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
381  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
382  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
383 
384  int myRank = comm->getRank();
385 
386  GO indexBase = 0;
387  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
388  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
389  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, 1);
390 
391  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
392 
393  if (myRank == 0) {
394  Teuchos::Array<GlobalOrdinal> inds;
395  Teuchos::Array<Scalar> vals;
396  for (int i = 0; i < m; i++) {
397  int row, rownnz;
398  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
399  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
400  inds.resize(rownnz);
401  vals.resize(rownnz);
402  for (int j = 0; j < rownnz; j++) {
403  int index;
404  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
405  inds[j] = Teuchos::as<GlobalOrdinal>(index);
406  }
407  for (int j = 0; j < rownnz; j++) {
408  double value;
409  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
410  vals[j] = Teuchos::as<SC>(value);
411  }
412  A->insertGlobalValues(row, inds, vals);
413  }
414  }
415 
416  A->fillComplete(domainMap, rangeMap);
417 
418  return A;
419  }
420 
421  return Teuchos::null;
422 
423  } //Read()
424 
425 
430  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
431  Read(const std::string& filename,
432  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rowMap,
433  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
434  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
435  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
436  const bool callFillComplete = true,
437  const bool binary = false,
438  const bool tolerant = false,
439  const bool debug = false) {
440  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
441 
442  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
443  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
444 
445  const Xpetra::UnderlyingLib lib = rowMap->lib();
446  if (binary == false) {
447  if (lib == Xpetra::UseEpetra) {
448 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
449  Epetra_CrsMatrix *eA;
450  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
451  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
452  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
453  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
454  int rv;
455  if (colMap.is_null()) {
456  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
457 
458  } else {
459  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
460  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
461  }
462 
463  if (rv != 0)
464  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
465 
466  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
467  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A =
468  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
469 
470  return A;
471 #else
472  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
473 #endif
474  } else if (lib == Xpetra::UseTpetra) {
475 #ifdef HAVE_XPETRA_TPETRA
476  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
477  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
478  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
479 
480  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
481  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
482  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
483  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
484 
485  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
486  callFillComplete, tolerant, debug);
487  if (tA.is_null())
488  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
489 
490  RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tA));
491  RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
492  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tmpA2));
493 
494  return A;
495 #else
496  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
497 #endif
498  } else {
499  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
500  }
501  } else {
502  // Custom file format (binary)
503  std::ifstream ifs(filename.c_str(), std::ios::binary);
504  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
505  int m, n, nnz;
506  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
507  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
508  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
509 
510  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap, 1);
511 
512  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
513 
514  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
515  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
516 
517  Teuchos::Array<GlobalOrdinal> inds;
518  Teuchos::Array<Scalar> vals;
519  for (int i = 0; i < m; i++) {
520  int row, rownnz;
521  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
522  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
523  inds.resize(rownnz);
524  vals.resize(rownnz);
525  for (int j = 0; j < rownnz; j++) {
526  int index;
527  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
528  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
529  }
530  for (int j = 0; j < rownnz; j++) {
531  double value;
532  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
533  vals[j] = Teuchos::as<SC>(value);
534  }
535  A->insertGlobalValues(rowElements[row], inds, vals);
536  }
537  A->fillComplete(domainMap, rangeMap);
538  return A;
539  }
540 
541  return Teuchos::null;
542  }
544 
545 
546  static RCP<MultiVector> ReadMultiVector (const std::string& fileName, const RCP<const Map>& map) {
547  Xpetra::UnderlyingLib lib = map->lib();
548 
549  if (lib == Xpetra::UseEpetra) {
550  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
551 
552  } else if (lib == Xpetra::UseTpetra) {
553 #ifdef HAVE_XPETRA_TPETRA
554  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
555  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
556  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
557  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
558 
559  RCP<const map_type> temp = toTpetra(map);
560  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),map->getNode(),temp);
561  RCP<MultiVector> rmv = Xpetra::toXpetra(TMV);
562  return rmv;
563 #else
564  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
565 #endif
566  } else {
567  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
568  }
569 
570  return Teuchos::null;
571  }
572 
573  static RCP<const Map> ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
574  if (lib == Xpetra::UseEpetra) {
575  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
576  } else if (lib == Xpetra::UseTpetra) {
577 #ifdef HAVE_XPETRA_TPETRA
578  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
579  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
580 
581  RCP<Node> node = rcp(new Node());
582 
583  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm, node);
584  if (tMap.is_null())
585  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
586 
587  return Xpetra::toXpetra(tMap);
588 #else
589  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
590 #endif
591  } else {
592  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
593  }
594 
595  return Teuchos::null;
596 
597  }
598 
600  static RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ReadBlockedCrsMatrix (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
603  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
604  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
608 
609  size_t numBlocks = 2; // TODO user parameter?
610 
611  std::vector<RCP<const XpMap> > rgMapVec;
612  for(size_t r = 0; r < numBlocks; ++r) {
613  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
614  rgMapVec.push_back(map);
615  }
616  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
617 
618  std::vector<RCP<const XpMap> > doMapVec;
619  for(size_t c = 0; c < numBlocks; ++c) {
620  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
621  doMapVec.push_back(map);
622  }
623  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
624 
625  /*std::vector<RCP<const XpMap> > testRgMapVec;
626  for(size_t r = 0; r < numBlocks; ++r) {
627  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
628  testRgMapVec.push_back(map);
629  }
630  std::vector<RCP<const XpMap> > testDoMapVec;
631  for(size_t c = 0; c < numBlocks; ++c) {
632  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
633  testDoMapVec.push_back(map);
634  }*/
635 
636  // create map extractors
637 
638  // range map extractor
639  bool bRangeUseThyraStyleNumbering = false;
640  /*GlobalOrdinal gMinGids = 0;
641  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
642  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
643  }
644  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
645  */
646  RCP<const XpMapExtractor> rangeMapExtractor =
647  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
648 
649 
650  // domain map extractor
651  bool bDomainUseThyraStyleNumbering = false;
652  /*gMinGids = 0;
653  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
654  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
655  }
656  if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
657  */
658  RCP<const XpMapExtractor> domainMapExtractor =
659  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
660 
661  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
662 
663  // write all matrices with their maps
664  for (size_t r = 0; r < numBlocks; ++r) {
665  for (size_t c = 0; c < numBlocks; ++c) {
666  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
667  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
668  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
669  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
670  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
671  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
672  bOp->setMatrix(r, c, mat);
673  }
674  }
675 
676  bOp->fillComplete();
677 
678  return bOp;
679  } //ReadBlockedCrsMatrix
680 
681 
683  template<class T>
684  static std::string toString(const T& what) {
685  std::ostringstream buf;
686  buf << what;
687  return buf.str();
688  }
689  };
690 
691 
692 #ifdef HAVE_XPETRA_EPETRA
693 
702  template <class Scalar>
703  class IO<Scalar,int,int,EpetraNode> {
704  public:
705  typedef int LocalOrdinal;
706  typedef int GlobalOrdinal;
707  typedef EpetraNode Node;
708 
709 #ifdef HAVE_XPETRA_EPETRA
710  // @{
712  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
713  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
714  if (xeMap == Teuchos::null)
715  throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
716  return xeMap->getEpetra_Map();
717  }
718  // @}
719 #endif
720 
721 #ifdef HAVE_XPETRA_TPETRA
722  // @{
724  static const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
725  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
726  if (tmp_TMap == Teuchos::null)
727  throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
728  return tmp_TMap->getTpetra_Map();
729  }
730 #endif
731 
732 
734 
735 
736  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
737  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > tmp_Map = rcpFromRef(M);
738 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
739  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
740  if (tmp_EMap != Teuchos::null) {
741  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
742  if (rv != 0)
743  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
744  return;
745  }
746 #endif // HAVE_XPETRA_EPETRA
747 
748 #ifdef HAVE_XPETRA_TPETRA
749 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
750  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
751  // do nothing
752 # else
753  const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> > &tmp_TMap =
754  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
755  if (tmp_TMap != Teuchos::null) {
756  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
757  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
758  return;
759  }
760 # endif
761 #endif // HAVE_XPETRA_TPETRA
762  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
763  }
764 
766  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
767  std::string mapfile = "map_" + fileName;
768  Write(mapfile, *(vec.getMap()));
769 
770  RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_Vec = Teuchos::rcpFromRef(vec);
771 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
772  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
773  if (tmp_EVec != Teuchos::null) {
774  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
775  if (rv != 0)
776  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
777  return;
778  }
779 #endif // HAVE_XPETRA_EPETRAEXT
780 
781 #ifdef HAVE_XPETRA_TPETRA
782 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
783  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
784  // do nothin
785 # else
786  const RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &tmp_TVec =
787  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
788  if (tmp_TVec != Teuchos::null) {
789  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
790  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
791  return;
792  }
793 # endif
794 #endif // HAVE_XPETRA_TPETRA
795 
796  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
797 
798  } //Write
799 
800 
801 
803  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
804 
805  Write("rowmap_" + fileName, *(Op.getRowMap()));
806  Write("colmap_" + fileName, *(Op.getColMap()));
807  Write("domainmap_" + fileName, *(Op.getDomainMap()));
808  Write("rangemap_" + fileName, *(Op.getRangeMap()));
809 
812  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmp_CrsMtx = crsOp.getCrsMatrix();
813 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
814  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
815  if (tmp_ECrsMtx != Teuchos::null) {
816  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
817  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
818  if (rv != 0)
819  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
820  return;
821  }
822 #endif // endif HAVE_XPETRA_EPETRA
823 
824 #ifdef HAVE_XPETRA_TPETRA
825 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
826  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
827  // do nothin
828 # else
829  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& tmp_TCrsMtx =
830  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
831  if (tmp_TCrsMtx != Teuchos::null) {
832  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
833  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
834  return;
835  }
836 # endif
837 #endif // HAVE_XPETRA_TPETRA
838 
839  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
840  } //Write
841 
846  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
850 
851  // write all matrices with their maps
852  for (size_t r = 0; r < Op.Rows(); ++r) {
853  for (size_t c = 0; c < Op.Cols(); ++c) {
854  RCP<const XpMat > m = Op.getMatrix(r,c);
855  if(m != Teuchos::null) { // skip empty blocks
856  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const XpCrsMatWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
857  "Sub block matrix (" << r << "," << c << ") is not of type CrsMatrixWrap.");
858  XpIO::Write(fileName + toString(r) + toString(c) + ".m", *m);
859  }
860  }
861  }
862 
863  // write map information of map extractors
864  RCP<const XpMapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
865  RCP<const XpMapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
866 
867  for(size_t r = 0; r < rangeMapExtractor->NumMaps(); ++r) {
868  RCP<const XpMap> map = rangeMapExtractor->getMap(r);
869  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", *map);
870  }
871  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
872 
873  for(size_t c = 0; c < domainMapExtractor->NumMaps(); ++c) {
874  RCP<const XpMap> map = domainMapExtractor->getMap(c);
875  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", *map);
876  }
877  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
878  } //WriteBlockCrsMatrix
879 
881  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) {
882  if (binary == false) {
883  // Matrix Market file format (ASCII)
884  if (lib == Xpetra::UseEpetra) {
885 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
886  Epetra_CrsMatrix *eA;
887  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
888  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
889  if (rv != 0)
890  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
891 
892  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
893 
894  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A =
895  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
896  return A;
897 #else
898  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
899 #endif
900  } else if (lib == Xpetra::UseTpetra) {
901 #ifdef HAVE_XPETRA_TPETRA
902 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
903  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
904  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
905 # else
906  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
907 
908  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
909 
910  //RCP<Node> node = Xpetra::DefaultPlatform::getDefaultPlatform().getNode();
911  Teuchos::ParameterList pl = Teuchos::ParameterList();
912  RCP<Node> node = rcp(new Node(pl));
913  bool callFillComplete = true;
914 
915  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, node, callFillComplete);
916 
917  if (tA.is_null())
918  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
919 
920  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tA));
921  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
922  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA2));
923 
924  return A;
925 # endif
926 #else
927  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
928 #endif
929  } else {
930  throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
931  }
932  } else {
933  // Custom file format (binary)
934  std::ifstream ifs(fileName.c_str(), std::ios::binary);
935  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
936  int m, n, nnz;
937  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
938  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
939  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
940 
941  int myRank = comm->getRank();
942 
943  GlobalOrdinal indexBase = 0;
944  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
945  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
946  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, colMap, 1);
947 
948  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
949 
950  if (myRank == 0) {
951  Teuchos::Array<GlobalOrdinal> inds;
952  Teuchos::Array<Scalar> vals;
953  for (int i = 0; i < m; i++) {
954  int row, rownnz;
955  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
956  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
957  inds.resize(rownnz);
958  vals.resize(rownnz);
959  for (int j = 0; j < rownnz; j++) {
960  int index;
961  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
962  inds[j] = Teuchos::as<GlobalOrdinal>(index);
963  }
964  for (int j = 0; j < rownnz; j++) {
965  double value;
966  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
967  vals[j] = Teuchos::as<Scalar>(value);
968  }
969  A->insertGlobalValues(row, inds, vals);
970  }
971  }
972 
973  A->fillComplete(domainMap, rangeMap);
974 
975  return A;
976  }
977 
978  return Teuchos::null;
979 
980  } //Read()
981 
982 
987  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Read(const std::string& filename,
988  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rowMap,
989  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
990  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
991  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
992  const bool callFillComplete = true,
993  const bool binary = false,
994  const bool tolerant = false,
995  const bool debug = false) {
996  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
997 
998  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
999  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
1000 
1001  const Xpetra::UnderlyingLib lib = rowMap->lib();
1002  if (binary == false) {
1003  if (lib == Xpetra::UseEpetra) {
1004 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1005  Epetra_CrsMatrix *eA;
1006  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
1007  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
1008  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
1009  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
1010  int rv;
1011  if (colMap.is_null()) {
1012  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
1013 
1014  } else {
1015  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
1016  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
1017  }
1018 
1019  if (rv != 0)
1020  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1021 
1022  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1023  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A =
1024  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1025 
1026  return A;
1027 #else
1028  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1029 #endif
1030  } else if (lib == Xpetra::UseTpetra) {
1031 #ifdef HAVE_XPETRA_TPETRA
1032 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1033  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1034  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1035 # else
1036  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1037  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1038  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1039 
1040  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
1041  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
1042  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
1043  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
1044 
1045  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
1046  callFillComplete, tolerant, debug);
1047  if (tA.is_null())
1048  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1049 
1050  RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA1 = rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tA));
1051  RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
1052  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tmpA2));
1053 
1054  return A;
1055 # endif
1056 #else
1057  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1058 #endif
1059  } else {
1060  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1061  }
1062  } else {
1063  // Custom file format (binary)
1064  std::ifstream ifs(filename.c_str(), std::ios::binary);
1065  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
1066  int m, n, nnz;
1067  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1068  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1069  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1070 
1071  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap, 1);
1072 
1073  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1074 
1075  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
1076  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
1077 
1078  Teuchos::Array<GlobalOrdinal> inds;
1079  Teuchos::Array<Scalar> vals;
1080  for (int i = 0; i < m; i++) {
1081  int row, rownnz;
1082  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1083  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1084  inds.resize(rownnz);
1085  vals.resize(rownnz);
1086  for (int j = 0; j < rownnz; j++) {
1087  int index;
1088  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1089  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
1090  }
1091  for (int j = 0; j < rownnz; j++) {
1092  double value;
1093  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1094  vals[j] = Teuchos::as<Scalar>(value);
1095  }
1096  A->insertGlobalValues(rowElements[row], inds, vals);
1097  }
1098  A->fillComplete(domainMap, rangeMap);
1099  return A;
1100  }
1101 
1102  return Teuchos::null;
1103  }
1105 
1106 
1107  static RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ReadMultiVector (const std::string& fileName, const RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& map) {
1108  Xpetra::UnderlyingLib lib = map->lib();
1109 
1110  if (lib == Xpetra::UseEpetra) {
1111  // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
1112  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1113 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1114  Epetra_MultiVector * MV;
1115  EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
1116  return Xpetra::toXpetra<int,Node>(rcp(MV));
1117 #else
1118  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1119 #endif
1120  } else if (lib == Xpetra::UseTpetra) {
1121 #ifdef HAVE_XPETRA_TPETRA
1122 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1123  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1124  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1125 # else
1126  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1127  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1128  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1129  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
1130 
1131  RCP<const map_type> temp = toTpetra(map);
1132  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),map->getNode(),temp);
1133  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rmv = Xpetra::toXpetra(TMV);
1134  return rmv;
1135 # endif
1136 #else
1137  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1138 #endif
1139  } else {
1140  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1141  }
1142 
1143  return Teuchos::null;
1144 
1145  }
1146 
1147 
1148  static RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1149  if (lib == Xpetra::UseEpetra) {
1150  // do we need another specialization for <double,int,int> ??
1151  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1152 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1153  Epetra_Map *eMap;
1154  int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
1155  if (rv != 0)
1156  throw Exceptions::RuntimeError("Error reading map from file " + fileName + " with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
1157 
1158  RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
1159  return Xpetra::toXpetra<int,Node>(*eMap1);
1160 #else
1161  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1162 #endif
1163  } else if (lib == Xpetra::UseTpetra) {
1164 #ifdef HAVE_XPETRA_TPETRA
1165 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1166  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1167  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1168 # else
1169  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1170  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1171 
1172  RCP<Node> node = rcp(new Node());
1173 
1174  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm, node);
1175  if (tMap.is_null())
1176  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1177 
1178  return Xpetra::toXpetra(tMap);
1179 # endif
1180 #else
1181  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1182 #endif
1183  } else {
1184  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1185  }
1186 
1187  return Teuchos::null;
1188 
1189  }
1190 
1192  static RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ReadBlockedCrsMatrix (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1195  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
1196  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
1200 
1201 
1202  size_t numBlocks = 2; // TODO user parameter?
1203 
1204  std::vector<RCP<const XpMap> > rgMapVec;
1205  for(size_t r = 0; r < numBlocks; ++r) {
1206  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
1207  rgMapVec.push_back(map);
1208  }
1209  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
1210 
1211  std::vector<RCP<const XpMap> > doMapVec;
1212  for(size_t c = 0; c < numBlocks; ++c) {
1213  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
1214  doMapVec.push_back(map);
1215  }
1216  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
1217 
1218  /*std::vector<RCP<const XpMap> > testRgMapVec;
1219  for(size_t r = 0; r < numBlocks; ++r) {
1220  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
1221  testRgMapVec.push_back(map);
1222  }
1223  std::vector<RCP<const XpMap> > testDoMapVec;
1224  for(size_t c = 0; c < numBlocks; ++c) {
1225  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
1226  testDoMapVec.push_back(map);
1227  }*/
1228 
1229  // create map extractors
1230 
1231  // range map extractor
1232  bool bRangeUseThyraStyleNumbering = false;
1233  /*
1234  GlobalOrdinal gMinGids = 0;
1235  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
1236  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
1237  }
1238  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;*/
1239  RCP<const XpMapExtractor> rangeMapExtractor =
1240  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
1241 
1242  // domain map extractor
1243  bool bDomainUseThyraStyleNumbering = false;
1244  /*gMinGids = 0;
1245  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
1246  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
1247  }
1248  if ( gMinGids==0 && testDoMapVec.size() > 1) bDomainUseThyraStyleNumbering = true;*/
1249  RCP<const XpMapExtractor> domainMapExtractor =
1250  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
1251 
1252  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
1253 
1254  // write all matrices with their maps
1255  for (size_t r = 0; r < numBlocks; ++r) {
1256  for (size_t c = 0; c < numBlocks; ++c) {
1257  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1258  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1259  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1260  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1261  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1262  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
1263  bOp->setMatrix(r, c, mat);
1264  }
1265  }
1266 
1267  bOp->fillComplete();
1268 
1269  return bOp;
1270  } //ReadBlockedCrsMatrix
1271 
1273  template<class T>
1274  static std::string toString(const T& what) {
1275  std::ostringstream buf;
1276  buf << what;
1277  return buf.str();
1278  }
1279  };
1280 #endif // HAVE_XPETRA_EPETRA
1281 
1282 
1283 } // end namespace Xpetra
1284 
1285 #define XPETRA_IO_SHORT
1286 
1287 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
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, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
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:431
RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getTpetra_Map() const
Get the underlying Tpetra map.
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:987
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
GlobalOrdinal GO
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:198
Xpetra namespace
static void Write(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:803
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:1192
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:724
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
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:186
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.
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:843
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
Definition: Xpetra_IO.hpp:573
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:881
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
Definition: Xpetra_IO.hpp:546
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const=0
The Map describing the parallel distribution of this object.
RCP< CrsMatrix > getCrsMatrix() const
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:1274
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
Get the underlying Tpetra matrix.
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:1107
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)
Definition: Xpetra_IO.hpp:104
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:329
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:157
Concrete implementation of Xpetra::Matrix.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
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.
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:1148
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_MultiVector() const
Get the underlying Tpetra multivector.
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:712
virtual size_t Rows() const
number of row blocks
virtual size_t Cols() const
number of column blocks
static void Write(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:256
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
Definition: Xpetra_IO.hpp:134
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:736
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:225
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:766
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.
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:291
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:600
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:684