47#ifndef THYRA_XPETRA_LINEAR_OP_HPP
48#define THYRA_XPETRA_LINEAR_OP_HPP
51#include "Teuchos_ScalarTraits.hpp"
52#include "Teuchos_TypeNameTraits.hpp"
55#include "Xpetra_MapExtractor.hpp"
63template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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
75 initializeImpl(rangeSpace, domainSpace, xpetraOperator);
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
86 initializeImpl(rangeSpace, domainSpace, xpetraOperator);
90template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
91RCP<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
94 return xpetraOperator_.getNonconstObj();
98template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
99RCP<const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
102 return xpetraOperator_;
109template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110RCP<const Thyra::VectorSpaceBase<Scalar> >
117template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
118RCP<const Thyra::VectorSpaceBase<Scalar> >
126template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
128 Thyra::EOpTransp M_trans)
const
130 if (is_null(xpetraOperator_))
133 if (M_trans == NOTRANS)
136 if (M_trans == CONJ) {
139 return !Teuchos::ScalarTraits<Scalar>::isComplex;
142 return xpetraOperator_->hasTransposeApply();
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,
155 using Teuchos::rcpFromRef;
156 using Teuchos::rcpFromPtr;
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();
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;
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.");
173 xpetraOperator_->apply(*tX_in, *tY_inout, transp, alpha, beta);
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) {
186 RCP<MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal,Node> > mueXop =
187 Teuchos::rcp_dynamic_cast<MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal,Node> >(xpetraOperator_.getNonconstObj());
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));
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));
197 rgMapExtractor = bA->getRangeMapExtractor();
198 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgMapExtractor));
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
215 TEUCHOS_ASSERT(nonnull(rangeSpace));
216 TEUCHOS_ASSERT(nonnull(domainSpace));
217 TEUCHOS_ASSERT(nonnull(xpetraOperator));
219 rangeSpace_ = rangeSpace;
220 domainSpace_ = domainSpace;
221 xpetraOperator_ = xpetraOperator;
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.