46 #ifndef MUELU_UTILITIES_DEF_HPP
47 #define MUELU_UTILITIES_DEF_HPP
49 #include <Teuchos_DefaultComm.hpp>
54 #ifdef HAVE_MUELU_EPETRA
56 # include "Epetra_MpiComm.h"
60 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
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_TpetraCrsMatrix.hpp>
81 #ifdef HAVE_MUELU_EPETRA
103 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_ML)
104 #include <ml_operator.h>
105 #include <ml_epetra_utils.h>
110 #ifdef HAVE_MUELU_EPETRA
115 #ifdef HAVE_MUELU_EPETRA
116 template<
typename SC,
typename LO,
typename GO,
typename NO>
118 return Xpetra::Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO>(epAB);
122 #ifdef HAVE_MUELU_EPETRA
123 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 if (tmpVec == Teuchos::null)
127 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
128 return tmpVec->getEpetra_MultiVector();
131 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
134 if (tmpVec == Teuchos::null)
135 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
136 return tmpVec->getEpetra_MultiVector();
139 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
145 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
151 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
154 if (crsOp == Teuchos::null)
157 if (tmp_ECrsMtx == Teuchos::null)
158 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
159 return tmp_ECrsMtx->getEpetra_CrsMatrix();
162 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
165 if (crsOp == Teuchos::null)
168 if (tmp_ECrsMtx == Teuchos::null)
169 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
170 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
173 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
180 }
catch (std::bad_cast) {
181 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
183 }
catch (std::bad_cast) {
188 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
195 }
catch (std::bad_cast) {
196 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
198 }
catch (std::bad_cast) {
203 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
206 if (xeMap == Teuchos::null)
207 throw Exceptions::BadCast(
"Utilities::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
208 return xeMap->getEpetra_Map();
212 #ifdef HAVE_MUELU_TPETRA
213 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
217 if (tmpVec == Teuchos::null)
218 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
219 return tmpVec->getTpetra_MultiVector();
222 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
225 if (tmpVec == Teuchos::null)
226 throw Exceptions::BadCast(
"Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
227 return tmpVec->getTpetra_MultiVector();
230 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
236 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
242 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
243 const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
249 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
253 if (crsOp == Teuchos::null)
256 if (tmp_ECrsMtx == Teuchos::null)
257 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
258 return tmp_ECrsMtx->getTpetra_CrsMatrix();
261 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
264 if (crsOp == Teuchos::null)
267 if (tmp_ECrsMtx == Teuchos::null)
268 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
269 return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
272 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
279 }
catch (std::bad_cast) {
280 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
282 }
catch (std::bad_cast) {
287 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
294 }
catch (std::bad_cast) {
295 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
297 }
catch (std::bad_cast) {
302 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
305 if (crsOp == Teuchos::null)
312 return tmp_Crs->getTpetra_CrsMatrixNonConst();
316 if (tmp_BlockCrs.is_null())
317 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
318 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
322 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
325 if (crsOp == Teuchos::null)
332 return tmp_Crs->getTpetra_CrsMatrixNonConst();
336 if (tmp_BlockCrs.is_null())
337 throw Exceptions::BadCast(
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
338 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
343 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
346 if (tmp_TMap == Teuchos::null)
347 throw Exceptions::BadCast(
"Utilities::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
348 return tmp_TMap->getTpetra_Map();
352 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
355 bool doOptimizeStorage)
360 for (
int i = 0; i < scalingVector.
size(); ++i)
361 sv[i] = one / scalingVector[i];
363 for (
int i = 0; i < scalingVector.
size(); ++i)
364 sv[i] = scalingVector[i];
367 switch (Op.getRowMap()->lib()) {
369 MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
373 MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
381 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
386 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
389 bool doOptimizeStorage)
391 #ifdef HAVE_MUELU_TPETRA
393 Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
399 size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
400 if (maxRowSize == Teuchos::as<size_t>(-1))
403 std::vector<Scalar> scaledVals(maxRowSize);
404 if (tpOp.isFillComplete())
411 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
412 tpOp.getLocalRowView(i, cols, vals);
413 size_t nnz = tpOp.getNumEntriesInLocalRow(i);
414 if (nnz > maxRowSize) {
416 scaledVals.resize(maxRowSize);
418 for (
size_t j = 0; j < nnz; ++j)
419 scaledVals[j] = vals[j]*scalingVector[i];
423 tpOp.replaceLocalValues(i, cols, valview);
431 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
432 GlobalOrdinal gid = rowMap->getGlobalElement(i);
433 tpOp.getGlobalRowView(gid, cols, vals);
434 size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
435 if (nnz > maxRowSize) {
437 scaledVals.resize(maxRowSize);
440 for (
size_t j = 0; j < nnz; ++j)
441 scaledVals[j] = vals[j]*scalingVector[i];
445 tpOp.replaceGlobalValues(gid, cols, valview);
450 if (doFillComplete) {
451 if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
452 throw Exceptions::RuntimeError(
"In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
455 params->
set(
"Optimize Storage", doOptimizeStorage);
456 params->
set(
"No Nonlocal Changes",
true);
467 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
471 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
472 std::string TorE =
"epetra";
474 std::string TorE =
"tpetra";
477 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
486 #ifdef HAVE_MUELU_TPETRA
487 if (TorE ==
"tpetra") {
492 Tpetra::RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label);
493 A = transposer.createTranspose(params);
498 if (!AAAA->isFillComplete())
501 if (Op.IsView(
"stridedMaps"))
502 AAAA->CreateView(
"stridedMaps", Teuchos::rcpFromRef(Op),
true);
506 }
catch (std::exception& e) {
507 std::cout <<
"threw exception '" << e.what() <<
"'" << std::endl;
514 std::cout <<
"Utilities::Transpose() not implemented for Epetra" << std::endl;
515 return Teuchos::null;
520 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
525 #if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))
528 if ((
typeid(Scalar).name() ==
typeid(std::complex<double>).name()) ||
529 (
typeid(Scalar).name() ==
typeid(std::complex<float>).name())) {
531 size_t numVecs = X->getNumVectors();
532 for (
size_t j=0;j<numVecs;j++) {
535 for(
size_t i = 0; i < static_cast<size_t>(XVec.
size()); ++i)
536 XVecScalar[i]=XVec[i];
545 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
556 #if defined(HAVE_MUELU_TPETRA)
562 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
563 typedef Tpetra::MultiVector<float, LocalOrdinal, GlobalOrdinal, Node> tfMV;
570 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_DOUBLE)
571 typedef Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node> tdMV;
575 doubleCoords = paramList.
get<
RCP<tdMV> >(
"Coordinates");
576 paramList.
remove(
"Coordinates");
578 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
582 paramList.
remove(
"Coordinates");
583 doubleCoords =
rcp(
new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
588 if(doubleCoords != Teuchos::null) {
597 throw Exceptions::RuntimeError(
"ExtractCoordinatesFromParameterList: The coordinates vector in parameter list is expected to be a Tpetra multivector with SC=double or float.");
599 #endif // endif HAVE_TPETRA
602 if(paramList.
isType<decltype(coordinates)>(
"Coordinates")) {
603 coordinates = paramList.
get<decltype(coordinates)>(
"Coordinates");
611 #define MUELU_UTILITIES_SHORT
612 #endif // MUELU_UTILITIES_DEF_HPP
Exception indicating invalid cast attempted.
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.
RCP< CrsMatrix > getCrsMatrix() const
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 RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
T & get(const std::string &name, T def_value)
RCP< Epetra_CrsMatrix > getEpetra_CrsMatrixNonConst() const
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList ¶mList)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
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.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
bool isParameter(const std::string &name) const
bool remove(std::string const &name, bool throwIfNotExists=true)
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > ¶ms=null)=0
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > RealValuedToScalarMultiVector(RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > X)
void deep_copy(const DynRankView< DT, DP...> &dst, typename ViewTraits< DT, DP...>::const_value_type &value, typename std::enable_if< std::is_same< typename ViewTraits< DT, DP...>::specialize, void >::value >::type *=0)
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
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 Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)
static void MyOldScaleMatrix_Epetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
virtual bool isLocallyIndexed() const =0
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrixNonConst() const
bool isType(const std::string &name) const
RCP< Epetra_MultiVector > getEpetra_MultiVector() const
virtual Teuchos::RCP< const Map > getRangeMap() const =0
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_MultiVector() const
Exception throws to report errors in the internal logical of the program.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< 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)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
virtual Teuchos::RCP< const Map > getDomainMap() const =0