47 #ifndef THYRA_MUELU_PRECONDITIONER_FACTORY_DECL_HPP 48 #define THYRA_MUELU_PRECONDITIONER_FACTORY_DECL_HPP 52 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 55 #include "Thyra_DefaultPreconditioner.hpp" 56 #include "Thyra_BlockedLinearOpBase.hpp" 58 #ifdef HAVE_MUELU_TPETRA 59 #include "Thyra_TpetraLinearOp.hpp" 60 #include "Thyra_TpetraThyraWrappers.hpp" 62 #ifdef HAVE_MUELU_EPETRA 63 #include "Thyra_EpetraLinearOp.hpp" 66 #include "Teuchos_Ptr.hpp" 67 #include "Teuchos_TestForException.hpp" 68 #include "Teuchos_Assert.hpp" 69 #include "Teuchos_Time.hpp" 71 #include <Xpetra_CrsMatrixWrap.hpp> 72 #include <Xpetra_CrsMatrix.hpp> 73 #include <Xpetra_Matrix.hpp> 74 #include <Xpetra_ThyraUtils.hpp> 76 #include <MueLu_Hierarchy.hpp> 78 #include <MueLu_HierarchyUtils.hpp> 79 #include <MueLu_Utilities.hpp> 80 #include <MueLu_ParameterListInterpreter.hpp> 81 #include <MueLu_MLParameterListInterpreter.hpp> 85 #ifdef HAVE_MUELU_TPETRA 86 #include <MueLu_TpetraOperator.hpp> 88 #ifdef HAVE_MUELU_EPETRA 92 #include "Thyra_PreconditionerFactoryBase.hpp" 94 #include "Kokkos_DefaultNode.hpp" 107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node = KokkosClassic::DefaultNode::DefaultNodeType>
108 class MueLuPreconditionerFactory :
public PreconditionerFactoryBase<Scalar> {
115 MueLuPreconditionerFactory();
122 bool isCompatible(
const LinearOpSourceBase<Scalar>& fwdOp)
const;
124 Teuchos::RCP<PreconditionerBase<Scalar> > createPrec()
const;
126 void initializePrec(
const Teuchos::RCP<
const LinearOpSourceBase<Scalar> >& fwdOp,
127 PreconditionerBase<Scalar>* prec,
128 const ESupportSolveUse supportSolveUse
131 void uninitializePrec(PreconditionerBase<Scalar>* prec,
132 Teuchos::RCP<
const LinearOpSourceBase<Scalar> >* fwdOp,
133 ESupportSolveUse* supportSolveUse
142 void setParameterList(
const Teuchos::RCP<Teuchos::ParameterList>& paramList);
144 Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
146 Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList();
148 Teuchos::RCP<const Teuchos::ParameterList> getParameterList()
const;
157 std::string description()
const;
167 Teuchos::RCP<Teuchos::ParameterList> paramList_;
171 #ifdef HAVE_MUELU_EPETRA 180 class MueLuPreconditionerFactory<double,int,int,Xpetra::
EpetraNode> :
public PreconditionerFactoryBase<double> {
191 MueLuPreconditionerFactory() : paramList_(rcp(new ParameterList())) { }
199 bool isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
200 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
202 #ifdef HAVE_MUELU_TPETRA 203 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
206 #ifdef HAVE_MUELU_EPETRA 207 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
210 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp))
return true;
216 Teuchos::RCP<PreconditionerBase<Scalar> > createPrec()
const {
217 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
221 void initializePrec(
const Teuchos::RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc,
222 PreconditionerBase<Scalar>* prec,
223 const ESupportSolveUse
225 using Teuchos::rcp_dynamic_cast;
228 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
229 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
230 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
231 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
232 typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
233 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
234 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
236 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
237 #ifdef HAVE_MUELU_TPETRA 239 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 240 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 242 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpOp;
243 typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
246 #if defined(HAVE_MUELU_EPETRA) 247 typedef MueLu::EpetraOperator MueEpOp;
248 typedef Thyra::EpetraLinearOp ThyEpLinOp;
256 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
257 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
258 TEUCHOS_ASSERT(prec);
261 ParameterList paramList = *paramList_;
264 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
265 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
268 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
269 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
270 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
271 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
272 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked ==
false);
273 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked ==
true);
275 RCP<XpMat> A = Teuchos::null;
277 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
278 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
279 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
281 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==
false);
283 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
284 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
286 RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
287 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
290 RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
291 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
294 RCP<XpMat> A00 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
295 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
297 RCP<const XpMap> rowmap00 = A00->getRowMap();
298 RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
301 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(
new XpBlockedCrsMat(ThyBlockedOp, comm));
302 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
307 RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
308 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
311 RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
312 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
315 A = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
317 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
320 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
321 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
324 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
325 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
328 RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = Teuchos::null;
333 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"reuse: type") || paramList.get<std::string>(
"reuse: type") ==
"none");
335 if (startingOver ==
true) {
337 Teuchos::RCP<XpMultVecDouble> coordinates = Teuchos::null;
341 RCP<XpMultVec> nullspace = Teuchos::null;
342 #ifdef HAVE_MUELU_TPETRA 344 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 345 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 346 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
347 RCP<tMV> tpetra_nullspace = Teuchos::null;
348 if (paramList.isType<Teuchos::RCP<tMV> >(
"Nullspace")) {
349 tpetra_nullspace = paramList.get<RCP<tMV> >(
"Nullspace");
350 paramList.remove(
"Nullspace");
351 nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
352 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
356 "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
360 #ifdef HAVE_MUELU_EPETRA 362 RCP<Epetra_MultiVector> epetra_nullspace = Teuchos::null;
363 if (paramList.isType<RCP<Epetra_MultiVector> >(
"Nullspace")) {
364 epetra_nullspace = paramList.get<RCP<Epetra_MultiVector> >(
"Nullspace");
365 paramList.remove(
"Nullspace");
366 RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpNullspace = Teuchos::rcp(
new Xpetra::EpetraMultiVectorT<int,Node>(epetra_nullspace));
367 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,int,int,
Node> > xpEpNullspaceMult = rcp_dynamic_cast<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,int,int,
Node> >(xpEpNullspace);
368 nullspace = rcp_dynamic_cast<XpMultVec>(xpEpNullspaceMult);
369 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
374 const std::string userName =
"user data";
375 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
376 if(Teuchos::nonnull(coordinates)) {
377 userParamList.set<RCP<XpMultVecDouble> >(
"Coordinates", coordinates);
379 if(Teuchos::nonnull(nullspace)) {
380 userParamList.set<RCP<XpMultVec> >(
"Nullspace", nullspace);
388 #if defined(HAVE_MUELU_TPETRA) 390 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 391 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 392 RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
393 RCP<MueTpOp> muelu_precOp = rcp_dynamic_cast<MueTpOp>(tpetr_precOp->getTpetraOperator(),
true);
395 H = muelu_precOp->GetHierarchy();
398 "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
402 #if defined(HAVE_MUELU_EPETRA)// && defined(HAVE_MUELU_SERIAL) 404 RCP<ThyEpLinOp> epetr_precOp = rcp_dynamic_cast<ThyEpLinOp>(thyra_precOp);
405 RCP<MueEpOp> muelu_precOp = rcp_dynamic_cast<MueEpOp>(epetr_precOp->epetra_op(),
true);
413 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
415 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
416 RCP<MueLu::Level> level0 = H->
GetLevel(0);
417 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >(
"A");
418 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
424 A->SetFixedBlockSize(A0->GetFixedBlockSize());
434 RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
435 #if defined(HAVE_MUELU_TPETRA) 437 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 438 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 439 RCP<MueTpOp> muelu_tpetraOp = rcp(
new MueTpOp(H));
440 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp));
441 RCP<TpOp> tpOp = Teuchos::rcp_dynamic_cast<TpOp>(muelu_tpetraOp);
442 thyraPrecOp = Thyra::createLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpOp);
445 "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
450 #if defined(HAVE_MUELU_EPETRA) 452 RCP<MueLu::Hierarchy<double,int,int,Xpetra::EpetraNode> > epetraH =
455 "Thyra::MueLuPreconditionerFactory: Failed to cast Hierarchy to Hierarchy<double,int,int,Xpetra::EpetraNode>. Epetra runs only on the Serial node.");
456 RCP<MueEpOp> muelu_epetraOp = rcp(
new MueEpOp(epetraH));
457 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_epetraOp));
459 set_extra_data(fwdOp,
"IFPF::fwdOp", Teuchos::inOutArg(muelu_epetraOp), Teuchos::POST_DESTROY,
false);
460 RCP<ThyEpLinOp> thyra_epetraOp = Thyra::nonconstEpetraLinearOp(muelu_epetraOp, NOTRANS, EPETRA_OP_APPLY_APPLY_INVERSE, EPETRA_OP_ADJOINT_UNSUPPORTED);
461 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyra_epetraOp));
462 thyraPrecOp = rcp_dynamic_cast<ThyLinOpBase>(thyra_epetraOp);
467 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::nonnull(thyraPrecOp));
470 const RCP<MueXpOp> muelu_xpetraOp = rcp(
new MueXpOp(H));
472 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getRangeMap());
473 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getDomainMap());
475 RCP <Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpOp = Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_xpetraOp);
476 thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
479 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
481 defaultPrec->initializeUnspecified(thyraPrecOp);
485 void uninitializePrec(PreconditionerBase<Scalar>* prec,
486 Teuchos::RCP<
const LinearOpSourceBase<Scalar> >* fwdOp,
487 ESupportSolveUse* supportSolveUse
489 TEUCHOS_ASSERT(prec);
492 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
493 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
497 *fwdOp = Teuchos::null;
500 if (supportSolveUse) {
502 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
505 defaultPrec->uninitialize();
514 void setParameterList(
const Teuchos::RCP<Teuchos::ParameterList>& paramList) {
515 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
516 paramList_ = paramList;
519 Teuchos::RCP<Teuchos::ParameterList> unsetParameterList() {
520 RCP<ParameterList> savedParamList = paramList_;
521 paramList_ = Teuchos::null;
522 return savedParamList;
525 Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList() {
return paramList_; }
527 Teuchos::RCP<const Teuchos::ParameterList> getParameterList()
const {
return paramList_; }
530 static RCP<const ParameterList> validPL;
532 if (Teuchos::is_null(validPL))
533 validPL = rcp(
new ParameterList());
543 std::string description()
const {
return "Thyra::MueLuPreconditionerFactory"; }
550 Teuchos::RCP<Teuchos::ParameterList> paramList_;
553 #endif // HAVE_MUELU_EPETRA 557 #endif // #ifdef HAVE_MUELU_STRATIMIKOS 559 #endif // THYRA_MUELU_PRECONDITIONER_FACTORY_DECL_HPP Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
MueLu::DefaultLocalOrdinal LocalOrdinal
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList ¶mList)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix...
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
Exception throws to report errors in the internal logical of the program.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
void getValidParameters(Teuchos::ParameterList ¶ms)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.