Teko Version of the Day
Loading...
Searching...
No Matches
Teko_EpetraOperatorWrapper.cpp
1// @HEADER
2//
3// ***********************************************************************
4//
5// Teko: A package for block and physics based preconditioning
6// Copyright 2010 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
39//
40// ***********************************************************************
41//
42// @HEADER
43
44
45#include "Teko_EpetraOperatorWrapper.hpp"
46#include "Thyra_SpmdVectorBase.hpp"
47#include "Thyra_MultiVectorStdOps.hpp"
48#ifdef HAVE_MPI
49#include "Teuchos_DefaultMpiComm.hpp"
50#endif
51#include "Teuchos_DefaultSerialComm.hpp"
52#include "Thyra_EpetraLinearOp.hpp"
53#include "Epetra_SerialComm.h"
54#include "Epetra_Vector.h"
55#ifdef HAVE_MPI
56#include "Epetra_MpiComm.h"
57#endif
58#include "Thyra_EpetraThyraWrappers.hpp"
59
60// #include "Thyra_LinearOperator.hpp"
61#include "Thyra_BlockedLinearOpBase.hpp"
62#include "Thyra_ProductVectorSpaceBase.hpp"
63#include "Thyra_get_Epetra_Operator.hpp"
64
65#include "Teko_EpetraThyraConverter.hpp"
66#include "Teuchos_Ptr.hpp"
67
68
69namespace Teko {
70namespace Epetra {
71
72
73using namespace Teuchos;
74using namespace Thyra;
75
76DefaultMappingStrategy::DefaultMappingStrategy(const RCP<const Thyra::LinearOpBase<double> > & thyraOp,const Epetra_Comm & comm)
77{
78 RCP<Epetra_Comm> newComm = rcp(comm.Clone());
79
80 // extract vector spaces from linear operator
81 domainSpace_ = thyraOp->domain();
82 rangeSpace_ = thyraOp->range();
83
84 domainMap_ = Teko::Epetra::thyraVSToEpetraMap(*domainSpace_,newComm);
85 rangeMap_ = Teko::Epetra::thyraVSToEpetraMap(*rangeSpace_,newComm);
86}
87
88void DefaultMappingStrategy::copyEpetraIntoThyra(const Epetra_MultiVector& x, const Ptr<Thyra::MultiVectorBase<double> > & thyraVec) const
89{
90 Teko::Epetra::blockEpetraToThyra(x,thyraVec);
91}
92
93void DefaultMappingStrategy::copyThyraIntoEpetra(const RCP<const Thyra::MultiVectorBase<double> > & thyraVec, Epetra_MultiVector& v) const
94{
95 Teko::Epetra::blockThyraToEpetra(thyraVec,v);
96}
97
98EpetraOperatorWrapper::EpetraOperatorWrapper()
99{
100 useTranspose_ = false;
101 mapStrategy_ = Teuchos::null;
102 thyraOp_ = Teuchos::null;
103 comm_ = Teuchos::null;
104 label_ = Teuchos::null;
105}
106
107EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<double> > & thyraOp)
108{
109 SetOperator(thyraOp);
110}
111
112EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<double> > & thyraOp,const RCP<const MappingStrategy> & mapStrategy)
113 : mapStrategy_(mapStrategy)
114{
115 SetOperator(thyraOp);
116}
117
118EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const MappingStrategy> & mapStrategy)
119 : mapStrategy_(mapStrategy)
120{
121 useTranspose_ = false;
122 thyraOp_ = Teuchos::null;
123 comm_ = Teuchos::null;
124 label_ = Teuchos::null;
125}
126
127void EpetraOperatorWrapper::SetOperator(const RCP<const LinearOpBase<double> > & thyraOp,bool buildMap)
128{
129 useTranspose_ = false;
130 thyraOp_ = thyraOp;
131 comm_ = getEpetraComm(*thyraOp);
132 label_ = thyraOp_->description();
133 if(mapStrategy_==Teuchos::null && buildMap)
134 mapStrategy_ = Teuchos::rcp(new DefaultMappingStrategy(thyraOp,*comm_));
135}
136
137double EpetraOperatorWrapper::NormInf() const
138{
139 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
140 "EpetraOperatorWrapper::NormInf not implemated");
141 return 1.0;
142}
143
144int EpetraOperatorWrapper::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
145{
146 if (!useTranspose_)
147 {
148 // allocate space for each vector
149 RCP<Thyra::MultiVectorBase<double> > tX;
150 RCP<Thyra::MultiVectorBase<double> > tY;
151
152 tX = Thyra::createMembers(thyraOp_->domain(),X.NumVectors());
153 tY = Thyra::createMembers(thyraOp_->range(),X.NumVectors());
154
155 Thyra::assign(tX.ptr(),0.0);
156 Thyra::assign(tY.ptr(),0.0);
157
158 // copy epetra X into thyra X
159 mapStrategy_->copyEpetraIntoThyra(X, tX.ptr());
160 mapStrategy_->copyEpetraIntoThyra(Y, tY.ptr()); // if this matrix isn't block square, this probably won't work!
161
162 // perform matrix vector multiplication
163 thyraOp_->apply(Thyra::NOTRANS,*tX,tY.ptr(),1.0,0.0);
164
165 // copy thyra Y into epetra Y
166 mapStrategy_->copyThyraIntoEpetra(tY, Y);
167 }
168 else
169 {
170 TEUCHOS_ASSERT(false);
171 }
172
173 return 0;
174}
175
176
177int EpetraOperatorWrapper::ApplyInverse(const Epetra_MultiVector& /* X */,
178 Epetra_MultiVector& /* Y */) const
179{
180 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
181 "EpetraOperatorWrapper::ApplyInverse not implemented");
182 return 1;
183}
184
185
186RCP<const Epetra_Comm>
187EpetraOperatorWrapper::getEpetraComm(const Thyra::LinearOpBase<double>& inOp) const
188{
189 RCP<const VectorSpaceBase<double> > vs = inOp.domain();
190
191 RCP<const SpmdVectorSpaceBase<double> > spmd;
192 RCP<const VectorSpaceBase<double> > current = vs;
193 while(current!=Teuchos::null) {
194 // try to cast to a product vector space first
195 RCP<const ProductVectorSpaceBase<double> > prod
196 = rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(current);
197
198 // figure out what type it is
199 if(prod==Teuchos::null) {
200 // hopfully this is a SPMD vector space
201 spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(current);
202
203 break;
204 }
205 else // get first convenient vector space
206 current = prod->getBlock(0);
207 }
208
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");
212
213 return Thyra::get_Epetra_Comm(*spmd->getComm());
214/*
215 const Thyra::ConstLinearOperator<double> thyraOp = rcpFromRef(inOp);
216
217 RCP<Epetra_Comm> rtn;
218 // VectorSpace<double> vs = thyraOp.domain().getBlock(0);
219 RCP<const VectorSpaceBase<double> > vs = thyraOp.domain().getBlock(0).constPtr();
220
221 // search for an SpmdVectorSpaceBase object
222 RCP<const SpmdVectorSpaceBase<double> > spmd;
223 RCP<const VectorSpaceBase<double> > current = vs;
224 while(current!=Teuchos::null) {
225 // try to cast to a product vector space first
226 RCP<const ProductVectorSpaceBase<double> > prod
227 = rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(current);
228
229 // figure out what type it is
230 if(prod==Teuchos::null) {
231 // hopfully this is a SPMD vector space
232 spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(current);
233
234 break;
235 }
236 else {
237 // get first convenient vector space
238 current = prod->getBlock(0);
239 }
240 }
241
242 TEUCHOS_TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error,
243 "EpetraOperatorWrapper requires std::vector space "
244 "blocks to be SPMD std::vector spaces");
245
246 const SerialComm<Thyra::Ordinal>* serialComm
247 = dynamic_cast<const SerialComm<Thyra::Ordinal>*>(spmd->getComm().get());
248
249#ifdef HAVE_MPI
250 const MpiComm<Thyra::Ordinal>* mpiComm
251 = dynamic_cast<const MpiComm<Thyra::Ordinal>*>(spmd->getComm().get());
252
253 TEUCHOS_TEST_FOR_EXCEPTION(mpiComm==0 && serialComm==0, std::runtime_error,
254 "SPMD std::vector space has a communicator that is "
255 "neither a serial comm nor an MPI comm");
256
257 if (mpiComm != 0)
258 {
259 rtn = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
260 }
261 else
262 {
263 rtn = rcp(new Epetra_SerialComm());
264 }
265#else
266 TEUCHOS_TEST_FOR_EXCEPTION(serialComm==0, std::runtime_error,
267 "SPMD std::vector space has a communicator that is "
268 "neither a serial comm nor an MPI comm");
269 rtn = rcp(new Epetra_SerialComm());
270
271#endif
272
273 TEUCHOS_TEST_FOR_EXCEPTION(rtn.get()==0, std::runtime_error, "null communicator created");
274 return rtn;
275*/
276}
277
279{
280 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
281 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
282
283 return blkOp->productRange()->numBlocks();
284}
285
287{
288 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
289 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
290
291 return blkOp->productDomain()->numBlocks();
292}
293
294Teuchos::RCP<const Epetra_Operator> EpetraOperatorWrapper::GetBlock(int i,int j) const
295{
296 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
297 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
298
299 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i,j));
300}
301
302} // namespace Epetra
303} // namespace Teko
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.