MueLu Version of the Day
Loading...
Searching...
No Matches
Xpetra_ThyraLinearOp.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_XPETRA_THYRALINEAROP_HPP
47#define MUELU_XPETRA_THYRALINEAROP_HPP
48
49#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
50
51#include "MueLu_ConfigDefs.hpp"
52
53#include <Xpetra_Operator.hpp>
54#include <Xpetra_MultiVector.hpp>
55
56// Stratimikos
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>
65# endif
66
67
68namespace MueLu {
69
72 template <class Scalar = DefaultScalar,
75 class Node = DefaultNode>
76 class XpetraThyraLinearOp : public Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
77 protected:
78 XpetraThyraLinearOp() = default;
79 public:
80
82
83
85 XpetraThyraLinearOp(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
86 RCP<ParameterList> params) : A_(A) {
87 throw Exceptions::RuntimeError("Interface not supported");
88 };
89
91 ~XpetraThyraLinearOp() = default;
92
94
96 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap() const {
97 throw Exceptions::RuntimeError("Interface not supported");
98 }
99
100 // //! Returns the Tpetra::Map object associated with the range of this operator.
101 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap() const {
102 throw Exceptions::RuntimeError("Interface not supported");
103 }
104
106
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");
116 }
117
119 bool hasTransposeApply() const { return false; }
120
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");
126 }
127
128
129 private:
130 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A_;
131 };
132
133
134 // Partial specialization for Scalar == double.
135 // Allows to avoid issues with Stokhos instantiating Thyra objects.
136 template <class LocalOrdinal,
137 class GlobalOrdinal,
138 class Node>
139 class XpetraThyraLinearOp<double,LocalOrdinal,GlobalOrdinal,Node> : public Xpetra::Operator<double,LocalOrdinal,GlobalOrdinal,Node> {
140
141 using Scalar = double;
142
143 protected:
144 XpetraThyraLinearOp() = default;
145 public:
146
148
149
151 XpetraThyraLinearOp(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
152 RCP<ParameterList> params) : A_(A) {
153 // Build Thyra linear algebra objects
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());
155
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
161 // Register Ifpack2 as a Stratimikos preconditioner strategy.
162 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Impl;
163 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
164#endif
165
166 linearSolverBuilder.setParameterList(params);
167
168 // Build a new "solver factory" according to the previously specified parameter list.
169 // RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
170
171 auto precFactory = Thyra::createPreconditioningStrategy(linearSolverBuilder);
172 auto prec = precFactory->createPrec();
173
174 precFactory->initializePrec(Thyra::defaultLinearOpSource(thyraA), prec.get(), Thyra::SUPPORT_SOLVE_UNSPECIFIED);
175 prec_ = prec;
176 };
177
179 ~XpetraThyraLinearOp() = default;
180
182
184 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap() const {
185 return A_->getDomainMap();
186 }
187
188 // //! Returns the Tpetra::Map object associated with the range of this operator.
189 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap() const {
190 return A_->getRangeMap();
191 }
192
194
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 {
203
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);
206
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));
209
210 prec_->getUnspecifiedPrecOp()->apply(Thyra::NOTRANS, *thyraX, thyraY.ptr(), alpha, beta);
211 Y = *Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraY, Y.getMap()->getComm());
212 }
213
215 bool hasTransposeApply() const { return false; }
216
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());
224 }
225
226
227 private:
228 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A_;
229 RCP<Thyra::PreconditionerBase<Scalar> > prec_;
230 };
231
232} // namespace
233
234#endif // defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
235
236#endif // MUELU_XPETRA_THYRALINEAROP_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Namespace for MueLu classes and methods.
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
Tpetra::Details::DefaultTypes::scalar_type DefaultScalar