42 #ifndef THYRA_TPETRA_LINEAR_OP_HPP
43 #define THYRA_TPETRA_LINEAR_OP_HPP
45 #include "Thyra_TpetraLinearOp_decl.hpp"
46 #include "Thyra_TpetraVectorSpace.hpp"
47 #include "Teuchos_ScalarTraits.hpp"
48 #include "Teuchos_TypeNameTraits.hpp"
50 #include "Tpetra_CrsMatrix.hpp"
52 #ifdef HAVE_THYRA_TPETRA_EPETRA
59 #ifdef HAVE_THYRA_TPETRA_EPETRA
65 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal>
66 class GetTpetraEpetraRowMatrixWrapper {
68 template<
class TpetraMatrixType>
70 RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
71 get(
const RCP<TpetraMatrixType> &tpetraMatrix)
81 class GetTpetraEpetraRowMatrixWrapper<double, int, int> {
83 template<
class TpetraMatrixType>
85 RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
86 get(
const RCP<TpetraMatrixType> &tpetraMatrix)
89 new Tpetra::EpetraRowMatrix<TpetraMatrixType>(tpetraMatrix,
99 #endif // HAVE_THYRA_TPETRA_EPETRA
102 template <
class Scalar>
105 convertConjNoTransToTeuchosTransMode()
109 Exceptions::OpNotSupported,
111 ", Tpetra does not support conjugation without transposition."
117 template <
class Scalar>
124 case CONJ:
return convertConjNoTransToTeuchosTransMode<Scalar>();
137 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
142 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
149 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
153 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
160 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
164 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
168 return tpetraOperator_.getNonconstObj();
172 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
176 return tpetraOperator_;
183 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
191 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
202 #ifdef HAVE_THYRA_TPETRA_EPETRA
205 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
217 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
218 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getEpetraOpView(
219 const Ptr<RCP<const Epetra_Operator> > &epetraOp,
220 const Ptr<EOpTransp> &epetraOpTransp,
221 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
222 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
225 using Teuchos::rcp_dynamic_cast;
227 if (
nonnull(tpetraOperator_)) {
229 epetraOp_ = GetTpetraEpetraRowMatrixWrapper<Scalar,LocalOrdinal,GlobalOrdinal>::get(
230 rcp_dynamic_cast<const TpetraRowMatrix_t>(tpetraOperator_.getConstObj(),
true));
232 *epetraOp = epetraOp_;
235 *epetraOpAdjointSupport = ( tpetraOperator_->hasTransposeApply()
239 *epetraOp = Teuchos::null;
244 #endif // HAVE_THYRA_TPETRA_EPETRA
250 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
260 if (M_trans ==
CONJ) {
266 return tpetraOperator_->hasTransposeApply();
270 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
279 using Teuchos::rcpFromRef;
280 using Teuchos::rcpFromPtr;
289 ConverterT::getConstTpetraMultiVector(rcpFromRef(X_in));
292 ConverterT::getTpetraMultiVector(rcpFromPtr(Y_inout));
294 const Teuchos::ETransp tTransp = convertToTeuchosTransMode<Scalar>(M_trans);
298 tpetraOperator_->apply(*tX, *tY, tTransp, alpha, beta);
305 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
312 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
319 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
324 using Teuchos::rcpFromRef;
332 rowMatrix->leftScale(*row_scaling);
337 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
342 using Teuchos::rcpFromRef;
350 rowMatrix->rightScale(*col_scaling);
357 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
360 const RowStatLinearOpBaseUtils::ERowStat rowStat)
const
366 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
367 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
376 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
378 const RowStatLinearOpBaseUtils::ERowStat rowStat,
385 typedef typename STS::magnitudeType MT;
388 if ( (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) ||
389 (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM) ) {
410 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tCrsMatrix->getRangeMap()));
411 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tRowSumVec->getMap()));
413 size_t numMyRows = tCrsMatrix->getLocalNumRows();
416 typename crs_t::local_inds_host_view_type indices;
417 typename crs_t::values_host_view_type values;
420 for (
size_t row=0; row < numMyRows; ++row) {
421 MT sum = STM::zero ();
422 tCrsMatrix->getLocalRowView (row, indices, values);
424 for (
int col = 0; col < (int) values.size(); ++col) {
425 sum += STS::magnitude (values[col]);
428 if (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) {
429 if (sum < STM::sfmin ()) {
431 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Inverse row sum "
432 <<
"requested for a matrix where one of the rows has a zero row sum!");
433 sum = STM::one () / STM::sfmin ();
436 sum = STM::one () / sum;
440 tRowSumVec->replaceLocalValue (row, Scalar (sum));
446 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Column sum support not implemented!");
455 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
456 template<
class TpetraOperator_t>
469 rangeSpace_ = rangeSpace;
470 domainSpace_ = domainSpace;
471 tpetraOperator_ = tpetraOperator;
478 #endif // THYRA_TPETRA_LINEAR_OP_HPP
bool is_null(const boost::shared_ptr< T > &p)
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
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)
Use the non-transposed operator.
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< Scalar > > &rowStatVec) const
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
Abstract interface for objects that represent a space for vectors.
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
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
Use the transposed operator.
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Interface for a collection of column vectors called a multi-vector.
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.
Abstract interface for finite-dimensional dense vectors.
TpetraLinearOp()
Construct to uninitialized.
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...
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.
#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)