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,
78 const RCP<
const VectorSpaceBase<double> > &range_in,
79 const RCP<
const VectorSpaceBase<double> > &domain_in
83 using Teuchos::rcp_dynamic_cast;
84 typedef SpmdVectorSpaceBase<double> SPMDVSB;
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);
110 rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(
op_);
121 const RCP<
const VectorSpaceBase<double> > &range_in,
122 const RCP<
const VectorSpaceBase<double> > &domain_in,
123 const RCP<Epetra_Operator> &op,
130 using Teuchos::rcp_dynamic_cast;
131 typedef SpmdVectorSpaceBase<double> SPMDVSB;
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);
151 rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(
op_);
169 RCP<Epetra_Operator> *op,
173 RCP<
const VectorSpaceBase<double> > *range_out,
174 RCP<
const VectorSpaceBase<double> > *domain_out
182 if(range_out) *range_out =
range_;
183 if(domain_out) *domain_out =
domain_;
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
243 const Ptr<RCP<const Epetra_Operator> > &epetraOp,
244 const Ptr<EOpTransp> &epetraOpTransp,
245 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
246 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
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) ) {
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";
362 const EOpTransp M_trans,
363 const MultiVectorBase<double> &X_in,
364 const Ptr<MultiVectorBase<double> > &Y_inout,
370 THYRA_FUNC_TIME_MONITOR(
"Thyra::EpetraLinearOp::euclideanApply");
372 const EOpTransp real_M_trans = real_trans(M_trans);
376 THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
377 "EpetraLinearOp::euclideanApply(...)", *
this, M_trans, X_in, Y_inout
379 TEUCHOS_TEST_FOR_EXCEPTION(
381 Exceptions::OpNotSupported,
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);
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);
538 using Teuchos::rcpFromRef;
539 const RCP<const Epetra_Vector> row_scaling =
547 using Teuchos::rcpFromRef;
548 const RCP<const Epetra_Vector> col_scaling =
558 const RowStatLinearOpBaseUtils::ERowStat rowStat)
const 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,
576 const Ptr<VectorBase<double> > &rowStatVec_in
579 using Teuchos::rcpFromPtr;
580 const RCP<Epetra_Vector> rowStatVec =
583 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
586 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
591 TEUCHOS_TEST_FOR_EXCEPT(
true);
599 RCP< const SpmdVectorSpaceBase<double> >
601 const RCP<Epetra_Operator> &op,
605 return Teuchos::rcp_dynamic_cast<
const SpmdVectorSpaceBase<double> >(
612 RCP<const SpmdVectorSpaceBase<double> >
614 const RCP<Epetra_Operator> &op,
618 return Teuchos::rcp_dynamic_cast<
const SpmdVectorSpaceBase<double> >(
630 ?
op_->OperatorRangeMap() :
op_->OperatorDomainMap() );
638 ?
op_->OperatorDomainMap() :
op_->OperatorRangeMap() );
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...
const Epetra_Map & getRangeMap() const
virtual RCP< const SpmdVectorSpaceBase< double > > allocateRange(const RCP< Epetra_Operator > &op, EOpTransp op_trans) const
Allocate the range space of the operator.
const std::string toString(const EAdjointEpetraOp adjointEpetraOp)
virtual bool supportsScaleLeftImpl() 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.
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.
RCP< const SpmdVectorSpaceBase< double > > spmdDomain() const
Return a smart pointer to the SpmdVectorSpaceBase object for the domain.
void computeAbsRowSum(Epetra_Vector &rowStatVec_in) const
Compute the absolute row sum for this matrix.
void getEpetraOpView(const Ptr< RCP< const Epetra_Operator > > &epetraOp, const Ptr< EOpTransp > &epetraOpTransp, const Ptr< EApplyEpetraOpAs > &epetraOpApplyAs, const Ptr< EAdjointEpetraOp > &epetraOpAdjointSupport) const
RCP< const SpmdVectorSpaceBase< double > > domain_
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.
RCP< const SpmdVectorSpaceBase< double > > range_
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)
EAdjointEpetraOp adjointSupport_
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.
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.
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
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
Apply using Epetra_Operator::Apply(...)
virtual void scaleRightImpl(const VectorBase< double > &col_scaling)
virtual bool supportsScaleRightImpl() const
EApplyEpetraOpAs applyAs_
const Epetra_Map & getDomainMap() const
RCP< Epetra_Operator > op_
RCP< Epetra_RowMatrix > rowMatrix_
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()