44 #include "Thyra_AddedLinearOpBase.hpp" 45 #include "Thyra_ScaledAdjointLinearOpBase.hpp" 49 #include "Thyra_DiagonalLinearOpBase.hpp" 50 #include "Thyra_DefaultDiagonalLinearOp.hpp" 51 #include "Thyra_IdentityLinearOpBase.hpp" 52 #include "Thyra_VectorStdOps.hpp" 53 #include "Epetra_Map.h" 54 #include "Epetra_LocalMap.h" 55 #include "Epetra_SerialComm.h" 56 #include "Epetra_Vector.h" 57 #include "Epetra_CrsMatrix.h" 58 #include "Teuchos_Assert.hpp" 59 #include "EpetraExt_ConfigDefs.h" 60 #include "EpetraExt_MatrixMatrix.h" 61 #include "EpetraExt_MMHelpers.h" 62 #include "EpetraExt_Transpose_RowMatrix.h" 65 #include "EpetraExt_RowMatrixOut.h" 75 const LinearOpBase<double> &op_in)
const 77 TEUCHOS_TEST_FOR_EXCEPT(
true);
82 RCP<LinearOpBase<double> >
85 return nonconstEpetraLinearOp();
90 const LinearOpBase<double> &op_in,
91 const Ptr<LinearOpBase<double> > &op_inout)
const 94 using EpetraExt::MatrixMatrix;
96 using Teuchos::rcp_dynamic_cast;
97 using Teuchos::dyn_cast;
103 const AddedLinearOpBase<double> & add_op =
104 dyn_cast<
const AddedLinearOpBase<double> >(op_in);
107 TEUCHOS_ASSERT_EQUALITY( add_op.numOps(), 2 );
112 const RCP<const LinearOpBase<double> > op_A = add_op.getOp(0);
113 double A_scalar = 0.0;
115 RCP<const LinearOpBase<double> > A;
116 unwrap( op_A, &A_scalar, &A_transp, &A );
120 const RCP<const LinearOpBase<double> > op_B = add_op.getOp(1);
121 double B_scalar = 0.0;
123 RCP<const LinearOpBase<double> > B;
124 unwrap( op_B, &B_scalar, &B_transp, &B );
132 if(rcp_dynamic_cast<
const Thyra::IdentityLinearOpBase<double> >(A)!=Teuchos::null) {
133 RCP<Thyra::VectorBase<double> > d = Thyra::createMember(A->domain(),
"d");
134 Thyra::V_S( d.ptr(), 1.0 );
135 A = Thyra::diagonal(d);
137 if(rcp_dynamic_cast<
const Thyra::IdentityLinearOpBase<double> >(B)!=Teuchos::null) {
138 RCP<Thyra::VectorBase<double> > d = Thyra::createMember(B->domain(),
"d");
139 Thyra::V_S( d.ptr(), 1.0 );
140 B = Thyra::diagonal(d);
145 RCP<const DiagonalLinearOpBase<double> > dA
146 = rcp_dynamic_cast<
const DiagonalLinearOpBase<double> >(A);
147 RCP<const DiagonalLinearOpBase<double> > dB
148 = rcp_dynamic_cast<
const DiagonalLinearOpBase<double> >(B);
151 RCP<const Epetra_CrsMatrix> epetra_A;
152 RCP<const Epetra_CrsMatrix> epetra_B;
153 if(dA==Teuchos::null)
155 if(dB==Teuchos::null)
162 if(epetra_A!=Teuchos::null && epetra_B!=Teuchos::null) {
169 RCP<Epetra_CrsMatrix> epetra_op =
170 rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.
epetra_op());
171 Epetra_CrsMatrix * epetra_op_raw = epetra_op.get();
174 const int add_epetra_B_err
175 = EpetraExt::MatrixMatrix::Add(*epetra_A,A_transp==
CONJTRANS,A_scalar,*epetra_B,B_transp==
CONJTRANS,B_scalar,epetra_op_raw);
176 if(epetra_op==Teuchos::null)
177 epetra_op = Teuchos::rcp(epetra_op_raw);
179 TEUCHOS_ASSERT_EQUALITY( add_epetra_B_err, 0 );
181 epetra_op->FillComplete();
186 else if((dA!=Teuchos::null && epetra_B!=Teuchos::null) ||
187 (dB!=Teuchos::null && epetra_A!=Teuchos::null)) {
190 RCP<const Epetra_CrsMatrix> crsMat = (dA!=Teuchos::null) ? epetra_B : epetra_A;
191 double matScalar = (dA!=Teuchos::null) ? B_scalar : A_scalar;
192 RCP<const DiagonalLinearOpBase<double> > diag = (dA!=Teuchos::null) ? dA : dB;
193 double diagScalar = (dA!=Teuchos::null) ? A_scalar : B_scalar;
195 TEUCHOS_ASSERT(crsMat!=Teuchos::null);
196 TEUCHOS_ASSERT(diag!=Teuchos::null);
200 RCP<Epetra_CrsMatrix> epetra_op =
201 rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.
epetra_op());
203 if(epetra_op==Teuchos::null)
204 epetra_op = Teuchos::rcp(
new Epetra_CrsMatrix(*crsMat));
206 *epetra_op = *crsMat;
209 RCP<const Epetra_Vector> v =
get_Epetra_Vector(epetra_op->OperatorDomainMap(),diag->getDiag());
212 epetra_op->Scale(matScalar);
215 RCP<Epetra_Vector> diagonal = rcp(
new Epetra_Vector(epetra_op->OperatorDomainMap()));
216 TEUCHOS_TEST_FOR_EXCEPTION(epetra_op->ExtractDiagonalCopy(*diagonal),std::runtime_error,
217 "Thyra::EpetraExtractAddTransformer::transform ExtractDiagonalCopy failed!");;
218 diagonal->Update(diagScalar,*v,1.0);
219 TEUCHOS_TEST_FOR_EXCEPTION(epetra_op->ReplaceDiagonalValues(*diagonal),std::runtime_error,
220 "Thyra::EpetraExtractAddTransformer::transform ReplaceDiagonalValues failed!");;
226 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error,
227 "Your case of adding Epetra operators is not yet implemented! Contact the Thyra developers.");
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()