43 #include "Thyra_EpetraExtDiagScaledMatProdTransformer.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 "EpetraExt_MatrixMatrix.h" 55 #include "Teuchos_Assert.hpp" 67 TEUCHOS_TEST_FOR_EXCEPT(
true);
72 RCP<LinearOpBase<double> >
75 return nonconstEpetraLinearOp();
84 using EpetraExt::MatrixMatrix;
86 using Teuchos::rcp_dynamic_cast;
87 using Teuchos::dyn_cast;
96 bool haveDiagScaling = (multi_op.
numOps()==3);
99 const RCP<const LinearOpBase<double> > op_B = multi_op.
getOp(0);
100 double B_scalar = 0.0;
102 RCP<const LinearOpBase<double> > B;
103 unwrap( op_B, &B_scalar, &B_transp, &B );
107 RCP<const VectorBase<double> > d;
108 double D_scalar = 1.0;
109 if(haveDiagScaling) {
110 const RCP<const LinearOpBase<double> > op_D = multi_op.
getOp(1);
112 RCP<const LinearOpBase<double> > D;
113 unwrap( op_D, &D_scalar, &D_transp, &D );
118 const RCP<const LinearOpBase<double> > op_G = multi_op.
getOp(haveDiagScaling ? 2 : 1);
119 double G_scalar = 0.0;
121 RCP<const LinearOpBase<double> > G;
122 unwrap( op_G, &G_scalar, &G_transp, &G );
130 const RCP<const Epetra_CrsMatrix> epetra_B =
135 RCP<const Epetra_Vector> epetra_d;
136 if(haveDiagScaling) {
142 const RCP<const Epetra_CrsMatrix> epetra_G =
146 const Epetra_Map op_inout_row_map
147 = (B_transp==
CONJTRANS ? epetra_B->ColMap() : epetra_B->RowMap());
148 const Epetra_Map op_inout_col_map
149 = (G_transp==
CONJTRANS ? epetra_B->RowMap() : epetra_B->ColMap());
160 RCP<Epetra_CrsMatrix> epetra_op =
161 rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.
epetra_op());
162 if(is_null(epetra_op)) {
163 epetra_op = Teuchos::rcp(
164 new Epetra_CrsMatrix(::Copy, op_inout_row_map, 0));
168 RCP<const Epetra_CrsMatrix> epetra_BD;
169 if(haveDiagScaling) {
171 RCP<Epetra_CrsMatrix> epetra_BD_temp = rcp(
new Epetra_CrsMatrix(*epetra_B));
175 epetra_BD_temp->LeftScale(*epetra_d);
177 epetra_BD_temp->RightScale(*epetra_d);
179 epetra_BD = epetra_BD_temp;
182 epetra_BD = epetra_B;
185 int mm_error = MatrixMatrix::Multiply( *epetra_BD, B_transp==
CONJTRANS,
186 *epetra_G, G_transp==
CONJTRANS, *epetra_op);
187 TEUCHOS_TEST_FOR_EXCEPTION(mm_error!=0,std::invalid_argument,
188 "EpetraExt::MatrixMatrix::Multiply failed returning error code " << mm_error <<
".");
191 if(B_scalar*G_scalar*D_scalar!=1.0)
192 epetra_op->Scale(B_scalar*G_scalar*D_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()