10 #ifndef THYRA_TPETRA_LINEAR_OP_HPP
11 #define THYRA_TPETRA_LINEAR_OP_HPP
13 #include "Thyra_TpetraLinearOp_decl.hpp"
14 #include "Kokkos_Core.hpp"
15 #include "Thyra_TpetraVectorSpace.hpp"
16 #include "Teuchos_ScalarTraits.hpp"
17 #include "Teuchos_TypeNameTraits.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>
74 convertConjNoTransToTeuchosTransMode()
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>
118 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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>
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;
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()
208 *epetraOp = Teuchos::null;
213 #endif // HAVE_THYRA_TPETRA_EPETRA
219 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
229 if (M_trans ==
CONJ) {
235 return tpetraOperator_->hasTransposeApply();
239 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
248 using Teuchos::rcpFromRef;
249 using Teuchos::rcpFromPtr;
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;
303 rowMatrix->leftScale(*row_scaling);
308 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
313 using Teuchos::rcpFromRef;
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,
356 typedef typename STS::magnitudeType MT;
359 if ( (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) ||
360 (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM) ) {
381 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tCrsMatrix->getRangeMap()));
382 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tRowSumVec->getMap()));
384 size_t numMyRows = tCrsMatrix->getLocalNumRows();
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>
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)
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)