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"
87 
88 #include "Xpetra_Map.hpp"
89 #include "Xpetra_StridedMap.hpp"
91 #include "Xpetra_MapExtractor.hpp"
92 #include "Xpetra_MatrixFactory.hpp"
93 
94 
95 
96 namespace Xpetra {
97 
98 
99 #ifdef HAVE_XPETRA_EPETRA
100  // This non-member templated function exists so that the matrix-matrix multiply will compile if Epetra, Tpetra, and ML are enabled.
101  template<class SC,class LO,class GO,class NO>
102  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >
105  "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
106  return Teuchos::null;
107  }
108 
109  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
110  template<>
111  inline RCP<Xpetra::CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_CrsMatrix> &epAB) {
112  typedef double SC;
113  typedef int LO;
114  typedef int GO;
115  typedef Xpetra::EpetraNode NO;
116 
118  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC,LO,GO,NO> >(tmpC1);
120 
121  return tmpC3;
122  }
123 #endif
124 
129  template <class Scalar,
130  class LocalOrdinal = int,
131  class GlobalOrdinal = LocalOrdinal,
133  class IO {
134 
135  private:
136 #undef XPETRA_IO_SHORT
137 #include "Xpetra_UseShortNames.hpp"
138 
139  public:
140 
141 #ifdef HAVE_XPETRA_EPETRA
142  // @{
144  /*static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec);
145  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec);
146 
147  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec);
148  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec);
149 
150  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op);
151  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op);
152 
153  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op);
154  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op);*/
155 
156  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
157  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
158  if (xeMap == Teuchos::null)
159  throw Exceptions::BadCast("Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
160  return xeMap->getEpetra_Map();
161  }
162  // @}
163 #endif
164 
165 #ifdef HAVE_XPETRA_TPETRA
166  // @{
168  /*static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec);
169  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec);
170  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec);
171 
172  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec);
173  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec);
174 
175  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op);
176  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op);
177 
178  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op);
179  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op);
180 
181  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op);
182  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op);*/
183 
184 
186  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
187  if (tmp_TMap == Teuchos::null)
188  throw Exceptions::BadCast("Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
189  return tmp_TMap->getTpetra_Map();
190  }
191 #endif
192 
193 
195 
196 
197  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
199 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
200  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
201  if (tmp_EMap != Teuchos::null) {
202  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
203  if (rv != 0)
204  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
205  return;
206  }
207 #endif // HAVE_XPETRA_EPETRAEXT
208 
209 #ifdef HAVE_XPETRA_TPETRA
211  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
212  if (tmp_TMap != Teuchos::null) {
213  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
214  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
215  return;
216  }
217 #endif // HAVE_XPETRA_TPETRA
218 
219  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
220 
221  } //Write
222 
224  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
225  std::string mapfile = "map_" + fileName;
226  Write(mapfile, *(vec.getMap()));
227 
229 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
230  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
231  if (tmp_EVec != Teuchos::null) {
232  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
233  if (rv != 0)
234  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
235  return;
236  }
237 #endif // HAVE_XPETRA_EPETRA
238 
239 #ifdef HAVE_XPETRA_TPETRA
241  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
242  if (tmp_TVec != Teuchos::null) {
243  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
244  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
245  return;
246  }
247 #endif // HAVE_XPETRA_TPETRA
248 
249  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
250 
251  } //Write
252 
253 
255  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
256 
257  Write("rowmap_" + fileName, *(Op.getRowMap()));
258  Write("colmap_" + fileName, *(Op.getColMap()));
259  Write("domainmap_" + fileName, *(Op.getDomainMap()));
260  Write("rangemap_" + fileName, *(Op.getRangeMap()));
261 
265 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
266  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
267  if (tmp_ECrsMtx != Teuchos::null) {
268  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
269  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
270  if (rv != 0)
271  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
272  return;
273  }
274 #endif
275 
276 #ifdef HAVE_XPETRA_TPETRA
278  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
279  if (tmp_TCrsMtx != Teuchos::null) {
280  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
281  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
282  return;
283  }
284 #endif // HAVE_XPETRA_TPETRA
285 
286  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
287  } //Write
288 
290  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) {
291  if (binary == false) {
292  // Matrix Market file format (ASCII)
293  if (lib == Xpetra::UseEpetra) {
294 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
295  Epetra_CrsMatrix *eA;
296  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
297  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
298  if (rv != 0)
299  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
300 
301  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
302 
304  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
305  return A;
306 #else
307  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
308 #endif
309  } else if (lib == Xpetra::UseTpetra) {
310 #ifdef HAVE_XPETRA_TPETRA
311  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
312 
313  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
314 
315  //RCP<Node> node = Xpetra::DefaultPlatform::getDefaultPlatform().getNode();
317  RCP<Node> node = rcp(new Node(pl));
318  bool callFillComplete = true;
319 
320  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, node, callFillComplete);
321 
322  if (tA.is_null())
323  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
324 
328 
329  return A;
330 #else
331  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
332 #endif
333  } else {
334  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
335  }
336  } else {
337  // Custom file format (binary)
338  std::ifstream ifs(fileName.c_str(), std::ios::binary);
339  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
340  int m, n, nnz;
341  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
342  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
343  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
344 
345  int myRank = comm->getRank();
346 
347  GO indexBase = 0;
348  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
349  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
351 
352  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
353 
354  if (myRank == 0) {
357  for (int i = 0; i < m; i++) {
358  int row, rownnz;
359  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
360  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
361  inds.resize(rownnz);
362  vals.resize(rownnz);
363  for (int j = 0; j < rownnz; j++) {
364  int index;
365  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
366  inds[j] = Teuchos::as<GlobalOrdinal>(index);
367  }
368  for (int j = 0; j < rownnz; j++) {
369  double value;
370  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
371  vals[j] = Teuchos::as<SC>(value);
372  }
373  A->insertGlobalValues(row, inds, vals);
374  }
375  }
376 
377  A->fillComplete(domainMap, rangeMap);
378 
379  return A;
380  }
381 
382  return Teuchos::null;
383 
384  } //Read()
385 
386 
392  Read(const std::string& filename,
394  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
395  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
396  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
397  const bool callFillComplete = true,
398  const bool binary = false,
399  const bool tolerant = false,
400  const bool debug = false) {
401  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
402 
403  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
404  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
405 
406  const Xpetra::UnderlyingLib lib = rowMap->lib();
407  if (binary == false) {
408  if (lib == Xpetra::UseEpetra) {
409 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
410  Epetra_CrsMatrix *eA;
411  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
412  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
413  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
414  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
415  int rv;
416  if (colMap.is_null()) {
417  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
418 
419  } else {
420  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
421  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
422  }
423 
424  if (rv != 0)
425  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
426 
427  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
429  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
430 
431  return A;
432 #else
433  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
434 #endif
435  } else if (lib == Xpetra::UseTpetra) {
436 #ifdef HAVE_XPETRA_TPETRA
437  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
438  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
439  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
440 
441  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
442  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
443  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
444  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
445 
446  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
447  callFillComplete, tolerant, debug);
448  if (tA.is_null())
449  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
450 
454 
455  return A;
456 #else
457  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
458 #endif
459  } else {
460  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
461  }
462  } else {
463  // Custom file format (binary)
464  std::ifstream ifs(filename.c_str(), std::ios::binary);
465  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
466  int m, n, nnz;
467  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
468  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
469  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
470 
472 
473  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
474 
475  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
476  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
477 
480  for (int i = 0; i < m; i++) {
481  int row, rownnz;
482  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
483  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
484  inds.resize(rownnz);
485  vals.resize(rownnz);
486  for (int j = 0; j < rownnz; j++) {
487  int index;
488  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
489  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
490  }
491  for (int j = 0; j < rownnz; j++) {
492  double value;
493  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
494  vals[j] = Teuchos::as<SC>(value);
495  }
496  A->insertGlobalValues(rowElements[row], inds, vals);
497  }
498  A->fillComplete(domainMap, rangeMap);
499  return A;
500  }
501 
502  return Teuchos::null;
503  }
505 
506 
507  static RCP<MultiVector> ReadMultiVector (const std::string& fileName, const RCP<const Map>& map) {
508  Xpetra::UnderlyingLib lib = map->lib();
509 
510  if (lib == Xpetra::UseEpetra) {
511  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
512 
513  } else if (lib == Xpetra::UseTpetra) {
514 #ifdef HAVE_XPETRA_TPETRA
515  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
516  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
517  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
518  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
519 
520  RCP<const map_type> temp = toTpetra(map);
521  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),map->getNode(),temp);
523  return rmv;
524 #else
525  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
526 #endif
527  } else {
528  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
529  }
530 
531  return Teuchos::null;
532  }
533 
534  static RCP<const Map> ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
535  if (lib == Xpetra::UseEpetra) {
536  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
537  } else if (lib == Xpetra::UseTpetra) {
538 #ifdef HAVE_XPETRA_TPETRA
539  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
540  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
541 
542  RCP<Node> node = rcp(new Node());
543 
544  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm, node);
545  if (tMap.is_null())
546  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
547 
548  return Xpetra::toXpetra(tMap);
549 #else
550  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
551 #endif
552  } else {
553  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
554  }
555 
556  return Teuchos::null;
557 
558  }
559 
560 
561  };
562 
563 
564 #ifdef HAVE_XPETRA_EPETRA
565 
574  template <class Scalar>
575  class IO<Scalar,int,int,EpetraNode> {
576  public:
577  typedef int LocalOrdinal;
578  typedef int GlobalOrdinal;
579  typedef EpetraNode Node;
580 
581 #ifdef HAVE_XPETRA_EPETRA
582  // @{
584  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
585  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
586  if (xeMap == Teuchos::null)
587  throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
588  return xeMap->getEpetra_Map();
589  }
590  // @}
591 #endif
592 
593 #ifdef HAVE_XPETRA_TPETRA
594  // @{
597  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
598  if (tmp_TMap == Teuchos::null)
599  throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
600  return tmp_TMap->getTpetra_Map();
601  }
602 #endif
603 
604 
606 
607 
608  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
610 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
611  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
612  if (tmp_EMap != Teuchos::null) {
613  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
614  if (rv != 0)
615  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
616  return;
617  }
618 #endif // HAVE_XPETRA_EPETRA
619 
620 #ifdef HAVE_XPETRA_TPETRA
621 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
622  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
623  // do nothing
624 # else
626  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
627  if (tmp_TMap != Teuchos::null) {
628  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
629  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
630  return;
631  }
632 # endif
633 #endif // HAVE_XPETRA_TPETRA
634  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
635  }
636 
638  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
639  std::string mapfile = "map_" + fileName;
640  Write(mapfile, *(vec.getMap()));
641 
643 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
644  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
645  if (tmp_EVec != Teuchos::null) {
646  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
647  if (rv != 0)
648  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
649  return;
650  }
651 #endif // HAVE_XPETRA_EPETRAEXT
652 
653 #ifdef HAVE_XPETRA_TPETRA
654 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
655  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
656  // do nothin
657 # else
659  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
660  if (tmp_TVec != Teuchos::null) {
661  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
662  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
663  return;
664  }
665 # endif
666 #endif // HAVE_XPETRA_TPETRA
667 
668  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
669 
670  } //Write
671 
672 
673 
675  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
676 
677  Write("rowmap_" + fileName, *(Op.getRowMap()));
678  Write("colmap_" + fileName, *(Op.getColMap()));
679  Write("domainmap_" + fileName, *(Op.getDomainMap()));
680  Write("rangemap_" + fileName, *(Op.getRangeMap()));
681 
685 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
686  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
687  if (tmp_ECrsMtx != Teuchos::null) {
688  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
689  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
690  if (rv != 0)
691  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
692  return;
693  }
694 #endif // endif HAVE_XPETRA_EPETRA
695 
696 #ifdef HAVE_XPETRA_TPETRA
697 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
698  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
699  // do nothin
700 # else
702  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
703  if (tmp_TCrsMtx != Teuchos::null) {
704  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
705  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
706  return;
707  }
708 # endif
709 #endif // HAVE_XPETRA_TPETRA
710 
711  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
712  } //Write
713 
715  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) {
716  if (binary == false) {
717  // Matrix Market file format (ASCII)
718  if (lib == Xpetra::UseEpetra) {
719 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
720  Epetra_CrsMatrix *eA;
721  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
722  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
723  if (rv != 0)
724  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
725 
726  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
727 
729  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
730  return A;
731 #else
732  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
733 #endif
734  } else if (lib == Xpetra::UseTpetra) {
735 #ifdef HAVE_XPETRA_TPETRA
736 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
737  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
738  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
739 # else
740  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
741 
742  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
743 
744  //RCP<Node> node = Xpetra::DefaultPlatform::getDefaultPlatform().getNode();
746  RCP<Node> node = rcp(new Node(pl));
747  bool callFillComplete = true;
748 
749  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, node, callFillComplete);
750 
751  if (tA.is_null())
752  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
753 
757 
758  return A;
759 # endif
760 #else
761  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
762 #endif
763  } else {
764  throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
765  }
766  } else {
767  // Custom file format (binary)
768  std::ifstream ifs(fileName.c_str(), std::ios::binary);
769  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
770  int m, n, nnz;
771  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
772  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
773  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
774 
775  int myRank = comm->getRank();
776 
777  GlobalOrdinal indexBase = 0;
778  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
779  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
781 
782  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
783 
784  if (myRank == 0) {
787  for (int i = 0; i < m; i++) {
788  int row, rownnz;
789  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
790  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
791  inds.resize(rownnz);
792  vals.resize(rownnz);
793  for (int j = 0; j < rownnz; j++) {
794  int index;
795  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
796  inds[j] = Teuchos::as<GlobalOrdinal>(index);
797  }
798  for (int j = 0; j < rownnz; j++) {
799  double value;
800  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
801  vals[j] = Teuchos::as<Scalar>(value);
802  }
803  A->insertGlobalValues(row, inds, vals);
804  }
805  }
806 
807  A->fillComplete(domainMap, rangeMap);
808 
809  return A;
810  }
811 
812  return Teuchos::null;
813 
814  } //Read()
815 
816 
823  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
824  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
825  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
826  const bool callFillComplete = true,
827  const bool binary = false,
828  const bool tolerant = false,
829  const bool debug = false) {
830  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
831 
832  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
833  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
834 
835  const Xpetra::UnderlyingLib lib = rowMap->lib();
836  if (binary == false) {
837  if (lib == Xpetra::UseEpetra) {
838 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
839  Epetra_CrsMatrix *eA;
840  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
841  const Epetra_Map& epetraRowMap = Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rowMap);
842  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
843  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
844  int rv;
845  if (colMap.is_null()) {
846  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
847 
848  } else {
849  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
850  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
851  }
852 
853  if (rv != 0)
854  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
855 
856  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
858  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
859 
860  return A;
861 #else
862  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
863 #endif
864  } else if (lib == Xpetra::UseTpetra) {
865 #ifdef HAVE_XPETRA_TPETRA
866 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
867  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
868  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
869 # else
870  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
871  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
872  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
873 
874  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
875  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
876  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
877  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
878 
879  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
880  callFillComplete, tolerant, debug);
881  if (tA.is_null())
882  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
883 
887 
888  return A;
889 # endif
890 #else
891  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
892 #endif
893  } else {
894  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
895  }
896  } else {
897  // Custom file format (binary)
898  std::ifstream ifs(filename.c_str(), std::ios::binary);
899  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
900  int m, n, nnz;
901  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
902  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
903  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
904 
906 
907  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
908 
909  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
910  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
911 
914  for (int i = 0; i < m; i++) {
915  int row, rownnz;
916  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
917  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
918  inds.resize(rownnz);
919  vals.resize(rownnz);
920  for (int j = 0; j < rownnz; j++) {
921  int index;
922  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
923  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
924  }
925  for (int j = 0; j < rownnz; j++) {
926  double value;
927  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
928  vals[j] = Teuchos::as<Scalar>(value);
929  }
930  A->insertGlobalValues(rowElements[row], inds, vals);
931  }
932  A->fillComplete(domainMap, rangeMap);
933  return A;
934  }
935 
936  return Teuchos::null;
937  }
939 
940 
942  Xpetra::UnderlyingLib lib = map->lib();
943 
944  if (lib == Xpetra::UseEpetra) {
945  // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
946  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
947 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
948  Epetra_MultiVector * MV;
949  EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
950  return Xpetra::toXpetra<int,Node>(rcp(MV));
951 #else
952  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
953 #endif
954  } else if (lib == Xpetra::UseTpetra) {
955 #ifdef HAVE_XPETRA_TPETRA
956 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
957  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
958  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
959 # else
960  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
961  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
962  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
963  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
964 
965  RCP<const map_type> temp = toTpetra(map);
966  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),map->getNode(),temp);
968  return rmv;
969 # endif
970 #else
971  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
972 #endif
973  } else {
974  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
975  }
976 
977  return Teuchos::null;
978 
979  }
980 
981 
982  static RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
983  if (lib == Xpetra::UseEpetra) {
984  // do we need another specialization for <double,int,int> ??
985  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
986 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
987  Epetra_Map *eMap;
988  int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
989  if (rv != 0)
990  throw Exceptions::RuntimeError("Error reading matrix with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
991 
992  RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
993  return Xpetra::toXpetra<int,Node>(*eMap1);
994 #else
995  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
996 #endif
997  } else if (lib == Xpetra::UseTpetra) {
998 #ifdef HAVE_XPETRA_TPETRA
999 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1000  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1001  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1002 # else
1003  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1004  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1005 
1006  RCP<Node> node = rcp(new Node());
1007 
1008  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm, node);
1009  if (tMap.is_null())
1010  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1011 
1012  return Xpetra::toXpetra(tMap);
1013 # endif
1014 #else
1015  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1016 #endif
1017  } else {
1018  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1019  }
1020 
1021  return Teuchos::null;
1022 
1023  }
1024 
1025 
1026  };
1027 #endif // HAVE_XPETRA_EPETRA
1028 
1029 
1030 } // end namespace Xpetra
1031 
1032 #define XPETRA_IO_SHORT
1033 
1034 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
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:392
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:821
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:197
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:675
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 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:596
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:185
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
Definition: Xpetra_IO.hpp:534
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:715
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
Definition: Xpetra_IO.hpp:507
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
void resize(size_type new_size, const value_type &x=value_type())
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:941
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< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)
Definition: Xpetra_IO.hpp:103
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:290
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:156
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...
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:982
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:584
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:255
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
Definition: Xpetra_IO.hpp:133
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:608
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:224
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:638
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)
bool is_null() const