46#ifndef MUELU_XPETRA_THYRALINEAROP_HPP
47#define MUELU_XPETRA_THYRALINEAROP_HPP
49#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
53#include <Xpetra_Operator.hpp>
54#include <Xpetra_MultiVector.hpp>
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>
63# ifdef HAVE_MUELU_IFPACK2
64# include <Thyra_Ifpack2PreconditionerFactory.hpp>
76 class XpetraThyraLinearOp :
public Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
78 XpetraThyraLinearOp() =
default;
85 XpetraThyraLinearOp(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
86 RCP<ParameterList> params) : A_(A) {
87 throw Exceptions::RuntimeError(
"Interface not supported");
91 ~XpetraThyraLinearOp() =
default;
96 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap()
const {
97 throw Exceptions::RuntimeError(
"Interface not supported");
101 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap()
const {
102 throw Exceptions::RuntimeError(
"Interface not supported");
110 void apply(
const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
111 Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
112 Teuchos::ETransp mode = Teuchos::NO_TRANS,
113 Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
114 Scalar beta = Teuchos::ScalarTraits<Scalar>::one())
const {
115 throw Exceptions::RuntimeError(
"Interface not supported");
119 bool hasTransposeApply()
const {
return false; }
122 void residual(
const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & X,
123 const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & B,
124 Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & R)
const {
125 throw Exceptions::RuntimeError(
"Interface not supported");
130 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A_;
144 XpetraThyraLinearOp() =
default;
151 XpetraThyraLinearOp(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
152 RCP<ParameterList> params) : A_(A) {
154 RCP<const Thyra::LinearOpBase<Scalar> > thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(A)->getCrsMatrix());
156 Stratimikos::LinearSolverBuilder<Scalar> linearSolverBuilder;
157 typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
158 typedef Thyra::MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> ImplMueLu;
159 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, ImplMueLu>(),
"MueLu");
160#ifdef HAVE_MUELU_IFPACK2
162 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Impl;
163 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(),
"Ifpack2");
166 linearSolverBuilder.setParameterList(params);
171 auto precFactory = Thyra::createPreconditioningStrategy(linearSolverBuilder);
172 auto prec = precFactory->createPrec();
174 precFactory->initializePrec(Thyra::defaultLinearOpSource(thyraA), prec.get(), Thyra::SUPPORT_SOLVE_UNSPECIFIED);
179 ~XpetraThyraLinearOp() =
default;
184 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap()
const {
185 return A_->getDomainMap();
189 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap()
const {
190 return A_->getRangeMap();
198 void apply(
const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
199 Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
200 Teuchos::ETransp mode = Teuchos::NO_TRANS,
201 Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
202 Scalar beta = Teuchos::ScalarTraits<Scalar>::one())
const {
204 RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rcpX = Teuchos::rcpFromRef(X);
205 RCP<const Thyra::MultiVectorBase<Scalar> > thyraX = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpX);
207 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rcpY = Teuchos::rcpFromRef(Y);
208 RCP<Thyra::MultiVectorBase<Scalar> > thyraY = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpY));
210 prec_->getUnspecifiedPrecOp()->apply(Thyra::NOTRANS, *thyraX, thyraY.ptr(), alpha, beta);
211 Y = *Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraY, Y.getMap()->getComm());
215 bool hasTransposeApply()
const {
return false; }
218 void residual(
const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & X,
219 const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & B,
220 Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & R)
const {
221 using STS = Teuchos::ScalarTraits<Scalar>;
222 R.update(STS::one(),B,STS::zero());
223 this->apply (X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
228 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A_;
229 RCP<Thyra::PreconditionerBase<Scalar> > prec_;
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Namespace for MueLu classes and methods.
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
Tpetra::Details::DefaultTypes::scalar_type DefaultScalar