45#include "Teko_TpetraOperatorWrapper.hpp"
46#include "Thyra_SpmdVectorBase.hpp"
47#include "Thyra_MultiVectorStdOps.hpp"
48#include "Teuchos_DefaultSerialComm.hpp"
49#include "Thyra_TpetraLinearOp.hpp"
50#include "Tpetra_Vector.hpp"
51#include "Thyra_TpetraThyraWrappers.hpp"
54#include "Thyra_BlockedLinearOpBase.hpp"
55#include "Thyra_ProductVectorSpaceBase.hpp"
56#include "Thyra_TpetraLinearOp.hpp"
58#include "Teko_TpetraThyraConverter.hpp"
59#include "Teuchos_Ptr.hpp"
63namespace TpetraHelpers {
66using namespace Teuchos;
69DefaultMappingStrategy::DefaultMappingStrategy(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp,
const Teuchos::Comm<Thyra::Ordinal> & comm)
71 RCP<Teuchos::Comm<Thyra::Ordinal> > newComm = comm.duplicate();
83 Teko::TpetraHelpers::blockTpetraToThyra(x,thyraVec);
88 Teko::TpetraHelpers::blockThyraToTpetra(thyraVec,v);
91TpetraOperatorWrapper::TpetraOperatorWrapper()
93 useTranspose_ =
false;
94 mapStrategy_ = Teuchos::null;
95 thyraOp_ = Teuchos::null;
96 comm_ = Teuchos::null;
97 label_ = Teuchos::null;
100TpetraOperatorWrapper::TpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<ST> > & thyraOp)
102 SetOperator(thyraOp);
105TpetraOperatorWrapper::TpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<ST> > & thyraOp,
const RCP<const MappingStrategy> & mapStrategy)
106 : mapStrategy_(mapStrategy)
108 SetOperator(thyraOp);
111TpetraOperatorWrapper::TpetraOperatorWrapper(
const RCP<const MappingStrategy> & mapStrategy)
112 : mapStrategy_(mapStrategy)
114 useTranspose_ =
false;
115 thyraOp_ = Teuchos::null;
116 comm_ = Teuchos::null;
117 label_ = Teuchos::null;
120void TpetraOperatorWrapper::SetOperator(
const RCP<
const LinearOpBase<ST> > & thyraOp,
bool buildMap)
122 useTranspose_ =
false;
124 comm_ = getThyraComm(*thyraOp);
125 label_ = thyraOp_->description();
126 if(mapStrategy_==Teuchos::null && buildMap)
127 mapStrategy_ = Teuchos::rcp(
new DefaultMappingStrategy(thyraOp,*comm_));
130double TpetraOperatorWrapper::NormInf()
const
132 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
133 "TpetraOperatorWrapper::NormInf not implemated");
137void TpetraOperatorWrapper::apply(
const Tpetra::MultiVector<ST,LO,GO,NT>& X, Tpetra::MultiVector<ST,LO,GO,NT>& Y,Teuchos::ETransp ,ST alpha, ST beta)
const
142 RCP<Thyra::MultiVectorBase<ST> > tX;
143 RCP<Thyra::MultiVectorBase<ST> > tY;
145 tX = Thyra::createMembers(thyraOp_->domain(),X.getNumVectors());
146 tY = Thyra::createMembers(thyraOp_->range(),X.getNumVectors());
148 Thyra::assign(tX.ptr(),0.0);
149 Thyra::assign(tY.ptr(),0.0);
152 mapStrategy_->copyTpetraIntoThyra(X, tX.ptr());
153 mapStrategy_->copyTpetraIntoThyra(Y, tY.ptr());
156 thyraOp_->apply(Thyra::NOTRANS,*tX,tY.ptr(),alpha,beta);
159 mapStrategy_->copyThyraIntoTpetra(tY, Y);
163 TEUCHOS_ASSERT(
false);
168void TpetraOperatorWrapper::applyInverse(
const Tpetra::MultiVector<ST,LO,GO,NT>& ,
169 Tpetra::MultiVector<ST,LO,GO,NT>& , Teuchos::ETransp , ST , ST )
const
171 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
172 "TpetraOperatorWrapper::applyInverse not implemented");
176RCP<const Teuchos::Comm<Thyra::Ordinal> >
177TpetraOperatorWrapper::getThyraComm(
const Thyra::LinearOpBase<ST>& inOp)
const
179 RCP<const VectorSpaceBase<ST> > vs = inOp.domain();
181 RCP<const SpmdVectorSpaceBase<ST> > spmd;
182 RCP<const VectorSpaceBase<ST> > current = vs;
183 while(current!=Teuchos::null) {
185 RCP<const ProductVectorSpaceBase<ST> > prod
186 = rcp_dynamic_cast<const ProductVectorSpaceBase<ST> >(current);
189 if(prod==Teuchos::null) {
191 spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<ST> >(current);
196 current = prod->getBlock(0);
199 TEUCHOS_TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error,
200 "TpetraOperatorWrapper requires std::vector space "
201 "blocks to be SPMD std::vector spaces");
203 return spmd->getComm();
270 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
271 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(
getThyraOp());
273 return blkOp->productRange()->numBlocks();
278 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
279 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(
getThyraOp());
281 return blkOp->productDomain()->numBlocks();
286 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
287 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(
getThyraOp());
289 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blkOp->getBlock(i,j));
291 return tOp->getConstTpetraOperator();
virtual void copyThyraIntoTpetra(const RCP< const Thyra::MultiVectorBase< ST > > &thyraX, Tpetra::MultiVector< ST, LO, GO, NT > &tpetraX) const
Copy an Thyra::MultiVectorBase into a Epetra_MultiVector.
RCP< const Thyra::VectorSpaceBase< ST > > rangeSpace_
Range space object.
RCP< const Tpetra::Map< LO, GO, NT > > domainMap_
Pointer to the constructed domain map.
RCP< const Thyra::VectorSpaceBase< ST > > domainSpace_
Domain space object.
virtual void copyTpetraIntoThyra(const Tpetra::MultiVector< ST, LO, GO, NT > &tpetraX, const Teuchos::Ptr< Thyra::MultiVectorBase< ST > > &thyraX) const
Copy an Epetra_MultiVector into a Thyra::MultiVectorBase.
RCP< const Tpetra::Map< LO, GO, NT > > rangeMap_
Pointer to the constructed range map.
virtual int GetBlockRowCount()
Get the number of block rows in this operator.
Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > GetBlock(int i, int j) const
Grab the i,j block.
virtual int GetBlockColCount()
Get the number of block columns in this operator.
const RCP< const Thyra::LinearOpBase< ST > > getThyraOp() const
Return the thyra operator associated with this wrapper.