47 #include "Teko_EpetraHelpers.hpp"
50 #include "Thyra_EpetraLinearOp.hpp"
51 #include "Thyra_BlockedLinearOpBase.hpp"
52 #include "Thyra_DefaultMultipliedLinearOp.hpp"
53 #include "Thyra_DefaultDiagonalLinearOp.hpp"
54 #include "Thyra_DefaultZeroLinearOp.hpp"
55 #include "Thyra_DefaultBlockedLinearOp.hpp"
56 #include "Thyra_EpetraThyraWrappers.hpp"
57 #include "Thyra_SpmdVectorBase.hpp"
58 #include "Thyra_SpmdVectorSpaceBase.hpp"
59 #include "Thyra_ScalarProdVectorSpaceBase.hpp"
62 #include "Epetra_Vector.h"
65 #include "EpetraExt_ProductOperator.h"
66 #include "EpetraExt_MatrixMatrix.h"
69 #include "Teko_EpetraOperatorWrapper.hpp"
75 using Teuchos::rcp_dynamic_cast;
76 using Teuchos::rcpFromRef;
91 const Teuchos::RCP<const Thyra::LinearOpBase<double> > thyraDiagOp(
92 const RCP<const Epetra_Vector>& ev,
const Epetra_Map& map,
const std::string& lbl) {
93 const RCP<const Thyra::VectorBase<double> > thyraVec
94 = Thyra::create_Vector(ev, Thyra::create_VectorSpace(rcpFromRef(map)));
95 Teuchos::RCP<Thyra::LinearOpBase<double> > op =
96 Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<double>(thyraVec));
97 op->setObjectLabel(lbl);
111 const Teuchos::RCP<Thyra::LinearOpBase<double> > thyraDiagOp(
const RCP<Epetra_Vector>& ev,
112 const Epetra_Map& map,
113 const std::string& lbl) {
114 const RCP<Thyra::VectorBase<double> > thyraVec
115 = Thyra::create_Vector(ev, Thyra::create_VectorSpace(rcpFromRef(map)));
116 Teuchos::RCP<Thyra::LinearOpBase<double> > op =
117 Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<double>(thyraVec));
118 op->setObjectLabel(lbl);
131 void fillDefaultSpmdMultiVector(Teuchos::RCP<Thyra::DefaultSpmdMultiVector<double> >& spmdMV,
132 Teuchos::RCP<Epetra_MultiVector>& epetraMV) {
134 const RCP<const Thyra::SpmdVectorSpaceBase<double> > range = spmdMV->spmdSpace();
135 const RCP<const Thyra::ScalarProdVectorSpaceBase<double> > domain =
136 rcp_dynamic_cast<
const Thyra::ScalarProdVectorSpaceBase<double> >(spmdMV->domain());
138 TEUCHOS_ASSERT(domain->dim() == epetraMV->NumVectors());
143 if (epetraMV->ConstantStride())
144 epetraMV->ExtractView(&localValues, &leadingDim);
146 TEUCHOS_TEST_FOR_EXCEPT(
true);
149 spmdMV->initialize(range, domain,
150 Teuchos::arcp(localValues, 0, leadingDim * epetraMV->NumVectors(),
false),
154 Teuchos::set_extra_data<RCP<Epetra_MultiVector> >(epetraMV,
"Epetra_MultiVector",
155 Teuchos::outArg(spmdMV));
167 void identityRowIndices(
const Epetra_Map& rowMap,
const Epetra_CrsMatrix& mat,
168 std::vector<int>& outIndices) {
169 int maxSz = mat.GlobalMaxNumEntries();
170 std::vector<double> values(maxSz);
171 std::vector<int> indices(maxSz);
174 for (
int i = 0; i < rowMap.NumMyElements(); i++) {
175 bool rowIsIdentity =
true;
177 int rowGID = rowMap.GID(i);
178 mat.ExtractGlobalRowCopy(rowGID, maxSz, sz, &values[0], &indices[0]);
181 for (
int j = 0; j < sz; j++) {
182 int colGID = indices[j];
185 if (colGID == rowGID)
186 rowIsIdentity &= values[j] == 1.0;
188 rowIsIdentity &= values[j] == 0.0;
191 if (not rowIsIdentity)
break;
195 if (rowIsIdentity) outIndices.push_back(rowGID);
209 void zeroMultiVectorRowIndices(Epetra_MultiVector& mv,
const std::vector<int>& zeroIndices) {
210 int colCnt = mv.NumVectors();
211 std::vector<int>::const_iterator itr;
214 for (itr = zeroIndices.begin(); itr != zeroIndices.end(); ++itr) {
216 for (
int j = 0; j < colCnt; j++) TEUCHOS_TEST_FOR_EXCEPT(mv.ReplaceGlobalValue(*itr, j, 0.0));
230 const Teuchos::RCP<const Epetra_Operator>& op)
231 : zeroIndices_(zeroIndices), epetraOp_(op) {
232 label_ =
"zeroed( " + std::string(epetraOp_->Label()) +
" )";
243 int result = epetraOp_->Apply(X, Y);
246 zeroMultiVectorRowIndices(Y, zeroIndices_);
ZeroedOperator(const std::vector< int > &zeroIndices, const Teuchos::RCP< const Epetra_Operator > &op)
Constructor for a ZeroedOperator.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Perform a matrix-vector product with certain rows zeroed out.