10 #ifndef THYRA_TPETRA_LINEAR_OP_HPP
11 #define THYRA_TPETRA_LINEAR_OP_HPP
18 #include "Tpetra_CrsMatrix.hpp"
20 #ifdef HAVE_THYRA_TPETRA_EPETRA
27 #ifdef HAVE_THYRA_TPETRA_EPETRA
33 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal>
34 class GetTpetraEpetraRowMatrixWrapper {
36 template<
class TpetraMatrixType>
38 RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
39 get(
const RCP<TpetraMatrixType> &tpetraMatrix)
49 class GetTpetraEpetraRowMatrixWrapper<double, int, int> {
51 template<
class TpetraMatrixType>
53 RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
54 get(
const RCP<TpetraMatrixType> &tpetraMatrix)
57 new Tpetra::EpetraRowMatrix<TpetraMatrixType>(tpetraMatrix,
67 #endif // HAVE_THYRA_TPETRA_EPETRA
70 template <
class Scalar>
77 Exceptions::OpNotSupported,
79 ", Tpetra does not support conjugation without transposition."
85 template <
class Scalar>
92 case CONJ:
return convertConjNoTransToTeuchosTransMode<Scalar>();
105 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 const RCP<
const VectorSpaceBase<Scalar> > &rangeSpace,
113 const RCP<
const VectorSpaceBase<Scalar> > &domainSpace,
114 const RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
117 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
121 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
123 const RCP<
const VectorSpaceBase<Scalar> > &rangeSpace,
124 const RCP<
const VectorSpaceBase<Scalar> > &domainSpace,
125 const RCP<
const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
128 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
136 return tpetraOperator_.getNonconstObj();
140 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
144 return tpetraOperator_;
151 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
159 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
170 #ifdef HAVE_THYRA_TPETRA_EPETRA
173 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
176 const Ptr<EOpTransp> &epetraOpTransp,
177 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
178 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
185 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
186 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getEpetraOpView(
187 const Ptr<RCP<const Epetra_Operator> > &epetraOp,
188 const Ptr<EOpTransp> &epetraOpTransp,
189 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
190 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
193 using Teuchos::rcp_dynamic_cast;
194 typedef Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpetraRowMatrix_t;
195 if (
nonnull(tpetraOperator_)) {
197 epetraOp_ = GetTpetraEpetraRowMatrixWrapper<Scalar,LocalOrdinal,GlobalOrdinal>::get(
198 rcp_dynamic_cast<const TpetraRowMatrix_t>(tpetraOperator_.getConstObj(),
true));
200 *epetraOp = epetraOp_;
203 *epetraOpAdjointSupport = ( tpetraOperator_->hasTransposeApply()
212 #endif // HAVE_THYRA_TPETRA_EPETRA
218 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
220 Thyra::EOpTransp M_trans)
const
228 if (M_trans == CONJ) {
234 return tpetraOperator_->hasTransposeApply();
238 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
240 const Thyra::EOpTransp M_trans,
241 const Thyra::MultiVectorBase<Scalar> &X_in,
242 const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
247 using Teuchos::rcpFromRef;
248 using Teuchos::rcpFromPtr;
251 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
257 ConverterT::getConstTpetraMultiVector(rcpFromRef(X_in));
260 ConverterT::getTpetraMultiVector(rcpFromPtr(Y_inout));
262 const Teuchos::ETransp tTransp = convertToTeuchosTransMode<Scalar>(M_trans);
266 tpetraOperator_->apply(*tX, *tY, tTransp, alpha, beta);
275 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
282 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
289 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
294 using Teuchos::rcpFromRef;
300 Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),
true);
302 rowMatrix->leftScale(*row_scaling);
307 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
312 using Teuchos::rcpFromRef;
318 Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),
true);
320 rowMatrix->rightScale(*col_scaling);
327 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
330 const RowStatLinearOpBaseUtils::ERowStat rowStat)
const
336 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
337 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
346 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
348 const RowStatLinearOpBaseUtils::ERowStat rowStat,
349 const Ptr<VectorBase<Scalar> > &rowStatVec_in
352 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
355 typedef typename STS::magnitudeType MT;
358 if ( (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) ||
359 (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM) ) {
375 Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getConstObj(),
true);
380 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tCrsMatrix->getRangeMap()));
381 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tRowSumVec->getMap()));
383 size_t numMyRows = tCrsMatrix->getLocalNumRows();
385 using crs_t = Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
386 typename crs_t::local_inds_host_view_type indices;
387 typename crs_t::values_host_view_type values;
390 for (
size_t row=0; row < numMyRows; ++row) {
391 MT sum = STM::zero ();
392 tCrsMatrix->getLocalRowView (row, indices, values);
394 for (
int col = 0; col < (int) values.size(); ++col) {
395 sum += STS::magnitude (values[col]);
398 if (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) {
399 if (sum < STM::sfmin ()) {
401 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Inverse row sum "
402 <<
"requested for a matrix where one of the rows has a zero row sum!");
403 sum = STM::one () / STM::sfmin ();
406 sum = STM::one () / sum;
410 tRowSumVec->replaceLocalValue (row, Scalar (sum));
416 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Column sum support not implemented!");
425 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
426 template<
class TpetraOperator_t>
428 const RCP<
const VectorSpaceBase<Scalar> > &rangeSpace,
429 const RCP<
const VectorSpaceBase<Scalar> > &domainSpace,
439 rangeSpace_ = rangeSpace;
440 domainSpace_ = domainSpace;
441 tpetraOperator_ = tpetraOperator;
448 #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)