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