42#include "Thyra_EpetraLinearOp.hpp"
44#include "Thyra_SpmdMultiVectorBase.hpp"
45#include "Thyra_MultiVectorStdOps.hpp"
46#include "Thyra_AssertOp.hpp"
47#include "Teuchos_dyn_cast.hpp"
48#include "Teuchos_Assert.hpp"
49#include "Teuchos_getConst.hpp"
50#include "Teuchos_as.hpp"
51#include "Teuchos_TimeMonitor.hpp"
53#include "Epetra_Map.h"
54#include "Epetra_Vector.h"
55#include "Epetra_Operator.h"
56#include "Epetra_CrsMatrix.h"
66 :isFullyInitialized_(false),
83 using Teuchos::rcp_dynamic_cast;
89 "Thyra::EpetraLinearOp::initialize(...): Error!" );
94 if(!is_null(range_in))
95 l_spmdRange = rcp_dynamic_cast<const SPMDVSB>(range_in,
true);
101 if(!is_null(domain_in))
102 l_spmdDomain = rcp_dynamic_cast<const SPMDVSB>(domain_in,
true);
108 isFullyInitialized_ =
true;
110 rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(op_);
113 adjointSupport_ = adjointSupport;
114 range_ = l_spmdRange;
115 domain_ = l_spmdDomain;
130 using Teuchos::rcp_dynamic_cast;
136 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
138 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
140 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
144 l_spmdRange = rcp_dynamic_cast<const SPMDVSB>(range_in,
true);
146 l_spmdDomain = rcp_dynamic_cast<const SPMDVSB>(domain_in,
true);
149 isFullyInitialized_ =
false;
151 rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(op_);
154 adjointSupport_ = adjointSupport;
155 range_ = l_spmdRange;
156 domain_ = l_spmdDomain;
164 isFullyInitialized_ =
true;
179 if(opTrans) *opTrans = opTrans_;
180 if(applyAs) *applyAs = applyAs_;
181 if(adjointSupport) *adjointSupport = adjointSupport_;
182 if(range_out) *range_out = range_;
183 if(domain_out) *domain_out = domain_;
185 isFullyInitialized_ =
false;
187 rowMatrix_ = Teuchos::null;
191 range_ = Teuchos::null;
192 domain_ = Teuchos::null;
236 *epetraOpTransp = opTrans_;
237 *epetraOpApplyAs = applyAs_;
238 *epetraOpAdjointSupport = adjointSupport_;
250 *epetraOpTransp = opTrans_;
251 *epetraOpApplyAs = applyAs_;
252 *epetraOpAdjointSupport = adjointSupport_;
277 return Teuchos::null;
286 std::ostringstream oss;
289 oss <<
"op=\'"<<
typeName(*op_)<<
"\'";
290 oss <<
",rangeDim="<<this->
range()->dim();
291 oss <<
",domainDim="<<this->
domain()->dim();
308 using Teuchos::rcp_dynamic_cast;
310 using Teuchos::describe;
319 <<
"rangeDim=" << this->
range()->dim()
320 <<
",domainDim=" << this->
domain()->dim()
325 out <<
"opTrans="<<
toString(opTrans_)<<
"\n";
326 out <<
"applyAs="<<
toString(applyAs_)<<
"\n";
327 out <<
"adjointSupport="<<
toString(adjointSupport_)<<
"\n";
333 csr_op = rcp_dynamic_cast<const Epetra_CrsMatrix>(op_);
334 if (!is_null(csr_op)) {
340 out <<
"op=NULL"<<
"\n";
354 if (!isFullyInitialized_)
370 THYRA_FUNC_TIME_MONITOR(
"Thyra::EpetraLinearOp::euclideanApply");
377 "EpetraLinearOp::euclideanApply(...)", *
this, M_trans, X_in, Y_inout
382 "EpetraLinearOp::apply(...): *this was informed that adjoints "
383 "are not supported when initialized."
388 const int numCols = XY_domain->dim();
401 THYRA_FUNC_TIME_MONITOR_DIFF(
402 "Thyra::EpetraLinearOp::euclideanApply: Convert MultiVectors", MultiVectors);
405 real_M_trans==
NOTRANS ? getDomainMap() : getRangeMap(), X_in );
409 real_M_trans==
NOTRANS ? getRangeMap() : getDomainMap(), *Y_inout );
421 bool oldState = op_->UseTranspose();
422 op_->SetUseTranspose(
429 THYRA_FUNC_TIME_MONITOR_DIFF(
430 "Thyra::EpetraLinearOp::euclideanApply: Apply", Apply);
434 THYRA_FUNC_TIME_MONITOR_DIFF(
435 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Apply",
437 op_->Apply( *X, *Y );
440 THYRA_FUNC_TIME_MONITOR_DIFF(
441 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): ApplyInverse",
443 op_->ApplyInverse( *X, *Y );
452 THYRA_FUNC_TIME_MONITOR_DIFF(
453 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Scale Y",
461 THYRA_FUNC_TIME_MONITOR_DIFF(
462 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Scale Y",
464 scale( beta, Y_inout );
467 THYRA_FUNC_TIME_MONITOR_DIFF(
468 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Y=0",
470 assign( Y_inout, 0.0 );
478 THYRA_FUNC_TIME_MONITOR_DIFF(
479 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Apply",
484 THYRA_FUNC_TIME_MONITOR_DIFF(
485 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): ApplyInverse",
487 op_->ApplyInverse( *X, T );
496 THYRA_FUNC_TIME_MONITOR_DIFF(
497 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Update Y",
513 op_->SetUseTranspose(oldState);
526 return nonnull(rowMatrix_);
532 return nonnull(rowMatrix_);
538 using Teuchos::rcpFromRef;
541 rowMatrix_->LeftScale(*row_scaling);
547 using Teuchos::rcpFromRef;
550 rowMatrix_->RightScale(*col_scaling);
558 const RowStatLinearOpBaseUtils::ERowStat rowStat)
const
560 if (is_null(rowMatrix_)) {
564 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
565 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
575 const RowStatLinearOpBaseUtils::ERowStat rowStat,
579 using Teuchos::rcpFromPtr;
583 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
584 rowMatrix_->InvRowSums(*rowStatVec);
586 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
588 computeAbsRowSum(*rowStatVec);
605 return Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
618 return Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
627const Epetra_Map& EpetraLinearOp::getRangeMap()
const
630 ? op_->OperatorRangeMap() : op_->OperatorDomainMap() );
635const Epetra_Map& EpetraLinearOp::getDomainMap()
const
638 ? op_->OperatorDomainMap() : op_->OperatorRangeMap() );
642void EpetraLinearOp::computeAbsRowSum(
Epetra_Vector & x)
const
646 RCP<Epetra_CrsMatrix> crsMatrix
647 = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(rowMatrix_);
650 Exceptions::OpNotSupported,
651 "EpetraLinearOp::computeAbsRowSum(...): wrapped matrix must be of type "
652 "Epetra_CrsMatrix for this method. Other operator types are not supported."
660 if (crsMatrix->Filled()) {
662 std::invalid_argument,
663 "EpetraLinearOp::computeAbsRowSum(...): Epetra_CrsMatrix must be filled"
668 double * xp = (
double*)x.
Values();
669 if (crsMatrix->Graph().RangeMap().SameAs(x.
Map()) && crsMatrix->Exporter() != 0) {
671 x_tmp.PutScalar(0.0);
672 double * x_tmp_p = (
double*)x_tmp.Values();
673 for (i=0; i < crsMatrix->NumMyRows(); i++) {
675 double * RowValues = 0;
676 crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
677 for (j=0; j < NumEntries; j++) x_tmp_p[i] += std::abs(RowValues[j]);
681 else if (crsMatrix->Graph().RowMap().SameAs(x.
Map())) {
682 for (i=0; i < crsMatrix->NumMyRows(); i++) {
684 double * RowValues = 0;
685 crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
687 for (j=0; j < NumEntries; j++) scale += std::abs(RowValues[j]);
703Thyra::nonconstEpetraLinearOp()
710Thyra::partialNonconstEpetraLinearOp(
711 const RCP<
const VectorSpaceBase<double> > &range,
712 const RCP<
const VectorSpaceBase<double> > &domain,
713 const RCP<Epetra_Operator> &op,
715 EApplyEpetraOpAs applyAs,
716 EAdjointEpetraOp adjointSupport
719 RCP<EpetraLinearOp> thyraEpetraOp =
Teuchos::rcp(
new EpetraLinearOp());
720 thyraEpetraOp->partiallyInitialize(
721 range, domain,op,opTrans, applyAs, adjointSupport
723 return thyraEpetraOp;
728Thyra::nonconstEpetraLinearOp(
729 const RCP<Epetra_Operator> &op,
731 EApplyEpetraOpAs applyAs,
732 EAdjointEpetraOp adjointSupport,
733 const RCP<
const VectorSpaceBase<double> > &range,
734 const RCP<
const VectorSpaceBase<double> > &domain
737 RCP<EpetraLinearOp> thyraEpetraOp =
Teuchos::rcp(
new EpetraLinearOp());
738 thyraEpetraOp->initialize(
739 op,opTrans, applyAs, adjointSupport, range, domain
741 return thyraEpetraOp;
746Thyra::epetraLinearOp(
747 const RCP<const Epetra_Operator> &op,
749 EApplyEpetraOpAs applyAs,
750 EAdjointEpetraOp adjointSupport,
751 const RCP<
const VectorSpaceBase<double> > &range,
752 const RCP<
const VectorSpaceBase<double> > &domain
755 RCP<EpetraLinearOp> thyraEpetraOp =
Teuchos::rcp(
new EpetraLinearOp());
756 thyraEpetraOp->initialize(
757 Teuchos::rcp_const_cast<Epetra_Operator>(op),
758 opTrans, applyAs, adjointSupport, range, domain
760 return thyraEpetraOp;
765Thyra::nonconstEpetraLinearOp(
766 const RCP<Epetra_Operator> &op,
767 const std::string &label,
769 EApplyEpetraOpAs applyAs,
770 EAdjointEpetraOp adjointSupport,
771 const RCP<
const VectorSpaceBase<double> > &range,
772 const RCP<
const VectorSpaceBase<double> > &domain
775 RCP<EpetraLinearOp> thyraEpetraOp =
Teuchos::rcp(
new EpetraLinearOp());
776 thyraEpetraOp->initialize(
777 op,opTrans, applyAs, adjointSupport, range, domain
779 thyraEpetraOp->setObjectLabel(label);
780 return thyraEpetraOp;
785Thyra::epetraLinearOp(
786 const RCP<const Epetra_Operator> &op,
787 const std::string &label,
789 EApplyEpetraOpAs applyAs,
790 EAdjointEpetraOp adjointSupport,
791 const RCP<
const SpmdVectorSpaceBase<double> > &range,
792 const RCP<
const SpmdVectorSpaceBase<double> > &domain
795 RCP<EpetraLinearOp> thyraEpetraOp =
Teuchos::rcp(
new EpetraLinearOp());
796 thyraEpetraOp->initialize(
797 Teuchos::rcp_const_cast<Epetra_Operator>(op),
798 opTrans, applyAs, adjointSupport, range, domain
800 thyraEpetraOp->setObjectLabel(label);
801 return thyraEpetraOp;
const Epetra_BlockMap & Map() const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int PutScalar(double ScalarConstant)
virtual std::string description() const
RCP< const SpmdVectorSpaceBase< double > > spmdDomain() const
Return a smart pointer to the SpmdVectorSpaceBase object for the domain.
std::string description() const
virtual void scaleLeftImpl(const VectorBase< double > &row_scaling)
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< double > &X, const Ptr< MultiVectorBase< double > > &Y, const double alpha, const double beta) const
virtual void scaleRightImpl(const VectorBase< double > &col_scaling)
void getNonconstEpetraOpView(const Ptr< RCP< Epetra_Operator > > &epetraOp, const Ptr< EOpTransp > &epetraOpTransp, const Ptr< EApplyEpetraOpAs > &epetraOpApplyAs, const Ptr< EAdjointEpetraOp > &epetraOpAdjointSupport)
void getEpetraOpView(const Ptr< RCP< const Epetra_Operator > > &epetraOp, const Ptr< EOpTransp > &epetraOpTransp, const Ptr< EApplyEpetraOpAs > &epetraOpApplyAs, const Ptr< EAdjointEpetraOp > &epetraOpAdjointSupport) const
void uninitialize(RCP< Epetra_Operator > *op=NULL, EOpTransp *opTrans=NULL, EApplyEpetraOpAs *applyAs=NULL, EAdjointEpetraOp *adjointSupport=NULL, RCP< const VectorSpaceBase< double > > *range=NULL, RCP< const VectorSpaceBase< double > > *domain=NULL)
Set to uninitialized and optionally return the current state.
EpetraLinearOp()
Construct to uninitialized.
virtual bool supportsScaleRightImpl() const
void setFullyInitialized(bool isFullyInitialized=true)
Set to fully initialized.
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< double > > &rowStatVec) const
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual RCP< const SpmdVectorSpaceBase< double > > allocateRange(const RCP< Epetra_Operator > &op, EOpTransp op_trans) const
Allocate the range space of the operator.
void initialize(const RCP< Epetra_Operator > &op, EOpTransp opTrans=NOTRANS, EApplyEpetraOpAs applyAs=EPETRA_OP_APPLY_APPLY, EAdjointEpetraOp adjointSupport=EPETRA_OP_ADJOINT_SUPPORTED, const RCP< const VectorSpaceBase< double > > &range=Teuchos::null, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Fully initialize.
RCP< const LinearOpBase< double > > clone() const
void partiallyInitialize(const RCP< const VectorSpaceBase< double > > &range, const RCP< const VectorSpaceBase< double > > &domain, const RCP< Epetra_Operator > &op, EOpTransp opTrans=NOTRANS, EApplyEpetraOpAs applyAs=EPETRA_OP_APPLY_APPLY, EAdjointEpetraOp adjointSupport=EPETRA_OP_ADJOINT_SUPPORTED)
Partially initialize.
RCP< const SpmdVectorSpaceBase< double > > spmdRange() const
Return a smart pointer to the SpmdVectorSpaceBase object for the range.
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
RCP< Epetra_Operator > epetra_op()
virtual RCP< const SpmdVectorSpaceBase< double > > allocateDomain(const RCP< Epetra_Operator > &op, EOpTransp op_trans) const
Allocate the domain space of the operator.
bool opSupportedImpl(EOpTransp M_trans) const
RCP< const VectorSpaceBase< double > > range() const
virtual bool supportsScaleLeftImpl() const
RCP< const VectorSpaceBase< double > > domain() const
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Interface for a collection of column vectors called a multi-vector.
Base abstract VectorSpaceBase class for all SPMD-based vector spaces.
Abstract interface for finite-dimensional dense vectors.
Abstract interface for objects that represent a space for vectors.
EApplyEpetraOpAs
Determine how the apply an Epetra_Operator as a linear operator.
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range=Teuchos::null, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Create a non-const MultiVectorBase object from a non-const Epetra_MultiVector object.
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
Get a non-const Epetra_Vector view from a non-const VectorBase object if possible.
RCP< Epetra_MultiVector > get_Epetra_MultiVector(const Epetra_Map &map, const RCP< MultiVectorBase< double > > &mv)
Get a non-const Epetra_MultiVector view from a non-const MultiVectorBase object if possible.
EAdjointEpetraOp
Determine if adjoints are supported on Epetra_Opeator or not.
@ EPETRA_OP_APPLY_APPLY
Apply using Epetra_Operator::Apply(...)
@ EPETRA_OP_APPLY_APPLY_INVERSE
Apply using Epetra_Operator::ApplyInverse(...)
@ 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)
bool is_null(const std::shared_ptr< T > &p)
#define THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
EOpTransp trans_trans(EOpTransp trans1, EOpTransp trans2)
Combine two transpose arguments.
@ TRANS
Use the transposed operator.
@ NOTRANS
Use the non-transposed operator.
TypeTo as(const TypeFrom &t)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
std::string typeName(const T &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)