MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
53 
55 #include "Xpetra_MapExtractor.hpp"
56 
57 namespace Thyra {
58 
59 
60 // Constructors/initializers
61 
62 
63 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65 {}
66 
67 
68 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
71  const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
73  )
74 {
75  initializeImpl(rangeSpace, domainSpace, xpetraOperator);
76 }
77 
78 
79 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
81  const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
82  const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
84  )
85 {
86  initializeImpl(rangeSpace, domainSpace, xpetraOperator);
87 }
88 
89 
90 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93 {
94  return xpetraOperator_.getNonconstObj();
95 }
96 
97 
98 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
101 {
102  return xpetraOperator_;
103 }
104 
105 
106 // Public Overridden functions from LinearOpBase
107 
108 
109 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112 {
113  return rangeSpace_;
114 }
115 
116 
117 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
120 {
121  return domainSpace_;
122 }
123 
124 // Protected Overridden functions from LinearOpBase
125 
126 template <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.
140  }
141 
142  return xpetraOperator_->hasTransposeApply();
143 }
144 
145 
146 template <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 
162  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(rcpFromRef(X_in), comm);
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
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
187  Teuchos::rcp_dynamic_cast<MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal,Node> >(xpetraOperator_.getNonconstObj());
188 
190  mueXop->GetHierarchy()->GetLevel(0)->template Get<RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal,Node> > >("A");
192 
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 
206 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
207 template<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
bool is_null(const boost::shared_ptr< T > &p)
RCP< const Thyra::VectorSpaceBase< Scalar > > range() const
RCP< const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstXpetraOperator() const
Get embedded const Xpetra::Operator.
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
void initializeImpl(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< XpetraOperator_t > &xpetraOperator)
XpetraLinearOp()
Construct to uninitialized.
Exception throws when you call an unimplemented method of MueLu.
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
RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getXpetraOperator()
Get embedded non-const Xpetra::Operator.
bool nonnull(const boost::shared_ptr< T > &p)
void initialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &xpetraOperator)
Initialize.
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
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.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.