43 #include "Thyra_EpetraExtDiagScalingTransformer.hpp" 44 #include "Thyra_MultipliedLinearOpBase.hpp" 45 #include "Thyra_DiagonalLinearOpBase.hpp" 46 #include "Thyra_ScaledAdjointLinearOpBase.hpp" 47 #include "Thyra_EpetraLinearOp.hpp" 48 #include "Thyra_get_Epetra_Operator.hpp" 50 #include "Epetra_Map.h" 51 #include "Epetra_LocalMap.h" 52 #include "Epetra_SerialComm.h" 53 #include "Epetra_CrsMatrix.h" 54 #include "Teuchos_Assert.hpp" 55 #include "Teuchos_RCP.hpp" 69 using Teuchos::rcp_dynamic_cast;
70 using Teuchos::dyn_cast;
83 RCP<const LinearOpBase<double> > A;
85 if(transp!=
NOTRANS)
return false;
88 RCP<const LinearOpBase<double> > B;
90 if(transp!=
NOTRANS)
return false;
93 RCP<const DiagonalLinearOpBase<double> > dA
95 RCP<const DiagonalLinearOpBase<double> > dB
98 if(dA==Teuchos::null && dB!=Teuchos::null)
101 if(dA!=Teuchos::null && dB==Teuchos::null)
108 RCP<LinearOpBase<double> >
111 return nonconstEpetraLinearOp();
121 using Teuchos::rcp_dynamic_cast;
122 using Teuchos::dyn_cast;
131 TEUCHOS_ASSERT(multi_op.
numOps()==2);
134 const RCP<const LinearOpBase<double> > op_A = multi_op.
getOp(0);
135 double A_scalar = 0.0;
137 RCP<const LinearOpBase<double> > A;
138 unwrap( op_A, &A_scalar, &A_transp, &A );
142 const RCP<const LinearOpBase<double> > op_B = multi_op.
getOp(1);
143 double B_scalar = 0.0;
145 RCP<const LinearOpBase<double> > B;
146 unwrap( op_B, &B_scalar, &B_transp, &B );
154 RCP<const DiagonalLinearOpBase<double> > dA
156 RCP<const DiagonalLinearOpBase<double> > dB
159 RCP<const Epetra_CrsMatrix> epetra_A;
160 RCP<const Epetra_CrsMatrix> epetra_B;
161 if(dA==Teuchos::null) {
162 bool exactly_one_op_must_be_diagonal__dB_neq_null = dB!=Teuchos::null;
163 TEUCHOS_ASSERT(exactly_one_op_must_be_diagonal__dB_neq_null);
168 else if(dB==Teuchos::null) {
169 bool exactly_one_op_must_be_diagonal__dA_neq_null = dA!=Teuchos::null;
170 TEUCHOS_ASSERT(exactly_one_op_must_be_diagonal__dA_neq_null);
176 bool exactly_one_op_must_be_diagonal=
false;
177 TEUCHOS_ASSERT(exactly_one_op_must_be_diagonal);
180 TEUCHOS_ASSERT( A_transp ==
NOTRANS );
181 TEUCHOS_ASSERT( B_transp ==
NOTRANS );
186 RCP<Epetra_CrsMatrix> epetra_op =
187 rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.
epetra_op());
189 bool rightScale = dB!=Teuchos::null;
192 RCP<const Epetra_Vector> v =
get_Epetra_Vector(epetra_A->OperatorDomainMap(),dB->getDiag());
196 if(epetra_op==Teuchos::null)
197 epetra_op = Teuchos::rcp(
new Epetra_CrsMatrix(*epetra_A));
199 *epetra_op = *epetra_A;
200 epetra_op->RightScale(*v);
203 RCP<const Epetra_Vector> v =
get_Epetra_Vector(epetra_B->OperatorRangeMap(),dA->getDiag());
207 if(epetra_op==Teuchos::null)
208 epetra_op = Teuchos::rcp(
new Epetra_CrsMatrix(*epetra_B));
210 *epetra_op = *epetra_B;
211 epetra_op->LeftScale(*v);
214 epetra_op->Scale(A_scalar*B_scalar);
Concrete LinearOpBase adapter subclass for Epetra_Operator object.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
void unwrap(const LinearOpBase< Scalar > &Op, Scalar *scalar, EOpTransp *transp, const LinearOpBase< Scalar > **origOp)
Extract the overallScalar, overallTransp and const origOp from a const LinearOpBase object...
Use the non-transposed operator.
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
virtual int numOps() const =0
Returns the number of constituent operators.
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.
Interface class for implicitly multiplied linear operators.
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.
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< Scalar > &)
Get smart pointer to non-const Epetra_Operator object from reference to a non-const EpetraLinearOp ac...
virtual Teuchos::RCP< const LinearOpBase< Scalar > > getOp(const int k) const =0
Return the kth constant constituent operator.
Interface class for for diagonal linear operators.
RCP< Epetra_Operator > epetra_op()