diff --git a/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_decl.hpp b/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_decl.hpp index 2ac74705d79a..5a2764b9f75c 100644 --- a/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_decl.hpp +++ b/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_decl.hpp @@ -54,6 +54,7 @@ // Stratimikos needs Thyra, so we don't need special guards for Thyra here #include "Thyra_DefaultPreconditioner.hpp" #include "Thyra_BlockedLinearOpBase.hpp" +#include "Thyra_DiagonalLinearOpBase.hpp" #include "Thyra_XpetraLinearOp.hpp" #ifdef HAVE_MUELU_TPETRA #include "Thyra_TpetraLinearOp.hpp" @@ -84,6 +85,7 @@ #include #ifdef HAVE_MUELU_TPETRA #include +#include #endif #ifdef HAVE_MUELU_EPETRA #include @@ -93,16 +95,20 @@ #include "Kokkos_DefaultNode.hpp" +#include namespace Thyra { + using Teuchos::RCP; + using Teuchos::rcp; + + template + static bool replaceWithXpetra(ParameterList& paramList, std::string parameterName); + /** @brief Concrete preconditioner factory subclass for Thyra based on MueLu. @ingroup MueLuAdapters Add support for MueLu preconditioners in Thyra. This class provides an interface both for Epetra and Tpetra. - - The general implementation only handles Tpetra. For Epetra there is a specialization - on SC=double, LO=int, GO=int and NO=EpetraNode. */ template class MueLuPreconditionerFactory : public PreconditionerFactoryBase { @@ -168,390 +174,6 @@ namespace Thyra { }; -#ifdef HAVE_MUELU_EPETRA - /** @brief Concrete preconditioner factory subclass for Thyra based on MueLu. - @ingroup MueLuAdapters - Add support for MueLu preconditioners in Thyra. This class provides an interface both - for Epetra and Tpetra. - - Specialization for Epetra - */ - template <> - class MueLuPreconditionerFactory : public PreconditionerFactoryBase { - public: - typedef double Scalar; - typedef int LocalOrdinal; - typedef int GlobalOrdinal; - typedef Xpetra::EpetraNode Node; - - /** @name Constructors/initializers/accessors */ - //@{ - - /** \brief . */ - MueLuPreconditionerFactory() : paramList_(rcp(new ParameterList())) { } - - //@} - - /** @name Overridden from PreconditionerFactoryBase */ - //@{ - - /** \brief . */ - bool isCompatible(const LinearOpSourceBase& fwdOpSrc) const { - const RCP > fwdOp = fwdOpSrc.getOp(); - -#ifdef HAVE_MUELU_TPETRA - if (Xpetra::ThyraUtils::isTpetra(fwdOp)) return true; -#endif - -#ifdef HAVE_MUELU_EPETRA - if (Xpetra::ThyraUtils::isEpetra(fwdOp)) return true; -#endif - - if (Xpetra::ThyraUtils::isBlockedOperator(fwdOp)) return true; - - return false; - } - - /** \brief . */ - Teuchos::RCP > createPrec() const { - return Teuchos::rcp(new DefaultPreconditioner); - } - - /** \brief . */ - void initializePrec(const Teuchos::RCP >& fwdOpSrc, - PreconditionerBase* prec, - const ESupportSolveUse /* supportSolveUse */ - ) const { - using Teuchos::rcp_dynamic_cast; - - // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...) - typedef Xpetra::Map XpMap; - typedef Xpetra::Operator XpOp; - typedef Xpetra::ThyraUtils XpThyUtils; - typedef Xpetra::CrsMatrix XpCrsMat; - typedef Xpetra::BlockedCrsMatrix XpBlockedCrsMat; - typedef Xpetra::Matrix XpMat; - typedef Xpetra::MultiVector XpMultVec; - typedef Xpetra::MultiVector::magnitudeType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble; - typedef Thyra::LinearOpBase ThyLinOpBase; -#ifdef HAVE_MUELU_TPETRA - // TAW 1/26/2016: We deal with Tpetra objects -#if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ - (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) - typedef MueLu::TpetraOperator MueTpOp; - typedef Tpetra::Operator TpOp; - typedef Thyra::TpetraLinearOp ThyTpLinOp; -#endif -#endif -#if defined(HAVE_MUELU_EPETRA) - typedef MueLu::EpetraOperator MueEpOp; - typedef Thyra::EpetraLinearOp ThyEpLinOp; -#endif - - //std::cout << "-======---------------------------------" << std::endl; - //std::cout << *paramList_ << std::endl; - //std::cout << "-======---------------------------------" << std::endl; - - // Check precondition - TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc)); - TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc)); - TEUCHOS_ASSERT(prec); - - // Create a copy, as we may remove some things from the list - ParameterList paramList = *paramList_; - - // Retrieve wrapped concrete Xpetra matrix from FwdOp - const RCP fwdOp = fwdOpSrc->getOp(); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp)); - - // Check whether it is Epetra/Tpetra - bool bIsEpetra = XpThyUtils::isEpetra(fwdOp); - bool bIsTpetra = XpThyUtils::isTpetra(fwdOp); - bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp); - TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true)); - TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false); - TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true); - - RCP A = Teuchos::null; - if(bIsBlocked) { - Teuchos::RCP > ThyBlockedOp = - Teuchos::rcp_dynamic_cast >(fwdOp); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp)); - - TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false); - - Teuchos::RCP > b00 = ThyBlockedOp->getBlock(0,0); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00)); - - RCP xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00)); - - // MueLu needs a non-const object as input - RCP xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast(xpetraFwdCrsMat00); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00)); - - // wrap the forward operator as an Xpetra::Matrix that MueLu can work with - RCP A00 = rcp(new Xpetra::CrsMatrixWrap(xpetraFwdCrsMatNonConst00)); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00)); - - RCP rowmap00 = A00->getRowMap(); - RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm(); - - // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with - RCP bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm)); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat)); - - // save blocked matrix - A = bMat; - } else { - RCP xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat)); - - // MueLu needs a non-const object as input - RCP xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast(xpetraFwdCrsMat); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst)); - - // wrap the forward operator as an Xpetra::Matrix that MueLu can work with - A = rcp(new Xpetra::CrsMatrixWrap(xpetraFwdCrsMatNonConst)); - } - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A)); - - // Retrieve concrete preconditioner object - const Teuchos::Ptr > defaultPrec = Teuchos::ptr(dynamic_cast *>(prec)); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec)); - - // extract preconditioner operator - RCP thyra_precOp = Teuchos::null; - thyra_precOp = rcp_dynamic_cast >(defaultPrec->getNonconstUnspecifiedPrecOp(), true); - - // Variable for multigrid hierarchy: either build a new one or reuse the existing hierarchy - RCP > H = Teuchos::null; - - // make a decision whether to (re)build the multigrid preconditioner or reuse the old one - // rebuild preconditioner if startingOver == true - // reuse preconditioner if startingOver == false - const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get("reuse: type") == "none"); - - if (startingOver == true) { - // extract coordinates from parameter list - Teuchos::RCP coordinates = Teuchos::null; - coordinates = MueLu::Utilities::ExtractCoordinatesFromParameterList(paramList); - - // TODO check for Xpetra or Thyra vectors? - RCP nullspace = Teuchos::null; -#ifdef HAVE_MUELU_TPETRA - if (bIsTpetra) { -#if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ - (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) - typedef Tpetra::MultiVector tMV; - RCP tpetra_nullspace = Teuchos::null; - if (paramList.isType >("Nullspace")) { - tpetra_nullspace = paramList.get >("Nullspace"); - paramList.remove("Nullspace"); - nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_nullspace); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace)); - } -#else - TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode."); -#endif - } -#endif -#ifdef HAVE_MUELU_EPETRA - if (bIsEpetra) { - RCP epetra_nullspace = Teuchos::null; - if (paramList.isType >("Nullspace")) { - epetra_nullspace = paramList.get >("Nullspace"); - paramList.remove("Nullspace"); - RCP > xpEpNullspace = Teuchos::rcp(new Xpetra::EpetraMultiVectorT(epetra_nullspace)); - RCP::magnitudeType,int,int,Node> > xpEpNullspaceMult = rcp_dynamic_cast::magnitudeType,int,int,Node> >(xpEpNullspace); - nullspace = rcp_dynamic_cast(xpEpNullspaceMult); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace)); - } - } -#endif - // build a new MueLu hierarchy - const std::string userName = "user data"; - Teuchos::ParameterList& userParamList = paramList.sublist(userName); - if(Teuchos::nonnull(coordinates)) { - userParamList.set >("Coordinates", coordinates); - } - if(Teuchos::nonnull(nullspace)) { - userParamList.set >("Nullspace", nullspace); - } - H = MueLu::CreateXpetraPreconditioner(A, paramList); - - } else { - // reuse old MueLu hierarchy stored in MueLu Tpetra/Epetra operator and put in new matrix - - // get old MueLu hierarchy -#if defined(HAVE_MUELU_TPETRA) - if (bIsTpetra) { -#if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ - (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) - RCP tpetr_precOp = rcp_dynamic_cast(thyra_precOp); - RCP muelu_precOp = rcp_dynamic_cast(tpetr_precOp->getTpetraOperator(),true); - - H = muelu_precOp->GetHierarchy(); -#else - TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode."); -#endif - } -#endif -#if defined(HAVE_MUELU_EPETRA)// && defined(HAVE_MUELU_SERIAL) - if (bIsEpetra) { - RCP epetr_precOp = rcp_dynamic_cast(thyra_precOp); - RCP muelu_precOp = rcp_dynamic_cast(epetr_precOp->epetra_op(),true); - - H = rcp_dynamic_cast >(muelu_precOp->GetHierarchy()); - } -#endif - // TODO add the blocked matrix case here... - - TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it"); - TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator"); - RCP level0 = H->GetLevel(0); - RCP O0 = level0->Get >("A"); - RCP A0 = rcp_dynamic_cast(O0); - - if (!A0.is_null()) { - // If a user provided a "number of equations" argument in a parameter list - // during the initial setup, we must honor that settings and reuse it for - // all consequent setups. - A->SetFixedBlockSize(A0->GetFixedBlockSize()); - } - - // set new matrix - level0->Set("A", A); - - H->SetupRe(); - } - - // wrap hierarchy H in thyraPrecOp - RCP thyraPrecOp = Teuchos::null; -#if defined(HAVE_MUELU_TPETRA) - if (bIsTpetra) { -#if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ - (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) - RCP muelu_tpetraOp = rcp(new MueTpOp(H)); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp)); - RCP tpOp = Teuchos::rcp_dynamic_cast(muelu_tpetraOp); - thyraPrecOp = Thyra::createLinearOp(tpOp); -#else - TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode."); -#endif - } -#endif - -#if defined(HAVE_MUELU_EPETRA) - if (bIsEpetra) { - RCP > epetraH = - rcp_dynamic_cast >(H); - TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(epetraH), MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Failed to cast Hierarchy to Hierarchy. Epetra runs only on the Serial node."); - RCP muelu_epetraOp = rcp(new MueEpOp(epetraH)); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_epetraOp)); - // attach fwdOp to muelu_epetraOp to guarantee that it will not go away - set_extra_data(fwdOp,"IFPF::fwdOp", Teuchos::inOutArg(muelu_epetraOp), Teuchos::POST_DESTROY,false); - RCP thyra_epetraOp = Thyra::nonconstEpetraLinearOp(muelu_epetraOp, NOTRANS, EPETRA_OP_APPLY_APPLY_INVERSE, EPETRA_OP_ADJOINT_UNSUPPORTED); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyra_epetraOp)); - thyraPrecOp = rcp_dynamic_cast(thyra_epetraOp); - } -#endif - - if(bIsBlocked) { - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::nonnull(thyraPrecOp)); - - typedef MueLu::XpetraOperator MueXpOp; - const RCP muelu_xpetraOp = rcp(new MueXpOp(H)); - - RCP > thyraRangeSpace = Xpetra::ThyraUtils::toThyra(muelu_xpetraOp->getRangeMap()); - RCP > thyraDomainSpace = Xpetra::ThyraUtils::toThyra(muelu_xpetraOp->getDomainMap()); - - RCP > xpOp = Teuchos::rcp_dynamic_cast >(muelu_xpetraOp); - thyraPrecOp = Thyra::xpetraLinearOp(thyraRangeSpace, thyraDomainSpace,xpOp); - } - - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp)); - - defaultPrec->initializeUnspecified(thyraPrecOp); - } - - /** \brief . */ - void uninitializePrec(PreconditionerBase* prec, - Teuchos::RCP >* fwdOp, - ESupportSolveUse* supportSolveUse - ) const { - TEUCHOS_ASSERT(prec); - - // Retrieve concrete preconditioner object - const Teuchos::Ptr > defaultPrec = Teuchos::ptr(dynamic_cast *>(prec)); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec)); - - if (fwdOp) { - // TODO: Implement properly instead of returning default value - *fwdOp = Teuchos::null; - } - - if (supportSolveUse) { - // TODO: Implement properly instead of returning default value - *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED; - } - - defaultPrec->uninitialize(); - } - - //@} - - /** @name Overridden from Teuchos::ParameterListAcceptor */ - //@{ - - /** \brief . */ - void setParameterList(const Teuchos::RCP& paramList) { - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList)); - paramList_ = paramList; - } - /** \brief . */ - Teuchos::RCP unsetParameterList() { - RCP savedParamList = paramList_; - paramList_ = Teuchos::null; - return savedParamList; - } - /** \brief . */ - Teuchos::RCP getNonconstParameterList() { return paramList_; } - /** \brief . */ - Teuchos::RCP getParameterList() const { return paramList_; } - /** \brief . */ - Teuchos::RCP getValidParameters() const { - static RCP validPL; - - if (Teuchos::is_null(validPL)) - validPL = rcp(new ParameterList()); - - return validPL; - } - //@} - - /** \name Public functions overridden from Describable. */ - //@{ - - /** \brief . */ - std::string description() const { return "Thyra::MueLuPreconditionerFactory"; } - - // ToDo: Add an override of describe(...) to give more detail! - - //@} - - private: - Teuchos::RCP paramList_; - }; // end specialization for Epetra - -#endif // HAVE_MUELU_EPETRA - } // namespace Thyra #endif // #ifdef HAVE_MUELU_STRATIMIKOS diff --git a/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp b/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp index f2db6fa2695b..f136a8975441 100644 --- a/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp +++ b/packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp @@ -55,8 +55,176 @@ namespace Thyra { using Teuchos::RCP; using Teuchos::rcp; using Teuchos::ParameterList; + using Teuchos::rcp_dynamic_cast; + using Teuchos::rcp_const_cast; + template + static bool replaceWithXpetra(ParameterList& paramList, std::string parameterName) { + typedef typename Teuchos::ScalarTraits::magnitudeType Magnitude; + typedef Xpetra::Operator XpOp; + typedef Xpetra::ThyraUtils XpThyUtils; + typedef Xpetra::CrsMatrixWrap XpCrsMatWrap; + typedef Xpetra::CrsMatrix XpCrsMat; + typedef Xpetra::Matrix XpMat; + typedef Xpetra::MultiVector XpMultVec; + typedef Xpetra::MultiVector XpMagMultVec; + typedef Xpetra::Vector XpVec; + + typedef Thyra::LinearOpBase ThyLinOpBase; + typedef Thyra::DiagonalLinearOpBase ThyDiagLinOpBase; + typedef Thyra::XpetraLinearOp ThyXpOp; + typedef Thyra::SpmdVectorSpaceBase ThyVSBase; + +#ifdef HAVE_MUELU_TPETRA + typedef Tpetra::CrsMatrix TpCrsMat; + typedef Tpetra::Vector tV; + typedef Thyra::TpetraVector thyTpV; + typedef Tpetra::MultiVector tMV; + typedef Tpetra::MultiVector tMagMV; +# if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) + typedef typename Teuchos::ScalarTraits::halfPrecision HalfMagnitude; + typedef Tpetra::MultiVector tHalfMagMV; +# endif +#endif +#if defined(HAVE_MUELU_EPETRA) + typedef Xpetra::EpetraCrsMatrixT XpEpCrsMat; +#endif + + if (paramList.isParameter(parameterName)) { + if (paramList.isType >(parameterName)) + return true; + else if (paramList.isType >(parameterName)) { + RCP constM = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP M = rcp_const_cast(constM); + paramList.set >(parameterName, M); + return true; + } + else if (paramList.isType >(parameterName)) + return true; + else if (paramList.isType >(parameterName)) { + RCP constX = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP X = rcp_const_cast(constX); + paramList.set >(parameterName, X); + return true; + } + else if (paramList.isType >(parameterName)) + return true; + else if (paramList.isType >(parameterName)) { + RCP constX = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP X = rcp_const_cast(constX); + paramList.set >(parameterName, X); + return true; + } +#ifdef HAVE_MUELU_TPETRA + else if (paramList.isType >(parameterName)) { + RCP tM = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP xM = rcp_dynamic_cast(tM, true); + paramList.set >(parameterName, xM); + return true; + } else if (paramList.isType >(parameterName)) { + RCP tpetra_X = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP X = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_X); + paramList.set >(parameterName, X); + TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X)); + return true; + } else if (paramList.isType >(parameterName)) { + RCP tpetra_X = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP X = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_X); + paramList.set >(parameterName, X); + TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X)); + return true; + } +# if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) + else if (paramList.isType >(parameterName)) { + RCP tpetra_hX = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP tpetra_X = rcp(new tMagMV(tpetra_hX->getMap(),tpetra_hX->getNumVectors())); + Tpetra::deep_copy(*tpetra_X,*tpetra_hX); + RCP X = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_X); + paramList.set >(parameterName, X); + TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X)); + return true; + } +# endif +#endif +#ifdef HAVE_MUELU_EPETRA + else if (paramList.isType >(parameterName)) { + RCP eM = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP xeM = rcp(new XpEpCrsMat(eM)); + RCP xCrsM = rcp_dynamic_cast(xeM, true); + RCP xwM = rcp(new XpCrsMatWrap(xCrsM)); + RCP xM = rcp_dynamic_cast(xwM); + paramList.set >(parameterName, xM); + return true; + } else if (paramList.isType >(parameterName)) { + RCP epetra_X = Teuchos::null; + epetra_X = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP > xpEpX = rcp(new Xpetra::EpetraMultiVectorT(epetra_X)); + RCP > xpEpXMult = rcp_dynamic_cast >(xpEpX, true); + RCP X = rcp_dynamic_cast(xpEpXMult, true); + paramList.set >(parameterName, X); + return true; + } +#endif + else if (paramList.isType >(parameterName)) { + RCP thyM = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP crsM = XpThyUtils::toXpetra(thyM); + TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM)); + // MueLu needs a non-const object as input + RCP crsMNonConst = rcp_const_cast(crsM); + TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMNonConst)); + // wrap as an Xpetra::Matrix that MueLu can work with + RCP M = rcp(new Xpetra::CrsMatrixWrap(crsMNonConst)); + paramList.set >(parameterName, M); + return true; + } else if (paramList.isType >(parameterName)) { + RCP thyM = paramList.get >(parameterName); + paramList.remove(parameterName); + RCP > diag = thyM->getDiag(); + + RCP xpDiag; +#ifdef HAVE_MUELU_TPETRA + if (!rcp_dynamic_cast(diag).is_null()) { + RCP tDiag = Thyra::TpetraOperatorVectorExtraction::getConstTpetraVector(diag); + if (!tDiag.is_null()) + xpDiag = Xpetra::toXpetra(tDiag); + } +#endif +#ifdef HAVE_MUELU_EPETRA + if (xpDiag.is_null()) { + RCP comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast(thyM->range())->getComm()); + RCP map = Thyra::get_Epetra_Map(*(thyM->range()), comm); + if (!map.is_null()) { + RCP eDiag = Thyra::get_Epetra_Vector(*map, diag); + RCP nceDiag = rcp_const_cast(eDiag); + RCP > xpEpDiag = rcp(new Xpetra::EpetraVectorT(nceDiag)); + xpDiag = rcp_dynamic_cast(xpEpDiag, true); + } + } +#endif + TEUCHOS_ASSERT(!xpDiag.is_null()); + RCP M = Xpetra::MatrixFactory::Build(xpDiag); + paramList.set >(parameterName, M); + return true; + } + else { + TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter " << parameterName << " has wrong type."); + return false; + } + } else + return false; + } + // Constructors/initializers/accessors template @@ -74,6 +242,11 @@ namespace Thyra { if (Xpetra::ThyraUtils::isTpetra(fwdOp)) return true; #endif +#ifdef HAVE_MUELU_EPETRA + if (Xpetra::ThyraUtils::isEpetra(fwdOp)) return true; +#endif + + if (Xpetra::ThyraUtils::isBlockedOperator(fwdOp)) return true; return false; @@ -93,6 +266,7 @@ namespace Thyra { // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...) typedef Xpetra::Map XpMap; typedef Xpetra::Operator XpOp; + typedef MueLu::XpetraOperator MueLuXpOp; typedef Xpetra::ThyraUtils XpThyUtils; typedef Xpetra::CrsMatrix XpCrsMat; typedef Xpetra::BlockedCrsMatrix XpBlockedCrsMat; @@ -100,12 +274,22 @@ namespace Thyra { typedef Xpetra::MultiVector XpMultVec; typedef Xpetra::MultiVector::coordinateType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble; typedef Thyra::LinearOpBase ThyLinOpBase; -#ifdef HAVE_MUELU_TPETRA - typedef MueLu::TpetraOperator MueTpOp; - typedef Tpetra::Operator TpOp; - typedef Thyra::TpetraLinearOp ThyTpLinOp; + typedef Thyra::XpetraLinearOp ThyXpOp; + typedef Xpetra::MultiVector XpMV; + typedef typename Teuchos::ScalarTraits::magnitudeType Magnitude; + typedef Xpetra::MultiVector XpmMV; +#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) + typedef Xpetra::TpetraHalfPrecisionOperator XpHalfPrecOp; + typedef typename XpHalfPrecOp::HalfScalar HalfScalar; + typedef Xpetra::Operator XpHalfOp; + typedef MueLu::XpetraOperator MueLuHalfXpOp; + typedef typename Teuchos::ScalarTraits::halfPrecision HalfMagnitude; + typedef Xpetra::MultiVector XphMV; + typedef Xpetra::MultiVector XphmMV; + typedef Xpetra::Matrix XphMat; #endif + // Check precondition TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc)); TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc)); @@ -178,104 +362,144 @@ namespace Thyra { RCP thyra_precOp = Teuchos::null; thyra_precOp = rcp_dynamic_cast >(defaultPrec->getNonconstUnspecifiedPrecOp(), true); - // Variable for multigrid hierarchy: either build a new one or reuse the existing hierarchy - RCP > H = Teuchos::null; - // make a decision whether to (re)build the multigrid preconditioner or reuse the old one // rebuild preconditioner if startingOver == true // reuse preconditioner if startingOver == false const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get("reuse: type") == "none"); + const bool useHalfPrecision = paramList.get("half precision", false); + if (useHalfPrecision) + TEUCHOS_TEST_FOR_EXCEPTION(!bIsTpetra, MueLu::Exceptions::RuntimeError, "The only scalar type Epetra knows is double, so a half precision preconditioner cannot be constructed."); + RCP xpPrecOp; if (startingOver == true) { - // extract coordinates from parameter list - RCP coordinates = Teuchos::null; - coordinates = MueLu::Utilities::ExtractCoordinatesFromParameterList(paramList); - - // TODO check for Xpetra or Thyra vectors? - RCP nullspace = Teuchos::null; -#ifdef HAVE_MUELU_TPETRA - if (bIsTpetra) { - typedef Tpetra::MultiVector tMV; - RCP tpetra_nullspace = Teuchos::null; - if (paramList.isType >("Nullspace")) { - tpetra_nullspace = paramList.get >("Nullspace"); + // Convert to Xpetra + std::list convertXpetra = {"Coordinates", "Nullspace"}; + for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it) + replaceWithXpetra(paramList,*it); + + if (useHalfPrecision) { +#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) + + // convert to half precision + RCP halfA = Xpetra::convertToHalfPrecision(A); + const std::string userName = "user data"; + Teuchos::ParameterList& userParamList = paramList.sublist(userName); + if (userParamList.isType >("Coordinates")) { + RCP coords = userParamList.get >("Coordinates"); + userParamList.remove("Coordinates"); + RCP halfCoords = Xpetra::convertToHalfPrecision(coords); + userParamList.set("Coordinates",halfCoords); + } + if (userParamList.isType >("Nullspace")) { + RCP nullspace = userParamList.get >("Nullspace"); + userParamList.remove("Nullspace"); + RCP halfNullspace = Xpetra::convertToHalfPrecision(nullspace); + userParamList.set("Nullspace",halfNullspace); + } + if (paramList.isType >("Coordinates")) { + RCP coords = paramList.get >("Coordinates"); + paramList.remove("Coordinates"); + RCP halfCoords = Xpetra::convertToHalfPrecision(coords); + userParamList.set("Coordinates",halfCoords); + } + if (paramList.isType >("Nullspace")) { + RCP nullspace = paramList.get >("Nullspace"); paramList.remove("Nullspace"); - nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_nullspace); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace)); + RCP halfNullspace = Xpetra::convertToHalfPrecision(nullspace); + userParamList.set("Nullspace",halfNullspace); } - } -#endif - // build a new MueLu hierarchy - ParameterList& userParamList = paramList.sublist("user data"); - if(Teuchos::nonnull(coordinates)) { - userParamList.set >("Coordinates", coordinates); - } - if(Teuchos::nonnull(nullspace)) { - userParamList.set >("Nullspace", nullspace); - } - H = MueLu::CreateXpetraPreconditioner(A, paramList); - } else { - // reuse old MueLu hierarchy stored in MueLu Tpetra/Epetra operator and put in new matrix - // get old MueLu hierarchy -#if defined(HAVE_MUELU_TPETRA) - if (bIsTpetra) { + // build a new half-precision MueLu preconditioner - RCP tpetr_precOp = rcp_dynamic_cast(thyra_precOp); - RCP muelu_precOp = rcp_dynamic_cast(tpetr_precOp->getTpetraOperator(),true); - - H = muelu_precOp->GetHierarchy(); - } + RCP > H = MueLu::CreateXpetraPreconditioner(halfA, paramList); + RCP xpOp = rcp(new MueLuHalfXpOp(H)); + xpPrecOp = rcp(new XpHalfPrecOp(xpOp)); +#else + TEUCHOS_TEST_FOR_EXCEPT(true); #endif - // TODO add the blocked matrix case here... - - TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it"); - TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError, - "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator"); - RCP level0 = H->GetLevel(0); - RCP O0 = level0->Get >("A"); - RCP A0 = rcp_dynamic_cast(O0); - - if (!A0.is_null()) { - // If a user provided a "number of equations" argument in a parameter list - // during the initial setup, we must honor that settings and reuse it for - // all consequent setups. - A->SetFixedBlockSize(A0->GetFixedBlockSize()); - } - - // set new matrix - level0->Set("A", A); + } else + { + const std::string userName = "user data"; + Teuchos::ParameterList& userParamList = paramList.sublist(userName); + if (paramList.isType >("Coordinates")) { + RCP coords = paramList.get >("Coordinates"); + paramList.remove("Coordinates"); + userParamList.set("Coordinates",coords); + } + if (paramList.isType >("Nullspace")) { + RCP nullspace = paramList.get >("Nullspace"); + paramList.remove("Nullspace"); + userParamList.set("Nullspace",nullspace); + } + + // build a new MueLu RefMaxwell preconditioner + RCP > H = MueLu::CreateXpetraPreconditioner(A, paramList); + xpPrecOp = rcp(new MueLuXpOp(H)); + } + } else { + // reuse old MueLu hierarchy stored in MueLu Xpetra operator and put in new matrix + RCP thyXpOp = rcp_dynamic_cast(thyra_precOp, true); + xpPrecOp = rcp_dynamic_cast(thyXpOp->getXpetraOperator(), true); +#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) + RCP xpHalfPrecOp = rcp_dynamic_cast(xpPrecOp); + if (!xpHalfPrecOp.is_null()) { + RCP > H = rcp_dynamic_cast(xpHalfPrecOp->GetHalfPrecisionOperator(), true)->GetHierarchy(); + RCP halfA = Xpetra::convertToHalfPrecision(A); + + TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError, + "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it"); + TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError, + "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator"); + RCP level0 = H->GetLevel(0); + RCP O0 = level0->Get >("A"); + RCP A0 = rcp_dynamic_cast(O0, true); + + if (!A0.is_null()) { + // If a user provided a "number of equations" argument in a parameter list + // during the initial setup, we must honor that settings and reuse it for + // all consequent setups. + halfA->SetFixedBlockSize(A0->GetFixedBlockSize()); + } - H->SetupRe(); - } + // set new matrix + level0->Set("A", halfA); - // wrap hierarchy H in thyraPrecOp - RCP thyraPrecOp = Teuchos::null; -#if defined(HAVE_MUELU_TPETRA) - if (bIsTpetra) { - RCP muelu_tpetraOp = rcp(new MueTpOp(H)); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp)); - RCP tpOp = Teuchos::rcp_dynamic_cast(muelu_tpetraOp); - thyraPrecOp = Thyra::createLinearOp(tpOp); - } + H->SetupRe(); + } else #endif + { + // get old MueLu hierarchy + RCP xpOp = rcp_dynamic_cast(thyXpOp->getXpetraOperator(), true); + RCP > H = xpOp->GetHierarchy();; + + TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError, + "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it"); + TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError, + "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator"); + RCP level0 = H->GetLevel(0); + RCP O0 = level0->Get >("A"); + RCP A0 = rcp_dynamic_cast(O0); + + if (!A0.is_null()) { + // If a user provided a "number of equations" argument in a parameter list + // during the initial setup, we must honor that settings and reuse it for + // all consequent setups. + A->SetFixedBlockSize(A0->GetFixedBlockSize()); + } - if(bIsBlocked) { - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::nonnull(thyraPrecOp)); - - typedef MueLu::XpetraOperator MueXpOp; - //typedef Thyra::XpetraLinearOp ThyXpLinOp; // unused - const RCP muelu_xpetraOp = rcp(new MueXpOp(H)); - - RCP > thyraRangeSpace = Xpetra::ThyraUtils::toThyra(muelu_xpetraOp->getRangeMap()); - RCP > thyraDomainSpace = Xpetra::ThyraUtils::toThyra(muelu_xpetraOp->getDomainMap()); + // set new matrix + level0->Set("A", A); - RCP > xpOp = Teuchos::rcp_dynamic_cast >(muelu_xpetraOp); - thyraPrecOp = Thyra::xpetraLinearOp(thyraRangeSpace, thyraDomainSpace,xpOp); + H->SetupRe(); + } } + // wrap preconditioner in thyraPrecOp + RCP > thyraRangeSpace = Xpetra::ThyraUtils::toThyra(xpPrecOp->getRangeMap()); + RCP > thyraDomainSpace = Xpetra::ThyraUtils::toThyra(xpPrecOp->getDomainMap()); + + RCP thyraPrecOp = Thyra::xpetraLinearOp(thyraRangeSpace, thyraDomainSpace, xpPrecOp); TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp)); defaultPrec->initializeUnspecified(thyraPrecOp); diff --git a/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_decl.hpp b/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_decl.hpp index aaf801cbb003..3af40c908e52 100644 --- a/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_decl.hpp +++ b/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_decl.hpp @@ -56,6 +56,7 @@ #include "Thyra_BlockedLinearOpBase.hpp" #include "Thyra_DiagonalLinearOpBase.hpp" #include "Thyra_XpetraLinearOp.hpp" +#include #ifdef HAVE_MUELU_TPETRA #include "Thyra_TpetraLinearOp.hpp" #include "Thyra_TpetraThyraWrappers.hpp" @@ -95,9 +96,6 @@ namespace Thyra { - template - bool replaceWithXpetra(ParameterList& paramList, std::string parameterName); - /** @brief Concrete preconditioner factory subclass for Thyra based on MueLu. @ingroup MueLuAdapters Add support for MueLu preconditioners in Thyra. This class provides an interface both @@ -283,7 +281,7 @@ namespace Thyra { // rebuild preconditioner if startingOver == true // reuse preconditioner if startingOver == false const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("refmaxwell: enable reuse") || !paramList.get("refmaxwell: enable reuse")); - const bool useHalfPrecision = paramList.get("refmaxwell: half precision", false) && bIsTpetra; + const bool useHalfPrecision = paramList.get("half precision", false) && bIsTpetra; RCP xpPrecOp; if (startingOver == true) { diff --git a/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_def.hpp b/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_def.hpp index df6aba25f482..a098d8065b99 100644 --- a/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_def.hpp +++ b/packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_def.hpp @@ -59,173 +59,6 @@ namespace Thyra { using Teuchos::rcp_dynamic_cast; using Teuchos::rcp_const_cast; - template - bool replaceWithXpetra(ParameterList& paramList, std::string parameterName) { - typedef typename Teuchos::ScalarTraits::magnitudeType Magnitude; - typedef Xpetra::Operator XpOp; - typedef Xpetra::ThyraUtils XpThyUtils; - typedef Xpetra::CrsMatrixWrap XpCrsMatWrap; - typedef Xpetra::CrsMatrix XpCrsMat; - typedef Xpetra::Matrix XpMat; - typedef Xpetra::MultiVector XpMultVec; - typedef Xpetra::MultiVector XpMagMultVec; - typedef Xpetra::Vector XpVec; - - typedef Thyra::LinearOpBase ThyLinOpBase; - typedef Thyra::DiagonalLinearOpBase ThyDiagLinOpBase; - typedef Thyra::XpetraLinearOp ThyXpOp; - typedef Thyra::SpmdVectorSpaceBase ThyVSBase; - -#ifdef HAVE_MUELU_TPETRA - typedef Tpetra::CrsMatrix TpCrsMat; - typedef Tpetra::Vector tV; - typedef Thyra::TpetraVector thyTpV; - typedef Tpetra::MultiVector tMV; - typedef Tpetra::MultiVector tMagMV; -# if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) - typedef typename Teuchos::ScalarTraits::halfPrecision HalfMagnitude; - typedef Tpetra::MultiVector tHalfMagMV; -# endif -#endif -#if defined(HAVE_MUELU_EPETRA) - typedef Xpetra::EpetraCrsMatrixT XpEpCrsMat; -#endif - - if (paramList.isParameter(parameterName)) { - if (paramList.isType >(parameterName)) - return true; - else if (paramList.isType >(parameterName)) { - RCP constM = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP M = rcp_const_cast(constM); - paramList.set >(parameterName, M); - return true; - } - else if (paramList.isType >(parameterName)) - return true; - else if (paramList.isType >(parameterName)) { - RCP constX = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP X = rcp_const_cast(constX); - paramList.set >(parameterName, X); - return true; - } - else if (paramList.isType >(parameterName)) - return true; - else if (paramList.isType >(parameterName)) { - RCP constX = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP X = rcp_const_cast(constX); - paramList.set >(parameterName, X); - return true; - } -#ifdef HAVE_MUELU_TPETRA - else if (paramList.isType >(parameterName)) { - RCP tM = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP xM = rcp_dynamic_cast(tM, true); - paramList.set >(parameterName, xM); - return true; - } else if (paramList.isType >(parameterName)) { - RCP tpetra_X = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP X = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_X); - paramList.set >(parameterName, X); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X)); - return true; - } else if (paramList.isType >(parameterName)) { - RCP tpetra_X = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP X = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_X); - paramList.set >(parameterName, X); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X)); - return true; - } -# if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) - else if (paramList.isType >(parameterName)) { - RCP tpetra_hX = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP tpetra_X = rcp(new tMagMV(tpetra_hX->getMap(),tpetra_hX->getNumVectors())); - Tpetra::deep_copy(*tpetra_X,*tpetra_hX); - RCP X = MueLu::TpetraMultiVector_To_XpetraMultiVector(tpetra_X); - paramList.set >(parameterName, X); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X)); - return true; - } -# endif -#endif -#ifdef HAVE_MUELU_EPETRA - else if (paramList.isType >(parameterName)) { - RCP eM = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP xeM = rcp(new XpEpCrsMat(eM)); - RCP xCrsM = rcp_dynamic_cast(xeM, true); - RCP xwM = rcp(new XpCrsMatWrap(xCrsM)); - RCP xM = rcp_dynamic_cast(xwM); - paramList.set >(parameterName, xM); - return true; - } else if (paramList.isType >(parameterName)) { - RCP epetra_X = Teuchos::null; - epetra_X = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP > xpEpX = rcp(new Xpetra::EpetraMultiVectorT(epetra_X)); - RCP > xpEpXMult = rcp_dynamic_cast >(xpEpX, true); - RCP X = rcp_dynamic_cast(xpEpXMult, true); - paramList.set >(parameterName, X); - return true; - } -#endif - else if (paramList.isType >(parameterName)) { - RCP thyM = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP crsM = XpThyUtils::toXpetra(thyM); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM)); - // MueLu needs a non-const object as input - RCP crsMNonConst = rcp_const_cast(crsM); - TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMNonConst)); - // wrap as an Xpetra::Matrix that MueLu can work with - RCP M = rcp(new Xpetra::CrsMatrixWrap(crsMNonConst)); - paramList.set >(parameterName, M); - return true; - } else if (paramList.isType >(parameterName)) { - RCP thyM = paramList.get >(parameterName); - paramList.remove(parameterName); - RCP > diag = thyM->getDiag(); - - RCP xpDiag; -#ifdef HAVE_MUELU_TPETRA - if (!rcp_dynamic_cast(diag).is_null()) { - RCP tDiag = Thyra::TpetraOperatorVectorExtraction::getConstTpetraVector(diag); - if (!tDiag.is_null()) - xpDiag = Xpetra::toXpetra(tDiag); - } -#endif -#ifdef HAVE_MUELU_EPETRA - if (xpDiag.is_null()) { - RCP comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast(thyM->range())->getComm()); - RCP map = Thyra::get_Epetra_Map(*(thyM->range()), comm); - if (!map.is_null()) { - RCP eDiag = Thyra::get_Epetra_Vector(*map, diag); - RCP nceDiag = rcp_const_cast(eDiag); - RCP > xpEpDiag = rcp(new Xpetra::EpetraVectorT(nceDiag)); - xpDiag = rcp_dynamic_cast(xpEpDiag, true); - } - } -#endif - TEUCHOS_ASSERT(!xpDiag.is_null()); - RCP M = Xpetra::MatrixFactory::Build(xpDiag); - paramList.set >(parameterName, M); - return true; - } - else { - TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter " << parameterName << " has wrong type."); - return false; - } - } else - return false; - } - - // Constructors/initializers/accessors template @@ -317,7 +150,7 @@ namespace Thyra { // rebuild preconditioner if startingOver == true // reuse preconditioner if startingOver == false const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("refmaxwell: enable reuse") || !paramList.get("refmaxwell: enable reuse")); - const bool useHalfPrecision = paramList.get("refmaxwell: half precision", false) && bIsTpetra; + const bool useHalfPrecision = paramList.get("half precision", false) && bIsTpetra; RCP xpPrecOp; if (startingOver == true) { diff --git a/packages/muelu/doc/UsersGuide/masterList.xml b/packages/muelu/doc/UsersGuide/masterList.xml index efdcf0a2371e..0084a47ec00d 100644 --- a/packages/muelu/doc/UsersGuide/masterList.xml +++ b/packages/muelu/doc/UsersGuide/masterList.xml @@ -146,6 +146,15 @@ Pass parameters to the underlying linear algebra not supported by ML false + + + + half precision + bool + Build half precision preconditioner + not supported by ML + false + false @@ -942,7 +951,7 @@ both (``both''). not supported by ML - + diff --git a/packages/muelu/doc/UsersGuide/paramlist_hidden.tex b/packages/muelu/doc/UsersGuide/paramlist_hidden.tex index 5b5c5c8be0bb..03e394f9e8e7 100644 --- a/packages/muelu/doc/UsersGuide/paramlist_hidden.tex +++ b/packages/muelu/doc/UsersGuide/paramlist_hidden.tex @@ -32,6 +32,8 @@ \cba{matvec params}{\parameterlist}{Pass parameters to the underlying linear algebra} +\cbb{half precision}{bool}{false}{Build half precision preconditioner} + \cbb{smoother: pre or post}{string}{"both"}{Pre- and post-smoother combination. Possible values: "pre" (only pre-smoother), "post" (only post-smoother), "both" (both pre-and post-smoothers), "none" (no smoothing).} \cbb{smoother: type}{string}{"RELAXATION"}{Smoother type. Possible values: see Table~\ref{tab:smoothers}.} diff --git a/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp b/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp index e9211c47efae..5599bc480fe2 100644 --- a/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp +++ b/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp @@ -1301,6 +1301,8 @@ namespace MueLu { oss << "\n--------------------------------------------------------------------------------\n"; oss << "--- Multigrid Summary " << std::setw(28) << std::left << label << "---\n"; oss << "--------------------------------------------------------------------------------" << std::endl; + if (verbLevel & Parameters1) + oss << "Scalar = " << Teuchos::ScalarTraits::name() << std::endl; oss << "Number of levels = " << numLevels << std::endl; oss << "Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed) << GetOperatorComplexity() << std::endl; diff --git a/packages/muelu/src/MueCentral/MueLu_MasterList.cpp b/packages/muelu/src/MueCentral/MueLu_MasterList.cpp index f1bdcba59603..61ba5bdf1062 100644 --- a/packages/muelu/src/MueCentral/MueLu_MasterList.cpp +++ b/packages/muelu/src/MueCentral/MueLu_MasterList.cpp @@ -184,6 +184,7 @@ namespace MueLu { "" "" "" + "" "" "" "" @@ -570,6 +571,8 @@ namespace MueLu { ("matvec params","matvec params") + ("half precision","half precision") + ("smoother: pre or post","smoother: pre or post") ("smoother: type","smoother: type")