44 #include "Teko_EpetraOperatorWrapper.hpp"
45 #include "Thyra_SpmdVectorBase.hpp"
46 #include "Thyra_MultiVectorStdOps.hpp"
48 #include "Teuchos_DefaultMpiComm.hpp"
50 #include "Teuchos_DefaultSerialComm.hpp"
51 #include "Thyra_EpetraLinearOp.hpp"
52 #include "Epetra_SerialComm.h"
53 #include "Epetra_Vector.h"
55 #include "Epetra_MpiComm.h"
57 #include "Thyra_EpetraThyraWrappers.hpp"
60 #include "Thyra_BlockedLinearOpBase.hpp"
61 #include "Thyra_ProductVectorSpaceBase.hpp"
62 #include "Thyra_get_Epetra_Operator.hpp"
64 #include "Teko_EpetraThyraConverter.hpp"
65 #include "Teuchos_Ptr.hpp"
70 using namespace Teuchos;
71 using namespace Thyra;
73 DefaultMappingStrategy::DefaultMappingStrategy(
74 const RCP<
const Thyra::LinearOpBase<double> >& thyraOp,
const Epetra_Comm& comm) {
75 RCP<Epetra_Comm> newComm = rcp(comm.Clone());
78 domainSpace_ = thyraOp->domain();
79 rangeSpace_ = thyraOp->range();
81 domainMap_ = Teko::Epetra::thyraVSToEpetraMap(*domainSpace_, newComm);
82 rangeMap_ = Teko::Epetra::thyraVSToEpetraMap(*rangeSpace_, newComm);
86 const Epetra_MultiVector& x,
const Ptr<Thyra::MultiVectorBase<double> >& thyraVec)
const {
87 Teko::Epetra::blockEpetraToThyra(x, thyraVec);
91 const RCP<
const Thyra::MultiVectorBase<double> >& thyraVec, Epetra_MultiVector& v)
const {
92 Teko::Epetra::blockThyraToEpetra(thyraVec, v);
95 EpetraOperatorWrapper::EpetraOperatorWrapper() {
96 useTranspose_ =
false;
97 mapStrategy_ = Teuchos::null;
98 thyraOp_ = Teuchos::null;
99 comm_ = Teuchos::null;
100 label_ = Teuchos::null;
103 EpetraOperatorWrapper::EpetraOperatorWrapper(
104 const RCP<
const Thyra::LinearOpBase<double> >& thyraOp) {
105 SetOperator(thyraOp);
108 EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<
const Thyra::LinearOpBase<double> >& thyraOp,
109 const RCP<const MappingStrategy>& mapStrategy)
110 : mapStrategy_(mapStrategy) {
111 SetOperator(thyraOp);
114 EpetraOperatorWrapper::EpetraOperatorWrapper(
const RCP<const MappingStrategy>& mapStrategy)
115 : mapStrategy_(mapStrategy) {
116 useTranspose_ =
false;
117 thyraOp_ = Teuchos::null;
118 comm_ = Teuchos::null;
119 label_ = Teuchos::null;
122 void EpetraOperatorWrapper::SetOperator(
const RCP<
const LinearOpBase<double> >& thyraOp,
124 useTranspose_ =
false;
126 comm_ = getEpetraComm(*thyraOp);
127 label_ = thyraOp_->description();
128 if (mapStrategy_ == Teuchos::null && buildMap)
129 mapStrategy_ = Teuchos::rcp(
new DefaultMappingStrategy(thyraOp, *comm_));
132 double EpetraOperatorWrapper::NormInf()
const {
133 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
134 "EpetraOperatorWrapper::NormInf not implemated");
138 int EpetraOperatorWrapper::Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const {
139 if (!useTranspose_) {
141 RCP<Thyra::MultiVectorBase<double> > tX;
142 RCP<Thyra::MultiVectorBase<double> > tY;
144 tX = Thyra::createMembers(thyraOp_->domain(), X.NumVectors());
145 tY = Thyra::createMembers(thyraOp_->range(), X.NumVectors());
147 Thyra::assign(tX.ptr(), 0.0);
148 Thyra::assign(tY.ptr(), 0.0);
151 mapStrategy_->copyEpetraIntoThyra(X, tX.ptr());
152 mapStrategy_->copyEpetraIntoThyra(
156 thyraOp_->apply(Thyra::NOTRANS, *tX, tY.ptr(), 1.0, 0.0);
159 mapStrategy_->copyThyraIntoEpetra(tY, Y);
161 TEUCHOS_ASSERT(
false);
167 int EpetraOperatorWrapper::ApplyInverse(
const Epetra_MultiVector& ,
168 Epetra_MultiVector& )
const {
169 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
170 "EpetraOperatorWrapper::ApplyInverse not implemented");
174 RCP<const Epetra_Comm> EpetraOperatorWrapper::getEpetraComm(
175 const Thyra::LinearOpBase<double>& inOp)
const {
176 RCP<const VectorSpaceBase<double> > vs = inOp.domain();
178 RCP<const SpmdVectorSpaceBase<double> > spmd;
179 RCP<const VectorSpaceBase<double> > current = vs;
180 while (current != Teuchos::null) {
182 RCP<const ProductVectorSpaceBase<double> > prod =
183 rcp_dynamic_cast<
const ProductVectorSpaceBase<double> >(current);
186 if (prod == Teuchos::null) {
188 spmd = rcp_dynamic_cast<
const SpmdVectorSpaceBase<double> >(current);
192 current = prod->getBlock(0);
195 TEUCHOS_TEST_FOR_EXCEPTION(spmd == Teuchos::null, std::runtime_error,
196 "EpetraOperatorWrapper requires std::vector space "
197 "blocks to be SPMD std::vector spaces");
199 return Thyra::get_Epetra_Comm(*spmd->getComm());
265 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp =
266 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
268 return blkOp->productRange()->numBlocks();
272 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp =
273 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
275 return blkOp->productDomain()->numBlocks();
279 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp =
280 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
282 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i, j));
virtual int GetBlockColCount()
Get the number of block columns 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.
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 copyThyraIntoEpetra(const RCP< const Thyra::MultiVectorBase< double > > &thyraX, Epetra_MultiVector &epetraX) const
Copy an Thyra::MultiVectorBase into a Epetra_MultiVector.
const RCP< const Thyra::LinearOpBase< double > > getThyraOp() const
Return the thyra operator associated with this wrapper.