45 #include "Teko_EpetraOperatorWrapper.hpp" 46 #include "Thyra_SpmdVectorBase.hpp" 47 #include "Thyra_MultiVectorStdOps.hpp" 49 #include "Teuchos_DefaultMpiComm.hpp" 51 #include "Teuchos_DefaultSerialComm.hpp" 52 #include "Thyra_EpetraLinearOp.hpp" 53 #include "Epetra_SerialComm.h" 54 #include "Epetra_Vector.h" 56 #include "Epetra_MpiComm.h" 58 #include "Thyra_EpetraThyraWrappers.hpp" 61 #include "Thyra_BlockedLinearOpBase.hpp" 62 #include "Thyra_ProductVectorSpaceBase.hpp" 63 #include "Thyra_get_Epetra_Operator.hpp" 65 #include "Teko_EpetraThyraConverter.hpp" 66 #include "Teuchos_Ptr.hpp" 74 using namespace Thyra;
76 DefaultMappingStrategy::DefaultMappingStrategy(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp,
const Epetra_Comm & comm)
78 RCP<Epetra_Comm> newComm = rcp(comm.Clone());
81 domainSpace_ = thyraOp->domain();
82 rangeSpace_ = thyraOp->range();
84 domainMap_ = Teko::Epetra::thyraVSToEpetraMap(*domainSpace_,newComm);
85 rangeMap_ = Teko::Epetra::thyraVSToEpetraMap(*rangeSpace_,newComm);
90 Teko::Epetra::blockEpetraToThyra(x,thyraVec);
95 Teko::Epetra::blockThyraToEpetra(thyraVec,v);
98 EpetraOperatorWrapper::EpetraOperatorWrapper()
100 useTranspose_ =
false;
101 mapStrategy_ = Teuchos::null;
102 thyraOp_ = Teuchos::null;
103 comm_ = Teuchos::null;
104 label_ = Teuchos::null;
107 EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp)
109 SetOperator(thyraOp);
112 EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp,
const RCP<const MappingStrategy> & mapStrategy)
113 : mapStrategy_(mapStrategy)
115 SetOperator(thyraOp);
118 EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<const MappingStrategy> & mapStrategy)
119 : mapStrategy_(mapStrategy)
121 useTranspose_ =
false;
122 thyraOp_ = Teuchos::null;
123 comm_ = Teuchos::null;
124 label_ = Teuchos::null;
127 void EpetraOperatorWrapper::SetOperator(
const RCP<
const LinearOpBase<double> > & thyraOp,
bool buildMap)
129 useTranspose_ =
false;
131 comm_ = getEpetraComm(*thyraOp);
132 label_ = thyraOp_->description();
133 if(mapStrategy_==Teuchos::null && buildMap)
134 mapStrategy_ = Teuchos::rcp(
new DefaultMappingStrategy(thyraOp,*comm_));
137 double EpetraOperatorWrapper::NormInf()
const 139 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
140 "EpetraOperatorWrapper::NormInf not implemated");
144 int EpetraOperatorWrapper::Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 149 RCP<Thyra::MultiVectorBase<double> > tX;
150 RCP<Thyra::MultiVectorBase<double> > tY;
152 tX = Thyra::createMembers(thyraOp_->domain(),X.NumVectors());
153 tY = Thyra::createMembers(thyraOp_->range(),X.NumVectors());
155 Thyra::assign(tX.ptr(),0.0);
156 Thyra::assign(tY.ptr(),0.0);
159 mapStrategy_->copyEpetraIntoThyra(X, tX.ptr());
160 mapStrategy_->copyEpetraIntoThyra(Y, tY.ptr());
163 thyraOp_->apply(Thyra::NOTRANS,*tX,tY.ptr(),1.0,0.0);
166 mapStrategy_->copyThyraIntoEpetra(tY, Y);
170 TEUCHOS_ASSERT(
false);
177 int EpetraOperatorWrapper::ApplyInverse(
const Epetra_MultiVector& X,
178 Epetra_MultiVector& Y)
const 180 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
181 "EpetraOperatorWrapper::ApplyInverse not implemented");
186 RCP<const Epetra_Comm>
187 EpetraOperatorWrapper::getEpetraComm(
const Thyra::LinearOpBase<double>& inOp)
const 189 RCP<const VectorSpaceBase<double> > vs = inOp.domain();
191 RCP<const SpmdVectorSpaceBase<double> > spmd;
192 RCP<const VectorSpaceBase<double> > current = vs;
193 while(current!=Teuchos::null) {
195 RCP<const ProductVectorSpaceBase<double> > prod
196 = rcp_dynamic_cast<
const ProductVectorSpaceBase<double> >(current);
199 if(prod==Teuchos::null) {
201 spmd = rcp_dynamic_cast<
const SpmdVectorSpaceBase<double> >(current);
206 current = prod->getBlock(0);
209 TEUCHOS_TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error,
210 "EpetraOperatorWrapper requires std::vector space " 211 "blocks to be SPMD std::vector spaces");
213 return Thyra::get_Epetra_Comm(*spmd->getComm());
280 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
281 = Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
283 return blkOp->productRange()->numBlocks();
288 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
289 = Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
291 return blkOp->productDomain()->numBlocks();
296 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
297 = Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
299 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i,j));
virtual int GetBlockColCount()
Get the number of block columns in this operator.
const RCP< const Thyra::LinearOpBase< double > > getThyraOp() const
Return the thyra operator associated with this wrapper.
virtual void copyThyraIntoEpetra(const RCP< const Thyra::MultiVectorBase< double > > &thyraX, Epetra_MultiVector &epetraX) const
Copy an Thyra::MultiVectorBase into a Epetra_MultiVector.
Teuchos::RCP< const Epetra_Operator > GetBlock(int i, int j) const
Grab the i,j block.
virtual int GetBlockRowCount()
Get the number of block rows in this operator.
virtual void copyEpetraIntoThyra(const Epetra_MultiVector &epetraX, const Teuchos::Ptr< Thyra::MultiVectorBase< double > > &thyraX) const
Copy an Epetra_MultiVector into a Thyra::MultiVectorBase.