47 #include "Teko_TpetraHelpers.hpp"
48 #include "Teko_ConfigDefs.hpp"
50 #ifdef TEKO_HAVE_EPETRA
51 #include "Thyra_EpetraLinearOp.hpp"
52 #include "Thyra_EpetraThyraWrappers.hpp"
55 #include "Epetra_Vector.h"
58 #include "EpetraExt_ProductOperator.h"
59 #include "EpetraExt_MatrixMatrix.h"
61 #include "Teko_EpetraOperatorWrapper.hpp"
65 #include "Thyra_BlockedLinearOpBase.hpp"
66 #include "Thyra_DefaultMultipliedLinearOp.hpp"
67 #include "Thyra_DefaultDiagonalLinearOp.hpp"
68 #include "Thyra_DefaultZeroLinearOp.hpp"
69 #include "Thyra_DefaultBlockedLinearOp.hpp"
71 #include "Thyra_SpmdVectorBase.hpp"
72 #include "Thyra_SpmdVectorSpaceBase.hpp"
73 #include "Thyra_ScalarProdVectorSpaceBase.hpp"
79 #include "Thyra_TpetraLinearOp.hpp"
80 #include "Thyra_TpetraMultiVector.hpp"
81 #include "Tpetra_CrsMatrix.hpp"
82 #include "Tpetra_Vector.hpp"
83 #include "Thyra_TpetraThyraWrappers.hpp"
84 #include "TpetraExt_MatrixMatrix.hpp"
89 using Teuchos::rcp_dynamic_cast;
90 using Teuchos::rcpFromRef;
93 namespace TpetraHelpers {
105 const Teuchos::RCP<const Thyra::LinearOpBase<ST> > thyraDiagOp(
106 const RCP<
const Tpetra::Vector<ST, LO, GO, NT> >& tv,
const Tpetra::Map<LO, GO, NT>& map,
107 const std::string& lbl) {
108 const RCP<const Thyra::VectorBase<ST> > thyraVec
109 = Thyra::createConstVector<ST, LO, GO, NT>(
110 tv, Thyra::createVectorSpace<ST, LO, GO, NT>(rcpFromRef(map)));
111 Teuchos::RCP<Thyra::LinearOpBase<ST> > op =
112 Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
113 op->setObjectLabel(lbl);
127 const Teuchos::RCP<Thyra::LinearOpBase<ST> > thyraDiagOp(
128 const RCP<Tpetra::Vector<ST, LO, GO, NT> >& tv,
const Tpetra::Map<LO, GO, NT>& map,
129 const std::string& lbl) {
130 const RCP<Thyra::VectorBase<ST> > thyraVec
131 = Thyra::createVector<ST, LO, GO, NT>(
132 tv, Thyra::createVectorSpace<ST, LO, GO, NT>(rcpFromRef(map)));
133 Teuchos::RCP<Thyra::LinearOpBase<ST> > op =
134 Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
135 op->setObjectLabel(lbl);
148 void fillDefaultSpmdMultiVector(Teuchos::RCP<Thyra::TpetraMultiVector<ST, LO, GO, NT> >& spmdMV,
149 Teuchos::RCP<Tpetra::MultiVector<ST, LO, GO, NT> >& tpetraMV) {
152 const RCP<Thyra::TpetraVectorSpace<ST, LO, GO, NT> > range =
153 Thyra::tpetraVectorSpace<ST, LO, GO, NT>(tpetraMV->getMap());
154 const RCP<const Thyra::ScalarProdVectorSpaceBase<ST> > domain =
155 rcp_dynamic_cast<
const Thyra::ScalarProdVectorSpaceBase<ST> >(spmdMV->domain());
157 TEUCHOS_ASSERT((
size_t)domain->dim() == tpetraMV->getNumVectors());
160 if (!tpetraMV->isConstantStride())
161 TEUCHOS_TEST_FOR_EXCEPT(
true);
164 spmdMV->initialize(range, domain, tpetraMV);
167 Teuchos::set_extra_data<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >(
168 tpetraMV,
"Tpetra::MultiVector", Teuchos::outArg(spmdMV));
180 void identityRowIndices(
const Tpetra::Map<LO, GO, NT>& rowMap,
181 const Tpetra::CrsMatrix<ST, LO, GO, NT>& mat, std::vector<GO>& outIndices) {
183 for (
size_t i = 0; i < rowMap.getLocalNumElements(); i++) {
184 bool rowIsIdentity =
true;
185 GO rowGID = rowMap.getGlobalElement(i);
187 size_t numEntries = mat.getNumEntriesInGlobalRow(i);
188 auto indices =
typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_global_inds_host_view_type(
189 Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"), numEntries);
190 auto values =
typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
191 Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"), numEntries);
193 mat.getGlobalRowCopy(rowGID, indices, values, numEntries);
196 for (
size_t j = 0; j < numEntries; j++) {
197 GO colGID = indices(j);
200 if (colGID == rowGID)
201 rowIsIdentity &= values(j) == 1.0;
203 rowIsIdentity &= values(j) == 0.0;
206 if (not rowIsIdentity)
break;
210 if (rowIsIdentity) outIndices.push_back(rowGID);
224 void zeroMultiVectorRowIndices(Tpetra::MultiVector<ST, LO, GO, NT>& mv,
225 const std::vector<GO>& zeroIndices) {
226 LO colCnt = mv.getNumVectors();
227 std::vector<GO>::const_iterator itr;
230 for (itr = zeroIndices.begin(); itr != zeroIndices.end(); ++itr) {
232 for (
int j = 0; j < colCnt; j++) mv.replaceGlobalValue(*itr, j, 0.0);
246 const Teuchos::RCP<
const Tpetra::Operator<ST, LO, GO, NT> >& op)
247 : zeroIndices_(zeroIndices), tpetraOp_(op) {}
251 Tpetra::MultiVector<ST, LO, GO, NT>& Y, Teuchos::ETransp mode, ST alpha,
259 tpetraOp_->apply(X, Y, mode, alpha, beta);
262 zeroMultiVectorRowIndices(Y, zeroIndices_);
265 bool isTpetraLinearOp(
const LinearOp& op) {
267 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
268 rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(op);
269 if (!tOp.is_null())
return true;
273 Thyra::EOpTransp transp = Thyra::NOTRANS;
274 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
275 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
276 tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(wrapped_op);
277 if (!tOp.is_null())
return true;
282 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > getTpetraCrsMatrix(
const LinearOp& op, ST* scalar,
285 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
286 rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(op);
287 if (!tOp.is_null()) {
288 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > matrix =
289 rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST, LO, GO, NT> >(tOp->getConstTpetraOperator(),
297 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
298 Thyra::EOpTransp eTransp = Thyra::NOTRANS;
299 Thyra::unwrap(op, scalar, &eTransp, &wrapped_op);
300 tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(wrapped_op,
true);
301 if (!tOp.is_null()) {
302 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > matrix =
303 rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST, LO, GO, NT> >(tOp->getConstTpetraOperator(),
306 if (eTransp == Thyra::NOTRANS) *transp =
false;
310 return Teuchos::null;
313 #ifdef TEKO_HAVE_EPETRA
314 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > epetraCrsMatrixToTpetra(
315 const RCP<const Epetra_CrsMatrix> A_e,
const RCP<
const Teuchos::Comm<int> > comm) {
320 int info = A_e->ExtractCrsDataPointers(ptr, ind, val);
321 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
322 "Could not extract data from Epetra_CrsMatrix");
323 const LO numRows = A_e->Graph().NumMyRows();
324 const LO nnz = A_e->Graph().NumMyEntries();
326 Teuchos::ArrayRCP<size_t> ptr2(numRows + 1);
327 Teuchos::ArrayRCP<int> ind2(nnz);
328 Teuchos::ArrayRCP<double> val2(nnz);
330 std::copy(ptr, ptr + numRows + 1, ptr2.begin());
331 std::copy(ind, ind + nnz, ind2.begin());
332 std::copy(val, val + nnz, val2.begin());
334 RCP<const Tpetra::Map<LO, GO, NT> > rowMap = epetraMapToTpetra(A_e->RowMap(), comm);
335 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > A_t =
336 Tpetra::createCrsMatrix<ST, LO, GO, NT>(rowMap, A_e->GlobalMaxNumEntries());
338 RCP<const Tpetra::Map<LO, GO, NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(), comm);
339 RCP<const Tpetra::Map<LO, GO, NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(), comm);
340 RCP<const Tpetra::Map<LO, GO, NT> > colMap = epetraMapToTpetra(A_e->ColMap(), comm);
342 A_t->replaceColMap(colMap);
343 A_t->setAllValues(ptr2, ind2, val2);
344 A_t->fillComplete(domainMap, rangeMap);
348 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > nonConstEpetraCrsMatrixToTpetra(
349 const RCP<Epetra_CrsMatrix> A_e,
const RCP<
const Teuchos::Comm<int> > comm) {
354 int info = A_e->ExtractCrsDataPointers(ptr, ind, val);
355 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
356 "Could not extract data from Epetra_CrsMatrix");
357 const LO numRows = A_e->Graph().NumMyRows();
358 const LO nnz = A_e->Graph().NumMyEntries();
360 Teuchos::ArrayRCP<size_t> ptr2(numRows + 1);
361 Teuchos::ArrayRCP<int> ind2(nnz);
362 Teuchos::ArrayRCP<double> val2(nnz);
364 std::copy(ptr, ptr + numRows + 1, ptr2.begin());
365 std::copy(ind, ind + nnz, ind2.begin());
366 std::copy(val, val + nnz, val2.begin());
368 RCP<const Tpetra::Map<LO, GO, NT> > rowMap = epetraMapToTpetra(A_e->RowMap(), comm);
369 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > A_t =
370 Tpetra::createCrsMatrix<ST, LO, GO, NT>(rowMap, A_e->GlobalMaxNumEntries());
372 RCP<const Tpetra::Map<LO, GO, NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(), comm);
373 RCP<const Tpetra::Map<LO, GO, NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(), comm);
374 RCP<const Tpetra::Map<LO, GO, NT> > colMap = epetraMapToTpetra(A_e->ColMap(), comm);
376 A_t->replaceColMap(colMap);
377 A_t->setAllValues(ptr2, ind2, val2);
378 A_t->fillComplete(domainMap, rangeMap);
382 RCP<const Tpetra::Map<LO, GO, NT> > epetraMapToTpetra(
const Epetra_Map eMap,
383 const RCP<
const Teuchos::Comm<int> > comm) {
384 std::vector<int> intGIDs(eMap.NumMyElements());
385 eMap.MyGlobalElements(&intGIDs[0]);
387 std::vector<GO> myGIDs(eMap.NumMyElements());
388 for (
int k = 0; k < eMap.NumMyElements(); k++) myGIDs[k] = (GO)intGIDs[k];
391 new const Tpetra::Map<LO, GO, NT>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
392 Teuchos::ArrayView<GO>(myGIDs), 0, comm));
394 #endif // TEKO_HAVE_EPETRA
void apply(const Tpetra::MultiVector< ST, LO, GO, NT > &X, Tpetra::MultiVector< ST, LO, GO, NT > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, ST alpha=Teuchos::ScalarTraits< ST >::one(), ST beta=Teuchos::ScalarTraits< ST >::zero()) const
Perform a matrix-vector product with certain rows zeroed out.
ZeroedOperator(const std::vector< GO > &zeroIndices, const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &op)
Constructor for a ZeroedOperator.