10 #ifndef THYRA_TPETRA_LINEAR_OP_HPP
11 #define THYRA_TPETRA_LINEAR_OP_HPP
14 #include "Kokkos_Core.hpp"
19 #include "Tpetra_CrsMatrix.hpp"
21 #ifdef HAVE_THYRA_TPETRA_EPETRA
28 #ifdef HAVE_THYRA_TPETRA_EPETRA
34 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal>
35 class GetTpetraEpetraRowMatrixWrapper {
37 template<
class TpetraMatrixType>
39 RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
40 get(
const RCP<TpetraMatrixType> &tpetraMatrix)
50 class GetTpetraEpetraRowMatrixWrapper<double, int, int> {
52 template<
class TpetraMatrixType>
54 RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
55 get(
const RCP<TpetraMatrixType> &tpetraMatrix)
58 new Tpetra::EpetraRowMatrix<TpetraMatrixType>(tpetraMatrix,
68 #endif // HAVE_THYRA_TPETRA_EPETRA
71 template <
class Scalar>
78 Exceptions::OpNotSupported,
80 ", Tpetra does not support conjugation without transposition."
86 template <
class Scalar>
93 case CONJ:
return convertConjNoTransToTeuchosTransMode<Scalar>();
106 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 const RCP<
const VectorSpaceBase<Scalar> > &rangeSpace,
114 const RCP<
const VectorSpaceBase<Scalar> > &domainSpace,
115 const RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
118 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
124 const RCP<
const VectorSpaceBase<Scalar> > &rangeSpace,
125 const RCP<
const VectorSpaceBase<Scalar> > &domainSpace,
126 const RCP<
const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
129 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
133 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
137 return tpetraOperator_.getNonconstObj();
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
145 return tpetraOperator_;
152 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
160 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 #ifdef HAVE_THYRA_TPETRA_EPETRA
174 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
177 const Ptr<EOpTransp> &epetraOpTransp,
178 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
179 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
186 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
187 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getEpetraOpView(
188 const Ptr<RCP<const Epetra_Operator> > &epetraOp,
189 const Ptr<EOpTransp> &epetraOpTransp,
190 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
191 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
194 using Teuchos::rcp_dynamic_cast;
195 typedef Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpetraRowMatrix_t;
196 if (
nonnull(tpetraOperator_)) {
198 epetraOp_ = GetTpetraEpetraRowMatrixWrapper<Scalar,LocalOrdinal,GlobalOrdinal>::get(
199 rcp_dynamic_cast<const TpetraRowMatrix_t>(tpetraOperator_.getConstObj(),
true));
201 *epetraOp = epetraOp_;
204 *epetraOpAdjointSupport = ( tpetraOperator_->hasTransposeApply()
213 #endif // HAVE_THYRA_TPETRA_EPETRA
219 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
221 Thyra::EOpTransp M_trans)
const
229 if (M_trans == CONJ) {
235 return tpetraOperator_->hasTransposeApply();
239 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
241 const Thyra::EOpTransp M_trans,
242 const Thyra::MultiVectorBase<Scalar> &X_in,
243 const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
248 using Teuchos::rcpFromRef;
249 using Teuchos::rcpFromPtr;
252 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
258 ConverterT::getConstTpetraMultiVector(rcpFromRef(X_in));
261 ConverterT::getTpetraMultiVector(rcpFromPtr(Y_inout));
263 const Teuchos::ETransp tTransp = convertToTeuchosTransMode<Scalar>(M_trans);
267 tpetraOperator_->apply(*tX, *tY, tTransp, alpha, beta);
276 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
283 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
290 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
295 using Teuchos::rcpFromRef;
301 Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),
true);
303 rowMatrix->leftScale(*row_scaling);
308 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
313 using Teuchos::rcpFromRef;
319 Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),
true);
321 rowMatrix->rightScale(*col_scaling);
328 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
331 const RowStatLinearOpBaseUtils::ERowStat rowStat)
const
337 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
338 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
347 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
349 const RowStatLinearOpBaseUtils::ERowStat rowStat,
350 const Ptr<VectorBase<Scalar> > &rowStatVec_in
353 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
356 typedef typename STS::magnitudeType MT;
359 if ( (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) ||
360 (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM) ) {
376 Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getConstObj(),
true);
381 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tCrsMatrix->getRangeMap()));
382 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tRowSumVec->getMap()));
384 size_t numMyRows = tCrsMatrix->getLocalNumRows();
386 using crs_t = Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
387 typename crs_t::local_inds_host_view_type indices;
388 typename crs_t::values_host_view_type values;
391 for (
size_t row=0; row < numMyRows; ++row) {
392 MT sum = STM::zero ();
393 tCrsMatrix->getLocalRowView (row, indices, values);
395 for (
int col = 0; col < (int) values.size(); ++col) {
396 sum += STS::magnitude (values[col]);
399 if (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) {
400 if (sum < STM::sfmin ()) {
402 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Inverse row sum "
403 <<
"requested for a matrix where one of the rows has a zero row sum!");
404 sum = STM::one () / STM::sfmin ();
407 sum = STM::one () / sum;
411 tRowSumVec->replaceLocalValue (row, Scalar (sum));
417 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Column sum support not implemented!");
426 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
427 template<
class TpetraOperator_t>
429 const RCP<
const VectorSpaceBase<Scalar> > &rangeSpace,
430 const RCP<
const VectorSpaceBase<Scalar> > &domainSpace,
440 rangeSpace_ = rangeSpace;
441 domainSpace_ = domainSpace;
442 tpetraOperator_ = tpetraOperator;
449 #endif // THYRA_TPETRA_LINEAR_OP_HPP
bool is_null(const boost::shared_ptr< T > &p)
void initializeImpl(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< TpetraOperator_t > &tpetraOperator)
void constInitialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< const Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraOperator)
Initialize.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< Scalar > > &rowStatVec) const
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
RCP< const Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraOperator() const
Get embedded const Tpetra::Operator.
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
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual bool supportsScaleLeftImpl() const
RCP< const Thyra::VectorSpaceBase< Scalar > > range() const
RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraOperator()
Get embedded non-const Tpetra::Operator.
Teuchos::ETransp convertConjNoTransToTeuchosTransMode()
TpetraLinearOp()
Construct to uninitialized.
Teuchos::ETransp convertToTeuchosTransMode(const Thyra::EOpTransp transp)
RCP< const Epetra_Comm > get_Epetra_Comm(const Teuchos::Comm< Ordinal > &comm)
Get (or create) and Epetra_Comm given a Teuchos::Comm object.
virtual void scaleLeftImpl(const VectorBase< Scalar > &row_scaling)
bool nonnull(const boost::shared_ptr< T > &p)
virtual void scaleRightImpl(const VectorBase< Scalar > &col_scaling)
virtual bool supportsScaleRightImpl() const
void initialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraOperator)
Initialize.
RCP< const Teuchos::Comm< Ordinal > > convertTpetraToThyraComm(const RCP< const Teuchos::Comm< int > > &tpetraComm)
Given an Tpetra Teuchos::Comm<int> object, return an equivalent Teuchos::Comm<Ordinal> object...
#define TEUCHOS_ASSERT(assertion_test)
RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
Concrete Thyra::LinearOpBase subclass for Tpetra::Operator.
Apply using Epetra_Operator::Apply(...)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)