Xpetra Version of the Day
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>
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_TpetraCrsGraph.hpp>
78#include <Xpetra_TpetraCrsMatrix.hpp>
79#include <Xpetra_TpetraBlockCrsMatrix.hpp>
80#endif
81
82#ifdef HAVE_XPETRA_EPETRA
83#include <Xpetra_EpetraMap.hpp>
84#endif
85
86#include "Xpetra_Matrix.hpp"
88#include "Xpetra_CrsGraph.hpp"
89#include "Xpetra_CrsMatrixWrap.hpp"
91
92#include "Xpetra_Map.hpp"
93#include "Xpetra_StridedMap.hpp"
94#include "Xpetra_StridedMapFactory.hpp"
95#include "Xpetra_MapExtractor.hpp"
97
98#include <Teuchos_MatrixMarket_Raw_Writer.hpp>
99#include <string>
100
101
102namespace Xpetra {
103
104
105#ifdef HAVE_XPETRA_EPETRA
106 // This non-member templated function exists so that the matrix-matrix multiply will compile if Epetra, Tpetra, and ML are enabled.
107 template<class SC,class LO,class GO,class NO>
108 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >
111 "Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
112 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
113 }
114
115 // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
116 template<>
117 inline RCP<Xpetra::CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_CrsMatrix> &epAB) {
118 typedef double SC;
119 typedef int LO;
120 typedef int GO;
121 typedef Xpetra::EpetraNode NO;
122
123 RCP<Xpetra::EpetraCrsMatrixT<GO,NO> > tmpC1 = rcp(new Xpetra::EpetraCrsMatrixT<GO,NO>(epAB));
124 RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > tmpC2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<SC,LO,GO,NO> >(tmpC1);
125 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > tmpC3 = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(tmpC2));
126
127 return tmpC3;
128 }
129
130
131 template<class SC,class LO,class GO,class NO>
132 RCP<Xpetra::MultiVector<SC,LO,GO,NO> >
135 "Convert_Epetra_MultiVector_ToXpetra_MultiVector cannot be used with Scalar != double, LocalOrdinal != int, GlobalOrdinal != int");
136 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
137 }
138
139 // specialization for the case of ScalarType=double and LocalOrdinal=GlobalOrdinal=int
140 template<>
141 inline RCP<Xpetra::MultiVector<double,int,int,Xpetra::EpetraNode> > Convert_Epetra_MultiVector_ToXpetra_MultiVector<double,int,int,Xpetra::EpetraNode> (RCP<Epetra_MultiVector> &epX) {
142 typedef double SC;
143 typedef int LO;
144 typedef int GO;
145 typedef Xpetra::EpetraNode NO;
146
147 RCP<Xpetra::MultiVector<SC,LO,GO,NO >> tmp = Xpetra::toXpetra<GO,NO>(epX);
148 return tmp;
149 }
150
151#endif
152
157 template <class Scalar,
158 class LocalOrdinal = int,
159 class GlobalOrdinal = LocalOrdinal,
161 class IO {
162
163 private:
164#undef XPETRA_IO_SHORT
166
167 public:
168
169#ifdef HAVE_XPETRA_EPETRA
171 // @{
172 /*static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec);
173 static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec);
174
175 static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec);
176 static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec);
177
178 static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op);
179 static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op);
180
181 static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op);
182 static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op);*/
183
185 RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
186 if (xeMap == Teuchos::null)
187 throw Exceptions::BadCast("Utils::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
188 return xeMap->getEpetra_Map();
189 }
190 // @}
191#endif
192
193#ifdef HAVE_XPETRA_TPETRA
195 // @{
196 /*static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec);
197 static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec);
198 static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec);
199
200 static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec);
201 static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec);
202
203 static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op);
204 static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op);
205
206 static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op);
207 static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op);
208
209 static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op);
210 static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op);*/
211
212
214 const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
215 if (tmp_TMap == Teuchos::null)
216 throw Exceptions::BadCast("Utils::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
217 return tmp_TMap->getTpetra_Map();
218 }
219#endif
220
221
223
224
225 static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
227#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
228 const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
229 if (tmp_EMap != Teuchos::null) {
230 int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
231 if (rv != 0)
232 throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
233 return;
234 }
235#endif // HAVE_XPETRA_EPETRAEXT
236
237#ifdef HAVE_XPETRA_TPETRA
239 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
240 if (tmp_TMap != Teuchos::null) {
241 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
242 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
243 return;
244 }
245#endif // HAVE_XPETRA_TPETRA
246
247 throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
248
249 } //Write
250
252 static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
253 std::string mapfile = "map_" + fileName;
254 Write(mapfile, *(vec.getMap()));
255
257#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
258 const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
259 if (tmp_EVec != Teuchos::null) {
260 int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
261 if (rv != 0)
262 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
263 return;
264 }
265#endif // HAVE_XPETRA_EPETRA
266
267#ifdef HAVE_XPETRA_TPETRA
269 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
270 if (tmp_TVec != Teuchos::null) {
271 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
272 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
273 return;
274 }
275#endif // HAVE_XPETRA_TPETRA
276
277 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
278
279 } //Write
280
281
283 static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
284
285 Write("rowmap_" + fileName, *(Op.getRowMap()));
286 if ( !Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
287 Write("domainmap_" + fileName, *(Op.getDomainMap()));
288 if ( !Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
289 Write("rangemap_" + fileName, *(Op.getRangeMap()));
290 if ( !Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps )
291 Write("colmap_" + fileName, *(Op.getColMap()));
292
296#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
297 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
298 if (tmp_ECrsMtx != Teuchos::null) {
299 RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
300 int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
301 if (rv != 0)
302 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
303 return;
304 }
305#endif
306
307#ifdef HAVE_XPETRA_TPETRA
309 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
310 if (tmp_TCrsMtx != Teuchos::null) {
311 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
312 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
313 return;
314 }
315#endif // HAVE_XPETRA_TPETRA
316
317 throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
318 } //Write
319
320
322 static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
326
327 ArrayRCP<const size_t> rowptr_RCP;
328 ArrayRCP<LocalOrdinal> rowptr2_RCP;
330 ArrayRCP<const Scalar> vals_RCP;
331 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
332
333 ArrayView<const size_t> rowptr = rowptr_RCP();
334 ArrayView<const LocalOrdinal> colind = colind_RCP();
335 ArrayView<const Scalar> vals = vals_RCP();
336
337 rowptr2_RCP.resize(rowptr.size());
338 ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
339 for (size_t j = 0; j<rowptr.size(); j++)
340 rowptr2[j] = rowptr[j];
341
343 writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
344 rowptr2,colind,vals,
345 rowptr.size()-1,Op.getColMap()->getNodeNumElements());
346 } //WriteLocal
347
348
363 static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
365
366 // Write all matrix blocks with their maps
367 for (size_t row = 0; row < Op.Rows(); ++row) {
368 for (size_t col = 0; col < Op.Cols(); ++col) {
369 RCP<const Matrix > m = Op.getMatrix(row,col);
370 if(m != Teuchos::null) { // skip empty blocks
371 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
372 "Sub block matrix (" << row << "," << col << ") is not of type CrsMatrixWrap.");
373 XpIO::Write(fileName + toString(row) + toString(col) + ".m", *m, writeAllMaps);
374 }
375 }
376 }
377
378 // write map information of map extractors
379 RCP<const MapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
380 RCP<const MapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
381
382 for(size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
383 RCP<const Map> map = rangeMapExtractor->getMap(row);
384 XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", *map);
385 }
386 XpIO::Write("fullRangeMap_" + fileName +".m",*(rangeMapExtractor->getFullMap()));
387
388 for(size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
389 RCP<const Map> map = domainMapExtractor->getMap(col);
390 XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", *map);
391 }
392 XpIO::Write("fullDomainMap_" + fileName+ ".m",*(domainMapExtractor->getFullMap()));
393 } // WriteBlockedCrsMatrix
394
396 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) {
397 if (binary == false) {
398 // Matrix Market file format (ASCII)
399 if (lib == Xpetra::UseEpetra) {
400#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
402 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
403 int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
404 if (rv != 0)
405 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
406
407 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
408
410 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
411 return A;
412#else
413 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
414#endif
415 } else if (lib == Xpetra::UseTpetra) {
416#ifdef HAVE_XPETRA_TPETRA
417 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
418
419 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
420
421 bool callFillComplete = true;
422
423 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
424
425 if (tA.is_null())
426 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
427
429 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
431
432 return A;
433#else
434 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
435#endif
436 } else {
437 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
438 }
439 } else {
440 // Custom file format (binary)
441 std::ifstream ifs(fileName.c_str(), std::ios::binary);
442 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
443 int m, n, nnz;
444 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
445 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
446 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
447
448 int myRank = comm->getRank();
449
450 GO indexBase = 0;
451 RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
452 RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
454
455 if (myRank == 0) {
458 // Scan matrix to determine the exact nnz per row.
459 Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
460 for (int i = 0; i < m; i++) {
461 int row, rownnz;
462 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
463 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
464 numEntriesPerRow[i] = rownnz;
465 for (int j = 0; j < rownnz; j++) {
466 int index;
467 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
468 }
469 for (int j = 0; j < rownnz; j++) {
470 double value;
471 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
472 }
473 }
474
476
477 // Now that nnz per row are known, reread and store the matrix.
478 ifs.seekg(0, ifs.beg); //rewind to beginning of file
479 int junk; //skip header info
480 ifs.read(reinterpret_cast<char*>(&m), sizeof(junk));
481 ifs.read(reinterpret_cast<char*>(&n), sizeof(junk));
482 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(junk));
483 for (int i = 0; i < m; i++) {
484 int row, rownnz;
485 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
486 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
487 inds.resize(rownnz);
488 vals.resize(rownnz);
489 for (int j = 0; j < rownnz; j++) {
490 int index;
491 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
492 inds[j] = Teuchos::as<GlobalOrdinal>(index);
493 }
494 for (int j = 0; j < rownnz; j++) {
495 double value;
496 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
497 vals[j] = Teuchos::as<Scalar>(value);
498 }
499 A->insertGlobalValues(row, inds, vals);
500 }
501 } //if (myRank == 0)
502
503 A->fillComplete(domainMap, rangeMap);
504
505 return A;
506 } //if (binary == false) ... else
507
508 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
509
510 } //Read()
511
512
518 Read(const std::string& filename,
520 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
521 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
522 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
523 const bool callFillComplete = true,
524 const bool binary = false,
525 const bool tolerant = false,
526 const bool debug = false) {
527 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
528
529 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
530 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
531
532 const Xpetra::UnderlyingLib lib = rowMap->lib();
533 if (binary == false) {
534 if (lib == Xpetra::UseEpetra) {
535#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
537 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
539 const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
540 const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
541 int rv;
542 if (colMap.is_null()) {
543 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
544
545 } else {
546 const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
547 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
548 }
549
550 if (rv != 0)
551 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
552
553 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
555 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
556
557 return A;
558#else
559 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
560#endif
561 } else if (lib == Xpetra::UseTpetra) {
562#ifdef HAVE_XPETRA_TPETRA
563 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
564 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
565 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
566
567 const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
568 RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
569 const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
570 const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
571
572 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
573 callFillComplete, tolerant, debug);
574 if (tA.is_null())
575 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
576
578 RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
580
581 return A;
582#else
583 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
584#endif
585 } else {
586 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
587 }
588 } else {
589 // Custom file format (binary)
590 std::ifstream ifs(filename.c_str(), std::ios::binary);
591 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
592 int m, n, nnz;
593 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
594 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
595 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
596
597 //2020-June-05 JHU : for Tpetra, this will probably fail because Tpetra now requires staticly-sized matrix graphs.
599
600 //2019-06-07 JHU I don't see why this should matter.
601 //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GO), Exceptions::RuntimeError, "Incompatible sizes");
602
603 Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
604 Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
605
608 for (int i = 0; i < m; i++) {
609 int row, rownnz;
610 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
611 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
612 inds.resize(rownnz);
613 vals.resize(rownnz);
614 for (int j = 0; j < rownnz; j++) {
615 int index;
616 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
617 inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
618 }
619 for (int j = 0; j < rownnz; j++) {
620 double value;
621 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
622 vals[j] = Teuchos::as<SC>(value);
623 }
624 //This implies that row is not a global index.
625 A->insertGlobalValues(rowElements[row], inds, vals);
626 }
627 A->fillComplete(domainMap, rangeMap);
628 return A;
629 }
630
631 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
632 }
634
635
636 static RCP<MultiVector> ReadMultiVector (const std::string& fileName, const RCP<const Map>& map) {
637 Xpetra::UnderlyingLib lib = map->lib();
638
639 if (lib == Xpetra::UseEpetra) {
640 TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
641
642 } else if (lib == Xpetra::UseTpetra) {
643#ifdef HAVE_XPETRA_TPETRA
644 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
645 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
646 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
647 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
648
649 RCP<const map_type> temp = toTpetra(map);
650 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
652 return rmv;
653#else
654 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
655#endif
656 } else {
657 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
658 }
659
660 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
661 }
662
663 static RCP<const Map> ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
664 if (lib == Xpetra::UseEpetra) {
665 TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
666 } else if (lib == Xpetra::UseTpetra) {
667#ifdef HAVE_XPETRA_TPETRA
668 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
669 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
670
671 RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
672 if (tMap.is_null())
673 throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
674
675 return Xpetra::toXpetra(tMap);
676#else
677 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
678#endif
679 } else {
680 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
681 }
682
683 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
684
685 }
686
702
703 size_t numBlocks = 2; // TODO user parameter?
704
705 std::vector<RCP<const Map> > rangeMapVec;
706 for(size_t row = 0; row < numBlocks; ++row) {
707 RCP<const Map> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", lib, comm);
708 rangeMapVec.push_back(map);
709 }
710 RCP<const Map> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
711
712 std::vector<RCP<const Map> > domainMapVec;
713 for(size_t col = 0; col < numBlocks; ++col) {
714 RCP<const Map> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", lib, comm);
715 domainMapVec.push_back(map);
716 }
717 RCP<const Map> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
718
719 /*std::vector<RCP<const XpMap> > testRgMapVec;
720 for(size_t r = 0; r < numBlocks; ++r) {
721 RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
722 testRgMapVec.push_back(map);
723 }
724 std::vector<RCP<const XpMap> > testDoMapVec;
725 for(size_t c = 0; c < numBlocks; ++c) {
726 RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
727 testDoMapVec.push_back(map);
728 }*/
729
730 // create map extractors
731
732 // range map extractor
733 bool bRangeUseThyraStyleNumbering = false;
734 /*GlobalOrdinal gMinGids = 0;
735 for(size_t v = 0; v < testRgMapVec.size(); ++v) {
736 gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
737 }
738 if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;
739 */
740 RCP<const MapExtractor> rangeMapExtractor =
741 Teuchos::rcp(new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
742
743
744 // domain map extractor
745 bool bDomainUseThyraStyleNumbering = false;
746 /*gMinGids = 0;
747 for(size_t v = 0; v < testDoMapVec.size(); ++v) {
748 gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
749 }
750 if ( gMinGids==0 && testDoMapVec.size() > 1 ) bDomainUseThyraStyleNumbering = true;
751 */
752 RCP<const MapExtractor> domainMapExtractor =
753 Teuchos::rcp(new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
754
755 RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor,33));
756
757 // Read all sub-matrices with their maps and place into blocked operator
758 for (size_t row = 0; row < numBlocks; ++row) {
759 for (size_t col = 0; col < numBlocks; ++col) {
760 RCP<const Map> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
761 RCP<const Map> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
762 RCP<const Map> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
763 RCP<const Map> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
764 RCP<Matrix> mat = XpIO::Read(fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
765 bOp->setMatrix(row, col, mat);
766 }
767 }
768
769 bOp->fillComplete();
770
771 return bOp;
772 } // ReadBlockedCrsMatrix
773
774
776 template<class T>
777 static std::string toString(const T& what) {
778 std::ostringstream buf;
779 buf << what;
780 return buf.str();
781 }
782 };
783
784
785#ifdef HAVE_XPETRA_EPETRA
795 template <class Scalar>
796 class IO<Scalar,int,int,EpetraNode> {
797 public:
798 typedef int LocalOrdinal;
799 typedef int GlobalOrdinal;
801
802#ifdef HAVE_XPETRA_EPETRA
804 // @{
806 RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(Teuchos::rcpFromRef(map));
807 if (xeMap == Teuchos::null)
808 throw Exceptions::BadCast("IO::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
809 return xeMap->getEpetra_Map();
810 }
811 // @}
812#endif
813
814#ifdef HAVE_XPETRA_TPETRA
816 // @{
818 const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >& tmp_TMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
819 if (tmp_TMap == Teuchos::null)
820 throw Exceptions::BadCast("IO::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
821 return tmp_TMap->getTpetra_Map();
822 }
823#endif
824
825
827
828
829 static void Write(const std::string& fileName, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> & M) {
831#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
832 const RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >& tmp_EMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(tmp_Map);
833 if (tmp_EMap != Teuchos::null) {
834 int rv = EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), tmp_EMap->getEpetra_Map());
835 if (rv != 0)
836 throw Exceptions::RuntimeError("EpetraExt::BlockMapToMatrixMarketFile() return value of " + Teuchos::toString(rv));
837 return;
838 }
839#endif // HAVE_XPETRA_EPETRA
840
841#ifdef HAVE_XPETRA_TPETRA
842# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
843 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
844 // do nothing
845# else
847 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Map);
848 if (tmp_TMap != Teuchos::null) {
849 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > TMap = tmp_TMap->getTpetra_Map();
850 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeMapFile(fileName, *TMap);
851 return;
852 }
853# endif
854#endif // HAVE_XPETRA_TPETRA
855 throw Exceptions::BadCast("Could not cast to EpetraMap or TpetraMap in map writing");
856 }
857
859 static void Write(const std::string& fileName, const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & vec) {
860 std::string mapfile = "map_" + fileName;
861 Write(mapfile, *(vec.getMap()));
862
864#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
865 const RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >& tmp_EVec = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(tmp_Vec);
866 if (tmp_EVec != Teuchos::null) {
867 int rv = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(), *(tmp_EVec->getEpetra_MultiVector()));
868 if (rv != 0)
869 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
870 return;
871 }
872#endif // HAVE_XPETRA_EPETRAEXT
873
874#ifdef HAVE_XPETRA_TPETRA
875# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
876 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
877 // do nothin
878# else
880 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_Vec);
881 if (tmp_TVec != Teuchos::null) {
882 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > TVec = tmp_TVec->getTpetra_MultiVector();
883 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeDenseFile(fileName, TVec);
884 return;
885 }
886# endif
887#endif // HAVE_XPETRA_TPETRA
888
889 throw Exceptions::BadCast("Could not cast to EpetraMultiVector or TpetraMultiVector in multivector writing");
890
891 } //Write
892
893
894
896 static void Write(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
897
898 Write("rowmap_" + fileName, *(Op.getRowMap()));
899 if ( !Op.getDomainMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
900 Write("domainmap_" + fileName, *(Op.getDomainMap()));
901 if ( !Op.getRangeMap()->isSameAs(*(Op.getRowMap())) || writeAllMaps )
902 Write("rangemap_" + fileName, *(Op.getRangeMap()));
903 if ( !Op.getColMap()->isSameAs(*(Op.getDomainMap())) || writeAllMaps )
904 Write("colmap_" + fileName, *(Op.getColMap()));
905
909#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
910 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
911 if (tmp_ECrsMtx != Teuchos::null) {
912 RCP<const Epetra_CrsMatrix> A = tmp_ECrsMtx->getEpetra_CrsMatrix();
913 int rv = EpetraExt::RowMatrixToMatrixMarketFile(fileName.c_str(), *A);
914 if (rv != 0)
915 throw Exceptions::RuntimeError("EpetraExt::RowMatrixToMatrixMarketFile return value of " + Teuchos::toString(rv));
916 return;
917 }
918#endif // endif HAVE_XPETRA_EPETRA
919
920#ifdef HAVE_XPETRA_TPETRA
921# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
922 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
923 // do nothin
924# else
926 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_CrsMtx);
927 if (tmp_TCrsMtx != Teuchos::null) {
928 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A = tmp_TCrsMtx->getTpetra_CrsMatrix();
929 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseFile(fileName, A);
930 return;
931 }
932# endif
933#endif // HAVE_XPETRA_TPETRA
934
935 throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
936 } //Write
937
938
940 static void Write(const std::string& fileName, const Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> & graph, const bool &writeAllMaps = false) {
941
942 Write("rowmap_" + fileName, *(graph.getRowMap()));
943 if ( !graph.getDomainMap()->isSameAs(*(graph.getRowMap())) || writeAllMaps )
944 Write("domainmap_" + fileName, *(graph.getDomainMap()));
945 if ( !graph.getRangeMap()->isSameAs(*(graph.getRowMap())) || writeAllMaps )
946 Write("rangemap_" + fileName, *(graph.getRangeMap()));
947 if ( !graph.getColMap()->isSameAs(*(graph.getDomainMap())) || writeAllMaps )
948 Write("colmap_" + fileName, *(graph.getColMap()));
949
951
952#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
953 const RCP<const Xpetra::EpetraCrsGraphT<GlobalOrdinal,Node> >& tmp_ECrsGraph = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsGraphT<GlobalOrdinal,Node> >(tmp_Graph);
954 if (tmp_ECrsGraph != Teuchos::null) {
955 throw Exceptions::BadCast("Writing not implemented for EpetraCrsGraphT");
956 }
957#endif // endif HAVE_XPETRA_EPETRA
958
959#ifdef HAVE_XPETRA_TPETRA
960# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
961 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
962 // do nothin
963# else
965 Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsGraph<LocalOrdinal, GlobalOrdinal, Node> >(tmp_Graph);
966 if (tmp_TCrsGraph != Teuchos::null) {
967 RCP<const Tpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > G = tmp_TCrsGraph->getTpetra_CrsGraph();
968 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >::writeSparseGraphFile(fileName, G);
969 return;
970 }
971# endif
972#endif // HAVE_XPETRA_TPETRA
973
974 throw Exceptions::BadCast("Could not cast to EpetraCrsMatrix or TpetraCrsMatrix in matrix writing");
975 } //Write
976
977
979 static void WriteLocal(const std::string& fileName, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
983
984 ArrayRCP<const size_t> rowptr_RCP;
985 ArrayRCP<LocalOrdinal> rowptr2_RCP;
987 ArrayRCP<const Scalar> vals_RCP;
988 tmp_CrsMtx->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
989
990 ArrayView<const size_t> rowptr = rowptr_RCP();
991 ArrayView<const LocalOrdinal> colind = colind_RCP();
992 ArrayView<const Scalar> vals = vals_RCP();
993
994 rowptr2_RCP.resize(rowptr.size());
995 ArrayView<LocalOrdinal> rowptr2 = rowptr2_RCP();
996 for (size_t j = 0; j<Teuchos::as<size_t>(rowptr.size()); j++)
997 rowptr2[j] = rowptr[j];
998
1000 writer.writeFile(fileName + "." + std::to_string(Op.getRowMap()->getComm()->getSize()) + "." + std::to_string(Op.getRowMap()->getComm()->getRank()),
1001 rowptr2,colind,vals,
1002 rowptr.size()-1,Op.getColMap()->getNodeNumElements());
1003 } //WriteLocal
1004
1021 static void WriteBlockedCrsMatrix(const std::string& fileName, const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op, const bool &writeAllMaps = false) {
1022#include "Xpetra_UseShortNames.hpp"
1024
1025 // write all matrices with their maps
1026 for (size_t row = 0; row < Op.Rows(); ++row) {
1027 for (size_t col = 0; col < Op.Cols(); ++col) {
1028 RCP<const Matrix> m = Op.getMatrix(row, col);
1029 if(m != Teuchos::null) { // skip empty blocks
1030 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(m) == Teuchos::null, Xpetra::Exceptions::BadCast,
1031 "Sub block matrix (" << row << "," << col << ") is not of type CrsMatrixWrap.");
1032 XpIO::Write(fileName + toString(row) + toString(col) + ".m", *m, writeAllMaps);
1033 }
1034 }
1035 }
1036
1037 // write map information of map extractors
1038 RCP<const MapExtractor> rangeMapExtractor = Op.getRangeMapExtractor();
1039 RCP<const MapExtractor> domainMapExtractor = Op.getDomainMapExtractor();
1040
1041 for(size_t row = 0; row < rangeMapExtractor->NumMaps(); ++row) {
1042 RCP<const Map> map = rangeMapExtractor->getMap(row);
1043 XpIO::Write("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", *map);
1044 }
1045 XpIO::Write("fullRangeMap_" + fileName +".m", *(rangeMapExtractor->getFullMap()));
1046
1047 for(size_t col = 0; col < domainMapExtractor->NumMaps(); ++col) {
1048 RCP<const Map> map = domainMapExtractor->getMap(col);
1049 XpIO::Write("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", *map);
1050 }
1051 XpIO::Write("fullDomainMap_" + fileName+ ".m", *(domainMapExtractor->getFullMap()));
1052 } // WriteBlockedCrsMatrix
1053
1055 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) {
1056 if (binary == false) {
1057 // Matrix Market file format (ASCII)
1058 if (lib == Xpetra::UseEpetra) {
1059#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1060 Epetra_CrsMatrix *eA;
1061 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(comm);
1062 int rv = EpetraExt::MatrixMarketFileToCrsMatrix(fileName.c_str(), *epcomm, eA);
1063 if (rv != 0)
1064 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1065
1066 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1067
1069 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1070 return A;
1071#else
1072 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1073#endif
1074 } else if (lib == Xpetra::UseTpetra) {
1075#ifdef HAVE_XPETRA_TPETRA
1076# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1077 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1078 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int enabled.");
1079# else
1080 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1081
1082 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1083
1084 bool callFillComplete = true;
1085
1086 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(fileName, comm, callFillComplete);
1087
1088 if (tA.is_null())
1089 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1090
1092 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmpA1);
1094
1095 return A;
1096# endif
1097#else
1098 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1099#endif
1100 } else {
1101 throw Exceptions::RuntimeError("Xpetra:IO: you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1102 }
1103 } else {
1104 // Custom file format (binary)
1105 std::ifstream ifs(fileName.c_str(), std::ios::binary);
1106 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << fileName << "\"");
1107 int m, n, nnz;
1108 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1109 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1110 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1111
1112 int myRank = comm->getRank();
1113
1114 GlobalOrdinal indexBase = 0;
1115 RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, m, (myRank == 0 ? m : 0), indexBase, comm), rangeMap = rowMap;
1116 RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, n, (myRank == 0 ? n : 0), indexBase, comm), domainMap = colMap;
1117
1119
1120 if (myRank == 0) {
1123 // Scan matrix to determine the exact nnz per row.
1124 Teuchos::ArrayRCP<size_t> numEntriesPerRow(m);
1125 for (int i = 0; i < m; i++) {
1126 int row, rownnz;
1127 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1128 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1129 numEntriesPerRow[i] = rownnz;
1130 for (int j = 0; j < rownnz; j++) {
1131 int index;
1132 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1133 }
1134 for (int j = 0; j < rownnz; j++) {
1135 double value;
1136 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1137 }
1138 }
1139
1141
1142 // Now that nnz per row are known, reread and store the matrix.
1143 ifs.seekg(0, ifs.beg); //rewind to beginning of file
1144 int junk; //skip header info
1145 ifs.read(reinterpret_cast<char*>(&m), sizeof(junk));
1146 ifs.read(reinterpret_cast<char*>(&n), sizeof(junk));
1147 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(junk));
1148 for (int i = 0; i < m; i++) {
1149 int row, rownnz;
1150 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1151 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1152 inds.resize(rownnz);
1153 vals.resize(rownnz);
1154 for (int j = 0; j < rownnz; j++) {
1155 int index;
1156 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1157 inds[j] = Teuchos::as<GlobalOrdinal>(index);
1158 }
1159 for (int j = 0; j < rownnz; j++) {
1160 double value;
1161 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1162 vals[j] = Teuchos::as<Scalar>(value);
1163 }
1164 A->insertGlobalValues(row, inds, vals);
1165 }
1166 } //if (myRank == 0)
1167
1168 A->fillComplete(domainMap, rangeMap);
1169
1170 return A;
1171 }
1172
1173 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1174
1175 } //Read()
1176
1177
1184 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > colMap = Teuchos::null,
1185 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::null,
1186 const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::null,
1187 const bool callFillComplete = true,
1188 const bool binary = false,
1189 const bool tolerant = false,
1190 const bool debug = false) {
1191 TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), Exceptions::RuntimeError, "Utils::Read() : rowMap cannot be null");
1192
1193 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domain = (domainMap.is_null() ? rowMap : domainMap);
1194 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > range = (rangeMap .is_null() ? rowMap : rangeMap);
1195
1196 const Xpetra::UnderlyingLib lib = rowMap->lib();
1197 if (binary == false) {
1198 if (lib == Xpetra::UseEpetra) {
1199#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1200 Epetra_CrsMatrix *eA;
1201 const RCP<const Epetra_Comm> epcomm = Xpetra::toEpetra(rowMap->getComm());
1203 const Epetra_Map& epetraDomainMap = (domainMap.is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*domainMap));
1204 const Epetra_Map& epetraRangeMap = (rangeMap .is_null() ? epetraRowMap : Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Map2EpetraMap(*rangeMap));
1205 int rv;
1206 if (colMap.is_null()) {
1207 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraRangeMap, epetraDomainMap, eA);
1208
1209 } else {
1210 const Epetra_Map& epetraColMap = Map2EpetraMap(*colMap);
1211 rv = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), epetraRowMap, epetraColMap, epetraRangeMap, epetraDomainMap, eA);
1212 }
1213
1214 if (rv != 0)
1215 throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToCrsMatrix return value of " + Teuchos::toString(rv));
1216
1217 RCP<Epetra_CrsMatrix> tmpA = rcp(eA);
1219 Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tmpA);
1220
1221 return A;
1222#else
1223 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1224#endif
1225 } else if (lib == Xpetra::UseTpetra) {
1226#ifdef HAVE_XPETRA_TPETRA
1227# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1228 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1229 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1230# else
1231 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1232 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1233 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1234
1235 const RCP<const map_type> tpetraRowMap = Map2TpetraMap(*rowMap);
1236 RCP<const map_type> tpetraColMap = (colMap.is_null() ? Teuchos::null : Map2TpetraMap(*colMap));
1237 const RCP<const map_type> tpetraRangeMap = (rangeMap.is_null() ? tpetraRowMap : Map2TpetraMap(*rangeMap));
1238 const RCP<const map_type> tpetraDomainMap = (domainMap.is_null() ? tpetraRowMap : Map2TpetraMap(*domainMap));
1239
1240 RCP<sparse_matrix_type> tA = reader_type::readSparseFile(filename, tpetraRowMap, tpetraColMap, tpetraDomainMap, tpetraRangeMap,
1241 callFillComplete, tolerant, debug);
1242 if (tA.is_null())
1243 throw Exceptions::RuntimeError("The Tpetra::CrsMatrix returned from readSparseFile() is null.");
1244
1246 RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpA2 = Teuchos::rcp_implicit_cast<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tmpA1);
1248
1249 return A;
1250# endif
1251#else
1252 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1253#endif
1254 } else {
1255 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1256 }
1257 } else {
1258 // Custom file format (binary)
1259 std::ifstream ifs(filename.c_str(), std::ios::binary);
1260 TEUCHOS_TEST_FOR_EXCEPTION(!ifs.good(), Exceptions::RuntimeError, "Can not read \"" << filename << "\"");
1261 int m, n, nnz;
1262 ifs.read(reinterpret_cast<char*>(&m), sizeof(m));
1263 ifs.read(reinterpret_cast<char*>(&n), sizeof(n));
1264 ifs.read(reinterpret_cast<char*>(&nnz), sizeof(nnz));
1265
1266 //2020-June-05 JHU : for Tpetra, this will probably fail because Tpetra now requires staticly-sized matrix graphs.
1268
1269 //2019-06-07 JHU I don't see why this should matter.
1270 //TEUCHOS_TEST_FOR_EXCEPTION(sizeof(int) != sizeof(GlobalOrdinal), Exceptions::RuntimeError, "Incompatible sizes");
1271
1272 Teuchos::ArrayView<const GlobalOrdinal> rowElements = rowMap->getNodeElementList();
1273 Teuchos::ArrayView<const GlobalOrdinal> colElements = colMap->getNodeElementList();
1274
1277 for (int i = 0; i < m; i++) {
1278 int row, rownnz;
1279 ifs.read(reinterpret_cast<char*>(&row), sizeof(row));
1280 ifs.read(reinterpret_cast<char*>(&rownnz), sizeof(rownnz));
1281 inds.resize(rownnz);
1282 vals.resize(rownnz);
1283 for (int j = 0; j < rownnz; j++) {
1284 int index;
1285 ifs.read(reinterpret_cast<char*>(&index), sizeof(index));
1286 inds[j] = colElements[Teuchos::as<LocalOrdinal>(index)];
1287 }
1288 for (int j = 0; j < rownnz; j++) {
1289 double value;
1290 ifs.read(reinterpret_cast<char*>(&value), sizeof(value));
1291 vals[j] = Teuchos::as<Scalar>(value);
1292 }
1293 //This implies that row is not a global index.
1294 A->insertGlobalValues(rowElements[row], inds, vals);
1295 }
1296 A->fillComplete(domainMap, rangeMap);
1297 return A;
1298 }
1299
1300 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1301 }
1303
1304
1306 Xpetra::UnderlyingLib lib = map->lib();
1307
1308 if (lib == Xpetra::UseEpetra) {
1309 // taw: Oct 9 2015: do we need a specialization for <double,int,int>??
1310 //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1311#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1312 Epetra_MultiVector * MV;
1313 int rv = EpetraExt::MatrixMarketFileToMultiVector(fileName.c_str(), toEpetra(map), MV);
1314 if(rv != 0) throw Exceptions::RuntimeError("EpetraExt::MatrixMarketFileToMultiVector failed");
1315 RCP<Epetra_MultiVector> MVrcp = rcp(MV);
1316 return Convert_Epetra_MultiVector_ToXpetra_MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(MVrcp);
1317#else
1318 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1319#endif
1320 } else if (lib == Xpetra::UseTpetra) {
1321#ifdef HAVE_XPETRA_TPETRA
1322# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1323 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1324 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1325# else
1326 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1327 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1328 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
1329 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type;
1330
1331 RCP<const map_type> temp = toTpetra(map);
1332 RCP<multivector_type> TMV = reader_type::readDenseFile(fileName,map->getComm(),temp);
1334 return rmv;
1335# endif
1336#else
1337 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1338#endif
1339 } else {
1340 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1341 }
1342
1343 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1344
1345 }
1346
1347
1348 static RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ReadMap (const std::string& fileName, Xpetra::UnderlyingLib lib, const RCP<const Teuchos::Comm<int> >& comm) {
1349 if (lib == Xpetra::UseEpetra) {
1350 // do we need another specialization for <double,int,int> ??
1351 //TEUCHOS_TEST_FOR_EXCEPTION(true, ::Xpetra::Exceptions::BadCast, "Epetra can only be used with Scalar=double and Ordinal=int");
1352#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1353 Epetra_Map *eMap;
1354 int rv = EpetraExt::MatrixMarketFileToMap(fileName.c_str(), *(Xpetra::toEpetra(comm)), eMap);
1355 if (rv != 0)
1356 throw Exceptions::RuntimeError("Error reading map from file " + fileName + " with EpetraExt::MatrixMarketToMap (returned " + Teuchos::toString(rv) + ")");
1357
1358 RCP<Epetra_Map> eMap1 = rcp(new Epetra_Map(*eMap));
1359 return Xpetra::toXpetra<int,Node>(*eMap1);
1360#else
1361 throw Exceptions::RuntimeError("Xpetra has not been compiled with Epetra and EpetraExt support.");
1362#endif
1363 } else if (lib == Xpetra::UseTpetra) {
1364#ifdef HAVE_XPETRA_TPETRA
1365# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1366 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1367 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra GO=int support.");
1368# else
1369 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> sparse_matrix_type;
1370 typedef Tpetra::MatrixMarket::Reader<sparse_matrix_type> reader_type;
1371
1372 RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > tMap = reader_type::readMapFile(fileName, comm);
1373 if (tMap.is_null())
1374 throw Exceptions::RuntimeError("The Tpetra::Map returned from readSparseFile() is null.");
1375
1376 return Xpetra::toXpetra(tMap);
1377# endif
1378#else
1379 throw Exceptions::RuntimeError("Xpetra has not been compiled with Tpetra support.");
1380#endif
1381 } else {
1382 throw Exceptions::RuntimeError("Utils::Read : you must specify Xpetra::UseEpetra or Xpetra::UseTpetra.");
1383 }
1384
1385 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1386
1387 }
1388
1405#include "Xpetra_UseShortNames.hpp"
1407
1408 size_t numBlocks = 2; // TODO user parameter?
1409
1410 std::vector<RCP<const Map> > rangeMapVec;
1411 for(size_t row = 0; row < numBlocks; ++row) {
1412 RCP<const Map> map = XpIO::ReadMap("subRangeMap_" + fileName + XpIO::toString<size_t>(row) + ".m", lib, comm);
1413 rangeMapVec.push_back(map);
1414 }
1415 RCP<const Map> fullRangeMap = XpIO::ReadMap("fullRangeMap_" + fileName + ".m", lib, comm);
1416
1417 std::vector<RCP<const Map> > domainMapVec;
1418 for(size_t col = 0; col < numBlocks; ++col) {
1419 RCP<const Map> map = XpIO::ReadMap("subDomainMap_" + fileName + XpIO::toString<size_t>(col) + ".m", lib, comm);
1420 domainMapVec.push_back(map);
1421 }
1422 RCP<const Map> fullDomainMap = XpIO::ReadMap("fullDomainMap_" + fileName + ".m", lib, comm);
1423
1424 /*std::vector<RCP<const XpMap> > testRgMapVec;
1425 for(size_t r = 0; r < numBlocks; ++r) {
1426 RCP<const XpMap> map = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(r) + "0.m", lib, comm);
1427 testRgMapVec.push_back(map);
1428 }
1429 std::vector<RCP<const XpMap> > testDoMapVec;
1430 for(size_t c = 0; c < numBlocks; ++c) {
1431 RCP<const XpMap> map = XpIO::ReadMap("domainmap_" + fileName + "0" + XpIO::toString<size_t>(c) + ".m", lib, comm);
1432 testDoMapVec.push_back(map);
1433 }*/
1434
1435 // create map extractors
1436
1437 // range map extractor
1438 bool bRangeUseThyraStyleNumbering = false;
1439 /*
1440 GlobalOrdinal gMinGids = 0;
1441 for(size_t v = 0; v < testRgMapVec.size(); ++v) {
1442 gMinGids += testRgMapVec[v]->getMinAllGlobalIndex();
1443 }
1444 if ( gMinGids==0 && testRgMapVec.size() > 1 ) bRangeUseThyraStyleNumbering = true;*/
1445 RCP<const MapExtractor> rangeMapExtractor =
1446 rcp(new MapExtractor(fullRangeMap, rangeMapVec, bRangeUseThyraStyleNumbering));
1447
1448 // domain map extractor
1449 bool bDomainUseThyraStyleNumbering = false;
1450 /*gMinGids = 0;
1451 for(size_t v = 0; v < testDoMapVec.size(); ++v) {
1452 gMinGids += testDoMapVec[v]->getMinAllGlobalIndex();
1453 }
1454 if ( gMinGids==0 && testDoMapVec.size() > 1) bDomainUseThyraStyleNumbering = true;*/
1455 RCP<const MapExtractor> domainMapExtractor =
1456 rcp(new MapExtractor(fullDomainMap, domainMapVec, bDomainUseThyraStyleNumbering));
1457
1458 RCP<BlockedCrsMatrix> bOp = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 33));
1459
1460 // Read all matrices with their maps and create the BlockedCrsMatrix
1461 for (size_t row = 0; row < numBlocks; ++row) {
1462 for (size_t col = 0; col < numBlocks; ++col) {
1463 RCP<const Map> rowSubMap = XpIO::ReadMap("rowmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1464 RCP<const Map> colSubMap = XpIO::ReadMap("colmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1465 RCP<const Map> domSubMap = XpIO::ReadMap("domainmap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1466 RCP<const Map> ranSubMap = XpIO::ReadMap("rangemap_" + fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", lib, comm);
1467 RCP<Matrix> mat = XpIO::Read(fileName + XpIO::toString<size_t>(row) + XpIO::toString<size_t>(col) + ".m", rowSubMap, colSubMap, domSubMap, ranSubMap);
1468 bOp->setMatrix(row, col, mat);
1469 }
1470 }
1471
1472 bOp->fillComplete();
1473
1474 return bOp;
1475 } // ReadBlockedCrsMatrix
1476
1478 template<class T>
1479 static std::string toString(const T& what) {
1480 std::ostringstream buf;
1481 buf << what;
1482 return buf.str();
1483 }
1484 };
1485#endif // HAVE_XPETRA_EPETRA
1486
1487
1488} // end namespace Xpetra
1489
1490#define XPETRA_IO_SHORT
1491
1492#endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_IO_HPP_ */
LocalOrdinal LO
GlobalOrdinal GO
void resize(const size_type n, const T &val=T())
size_type size() const
void resize(size_type new_size, const value_type &x=value_type())
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)
bool is_null() const
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual size_t Rows() const
number of row blocks
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
virtual size_t Cols() const
number of column blocks
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
Returns the Map associated with the domain of this graph.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
Returns the Map that describes the row distribution in this graph.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
Returns the Map associated with the domain of this graph.
virtual RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const =0
Returns the Map that describes the column distribution in this graph.
Concrete implementation of Xpetra::Matrix.
RCP< CrsMatrix > getCrsMatrix() const
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
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:1182
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:1055
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save block matrix to one file per block in Matrix Market format.
Definition: Xpetra_IO.hpp:1021
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:805
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:1305
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:979
static void Write(const std::string &fileName, const Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, const bool &writeAllMaps=false)
Save CrsGraph to file in Matrix Market format.
Definition: Xpetra_IO.hpp:940
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:1348
static void Write(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:896
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 block matrix from one file per block in Matrix Market format.
Definition: Xpetra_IO.hpp:1404
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:1479
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:817
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:829
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:859
Xpetra utility class containing IO routines to read/write vectors, matrices etc...
Definition: Xpetra_IO.hpp:161
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:213
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:396
static std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: Xpetra_IO.hpp:777
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:252
static RCP< const Map > ReadMap(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm)
Definition: Xpetra_IO.hpp:663
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
Definition: Xpetra_IO.hpp:636
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Read/Write methods.
Definition: Xpetra_IO.hpp:225
static void Write(const std::string &fileName, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save matrix to file in Matrix Market format.
Definition: Xpetra_IO.hpp:283
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:322
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:518
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:184
static void WriteBlockedCrsMatrix(const std::string &fileName, const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const bool &writeAllMaps=false)
Save block matrix to one file per block in Matrix Market format.
Definition: Xpetra_IO.hpp:363
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 block matrix from one file per block in Matrix Market format.
Definition: Xpetra_IO.hpp:700
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Xpetra-specific matrix class.
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this 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 Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X....
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
std::string toString(const T &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > Convert_Epetra_MultiVector_ToXpetra_MultiVector(RCP< Epetra_MultiVector > &epX)
Definition: Xpetra_IO.hpp:133
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &)
Definition: Xpetra_IO.hpp:109