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()
107 TEUCHOS_TEST_FOR_EXCEPTION(
108 Teuchos::ScalarTraits<Scalar>::isComplex,
109 Exceptions::OpNotSupported,
110 "For complex scalars such as " + Teuchos::TypeNameTraits<Scalar>::name() +
111 ", Tpetra does not support conjugation without transposition." 113 return Teuchos::NO_TRANS;
117 template <
class Scalar>
123 case NOTRANS:
return Teuchos::NO_TRANS;
124 case CONJ:
return convertConjNoTransToTeuchosTransMode<Scalar>();
126 case CONJTRANS:
return Teuchos::CONJ_TRANS;
130 TEUCHOS_TEST_FOR_EXCEPT(
true);
137 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
142 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
146 const RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
149 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
153 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
157 const RCP<
const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
160 initializeImpl(rangeSpace, domainSpace, tpetraOperator);
164 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
165 RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
168 return tpetraOperator_.getNonconstObj();
172 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
173 RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
176 return tpetraOperator_;
183 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
184 RCP<const Thyra::VectorSpaceBase<Scalar> >
191 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
192 RCP<const Thyra::VectorSpaceBase<Scalar> >
202 #ifdef HAVE_THYRA_TPETRA_EPETRA 205 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
207 const Ptr<RCP<Epetra_Operator> > &epetraOp,
208 const Ptr<EOpTransp> &epetraOpTransp,
209 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
210 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
213 TEUCHOS_TEST_FOR_EXCEPT(
true);
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;
226 typedef Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpetraRowMatrix_t;
227 if (nonnull(tpetraOperator_)) {
228 if (is_null(epetraOp_)) {
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>
254 if (is_null(tpetraOperator_))
260 if (M_trans ==
CONJ) {
263 return !Teuchos::ScalarTraits<Scalar>::isComplex;
266 return tpetraOperator_->hasTransposeApply();
270 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
279 using Teuchos::rcpFromRef;
280 using Teuchos::rcpFromPtr;
283 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
288 const RCP<const TpetraMultiVector_t> tX =
289 ConverterT::getConstTpetraMultiVector(rcpFromRef(X_in));
291 const RCP<TpetraMultiVector_t> tY =
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;
326 const RCP<const Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > row_scaling =
329 const RCP<typename Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rowMatrix =
330 Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),
true);
332 rowMatrix->leftScale(*row_scaling);
336 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
341 using Teuchos::rcpFromRef;
343 const RCP<const Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > col_scaling =
346 const RCP<typename Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rowMatrix =
347 Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),
true);
349 rowMatrix->rightScale(*col_scaling);
355 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
358 const RowStatLinearOpBaseUtils::ERowStat rowStat)
const 360 if (is_null(tpetraOperator_))
364 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
365 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
375 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
377 const RowStatLinearOpBaseUtils::ERowStat rowStat,
381 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
383 typedef Teuchos::ScalarTraits<Scalar> STS;
384 typedef typename STS::magnitudeType MT;
385 typedef Teuchos::ScalarTraits<MT> STM;
387 if ( (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) ||
388 (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM) ) {
390 TEUCHOS_ASSERT(nonnull(tpetraOperator_));
391 TEUCHOS_ASSERT(nonnull(rowStatVec_in));
400 const RCP<TpetraVector_t> tRowSumVec =
403 const RCP<const typename Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tCrsMatrix =
404 Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getConstObj(),
true);
409 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tCrsMatrix->getRangeMap()));
410 TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tRowSumVec->getMap()));
412 size_t numMyRows = tCrsMatrix->getNodeNumRows();
414 Teuchos::ArrayView<const LocalOrdinal> indices;
415 Teuchos::ArrayView<const Scalar> values;
417 for (
size_t row=0; row < numMyRows; ++row) {
418 MT sum = STM::zero ();
419 tCrsMatrix->getLocalRowView (row, indices, values);
421 for (
int col = 0; col < values.size(); ++col) {
422 sum += STS::magnitude (values[col]);
425 if (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) {
426 if (sum < STM::sfmin ()) {
427 TEUCHOS_TEST_FOR_EXCEPTION(sum == STM::zero (), std::runtime_error,
428 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Inverse row sum " 429 <<
"requested for a matrix where one of the rows has a zero row sum!");
430 sum = STM::one () / STM::sfmin ();
433 sum = STM::one () / sum;
437 tRowSumVec->replaceLocalValue (row, Scalar (sum));
442 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error,
443 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Column sum support not implemented!");
451 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
452 template<
class TpetraOperator_t>
456 const RCP<TpetraOperator_t> &tpetraOperator
460 TEUCHOS_ASSERT(nonnull(rangeSpace));
461 TEUCHOS_ASSERT(nonnull(domainSpace));
462 TEUCHOS_ASSERT(nonnull(tpetraOperator));
465 rangeSpace_ = rangeSpace;
466 domainSpace_ = domainSpace;
467 tpetraOperator_ = tpetraOperator;
474 #endif // THYRA_TPETRA_LINEAR_OP_HPP bool opSupportedImpl(Thyra::EOpTransp M_trans) const
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
RCP< const Thyra::VectorSpaceBase< Scalar > > range() const
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.
RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
Use the non-transposed operator.
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
virtual bool supportsScaleLeftImpl() const
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...
Use the transposed operator.
Interface for a collection of column vectors called a multi-vector.
RCP< const Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraOperator() const
Get embedded const Tpetra::Operator.
RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraOperator()
Get embedded non-const Tpetra::Operator.
Abstract interface for finite-dimensional dense vectors.
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
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...
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< Scalar > > &rowStatVec) const
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)
virtual void scaleRightImpl(const VectorBase< Scalar > &col_scaling)
void initialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraOperator)
Initialize.
virtual bool supportsScaleRightImpl() const
Concrete Thyra::LinearOpBase subclass for Tpetra::Operator.
Apply using Epetra_Operator::Apply(...)