All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Xpetra_IO.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 
48 #ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_
50 
51 #include <fstream>
52 #include "Xpetra_ConfigDefs.hpp"
53 
54 #ifdef HAVE_XPETRA_EPETRA
55 # ifdef HAVE_MPI
56 # include "Epetra_MpiComm.h"
57 # endif
58 #endif
59 
60 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
61 #include <EpetraExt_MatrixMatrix.h>
62 #include <EpetraExt_RowMatrixOut.h>
63 #include <EpetraExt_MultiVectorOut.h>
64 #include <EpetraExt_CrsMatrixIn.h>
65 #include <EpetraExt_MultiVectorIn.h>
66 #include <EpetraExt_BlockMapIn.h>
67 #include <Xpetra_EpetraUtils.hpp>
69 #include <EpetraExt_BlockMapOut.h>
70 #endif
71 
72 #ifdef HAVE_XPETRA_TPETRA
73 #include <MatrixMarket_Tpetra.hpp>
74 #include <Tpetra_RowMatrixTransposer.hpp>
75 #include <TpetraExt_MatrixMatrix.hpp>
76 #include <Xpetra_TpetraMultiVector.hpp>
77 #include <Xpetra_TpetraCrsMatrix.hpp>
79 #endif
80 
81 #ifdef HAVE_XPETRA_EPETRA
82 #include <Xpetra_EpetraMap.hpp>
83 #endif
84 
85 #include "Xpetra_Matrix.hpp"
86 #include "Xpetra_MatrixMatrix.hpp"
87 #include "Xpetra_CrsMatrixWrap.hpp"
89 
90 #include "Xpetra_Map.hpp"
91 #include "Xpetra_StridedMap.hpp"
93 #include "Xpetra_MapExtractor.hpp"
94 #include "Xpetra_MatrixFactory.hpp"
95 
96 #include <Teuchos_MatrixMarket_Raw_Writer.hpp>
97 #include <string>
98 
99 
100 namespace Xpetra {
101 
102 
103 #ifdef HAVE_XPETRA_EPETRA
104  // This non-member templated function exists so that the matrix-matrix multiply will compile if Epetra, Tpetra, and ML are enabled.
105  template<class SC,class LO,class GO,class NO>
106  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >
109  "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
110  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
111  }
112 
113  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
114  template<>
115  inline RCP<Xpetra::CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_CrsMatrix> &epAB) {
116  typedef double SC;
117  typedef int LO;
118  typedef int GO;
119  typedef Xpetra::EpetraNode NO;
120 
122  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC,LO,GO,NO> >(tmpC1);
124 
125  return tmpC3;
126  }
127 
128 
129  template<class SC,class LO,class GO,class NO>
130  RCP<Xpetra::MultiVector<SC,LO,GO,NO> >
133  "Convert_Epetra_MultiVector_ToXpetra_MultiVector cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
134  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
135  }
136 
137  // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
138  template<>
139  inline RCP<Xpetra::MultiVector<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_MultiVector_ToXpetra_MultiVector<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_MultiVector> &epX) {
140  typedef double SC;
141  typedef int LO;
142  typedef int GO;
143  typedef Xpetra::EpetraNode NO;
144 
145  RCP<Xpetra::MultiVector<SC,LO,GO,NO >> tmp = Xpetra::toXpetra<GO,NO>(epX);
146  return tmp;
147  }
148 
149 #endif
150 
155  template <class Scalar,
156  class LocalOrdinal = int,
157  class GlobalOrdinal = LocalOrdinal,
159  class IO {
160 
161  private:
162 #undef XPETRA_IO_SHORT
163 #include "Xpetra_UseShortNames.hpp"
164 
165  public:
166 
167 #ifdef HAVE_XPETRA_EPETRA
168  // @{
170  /*static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec);
171  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec);
172 
173  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec);
174  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec);
175 
176  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op);
177  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op);
178 
179  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op);
180  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op);*/
181 
183  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
184  if (xeMap == Teuchos::null)
185  throw Exceptions::BadCast("Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
186  return xeMap->getEpetra_Map();
187  }
188  // @}
189 #endif
190 
191 #ifdef HAVE_XPETRA_TPETRA
192  // @{
194  /*static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec);
195  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec);
196  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec);
197 
198  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec);
199  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec);
200 
201  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op);
202  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op);
203 
204  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op);
205  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op);
206 
207  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op);
208  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op);*/
209 
210 
212  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
213  if (tmp_TMap == Teuchos::null)
214  throw Exceptions::BadCast("Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
215  return tmp_TMap->getTpetra_Map();
216  }
217 #endif
218 
219 
221 
222 
223  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
225 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
226  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
227  if (tmp_EMap != Teuchos::null) {
228  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
229  if (rv != 0)
230  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
231  return;
232  }
233 #endif // HAVE_XPETRA_EPETRAEXT
234 
235 #ifdef HAVE_XPETRA_TPETRA
237  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
238  if (tmp_TMap != Teuchos::null) {
239  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
240  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
241  return;
242  }
243 #endif // HAVE_XPETRA_TPETRA
244 
245  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
246 
247  } //Write
248 
250  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
251  std::string mapfile = "map_" + fileName;
252  Write(mapfile, *(vec.getMap()));
253 
255 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
256  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
257  if (tmp_EVec != Teuchos::null) {
258  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
259  if (rv != 0)
260  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
261  return;
262  }
263 #endif // HAVE_XPETRA_EPETRA
264 
265 #ifdef HAVE_XPETRA_TPETRA
267  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
268  if (tmp_TVec != Teuchos::null) {
269  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
270  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
271  return;
272  }
273 #endif // HAVE_XPETRA_TPETRA
274 
275  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
276 
277  } //Write
278 
279 
281  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
282 
283  Write("rowmap_" + fileName, *(Op.getRowMap()));
284  Write("colmap_" + fileName, *(Op.getColMap()));
285  Write("domainmap_" + fileName, *(Op.getDomainMap()));
286  Write("rangemap_" + fileName, *(Op.getRangeMap()));
287 
291 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
292  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
293  if (tmp_ECrsMtx != Teuchos::null) {
294  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
295  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
296  if (rv != 0)
297  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
298  return;
299  }
300 #endif
301 
302 #ifdef HAVE_XPETRA_TPETRA
304  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
305  if (tmp_TCrsMtx != Teuchos::null) {
306  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
307  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
308  return;
309  }
310 #endif // HAVE_XPETRA_TPETRA
311 
312  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
313  } //Write
314 
315 
317  static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
321 
322  ArrayRCP<const size_t> rowptr_RCP;
323  ArrayRCP<LocalOrdinal> rowptr2_RCP;
324  ArrayRCP<const LocalOrdinal> colind_RCP;
325  ArrayRCP<const Scalar> vals_RCP;
326  tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
327 
328  ArrayView<const size_t> rowptr = rowptr_RCP();
329  ArrayView<const LocalOrdinal> colind = colind_RCP();
330  ArrayView<const Scalar> vals = vals_RCP();
331 
332  rowptr2_RCP.resize(rowptr.size());
333  ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
334  for (size_t j = 0; j<rowptr.size(); j++)
335  rowptr2[j] = rowptr[j];
336 
338  writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
339  rowptr2,colind,vals,
340  rowptr.size()-1,Op.getColMap()->getNodeNumElements());
341  } //WriteLocal
342 
343 
348  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
352 
353  // write all matrices with their maps
354  for (size_t r = 0; r < Op.Rows(); ++r) {
355  for (size_t c = 0; c < Op.Cols(); ++c) {
356  RCP<const XpMat > m = Op.getMatrix(r,c);
357  if(m != Teuchos::null) { // skip empty blocks
358  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const XpCrsMatWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
359  "Sub block matrix (" << r << "," << c << ") is not of type CrsMatrixWrap.");
360  XpIO::Write(fileName + toString(r) + toString(c) + ".m", *m);
361  }
362  }
363  }
364 
365  // write map information of map extractors
366  RCP<const XpMapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
367  RCP<const XpMapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
368 
369  for(size_t r = 0; r < rangeMapExtractor->NumMaps(); ++r) {
370  RCP<const XpMap> map = rangeMapExtractor->getMap(r);
371  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", *map);
372  }
373  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
374 
375  for(size_t c = 0; c < domainMapExtractor->NumMaps(); ++c) {
376  RCP<const XpMap> map = domainMapExtractor->getMap(c);
377  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", *map);
378  }
379  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
380  } //WriteBlockCrsMatrix
381 
383  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) {
384  if (binary == false) {
385  // Matrix Market file format (ASCII)
386  if (lib == Xpetra::UseEpetra) {
387 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
388  Epetra_CrsMatrix *eA;
389  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
390  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
391  if (rv != 0)
392  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
393 
394  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
395 
397  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
398  return A;
399 #else
400  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
401 #endif
402  } else if (lib == Xpetra::UseTpetra) {
403 #ifdef HAVE_XPETRA_TPETRA
404  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
405 
406  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
407 
408  bool callFillComplete = true;
409 
410  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
411 
412  if (tA.is_null())
413  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
414 
418 
419  return A;
420 #else
421  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
422 #endif
423  } else {
424  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
425  }
426  } else {
427  // Custom file format (binary)
428  std::ifstream ifs(fileName.c_str(), std::ios::binary);
429  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
430  int m, n, nnz;
431  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
432  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
433  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
434 
435  int myRank = comm->getRank();
436 
437  GO indexBase = 0;
438  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
439  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
441 
442  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
443 
444  if (myRank == 0) {
447  for (int i = 0; i < m; i++) {
448  int row, rownnz;
449  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
450  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
451  inds.resize(rownnz);
452  vals.resize(rownnz);
453  for (int j = 0; j < rownnz; j++) {
454  int index;
455  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
456  inds[j] = Teuchos::as<GlobalOrdinal>(index);
457  }
458  for (int j = 0; j < rownnz; j++) {
459  double value;
460  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
461  vals[j] = Teuchos::as<SC>(value);
462  }
463  A->insertGlobalValues(row, inds, vals);
464  }
465  }
466 
467  A->fillComplete(domainMap, rangeMap);
468 
469  return A;
470  }
471 
472  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
473 
474  } //Read()
475 
476 
482  Read(const std::string& filename,
484  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
485  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
486  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
487  const bool callFillComplete = true,
488  const bool binary = false,
489  const bool tolerant = false,
490  const bool debug = false) {
491  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
492 
493  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
494  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
495 
496  const Xpetra::UnderlyingLib lib = rowMap->lib();
497  if (binary == false) {
498  if (lib == Xpetra::UseEpetra) {
499 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
500  Epetra_CrsMatrix *eA;
501  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
503  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
504  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
505  int rv;
506  if (colMap.is_null()) {
507  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
508 
509  } else {
510  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
511  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
512  }
513 
514  if (rv != 0)
515  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
516 
517  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
519  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
520 
521  return A;
522 #else
523  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
524 #endif
525  } else if (lib == Xpetra::UseTpetra) {
526 #ifdef HAVE_XPETRA_TPETRA
527  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
528  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
529  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
530 
531  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
532  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
533  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
534  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
535 
536  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
537  callFillComplete, tolerant, debug);
538  if (tA.is_null())
539  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
540 
544 
545  return A;
546 #else
547  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
548 #endif
549  } else {
550  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
551  }
552  } else {
553  // Custom file format (binary)
554  std::ifstream ifs(filename.c_str(), std::ios::binary);
555  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
556  int m, n, nnz;
557  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
558  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
559  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
560 
562 
563  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
564 
565  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
566  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
567 
570  for (int i = 0; i < m; i++) {
571  int row, rownnz;
572  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
573  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
574  inds.resize(rownnz);
575  vals.resize(rownnz);
576  for (int j = 0; j < rownnz; j++) {
577  int index;
578  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
579  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
580  }
581  for (int j = 0; j < rownnz; j++) {
582  double value;
583  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
584  vals[j] = Teuchos::as<SC>(value);
585  }
586  A->insertGlobalValues(rowElements[row], inds, vals);
587  }
588  A->fillComplete(domainMap, rangeMap);
589  return A;
590  }
591 
592  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
593  }
595 
596 
597  static RCP<MultiVector> ReadMultiVector (const std::string& fileName, const RCP<const Map>& map) {
598  Xpetra::UnderlyingLib lib = map->lib();
599 
600  if (lib == Xpetra::UseEpetra) {
601  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
602 
603  } else if (lib == Xpetra::UseTpetra) {
604 #ifdef HAVE_XPETRA_TPETRA
605  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
606  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
607  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
608  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
609 
610  RCP<const map_type> temp = toTpetra(map);
611  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
613  return rmv;
614 #else
615  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
616 #endif
617  } else {
618  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
619  }
620 
621  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
622  }
623 
624  static RCP<const Map> ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
625  if (lib == Xpetra::UseEpetra) {
626  TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
627  } else if (lib == Xpetra::UseTpetra) {
628 #ifdef HAVE_XPETRA_TPETRA
629  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
630  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
631 
632  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
633  if (tMap.is_null())
634  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
635 
636  return Xpetra::toXpetra(tMap);
637 #else
638  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
639 #endif
640  } else {
641  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
642  }
643 
644  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
645 
646  }
647 
652  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
653  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
657 
658  size_t numBlocks = 2; // TODO user parameter?
659 
660  std::vector<RCP<const XpMap> > rgMapVec;
661  for(size_t r = 0; r < numBlocks; ++r) {
662  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
663  rgMapVec.push_back(map);
664  }
665  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
666 
667  std::vector<RCP<const XpMap> > doMapVec;
668  for(size_t c = 0; c < numBlocks; ++c) {
669  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
670  doMapVec.push_back(map);
671  }
672  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
673 
674  /*std::vector<RCP<const XpMap> > testRgMapVec;
675  for(size_t r = 0; r < numBlocks; ++r) {
676  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
677  testRgMapVec.push_back(map);
678  }
679  std::vector<RCP<const XpMap> > testDoMapVec;
680  for(size_t c = 0; c < numBlocks; ++c) {
681  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
682  testDoMapVec.push_back(map);
683  }*/
684 
685  // create map extractors
686 
687  // range map extractor
688  bool bRangeUseThyraStyleNumbering = false;
689  /*GlobalOrdinal gMinGids = 0;
690  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
691  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
692  }
693  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
694  */
695  RCP<const XpMapExtractor> rangeMapExtractor =
696  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
697 
698 
699  // domain map extractor
700  bool bDomainUseThyraStyleNumbering = false;
701  /*gMinGids = 0;
702  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
703  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
704  }
705  if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
706  */
707  RCP<const XpMapExtractor> domainMapExtractor =
708  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
709 
710  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
711 
712  // write all matrices with their maps
713  for (size_t r = 0; r < numBlocks; ++r) {
714  for (size_t c = 0; c < numBlocks; ++c) {
715  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
716  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
717  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
718  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
719  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
720  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
721  bOp->setMatrix(r, c, mat);
722  }
723  }
724 
725  bOp->fillComplete();
726 
727  return bOp;
728  } //ReadBlockedCrsMatrix
729 
730 
732  template<class T>
733  static std::string toString(const T& what) {
734  std::ostringstream buf;
735  buf << what;
736  return buf.str();
737  }
738  };
739 
740 
741 #ifdef HAVE_XPETRA_EPETRA
742 
751  template <class Scalar>
752  class IO<Scalar,int,int,EpetraNode> {
753  public:
754  typedef int LocalOrdinal;
755  typedef int GlobalOrdinal;
756  typedef EpetraNode Node;
757 
758 #ifdef HAVE_XPETRA_EPETRA
759  // @{
762  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
763  if (xeMap == Teuchos::null)
764  throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
765  return xeMap->getEpetra_Map();
766  }
767  // @}
768 #endif
769 
770 #ifdef HAVE_XPETRA_TPETRA
771  // @{
774  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
775  if (tmp_TMap == Teuchos::null)
776  throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
777  return tmp_TMap->getTpetra_Map();
778  }
779 #endif
780 
781 
783 
784 
785  static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
787 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
788  const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
789  if (tmp_EMap != Teuchos::null) {
790  int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
791  if (rv != 0)
792  throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
793  return;
794  }
795 #endif // HAVE_XPETRA_EPETRA
796 
797 #ifdef HAVE_XPETRA_TPETRA
798 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
799  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
800  // do nothing
801 # else
803  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
804  if (tmp_TMap != Teuchos::null) {
805  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
806  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
807  return;
808  }
809 # endif
810 #endif // HAVE_XPETRA_TPETRA
811  throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
812  }
813 
815  static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
816  std::string mapfile = "map_" + fileName;
817  Write(mapfile, *(vec.getMap()));
818 
820 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
821  const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
822  if (tmp_EVec != Teuchos::null) {
823  int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
824  if (rv != 0)
825  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
826  return;
827  }
828 #endif // HAVE_XPETRA_EPETRAEXT
829 
830 #ifdef HAVE_XPETRA_TPETRA
831 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
832  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
833  // do nothin
834 # else
836  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
837  if (tmp_TVec != Teuchos::null) {
838  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
839  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
840  return;
841  }
842 # endif
843 #endif // HAVE_XPETRA_TPETRA
844 
845  throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
846 
847  } //Write
848 
849 
850 
852  static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
853 
854  Write("rowmap_" + fileName, *(Op.getRowMap()));
855  Write("colmap_" + fileName, *(Op.getColMap()));
856  Write("domainmap_" + fileName, *(Op.getDomainMap()));
857  Write("rangemap_" + fileName, *(Op.getRangeMap()));
858 
862 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
863  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
864  if (tmp_ECrsMtx != Teuchos::null) {
865  RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
866  int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
867  if (rv != 0)
868  throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
869  return;
870  }
871 #endif // endif HAVE_XPETRA_EPETRA
872 
873 #ifdef HAVE_XPETRA_TPETRA
874 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
875  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
876  // do nothin
877 # else
879  Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
880  if (tmp_TCrsMtx != Teuchos::null) {
881  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
882  Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
883  return;
884  }
885 # endif
886 #endif // HAVE_XPETRA_TPETRA
887 
888  throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
889  } //Write
890 
891 
893  static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
897 
898  ArrayRCP<const size_t> rowptr_RCP;
899  ArrayRCP<LocalOrdinal> rowptr2_RCP;
900  ArrayRCP<const LocalOrdinal> colind_RCP;
901  ArrayRCP<const Scalar> vals_RCP;
902  tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
903 
904  ArrayView<const size_t> rowptr = rowptr_RCP();
905  ArrayView<const LocalOrdinal> colind = colind_RCP();
906  ArrayView<const Scalar> vals = vals_RCP();
907 
908  rowptr2_RCP.resize(rowptr.size());
909  ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
910  for (size_t j = 0; j<rowptr.size(); j++)
911  rowptr2[j] = rowptr[j];
912 
914  writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
915  rowptr2,colind,vals,
916  rowptr.size()-1,Op.getColMap()->getNodeNumElements());
917  } //WriteLocal
918 
923  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
927 
928  // write all matrices with their maps
929  for (size_t r = 0; r < Op.Rows(); ++r) {
930  for (size_t c = 0; c < Op.Cols(); ++c) {
931  RCP<const XpMat > m = Op.getMatrix(r,c);
932  if(m != Teuchos::null) { // skip empty blocks
933  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const XpCrsMatWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
934  "Sub block matrix (" << r << "," << c << ") is not of type CrsMatrixWrap.");
935  XpIO::Write(fileName + toString(r) + toString(c) + ".m", *m);
936  }
937  }
938  }
939 
940  // write map information of map extractors
941  RCP<const XpMapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
942  RCP<const XpMapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
943 
944  for(size_t r = 0; r < rangeMapExtractor->NumMaps(); ++r) {
945  RCP<const XpMap> map = rangeMapExtractor->getMap(r);
946  XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", *map);
947  }
948  XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
949 
950  for(size_t c = 0; c < domainMapExtractor->NumMaps(); ++c) {
951  RCP<const XpMap> map = domainMapExtractor->getMap(c);
952  XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", *map);
953  }
954  XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
955  } //WriteBlockCrsMatrix
956 
958  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) {
959  if (binary == false) {
960  // Matrix Market file format (ASCII)
961  if (lib == Xpetra::UseEpetra) {
962 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
963  Epetra_CrsMatrix *eA;
964  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
965  int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
966  if (rv != 0)
967  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
968 
969  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
970 
972  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
973  return A;
974 #else
975  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
976 #endif
977  } else if (lib == Xpetra::UseTpetra) {
978 #ifdef HAVE_XPETRA_TPETRA
979 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
980  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
981  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
982 # else
983  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
984 
985  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
986 
987  bool callFillComplete = true;
988 
989  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
990 
991  if (tA.is_null())
992  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
993 
997 
998  return A;
999 # endif
1000 #else
1001  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1002 #endif
1003  } else {
1004  throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1005  }
1006  } else {
1007  // Custom file format (binary)
1008  std::ifstream ifs(fileName.c_str(), std::ios::binary);
1009  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
1010  int m, n, nnz;
1011  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1012  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1013  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1014 
1015  int myRank = comm->getRank();
1016 
1017  GlobalOrdinal indexBase = 0;
1018  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
1019  RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
1021 
1022  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1023 
1024  if (myRank == 0) {
1027  for (int i = 0; i < m; i++) {
1028  int row, rownnz;
1029  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1030  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1031  inds.resize(rownnz);
1032  vals.resize(rownnz);
1033  for (int j = 0; j < rownnz; j++) {
1034  int index;
1035  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1036  inds[j] = Teuchos::as<GlobalOrdinal>(index);
1037  }
1038  for (int j = 0; j < rownnz; j++) {
1039  double value;
1040  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1041  vals[j] = Teuchos::as<Scalar>(value);
1042  }
1043  A->insertGlobalValues(row, inds, vals);
1044  }
1045  }
1046 
1047  A->fillComplete(domainMap, rangeMap);
1048 
1049  return A;
1050  }
1051 
1052  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1053 
1054  } //Read()
1055 
1056 
1063  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
1064  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
1065  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
1066  const bool callFillComplete = true,
1067  const bool binary = false,
1068  const bool tolerant = false,
1069  const bool debug = false) {
1070  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
1071 
1072  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
1073  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
1074 
1075  const Xpetra::UnderlyingLib lib = rowMap->lib();
1076  if (binary == false) {
1077  if (lib == Xpetra::UseEpetra) {
1078 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1079  Epetra_CrsMatrix *eA;
1080  const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
1082  const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
1083  const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
1084  int rv;
1085  if (colMap.is_null()) {
1086  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
1087 
1088  } else {
1089  const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
1090  rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
1091  }
1092 
1093  if (rv != 0)
1094  throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1095 
1096  RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1098  Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1099 
1100  return A;
1101 #else
1102  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1103 #endif
1104  } else if (lib == Xpetra::UseTpetra) {
1105 #ifdef HAVE_XPETRA_TPETRA
1106 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1107  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1108  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1109 # else
1110  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1111  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1112  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1113 
1114  const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
1115  RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
1116  const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
1117  const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
1118 
1119  RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
1120  callFillComplete, tolerant, debug);
1121  if (tA.is_null())
1122  throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1123 
1127 
1128  return A;
1129 # endif
1130 #else
1131  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1132 #endif
1133  } else {
1134  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1135  }
1136  } else {
1137  // Custom file format (binary)
1138  std::ifstream ifs(filename.c_str(), std::ios::binary);
1139  TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
1140  int m, n, nnz;
1141  ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1142  ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1143  ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1144 
1146 
1147  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1148 
1149  Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
1150  Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
1151 
1154  for (int i = 0; i < m; i++) {
1155  int row, rownnz;
1156  ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1157  ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1158  inds.resize(rownnz);
1159  vals.resize(rownnz);
1160  for (int j = 0; j < rownnz; j++) {
1161  int index;
1162  ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1163  inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
1164  }
1165  for (int j = 0; j < rownnz; j++) {
1166  double value;
1167  ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1168  vals[j] = Teuchos::as<Scalar>(value);
1169  }
1170  A->insertGlobalValues(rowElements[row], inds, vals);
1171  }
1172  A->fillComplete(domainMap, rangeMap);
1173  return A;
1174  }
1175 
1176  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1177  }
1179 
1180 
1182  Xpetra::UnderlyingLib lib = map->lib();
1183 
1184  if (lib == Xpetra::UseEpetra) {
1185  // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
1186  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1187 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1188  Epetra_MultiVector * MV;
1189  EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
1190  RCP<Epetra_MultiVector> MVrcp = rcp(MV);
1191  return Convert_Epetra_MultiVector_ToXpetra_MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(MVrcp);
1192 #else
1193  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1194 #endif
1195  } else if (lib == Xpetra::UseTpetra) {
1196 #ifdef HAVE_XPETRA_TPETRA
1197 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1198  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1199  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1200 # else
1201  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1202  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1203  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1204  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
1205 
1206  RCP<const map_type> temp = toTpetra(map);
1207  RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
1209  return rmv;
1210 # endif
1211 #else
1212  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1213 #endif
1214  } else {
1215  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1216  }
1217 
1218  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1219 
1220  }
1221 
1222 
1223  static RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1224  if (lib == Xpetra::UseEpetra) {
1225  // do we need another specialization for <double,int,int> ??
1226  //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1227 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1228  Epetra_Map *eMap;
1229  int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
1230  if (rv != 0)
1231  throw Exceptions::RuntimeError("Error reading map from file " + fileName + " with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
1232 
1233  RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
1234  return Xpetra::toXpetra<int,Node>(*eMap1);
1235 #else
1236  throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1237 #endif
1238  } else if (lib == Xpetra::UseTpetra) {
1239 #ifdef HAVE_XPETRA_TPETRA
1240 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1241  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1242  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1243 # else
1244  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1245  typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1246 
1247  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
1248  if (tMap.is_null())
1249  throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1250 
1251  return Xpetra::toXpetra(tMap);
1252 # endif
1253 #else
1254  throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1255 #endif
1256  } else {
1257  throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1258  }
1259 
1260  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1261 
1262  }
1263 
1268  //typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
1269  //typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
1273 
1274 
1275  size_t numBlocks = 2; // TODO user parameter?
1276 
1277  std::vector<RCP<const XpMap> > rgMapVec;
1278  for(size_t r = 0; r < numBlocks; ++r) {
1279  RCP<const XpMap> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(r) + ".m", lib, comm);
1280  rgMapVec.push_back(map);
1281  }
1282  RCP<const XpMap> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
1283 
1284  std::vector<RCP<const XpMap> > doMapVec;
1285  for(size_t c = 0; c < numBlocks; ++c) {
1286  RCP<const XpMap> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(c) + ".m", lib, comm);
1287  doMapVec.push_back(map);
1288  }
1289  RCP<const XpMap> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
1290 
1291  /*std::vector<RCP<const XpMap> > testRgMapVec;
1292  for(size_t r = 0; r < numBlocks; ++r) {
1293  RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
1294  testRgMapVec.push_back(map);
1295  }
1296  std::vector<RCP<const XpMap> > testDoMapVec;
1297  for(size_t c = 0; c < numBlocks; ++c) {
1298  RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
1299  testDoMapVec.push_back(map);
1300  }*/
1301 
1302  // create map extractors
1303 
1304  // range map extractor
1305  bool bRangeUseThyraStyleNumbering = false;
1306  /*
1307  GlobalOrdinal gMinGids = 0;
1308  for(size_t v = 0; v < testRgMapVec.size(); ++v) {
1309  gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
1310  }
1311  if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;*/
1312  RCP<const XpMapExtractor> rangeMapExtractor =
1313  Teuchos::rcp(new XpMapExtractor(fullRangeMap, rgMapVec, bRangeUseThyraStyleNumbering));
1314 
1315  // domain map extractor
1316  bool bDomainUseThyraStyleNumbering = false;
1317  /*gMinGids = 0;
1318  for(size_t v = 0; v < testDoMapVec.size(); ++v) {
1319  gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
1320  }
1321  if ( gMinGids==0 && testDoMapVec.size() > 1) bDomainUseThyraStyleNumbering = true;*/
1322  RCP<const XpMapExtractor> domainMapExtractor =
1323  Teuchos::rcp(new XpMapExtractor(fullDomainMap, doMapVec, bDomainUseThyraStyleNumbering));
1324 
1325  RCP<XpBlockedCrsMat> bOp = Teuchos::rcp(new XpBlockedCrsMat(rangeMapExtractor,domainMapExtractor,33));
1326 
1327  // write all matrices with their maps
1328  for (size_t r = 0; r < numBlocks; ++r) {
1329  for (size_t c = 0; c < numBlocks; ++c) {
1330  RCP<const XpMap> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1331  RCP<const XpMap> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1332  RCP<const XpMap> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1333  RCP<const XpMap> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", lib, comm);
1334  RCP<XpMat> mat = XpIO::Read(fileName + XpIO::toString<size_t>(r) + XpIO::toString<size_t>(c) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1335  //RCP<XpCrsMatWrap> cmat = Teuchos::rcp_dynamic_cast<XpCrsMatWrap>(mat);
1336  bOp->setMatrix(r, c, mat);
1337  }
1338  }
1339 
1340  bOp->fillComplete();
1341 
1342  return bOp;
1343  } //ReadBlockedCrsMatrix
1344 
1346  template<class T>
1347  static std::string toString(const T& what) {
1348  std::ostringstream buf;
1349  buf << what;
1350  return buf.str();
1351  }
1352  };
1353 #endif // HAVE_XPETRA_EPETRA
1354 
1355 
1356 } // end namespace Xpetra
1357 
1358 #define XPETRA_IO_SHORT
1359 
1360 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
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:482
RCP< CrsMatrix > getCrsMatrix() const
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &)
Definition: Xpetra_IO.hpp:107
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Read(const std::string &filename, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > rowMap, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > colMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > domainMap=Teuchos::null, const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > rangeMap=Teuchos::null, const bool callFillComplete=true, const bool binary=false, const bool tolerant=false, const bool debug=false)
Read matrix from file in Matrix Market or binary format.
Definition: Xpetra_IO.hpp:1061
GlobalOrdinal GO
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:223
size_type size() const
static void 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:852
static void WriteLocal(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save local parts of matrix to files in Matrix Market format.
Definition: Xpetra_IO.hpp:893
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:1265
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:773
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
Definition: Xpetra_IO.hpp:211
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
Exception indicating invalid cast attempted.
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:920
void resize(const size_type n, const T &val=T())
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:624
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Read(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm, bool binary=false)
Read matrix from file in Matrix Market or binary format.
Definition: Xpetra_IO.hpp:958
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
Definition: Xpetra_IO.hpp:597
virtual size_t Cols() const
number of column blocks
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > Convert_Epetra_MultiVector_ToXpetra_MultiVector(RCP< Epetra_MultiVector > &epX)
Definition: Xpetra_IO.hpp:131
void resize(size_type new_size, const value_type &x=value_type())
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:1347
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:1181
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
void writeFile(const std::string &filename, const ArrayView< const OrdinalType > &rowptr, const ArrayView< const OrdinalType > &colind, const ArrayView< const ScalarType > &values, const OrdinalType numRows, const OrdinalType numCols)
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
static void WriteLocal(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Save local parts of matrix to files in Matrix Market format.
Definition: Xpetra_IO.hpp:317
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:383
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
Definition: Xpetra_IO.hpp:182
Concrete implementation of Xpetra::Matrix.
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
virtual size_t Rows() const
number of row blocks
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
static RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
Definition: Xpetra_IO.hpp:1223
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
static 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:761
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:281
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
Definition: Xpetra_IO.hpp:159
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:785
static void Write(const std::string &fileName, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save vector to file in Matrix Market format.
Definition: Xpetra_IO.hpp:250
static void Write(const std::string &fileName, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Save vector to file in Matrix Market format.
Definition: Xpetra_IO.hpp:815
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra-specific matrix class.
std::string toString(const T &t)
static 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:345
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:649
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
bool is_null() const
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:733