46#ifndef MUELU_UTILITIES_DEF_HPP
47#define MUELU_UTILITIES_DEF_HPP
49#include <Teuchos_DefaultComm.hpp>
50#include <Teuchos_ParameterList.hpp>
54#ifdef HAVE_MUELU_EPETRA
56# include "Epetra_MpiComm.h"
60#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
67#include <Xpetra_EpetraUtils.hpp>
68#include <Xpetra_EpetraMultiVector.hpp>
72#ifdef HAVE_MUELU_TPETRA
73#include <MatrixMarket_Tpetra.hpp>
74#include <Tpetra_RowMatrixTransposer.hpp>
75#include <TpetraExt_MatrixMatrix.hpp>
76#include <Xpetra_TpetraMultiVector.hpp>
77#include <Xpetra_TpetraOperator.hpp>
78#include <Xpetra_TpetraCrsMatrix.hpp>
79#include <Xpetra_TpetraBlockCrsMatrix.hpp>
82#ifdef HAVE_MUELU_EPETRA
83#include <Xpetra_EpetraMap.hpp>
86#include <Xpetra_BlockedCrsMatrix.hpp>
88#include <Xpetra_IO.hpp>
89#include <Xpetra_Import.hpp>
90#include <Xpetra_ImportFactory.hpp>
91#include <Xpetra_Map.hpp>
92#include <Xpetra_MapFactory.hpp>
93#include <Xpetra_Matrix.hpp>
94#include <Xpetra_MatrixFactory.hpp>
95#include <Xpetra_MultiVector.hpp>
96#include <Xpetra_MultiVectorFactory.hpp>
97#include <Xpetra_Operator.hpp>
98#include <Xpetra_Vector.hpp>
99#include <Xpetra_VectorFactory.hpp>
101#include <Xpetra_MatrixMatrix.hpp>
104#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_ML)
105#include <ml_operator.h>
106#include <ml_epetra_utils.h>
111#ifdef HAVE_MUELU_EPETRA
116#ifdef HAVE_MUELU_EPETRA
117 template<
typename SC,
typename LO,
typename GO,
typename NO>
119 return Xpetra::Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO>(epAB);
123#ifdef HAVE_MUELU_EPETRA
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
127 if (tmpVec == Teuchos::null)
128 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
129 return tmpVec->getEpetra_MultiVector();
132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
134 RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
135 if (tmpVec == Teuchos::null)
136 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
137 return tmpVec->getEpetra_MultiVector();
140 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
142 const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec =
dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &
>(vec);
143 return *(tmpVec.getEpetra_MultiVector());
146 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
148 const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec =
dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &
>(vec);
149 return *(tmpVec.getEpetra_MultiVector());
152 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
154 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
155 if (crsOp == Teuchos::null)
157 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
158 if (tmp_ECrsMtx == Teuchos::null)
159 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
160 return tmp_ECrsMtx->getEpetra_CrsMatrix();
163 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
165 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
166 if (crsOp == Teuchos::null)
168 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
169 if (tmp_ECrsMtx == Teuchos::null)
170 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
171 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
174 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
177 const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp =
dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(Op);
179 const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>& tmp_ECrsMtx =
dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>&
>(*crsOp.getCrsMatrix());
180 return *tmp_ECrsMtx.getEpetra_CrsMatrix();
181 }
catch (std::bad_cast&) {
182 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
184 }
catch (std::bad_cast&) {
189 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
192 Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp =
dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(Op);
194 Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>& tmp_ECrsMtx =
dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>&
>(*crsOp.getCrsMatrix());
195 return *tmp_ECrsMtx.getEpetra_CrsMatrixNonConst();
196 }
catch (std::bad_cast&) {
197 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
199 }
catch (std::bad_cast&) {
204 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
206 RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(rcpFromRef(map));
207 if (xeMap == Teuchos::null)
208 throw Exceptions::BadCast(
"Utilities::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
209 return xeMap->getEpetra_Map();
213#ifdef HAVE_MUELU_TPETRA
214 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
215 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
217 RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
218 if (tmpVec == Teuchos::null)
219 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
220 return tmpVec->getTpetra_MultiVector();
223 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
225 RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
226 if (tmpVec == Teuchos::null)
227 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
228 return tmpVec->getTpetra_MultiVector();
231 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
233 const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec =
dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(vec);
234 return *(tmpVec.getTpetra_MultiVector());
237 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
239 const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec =
dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(vec);
240 return tmpVec.getTpetra_MultiVector();
243 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
244 const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
246 const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec =
dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(vec);
247 return *(tmpVec.getTpetra_MultiVector());
250 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
253 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
254 if (crsOp == Teuchos::null)
256 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
257 if (tmp_ECrsMtx == Teuchos::null)
258 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
259 return tmp_ECrsMtx->getTpetra_CrsMatrix();
262 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
264 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
265 if (crsOp == Teuchos::null)
267 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
268 if (tmp_ECrsMtx == Teuchos::null)
269 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
270 return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
273 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
276 const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp =
dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(Op);
278 const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx =
dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*crsOp.getCrsMatrix());
279 return *tmp_ECrsMtx.getTpetra_CrsMatrix();
280 }
catch (std::bad_cast&) {
281 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
283 }
catch (std::bad_cast&) {
288 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
291 Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp =
dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(Op);
293 Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx =
dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*crsOp.getCrsMatrix());
294 return *tmp_ECrsMtx.getTpetra_CrsMatrixNonConst();
295 }
catch (std::bad_cast&) {
296 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
298 }
catch (std::bad_cast&) {
304 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
306 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
308 RCP<const XCrsMatrixWrap> crsOp = rcp_dynamic_cast<const XCrsMatrixWrap>(Op);
309 if (crsOp == Teuchos::null)
311 const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
312 if (tmp_ECrsMtx == Teuchos::null)
313 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
314 return tmp_ECrsMtx->getTpetra_BlockCrsMatrix();
317 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
319 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
320 RCP<const XCrsMatrixWrap> crsOp = rcp_dynamic_cast<const XCrsMatrixWrap>(Op);
321 if (crsOp == Teuchos::null)
323 const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
324 if (tmp_ECrsMtx == Teuchos::null)
325 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
326 return tmp_ECrsMtx->getTpetra_BlockCrsMatrixNonConst();
329 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
332 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
333 const XCrsMatrixWrap& crsOp =
dynamic_cast<const XCrsMatrixWrap&
>(Op);
335 const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx =
dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*crsOp.getCrsMatrix());
336 return *tmp_ECrsMtx.getTpetra_BlockCrsMatrix();
337 }
catch (std::bad_cast&) {
338 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
340 }
catch (std::bad_cast&) {
345 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
348 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
349 XCrsMatrixWrap& crsOp =
dynamic_cast<XCrsMatrixWrap&
>(Op);
351 Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx =
dynamic_cast<Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*crsOp.getCrsMatrix());
352 return *tmp_ECrsMtx.getTpetra_BlockCrsMatrixNonConst();
353 }
catch (std::bad_cast&) {
354 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
356 }
catch (std::bad_cast&) {
362 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
364 RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mat = rcp_dynamic_cast<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
365 if (!mat.is_null()) {
366 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mat);
367 if (crsOp == Teuchos::null)
370 RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
371 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
372 RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
373 if(!tmp_Crs.is_null()) {
374 return tmp_Crs->getTpetra_CrsMatrixNonConst();
377 tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
378 if (tmp_BlockCrs.is_null())
379 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
380 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
383 RCP<const Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpOp = rcp_dynamic_cast<const Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op,
true);
384 RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> tOp = tpOp->getOperatorConst();
385 RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tRow = rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tOp,
true);
390 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
392 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mat = rcp_dynamic_cast<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
393 if (!mat.is_null()) {
394 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mat);
395 if (crsOp == Teuchos::null)
398 RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
399 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
400 RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
401 if(!tmp_Crs.is_null()) {
402 return tmp_Crs->getTpetra_CrsMatrixNonConst();
405 tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
406 if (tmp_BlockCrs.is_null())
407 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
408 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
411 RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op,
true);
412 RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> tOp = tpOp->getOperator();
413 RCP<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tRow = rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tOp,
true);
419 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
421 const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node>>& tmp_TMap = rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
422 if (tmp_TMap == Teuchos::null)
423 throw Exceptions::BadCast(
"Utilities::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
424 return tmp_TMap->getTpetra_Map();
428 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
431 bool doOptimizeStorage)
433 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
434 Teuchos::ArrayRCP<Scalar> sv(scalingVector.size());
436 for (
int i = 0; i < scalingVector.size(); ++i)
437 sv[i] = one / scalingVector[i];
439 for (
int i = 0; i < scalingVector.size(); ++i)
440 sv[i] = scalingVector[i];
443 switch (Op.getRowMap()->lib()) {
444 case Xpetra::UseTpetra:
445 MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
448 case Xpetra::UseEpetra:
449 MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
457 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
462 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
465 bool doOptimizeStorage)
467#ifdef HAVE_MUELU_TPETRA
469 Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
471 const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = tpOp.getRowMap();
472 const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domainMap = tpOp.getDomainMap();
473 const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rangeMap = tpOp.getRangeMap();
475 size_t maxRowSize = tpOp.getLocalMaxNumRowEntries();
476 if (maxRowSize == Teuchos::as<size_t>(-1))
480 if (tpOp.isFillComplete())
483 if (Op.isLocallyIndexed() ==
true) {
484 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_inds_host_view_type cols;
485 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
487 for (
size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
488 tpOp.getLocalRowView(i, cols, vals);
490 size_t nnz = tpOp.getNumEntriesInLocalRow(i);
491 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::nonconst_values_host_view_type scaledVals(
"ScaledVals",nnz);
493 for (
size_t j = 0; j < nnz; ++j) {
494 scaledVals[j] = scalingVector[i]*vals[j];
498 tpOp.replaceLocalValues(i, cols, scaledVals);
503 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_inds_host_view_type cols;
504 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
506 for (
size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
508 tpOp.getGlobalRowView(gid, cols, vals);
509 size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
510 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::nonconst_values_host_view_type scaledVals(
"ScaledVals",nnz);
513 for (
size_t j = 0; j < nnz; ++j)
514 scaledVals[j] = scalingVector[i]*vals[j];
517 tpOp.replaceGlobalValues(gid, cols, scaledVals);
522 if (doFillComplete) {
523 if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
524 throw Exceptions::RuntimeError(
"In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
526 RCP<Teuchos::ParameterList> params = rcp(
new Teuchos::ParameterList());
527 params->set(
"Optimize Storage", doOptimizeStorage);
528 params->set(
"No Nonlocal Changes",
true);
529 Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
539 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
540 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
542 Transpose (Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op,
bool ,
const std::string & label,
const Teuchos::RCP<Teuchos::ParameterList> ¶ms) {
543#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
544 std::string TorE =
"epetra";
546 std::string TorE =
"tpetra";
549#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
558#ifdef HAVE_MUELU_TPETRA
559 if (TorE ==
"tpetra") {
560 using Helpers = Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
562 if(Helpers::isTpetraCrs(Op)) {
565 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
566 Tpetra::RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label);
569 using Teuchos::ParameterList;
571 RCP<ParameterList> transposeParams = params.is_null () ?
572 rcp (
new ParameterList) :
573 rcp (
new ParameterList (*params));
574 transposeParams->set (
"sort",
false);
575 A = transposer.createTranspose (transposeParams);
578 RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AA = rcp(
new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A) );
579 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAA = rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(AA);
580 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAAA = rcp(
new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(AAA) );
581 if (!AAAA->isFillComplete())
582 AAAA->fillComplete(Op.getRangeMap(), Op.getDomainMap());
584 if (Op.IsView(
"stridedMaps"))
585 AAAA->CreateView(
"stridedMaps", Teuchos::rcpFromRef(Op),
true);
589 else if(Helpers::isTpetraBlockCrs(Op)) {
590 using XMatrix = Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
591 using XCrsMatrix = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
592 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
593 using BCRS = Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
594 using CRS = Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
599 Tpetra::BlockCrsMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label);
601 using Teuchos::ParameterList;
603 RCP<ParameterList> transposeParams = params.is_null () ?
604 rcp (
new ParameterList) :
605 rcp (
new ParameterList (*params));
606 transposeParams->set (
"sort",
false);
607 At = transposer.createTranspose(transposeParams);
610 RCP<Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AA = rcp(
new Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(At));
611 RCP<XCrsMatrix> AAA = rcp_implicit_cast<XCrsMatrix>(AA);
612 RCP<XMatrix> AAAA = rcp(
new XCrsMatrixWrap(AAA));
614 if (Op.IsView(
"stridedMaps"))
615 AAAA->CreateView(
"stridedMaps", Teuchos::rcpFromRef(Op),
true);
625 std::cout <<
"Utilities::Transpose() not implemented for Epetra" << std::endl;
626 return Teuchos::null;
631 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
632 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
635 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xscalar;
636#if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))
637 typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType real_type;
639 if ((
typeid(
Scalar).name() ==
typeid(std::complex<double>).name()) ||
640 (
typeid(
Scalar).name() ==
typeid(std::complex<float>).name())) {
641 Xscalar = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X->getMap(),X->getNumVectors());
642 size_t numVecs = X->getNumVectors();
643 for (
size_t j=0;j<numVecs;j++) {
644 Teuchos::ArrayRCP<const real_type> XVec = X->getData(j);
645 Teuchos::ArrayRCP<Scalar> XVecScalar = Xscalar->getDataNonConst(j);
646 for(
size_t i = 0; i < static_cast<size_t>(XVec.size()); ++i)
647 XVecScalar[i]=XVec[i];
651 Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(X);
656 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
661 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,
LocalOrdinal,
GlobalOrdinal,
Node> > coordinates = Teuchos::null;
664 if(paramList.isParameter (
"Coordinates") ==
false)
667#if defined(HAVE_MUELU_TPETRA)
673#if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
674 typedef Tpetra::MultiVector<float, LocalOrdinal, GlobalOrdinal, Node> tfMV;
675 RCP<tfMV> floatCoords = Teuchos::null;
681#if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_DOUBLE)
683 RCP<tdMV> doubleCoords = Teuchos::null;
684 if (paramList.isType<RCP<tdMV> >(
"Coordinates")) {
686 doubleCoords = paramList.get<RCP<tdMV> >(
"Coordinates");
687 paramList.remove(
"Coordinates");
689#if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
690 else if (paramList.isType<RCP<tfMV> >(
"Coordinates")) {
692 floatCoords = paramList.get<RCP<tfMV> >(
"Coordinates");
693 paramList.remove(
"Coordinates");
694 doubleCoords = rcp(
new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
695 deep_copy(*doubleCoords, *floatCoords);
699 if(doubleCoords != Teuchos::null) {
701 coordinates = Teuchos::rcp_dynamic_cast<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,
LocalOrdinal,
GlobalOrdinal,
Node> >(MueLu::TpetraMultiVector_To_XpetraMultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,
LocalOrdinal,
GlobalOrdinal,
Node>(doubleCoords));
702 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(coordinates));
703 TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->getNumVectors() != coordinates->getNumVectors());
708 throw Exceptions::RuntimeError(
"ExtractCoordinatesFromParameterList: The coordinates vector in parameter list is expected to be a Tpetra multivector with SC=double or float.");
713 if(paramList.isType<
decltype(coordinates)>(
"Coordinates")) {
714 coordinates = paramList.get<
decltype(coordinates)>(
"Coordinates");
722#define MUELU_UTILITIES_SHORT
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraBlockCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static void MyOldScaleMatrix_Tpetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > RealValuedToScalarMultiVector(RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LocalOrdinal, GlobalOrdinal, Node > > X)
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
static void MyOldScaleMatrix_Epetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraBlockCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList ¶mList)
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
Namespace for MueLu classes and methods.
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)