42#ifndef THYRA_TPETRA_LINEAR_OP_HPP
43#define THYRA_TPETRA_LINEAR_OP_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>
66class GetTpetraEpetraRowMatrixWrapper {
68 template<
class TpetraMatrixType>
81class GetTpetraEpetraRowMatrixWrapper<double, int, int> {
83 template<
class TpetraMatrixType>
89 new Tpetra::EpetraRowMatrix<TpetraMatrixType>(
tpetraMatrix,
102template <
class Scalar>
109 Exceptions::OpNotSupported,
111 ", Tpetra does not support conjugation without transposition."
117template <
class Scalar>
137template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
142template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
146 const RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &
tpetraOperator
153template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
157 const RCP<
const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &
tpetraOperator
164template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
168 return tpetraOperator_.getNonconstObj();
172template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
176 return tpetraOperator_;
183template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
191template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
202#ifdef HAVE_THYRA_TPETRA_EPETRA
205template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
217template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
218void 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_)) {
229 epetraOp_ = GetTpetraEpetraRowMatrixWrapper<Scalar,LocalOrdinal,GlobalOrdinal>::get(
230 rcp_dynamic_cast<const TpetraRowMatrix_t>(tpetraOperator_.getConstObj(),
true));
250template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
266 return tpetraOperator_->hasTransposeApply();
270template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
279 using Teuchos::rcpFromRef;
280 using Teuchos::rcpFromPtr;
283 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
298 tpetraOperator_->apply(*
tX, *
tY,
tTransp, alpha, beta);
305template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
312template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
319template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
324 using Teuchos::rcpFromRef;
337template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
342 using Teuchos::rcpFromRef;
357template <
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:
376template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
378 const RowStatLinearOpBaseUtils::ERowStat
rowStat,
382 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
385 typedef typename STS::magnitudeType
MT;
388 if ( (
rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) ||
389 (
rowStat == RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM) ) {
415 using crs_t = Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
416 typename crs_t::local_inds_host_view_type indices;
417 typename crs_t::values_host_view_type values;
421 MT sum = STM::zero ();
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;
446 "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Column sum support not implemented!");
455template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
456template<
class TpetraOperator_t>
RCP< T > rcpFromPtr(const Ptr< T > &ptr)
RCP< T > rcpFromRef(T &r)
RCP(T *p, const RCPNodeHandle &node)
bool nonnull(const RCP< T > &p)
Concrete Thyra::LinearOpBase subclass for Tpetra::Operator.
TpetraLinearOp()
Construct to uninitialized.
RCP< const Epetra_Comm > get_Epetra_Comm(const Teuchos::Comm< Ordinal > &comm)
Get (or create) and Epetra_Comm given a Teuchos::Comm object.
@ EPETRA_OP_APPLY_APPLY
Apply using Epetra_Operator::Apply(...)
@ EPETRA_OP_ADJOINT_UNSUPPORTED
Adjoint not supported.
@ EPETRA_OP_ADJOINT_SUPPORTED
Adjoint supported.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::ETransp convertConjNoTransToTeuchosTransMode()
Teuchos::ETransp convertToTeuchosTransMode(const Thyra::EOpTransp transp)
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.