MueLu Version of the Day
Loading...
Searching...
No Matches
Thyra_XpetraLinearOp_def.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// Tobias Wiesner (tawiesn@sandia.gov)
43//
44// ***********************************************************************
45//
46// @HEADER
47#ifndef THYRA_XPETRA_LINEAR_OP_HPP
48#define THYRA_XPETRA_LINEAR_OP_HPP
49
51#include "Teuchos_ScalarTraits.hpp"
52#include "Teuchos_TypeNameTraits.hpp"
53
55#include "Xpetra_MapExtractor.hpp"
56
57namespace Thyra {
58
59
60// Constructors/initializers
61
62
63template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65{}
66
67
68template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70 const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
71 const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
72 const RCP<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &xpetraOperator
73 )
74{
75 initializeImpl(rangeSpace, domainSpace, xpetraOperator);
76}
77
78
79template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
81 const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
82 const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
83 const RCP<const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &xpetraOperator
84 )
85{
86 initializeImpl(rangeSpace, domainSpace, xpetraOperator);
87}
88
89
90template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
91RCP<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
93{
94 return xpetraOperator_.getNonconstObj();
95}
96
97
98template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
99RCP<const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
101{
102 return xpetraOperator_;
103}
104
105
106// Public Overridden functions from LinearOpBase
107
108
109template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
110RCP<const Thyra::VectorSpaceBase<Scalar> >
112{
113 return rangeSpace_;
114}
115
116
117template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118RCP<const Thyra::VectorSpaceBase<Scalar> >
120{
121 return domainSpace_;
122}
123
124// Protected Overridden functions from LinearOpBase
125
126template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
128 Thyra::EOpTransp M_trans) const
129{
130 if (is_null(xpetraOperator_))
131 return false;
132
133 if (M_trans == NOTRANS)
134 return true;
135
136 if (M_trans == CONJ) {
137 // For non-complex scalars, CONJ is always supported since it is equivalent to NO_TRANS.
138 // For complex scalars, Xpetra does not support conjugation without transposition.
139 return !Teuchos::ScalarTraits<Scalar>::isComplex;
140 }
141
142 return xpetraOperator_->hasTransposeApply();
143}
144
145
146template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
148 const Thyra::EOpTransp M_trans,
149 const Thyra::MultiVectorBase<Scalar> &X_in,
150 const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
151 const Scalar alpha,
152 const Scalar beta
153 ) const
154{
155 using Teuchos::rcpFromRef;
156 using Teuchos::rcpFromPtr;
157
158 TEUCHOS_TEST_FOR_EXCEPTION(getConstXpetraOperator() == Teuchos::null, MueLu::Exceptions::RuntimeError, "XpetraLinearOp::applyImpl: internal Xpetra::Operator is null.");
159 RCP< const Teuchos::Comm< int > > comm = getConstXpetraOperator()->getRangeMap()->getComm();
160
161 const RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tX_in =
162 Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(rcpFromRef(X_in), comm);
163 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tY_inout =
164 Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(rcpFromPtr(Y_inout), comm);
165 Teuchos::ETransp transp;
166 switch (M_trans) {
167 case NOTRANS: transp = Teuchos::NO_TRANS; break;
168 case TRANS: transp = Teuchos::TRANS; break;
169 case CONJTRANS: transp = Teuchos::CONJ_TRANS; break;
170 default: TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::NotImplemented, "Thyra::XpetraLinearOp::apply. Unknown value for M_trans. Only NOTRANS, TRANS and CONJTRANS are supported.");
171 }
172
173 xpetraOperator_->apply(*tX_in, *tY_inout, transp, alpha, beta);
174
175 // check whether Y is a product vector
176 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal,Node> > rgMapExtractor = Teuchos::null;
177 Teuchos::Ptr<Thyra::ProductMultiVectorBase<Scalar> > prodY_inout =
178 Teuchos::ptr_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(Y_inout);
179 if(prodY_inout != Teuchos::null) {
180 // If Y is a product vector we split up the data from tY and merge them
181 // into the product vector. The necessary Xpetra::MapExtractor is extracted
182 // from the fine level operator (not this!)
183
184 // get underlying fine level operator (BlockedCrsMatrix)
185 // to extract the range MapExtractor
186 RCP<MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal,Node> > mueXop =
187 Teuchos::rcp_dynamic_cast<MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal,Node> >(xpetraOperator_.getNonconstObj());
188
189 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal,Node> > A =
190 mueXop->GetHierarchy()->GetLevel(0)->template Get<RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal,Node> > >("A");
191 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
192
193 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal,Node> > bA =
194 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal,Node> >(A);
195 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bA));
196
197 rgMapExtractor = bA->getRangeMapExtractor();
198 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgMapExtractor));
199 }
200}
201
202
203// private
204
205
206template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
207template<class XpetraOperator_t>
209 const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
210 const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
211 const RCP<XpetraOperator_t> &xpetraOperator
212 )
213{
214#ifdef THYRA_DEBUG
215 TEUCHOS_ASSERT(nonnull(rangeSpace));
216 TEUCHOS_ASSERT(nonnull(domainSpace));
217 TEUCHOS_ASSERT(nonnull(xpetraOperator));
218#endif
219 rangeSpace_ = rangeSpace;
220 domainSpace_ = domainSpace;
221 xpetraOperator_ = xpetraOperator;
222}
223
224
225} // namespace Thyra
226
227
228#endif // THYRA_XPETRA_LINEAR_OP_HPP
MueLu::DefaultScalar Scalar
Exception throws when you call an unimplemented method of MueLu.
Exception throws to report errors in the internal logical of the program.
void applyImpl(const Thyra::EOpTransp M_trans, const Thyra::MultiVectorBase< Scalar > &X_in, const Teuchos::Ptr< Thyra::MultiVectorBase< Scalar > > &Y_inout, const Scalar alpha, const Scalar beta) const
void constInitialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &xpetraOperator)
Initialize.
RCP< const Thyra::VectorSpaceBase< Scalar > > range() const
RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
RCP< const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstXpetraOperator() const
Get embedded const Xpetra::Operator.
void initialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &xpetraOperator)
Initialize.
RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getXpetraOperator()
Get embedded non-const Xpetra::Operator.
void initializeImpl(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< XpetraOperator_t > &xpetraOperator)
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
XpetraLinearOp()
Construct to uninitialized.