Teko  Version of the Day
 All Classes Files Functions Variables Pages
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 
69 namespace Teko {
70 namespace Epetra {
71 
72 
73 using namespace Teuchos;
74 using namespace Thyra;
75 
76 DefaultMappingStrategy::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 
88 void DefaultMappingStrategy::copyEpetraIntoThyra(const Epetra_MultiVector& x, const Ptr<Thyra::MultiVectorBase<double> > & thyraVec) const
89 {
90  Teko::Epetra::blockEpetraToThyra(x,thyraVec);
91 }
92 
93 void DefaultMappingStrategy::copyThyraIntoEpetra(const RCP<const Thyra::MultiVectorBase<double> > & thyraVec, Epetra_MultiVector& v) const
94 {
95  Teko::Epetra::blockThyraToEpetra(thyraVec,v);
96 }
97 
98 EpetraOperatorWrapper::EpetraOperatorWrapper()
99 {
100  useTranspose_ = false;
101  mapStrategy_ = Teuchos::null;
102  thyraOp_ = Teuchos::null;
103  comm_ = Teuchos::null;
104  label_ = Teuchos::null;
105 }
106 
107 EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<double> > & thyraOp)
108 {
109  SetOperator(thyraOp);
110 }
111 
112 EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<double> > & thyraOp,const RCP<const MappingStrategy> & mapStrategy)
113  : mapStrategy_(mapStrategy)
114 {
115  SetOperator(thyraOp);
116 }
117 
118 EpetraOperatorWrapper::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 
127 void 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 
137 double EpetraOperatorWrapper::NormInf() const
138 {
139  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
140  "EpetraOperatorWrapper::NormInf not implemated");
141  return 1.0;
142 }
143 
144 int 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 
177 int 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 
186 RCP<const Epetra_Comm>
187 EpetraOperatorWrapper::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 
294 Teuchos::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
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.