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"
88using Teuchos::rcpFromRef;
89using Teuchos::rcp_dynamic_cast;
93namespace TpetraHelpers {
105const Teuchos::RCP<const Thyra::LinearOpBase<ST> > thyraDiagOp(
const RCP<
const Tpetra::Vector<ST,LO,GO,NT> > & tv,
const Tpetra::Map<LO,GO,NT> & map,
106 const std::string & lbl)
108 const RCP<const Thyra::VectorBase<ST> > thyraVec
109 = Thyra::createConstVector<ST,LO,GO,NT>(tv,Thyra::createVectorSpace<ST,LO,GO,NT>(rcpFromRef(map)));
110 Teuchos::RCP<Thyra::LinearOpBase<ST> > op
111 = Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
112 op->setObjectLabel(lbl);
126const Teuchos::RCP<Thyra::LinearOpBase<ST> > thyraDiagOp(
const RCP<Tpetra::Vector<ST,LO,GO,NT> > & tv,
const Tpetra::Map<LO,GO,NT> & map,
127 const std::string & lbl)
129 const RCP<Thyra::VectorBase<ST> > thyraVec
130 = Thyra::createVector<ST,LO,GO,NT>(tv,Thyra::createVectorSpace<ST,LO,GO,NT>(rcpFromRef(map)));
131 Teuchos::RCP<Thyra::LinearOpBase<ST> > op
132 = Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
133 op->setObjectLabel(lbl);
146void fillDefaultSpmdMultiVector(Teuchos::RCP<Thyra::TpetraMultiVector<ST,LO,GO,NT> > & spmdMV,
147 Teuchos::RCP<Tpetra::MultiVector<ST,LO,GO,NT> > & tpetraMV)
151 const RCP<Thyra::TpetraVectorSpace<ST,LO,GO,NT> > range = Thyra::tpetraVectorSpace<ST,LO,GO,NT>(tpetraMV->getMap());
152 const RCP<const Thyra::ScalarProdVectorSpaceBase<ST> > domain
153 = rcp_dynamic_cast<const Thyra::ScalarProdVectorSpaceBase<ST> >(spmdMV->domain());
155 TEUCHOS_ASSERT((
size_t) domain->dim()==tpetraMV->getNumVectors());
158 if(!tpetraMV->isConstantStride())
159 TEUCHOS_TEST_FOR_EXCEPT(
true);
162 spmdMV->initialize(range, domain, tpetraMV);
165 Teuchos::set_extra_data<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > >(tpetraMV,
"Tpetra::MultiVector",Teuchos::outArg(spmdMV));
177void identityRowIndices(
const Tpetra::Map<LO,GO,NT> & rowMap,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & mat,std::vector<GO> & outIndices)
180 for(
size_t i=0;i<rowMap.getLocalNumElements();i++) {
181 bool rowIsIdentity =
true;
182 GO rowGID = rowMap.getGlobalElement(i);
184 size_t numEntries = mat.getNumEntriesInGlobalRow (i);
185 auto indices =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),numEntries);
186 auto values =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),numEntries);
188 mat.getGlobalRowCopy(rowGID,indices,values,numEntries);
191 for(
size_t j=0;j<numEntries;j++) {
192 GO colGID = indices(j);
195 if(colGID==rowGID) rowIsIdentity &= values(j)==1.0;
196 else rowIsIdentity &= values(j)==0.0;
199 if(not rowIsIdentity)
205 outIndices.push_back(rowGID);
219void zeroMultiVectorRowIndices(Tpetra::MultiVector<ST,LO,GO,NT> & mv,
const std::vector<GO> & zeroIndices)
221 LO colCnt = mv.getNumVectors();
222 std::vector<GO>::const_iterator itr;
225 for(itr=zeroIndices.begin();itr!=zeroIndices.end();++itr) {
228 for(
int j=0;j<colCnt;j++)
229 mv.replaceGlobalValue(*itr,j,0.0);
243 const Teuchos::RCP<
const Tpetra::Operator<ST,LO,GO,NT> > & op)
244 : zeroIndices_(zeroIndices), tpetraOp_(op)
248void ZeroedOperator::apply(
const Tpetra::MultiVector<ST,LO,GO,NT> & X, Tpetra::MultiVector<ST,LO,GO,NT> & Y, Teuchos::ETransp mode, ST alpha, ST beta)
const
256 tpetraOp_->apply(X,Y,mode,alpha,beta);
259 zeroMultiVectorRowIndices(Y,zeroIndices_);
262bool isTpetraLinearOp(
const LinearOp & op)
265 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
271 Thyra::EOpTransp transp = Thyra::NOTRANS;
272 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
273 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
274 tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(wrapped_op);
281RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > getTpetraCrsMatrix(
const LinearOp & op, ST *scalar,
bool *transp)
284 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
286 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > matrix = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
293 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
294 Thyra::EOpTransp eTransp = Thyra::NOTRANS;
295 Thyra::unwrap(op, scalar, &eTransp, &wrapped_op);
296 tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(wrapped_op,
true);
298 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > matrix = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
300 if(eTransp == Thyra::NOTRANS)
305 return Teuchos::null;
308#ifdef TEKO_HAVE_EPETRA
309RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > epetraCrsMatrixToTpetra(
const RCP<const Epetra_CrsMatrix> A_e,
const RCP<
const Teuchos::Comm<int> > comm)
315 int info = A_e->ExtractCrsDataPointers (ptr, ind, val);
316 TEUCHOS_TEST_FOR_EXCEPTION(info!=0,std::logic_error,
"Could not extract data from Epetra_CrsMatrix");
317 const LO numRows = A_e->Graph ().NumMyRows ();
318 const LO nnz = A_e->Graph ().NumMyEntries ();
320 Teuchos::ArrayRCP<size_t> ptr2 (numRows+1);
321 Teuchos::ArrayRCP<int> ind2 (nnz);
322 Teuchos::ArrayRCP<double> val2 (nnz);
324 std::copy (ptr, ptr + numRows + 1, ptr2.begin ());
325 std::copy (ind, ind + nnz, ind2.begin ());
326 std::copy (val, val + nnz, val2.begin ());
328 RCP<const Tpetra::Map<LO,GO,NT> > rowMap = epetraMapToTpetra(A_e->RowMap(),comm);
329 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > A_t = Tpetra::createCrsMatrix<ST,LO,GO,NT>(rowMap, A_e->GlobalMaxNumEntries());
331 RCP<const Tpetra::Map<LO,GO,NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(),comm);
332 RCP<const Tpetra::Map<LO,GO,NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(),comm);
333 RCP<const Tpetra::Map<LO,GO,NT> > colMap = epetraMapToTpetra(A_e->ColMap(),comm);
335 A_t->replaceColMap(colMap);
336 A_t->setAllValues (ptr2, ind2, val2);
337 A_t->fillComplete(domainMap,rangeMap);
342RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > nonConstEpetraCrsMatrixToTpetra(
const RCP<Epetra_CrsMatrix> A_e,
const RCP<
const Teuchos::Comm<int> > comm)
348 int info = A_e->ExtractCrsDataPointers (ptr, ind, val);
349 TEUCHOS_TEST_FOR_EXCEPTION(info!=0,std::logic_error,
"Could not extract data from Epetra_CrsMatrix");
350 const LO numRows = A_e->Graph ().NumMyRows ();
351 const LO nnz = A_e->Graph ().NumMyEntries ();
353 Teuchos::ArrayRCP<size_t> ptr2 (numRows+1);
354 Teuchos::ArrayRCP<int> ind2 (nnz);
355 Teuchos::ArrayRCP<double> val2 (nnz);
357 std::copy (ptr, ptr + numRows + 1, ptr2.begin ());
358 std::copy (ind, ind + nnz, ind2.begin ());
359 std::copy (val, val + nnz, val2.begin ());
361 RCP<const Tpetra::Map<LO,GO,NT> > rowMap = epetraMapToTpetra(A_e->RowMap(),comm);
362 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > A_t = Tpetra::createCrsMatrix<ST,LO,GO,NT>(rowMap, A_e->GlobalMaxNumEntries());
364 RCP<const Tpetra::Map<LO,GO,NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(),comm);
365 RCP<const Tpetra::Map<LO,GO,NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(),comm);
366 RCP<const Tpetra::Map<LO,GO,NT> > colMap = epetraMapToTpetra(A_e->ColMap(),comm);
368 A_t->replaceColMap(colMap);
369 A_t->setAllValues (ptr2, ind2, val2);
370 A_t->fillComplete(domainMap,rangeMap);
375RCP<const Tpetra::Map<LO,GO,NT> > epetraMapToTpetra(
const Epetra_Map eMap,
const RCP<
const Teuchos::Comm<int> > comm)
377 std::vector<int> intGIDs(eMap.NumMyElements());
378 eMap.MyGlobalElements(&intGIDs[0]);
380 std::vector<GO> myGIDs(eMap.NumMyElements());
381 for(
int k = 0; k < eMap.NumMyElements(); k++)
382 myGIDs[k] = (GO) intGIDs[k];
384 return rcp(
new const Tpetra::Map<LO,GO,NT>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),Teuchos::ArrayView<GO>(myGIDs),0,comm));
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.