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),
74 const RCP<Epetra_Operator> &op,
83 using Teuchos::rcp_dynamic_cast;
88 TEUCHOS_TEST_FOR_EXCEPTION( is_null(op), std::invalid_argument,
89 "Thyra::EpetraLinearOp::initialize(...): Error!" );
93 RCP<const SPMDVSB> l_spmdRange;
94 if(!is_null(range_in))
95 l_spmdRange = rcp_dynamic_cast<
const SPMDVSB>(range_in,
true);
100 RCP<const SPMDVSB> l_spmdDomain;
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;
123 const RCP<Epetra_Operator> &op,
130 using Teuchos::rcp_dynamic_cast;
135 TEUCHOS_TEST_FOR_EXCEPTION( is_null(range_in), std::invalid_argument,
136 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
137 TEUCHOS_TEST_FOR_EXCEPTION( is_null(domain_in), std::invalid_argument,
138 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
139 TEUCHOS_TEST_FOR_EXCEPTION( is_null(op), std::invalid_argument,
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;
169 RCP<Epetra_Operator> *op,
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;
197 RCP< const SpmdVectorSpaceBase<double> >
204 RCP< const SpmdVectorSpaceBase<double> >
218 RCP<const Epetra_Operator>
229 const Ptr<RCP<Epetra_Operator> > &epetraOp,
230 const Ptr<EOpTransp> &epetraOpTransp,
231 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
232 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
236 *epetraOpTransp = opTrans_;
237 *epetraOpApplyAs = applyAs_;
238 *epetraOpAdjointSupport = adjointSupport_;
243 const Ptr<RCP<const Epetra_Operator> > &epetraOp,
244 const Ptr<EOpTransp> &epetraOpTransp,
245 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
246 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
250 *epetraOpTransp = opTrans_;
251 *epetraOpApplyAs = applyAs_;
252 *epetraOpAdjointSupport = adjointSupport_;
259 RCP<const VectorSpaceBase<double> >
266 RCP<const VectorSpaceBase<double> >
273 RCP<const LinearOpBase<double> >
277 return Teuchos::null;
286 std::ostringstream oss;
287 oss << Teuchos::Describable::description() <<
"{";
289 oss <<
"op=\'"<<typeName(*op_)<<
"\'";
290 oss <<
",rangeDim="<<this->
range()->dim();
291 oss <<
",domainDim="<<this->
domain()->dim();
303 const Teuchos::EVerbosityLevel verbLevel
306 using Teuchos::includesVerbLevel;
308 using Teuchos::rcp_dynamic_cast;
309 using Teuchos::OSTab;
310 using Teuchos::describe;
312 if ( as<int>(verbLevel) == as<int>(Teuchos::VERB_LOW) || is_null(op_)) {
315 else if (includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM)) {
317 << Teuchos::Describable::description()
319 <<
"rangeDim=" << this->
range()->dim()
320 <<
",domainDim=" << this->
domain()->dim()
324 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
325 out <<
"opTrans="<<
toString(opTrans_)<<
"\n";
326 out <<
"applyAs="<<
toString(applyAs_)<<
"\n";
327 out <<
"adjointSupport="<<
toString(adjointSupport_)<<
"\n";
328 out <<
"op="<<typeName(*op_)<<
"\n";
330 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
332 RCP<const Epetra_CrsMatrix>
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");
375 TEUCHOS_TEST_FOR_EXCEPT(!isFullyInitialized_);
377 "EpetraLinearOp::euclideanApply(...)", *
this, M_trans, X_in, Y_inout
379 TEUCHOS_TEST_FOR_EXCEPTION(
382 "EpetraLinearOp::apply(...): *this was informed that adjoints " 383 "are not supported when initialized." 387 const RCP<const VectorSpaceBase<double> > XY_domain = X_in.
domain();
388 const int numCols = XY_domain->dim();
398 RCP<const Epetra_MultiVector> X;
399 RCP<Epetra_MultiVector> Y;
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 );
447 TEUCHOS_TEST_FOR_EXCEPT(
true);
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 );
473 Epetra_MultiVector T(op_->OperatorRangeMap(), numCols,
false);
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 );
491 TEUCHOS_TEST_FOR_EXCEPT(
true);
496 THYRA_FUNC_TIME_MONITOR_DIFF(
497 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Update Y",
502 Teuchos::rcp(&Teuchos::getConst(T),
false),
513 op_->SetUseTranspose(oldState);
526 return nonnull(rowMatrix_);
532 return nonnull(rowMatrix_);
538 using Teuchos::rcpFromRef;
539 const RCP<const Epetra_Vector> row_scaling =
541 rowMatrix_->LeftScale(*row_scaling);
547 using Teuchos::rcpFromRef;
548 const RCP<const Epetra_Vector> col_scaling =
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:
568 TEUCHOS_TEST_FOR_EXCEPT(
true);
575 const RowStatLinearOpBaseUtils::ERowStat rowStat,
579 using Teuchos::rcpFromPtr;
580 const RCP<Epetra_Vector> rowStatVec =
583 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
584 rowMatrix_->InvRowSums(*rowStatVec);
586 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
588 computeAbsRowSum(*rowStatVec);
591 TEUCHOS_TEST_FOR_EXCEPT(
true);
599 RCP< const SpmdVectorSpaceBase<double> >
601 const RCP<Epetra_Operator> &op,
612 RCP<const SpmdVectorSpaceBase<double> >
614 const RCP<Epetra_Operator> &op,
627 const Epetra_Map& EpetraLinearOp::getRangeMap()
const 630 ? op_->OperatorRangeMap() : op_->OperatorDomainMap() );
635 const Epetra_Map& EpetraLinearOp::getDomainMap()
const 638 ? op_->OperatorDomainMap() : op_->OperatorRangeMap() );
642 void EpetraLinearOp::computeAbsRowSum(Epetra_Vector & x)
const 644 TEUCHOS_ASSERT(!is_null(rowMatrix_));
646 RCP<Epetra_CrsMatrix> crsMatrix
647 = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(rowMatrix_);
649 TEUCHOS_TEST_FOR_EXCEPTION(is_null(crsMatrix),
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()) {
661 TEUCHOS_TEST_FOR_EXCEPTION(is_null(crsMatrix),
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) {
670 Epetra_Vector x_tmp(crsMatrix->RowMap());
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]);
679 TEUCHOS_TEST_FOR_EXCEPT(0!=x.Export(x_tmp, *crsMatrix->Exporter(), Add));
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]);
692 TEUCHOS_TEST_FOR_EXCEPT(
true);
702 Teuchos::RCP<Thyra::EpetraLinearOp>
703 Thyra::nonconstEpetraLinearOp()
705 return Teuchos::rcp(
new EpetraLinearOp());
709 Teuchos::RCP<Thyra::EpetraLinearOp>
710 Thyra::partialNonconstEpetraLinearOp(
711 const RCP<
const VectorSpaceBase<double> > &range,
712 const RCP<
const VectorSpaceBase<double> > &domain,
713 const RCP<Epetra_Operator> &op,
719 RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(
new EpetraLinearOp());
720 thyraEpetraOp->partiallyInitialize(
721 range, domain,op,opTrans, applyAs, adjointSupport
723 return thyraEpetraOp;
727 Teuchos::RCP<Thyra::EpetraLinearOp>
728 Thyra::nonconstEpetraLinearOp(
729 const RCP<Epetra_Operator> &op,
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;
745 Teuchos::RCP<const Thyra::EpetraLinearOp>
746 Thyra::epetraLinearOp(
747 const RCP<const Epetra_Operator> &op,
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;
764 Teuchos::RCP<Thyra::EpetraLinearOp>
765 Thyra::nonconstEpetraLinearOp(
766 const RCP<Epetra_Operator> &op,
767 const std::string &label,
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;
784 Teuchos::RCP<const Thyra::EpetraLinearOp>
785 Thyra::epetraLinearOp(
786 const RCP<const Epetra_Operator> &op,
787 const std::string &label,
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;
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...
virtual RCP< const SpmdVectorSpaceBase< double > > allocateRange(const RCP< Epetra_Operator > &op, EOpTransp op_trans) const
Allocate the range space of the operator.
virtual bool supportsScaleLeftImpl() const
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
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.
#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...
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< double > > &rowStatVec) const
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
Use the non-transposed operator.
RCP< const SpmdVectorSpaceBase< double > > spmdDomain() const
Return a smart pointer to the SpmdVectorSpaceBase object for the domain.
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
void getEpetraOpView(const Ptr< RCP< const Epetra_Operator > > &epetraOp, const Ptr< EOpTransp > &epetraOpTransp, const Ptr< EApplyEpetraOpAs > &epetraOpApplyAs, const Ptr< EAdjointEpetraOp > &epetraOpAdjointSupport) const
Abstract interface for objects that represent a space for vectors.
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Create a non-const MultiVectorBase object from a non-const Epetra_MultiVector object.
Use the transposed operator.
Apply using Epetra_Operator::ApplyInverse(...)
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.
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void scaleLeftImpl(const VectorBase< double > &row_scaling)
Interface for a collection of column vectors called a multi-vector.
RCP< const LinearOpBase< double > > clone() const
EAdjointEpetraOp
Determine if adjoints are supported on Epetra_Opeator or not.
EpetraLinearOp()
Construct to uninitialized.
void setFullyInitialized(bool isFullyInitialized=true)
Set to fully initialized.
EOpTransp trans_trans(EOpTransp trans1, EOpTransp trans2)
Combine two transpose arguments.
void getNonconstEpetraOpView(const Ptr< RCP< Epetra_Operator > > &epetraOp, const Ptr< EOpTransp > &epetraOpTransp, const Ptr< EApplyEpetraOpAs > &epetraOpApplyAs, const Ptr< EAdjointEpetraOp > &epetraOpAdjointSupport)
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.
Abstract interface for finite-dimensional dense vectors.
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
virtual RCP< const SpmdVectorSpaceBase< double > > allocateDomain(const RCP< Epetra_Operator > &op, EOpTransp op_trans) const
Allocate the domain space of the operator.
RCP< const VectorSpaceBase< double > > range() const
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
RCP< const VectorSpaceBase< double > > domain() const
RCP< const SpmdVectorSpaceBase< double > > spmdRange() const
Return a smart pointer to the SpmdVectorSpaceBase object for the range.
EApplyEpetraOpAs
Determine how the apply an Epetra_Operator as a linear operator.
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.
bool opSupportedImpl(EOpTransp M_trans) const
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
Apply using Epetra_Operator::Apply(...)
Base abstract VectorSpaceBase class for all SPMD-based vector spaces.
virtual void scaleRightImpl(const VectorBase< double > &col_scaling)
virtual bool supportsScaleRightImpl() const
std::string description() const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< double > &X, const Ptr< MultiVectorBase< double > > &Y, const double alpha, const double beta) const
RCP< Epetra_Operator > epetra_op()