46 #ifndef MUELU_XPETRA_THYRALINEAROP_HPP
47 #define MUELU_XPETRA_THYRALINEAROP_HPP
49 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
57 #include <Thyra_VectorBase.hpp>
58 #include <Thyra_SolveSupportTypes.hpp>
59 #include <Thyra_LinearOpWithSolveBase.hpp>
60 #include <Teuchos_AbstractFactoryStd.hpp>
61 #include <Stratimikos_LinearSolverBuilder.hpp>
72 class XpetraThyraLinearOp :
public Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
74 XpetraThyraLinearOp() =
default;
82 RCP<ParameterList> params)
84 throw Exceptions::RuntimeError(
"Interface not supported");
88 ~XpetraThyraLinearOp() =
default;
94 throw Exceptions::RuntimeError(
"Interface not supported");
99 throw Exceptions::RuntimeError(
"Interface not supported");
112 throw Exceptions::RuntimeError(
"Interface not supported");
116 bool hasTransposeApply()
const {
return false; }
122 throw Exceptions::RuntimeError(
"Interface not supported");
126 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A_;
138 XpetraThyraLinearOp() =
default;
146 RCP<ParameterList> params)
151 Stratimikos::LinearSolverBuilder<Scalar> linearSolverBuilder;
152 typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
153 typedef Thyra::MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> ImplMueLu;
154 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, ImplMueLu>(),
"MueLu");
156 linearSolverBuilder.setParameterList(params);
161 auto precFactory = Thyra::createPreconditioningStrategy(linearSolverBuilder);
162 auto prec = precFactory->createPrec();
164 precFactory->initializePrec(Thyra::defaultLinearOpSource(thyraA), prec.get(), Thyra::SUPPORT_SOLVE_UNSPECIFIED);
169 ~XpetraThyraLinearOp() =
default;
175 return A_->getDomainMap();
180 return A_->getRangeMap();
193 RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> rcpX = Teuchos::rcpFromRef(X);
194 RCP<const Thyra::MultiVectorBase<Scalar>> thyraX = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(rcpX);
196 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> rcpY = Teuchos::rcpFromRef(Y);
197 RCP<Thyra::MultiVectorBase<Scalar>> thyraY = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar>>(Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(rcpY));
199 prec_->getUnspecifiedPrecOp()->apply(Thyra::NOTRANS, *thyraX, thyraY.ptr(), alpha, beta);
200 Y = *Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(thyraY, Y.
getMap()->getComm());
204 bool hasTransposeApply()
const {
return false; }
211 R.
update(STS::one(), B, STS::zero());
216 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A_;
217 RCP<Thyra::PreconditionerBase<Scalar>> prec_;
222 #endif // defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
224 #endif // MUELU_XPETRA_THYRALINEAROP_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
virtual void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)=0
MueLu::DefaultScalar Scalar
Tpetra::Details::DefaultTypes::scalar_type DefaultScalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)