44 #include "Thyra_MultipliedLinearOpBase.hpp" 45 #include "Thyra_DiagonalLinearOpBase.hpp" 46 #include "Thyra_ScaledAdjointLinearOpBase.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" 65 const LinearOpBase<double> &op_in)
const 69 using Teuchos::rcp_dynamic_cast;
70 using Teuchos::dyn_cast;
72 const MultipliedLinearOpBase<double> &multi_op =
73 dyn_cast<
const MultipliedLinearOpBase<double> >(op_in);
76 if(multi_op.numOps()!=2)
83 RCP<const LinearOpBase<double> > A;
84 unwrap( multi_op.getOp(0), &scalar, &transp, &A );
85 if(transp!=
NOTRANS)
return false;
88 RCP<const LinearOpBase<double> > B;
89 unwrap( multi_op.getOp(1), &scalar, &transp, &B );
90 if(transp!=
NOTRANS)
return false;
93 RCP<const DiagonalLinearOpBase<double> > dA
94 = rcp_dynamic_cast<
const DiagonalLinearOpBase<double> >(A);
95 RCP<const DiagonalLinearOpBase<double> > dB
96 = rcp_dynamic_cast<
const DiagonalLinearOpBase<double> >(B);
98 if(dA==Teuchos::null && dB!=Teuchos::null)
101 if(dA!=Teuchos::null && dB==Teuchos::null)
108 RCP<LinearOpBase<double> >
111 return nonconstEpetraLinearOp();
116 const LinearOpBase<double> &op_in,
117 const Ptr<LinearOpBase<double> > &op_inout)
const 121 using Teuchos::rcp_dynamic_cast;
122 using Teuchos::dyn_cast;
128 const MultipliedLinearOpBase<double> &multi_op =
129 dyn_cast<
const MultipliedLinearOpBase<double> >(op_in);
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
155 = rcp_dynamic_cast<
const DiagonalLinearOpBase<double> >(A);
156 RCP<const DiagonalLinearOpBase<double> > dB
157 = rcp_dynamic_cast<
const DiagonalLinearOpBase<double> >(B);
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.
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< 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< double > &op)
Full specialization for Scalar=double.
RCP< Epetra_Operator > epetra_op()