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" 63 namespace TpetraHelpers {
67 using namespace Thyra;
69 DefaultMappingStrategy::DefaultMappingStrategy(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp,
const Teuchos::Comm<Thyra::Ordinal> & comm)
71 RCP<Teuchos::Comm<Thyra::Ordinal> > newComm = comm.duplicate();
74 domainSpace_ = thyraOp->domain();
75 rangeSpace_ = thyraOp->range();
77 domainMap_ = Teko::TpetraHelpers::thyraVSToTpetraMap(*domainSpace_,newComm);
78 rangeMap_ = Teko::TpetraHelpers::thyraVSToTpetraMap(*rangeSpace_,newComm);
83 Teko::TpetraHelpers::blockTpetraToThyra(x,thyraVec);
88 Teko::TpetraHelpers::blockThyraToTpetra(thyraVec,v);
91 TpetraOperatorWrapper::TpetraOperatorWrapper()
93 useTranspose_ =
false;
94 mapStrategy_ = Teuchos::null;
95 thyraOp_ = Teuchos::null;
96 comm_ = Teuchos::null;
97 label_ = Teuchos::null;
100 TpetraOperatorWrapper::TpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<ST> > & thyraOp)
102 SetOperator(thyraOp);
105 TpetraOperatorWrapper::TpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<ST> > & thyraOp,
const RCP<const MappingStrategy> & mapStrategy)
106 : mapStrategy_(mapStrategy)
108 SetOperator(thyraOp);
111 TpetraOperatorWrapper::TpetraOperatorWrapper(
const RCP<const MappingStrategy> & mapStrategy)
112 : mapStrategy_(mapStrategy)
114 useTranspose_ =
false;
115 thyraOp_ = Teuchos::null;
116 comm_ = Teuchos::null;
117 label_ = Teuchos::null;
120 void 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)
130 double TpetraOperatorWrapper::NormInf()
const 132 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
133 "TpetraOperatorWrapper::NormInf not implemated");
137 void TpetraOperatorWrapper::apply(
const Tpetra::MultiVector<ST,LO,GO,NT>& X, Tpetra::MultiVector<ST,LO,GO,NT>& Y,Teuchos::ETransp mode,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);
168 void TpetraOperatorWrapper::applyInverse(
const Tpetra::MultiVector<ST,LO,GO,NT>& X,
169 Tpetra::MultiVector<ST,LO,GO,NT>& Y, Teuchos::ETransp mode, ST alpha, ST beta)
const 171 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
172 "TpetraOperatorWrapper::applyInverse not implemented");
176 RCP<const Teuchos::Comm<Thyra::Ordinal> >
177 TpetraOperatorWrapper::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();
default mapping strategy for the basic TpetraOperatorWrapper
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.
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.
Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > GetBlock(int i, int j) const
Grab the i,j block.
virtual int GetBlockRowCount()
Get the number of block rows in this operator.
virtual int GetBlockColCount()
Get the number of block columns in this operator.