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"
73using namespace Teuchos;
76DefaultMappingStrategy::DefaultMappingStrategy(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp,
const Epetra_Comm & comm)
78 RCP<Epetra_Comm> newComm = rcp(comm.Clone());
90 Teko::Epetra::blockEpetraToThyra(x,thyraVec);
95 Teko::Epetra::blockThyraToEpetra(thyraVec,v);
98EpetraOperatorWrapper::EpetraOperatorWrapper()
100 useTranspose_ =
false;
101 mapStrategy_ = Teuchos::null;
102 thyraOp_ = Teuchos::null;
103 comm_ = Teuchos::null;
104 label_ = Teuchos::null;
107EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp)
109 SetOperator(thyraOp);
112EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<double> > & thyraOp,
const RCP<const MappingStrategy> & mapStrategy)
113 : mapStrategy_(mapStrategy)
115 SetOperator(thyraOp);
118EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<const MappingStrategy> & mapStrategy)
119 : mapStrategy_(mapStrategy)
121 useTranspose_ =
false;
122 thyraOp_ = Teuchos::null;
123 comm_ = Teuchos::null;
124 label_ = Teuchos::null;
127void 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_));
137double EpetraOperatorWrapper::NormInf()
const
139 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
140 "EpetraOperatorWrapper::NormInf not implemated");
144int 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);
177int EpetraOperatorWrapper::ApplyInverse(
const Epetra_MultiVector& ,
178 Epetra_MultiVector& )
const
180 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
181 "EpetraOperatorWrapper::ApplyInverse not implemented");
186RCP<const Epetra_Comm>
187EpetraOperatorWrapper::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));
RCP< const Thyra::VectorSpaceBase< double > > rangeSpace_
Range space object.
RCP< const Epetra_Map > domainMap_
Pointer to the constructed domain map.
RCP< const Thyra::VectorSpaceBase< double > > domainSpace_
Domain space object.
RCP< const Epetra_Map > rangeMap_
Pointer to the constructed range map.
virtual void copyEpetraIntoThyra(const Epetra_MultiVector &epetraX, const Teuchos::Ptr< Thyra::MultiVectorBase< double > > &thyraX) const
Copy an Epetra_MultiVector into a Thyra::MultiVectorBase.
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 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.