Teko  Version of the Day
Teko_TpetraOperatorWrapper.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_TpetraOperatorWrapper.hpp"
46 #include "Thyra_SpmdVectorBase.hpp"
47 #include "Thyra_MultiVectorStdOps.hpp"
48 #include "Teuchos_DefaultSerialComm.hpp"
49 #include "Thyra_TpetraLinearOp.hpp"
50 #include "Tpetra_Vector.hpp"
51 #include "Thyra_TpetraThyraWrappers.hpp"
52 
53 // #include "Thyra_LinearOperator.hpp"
54 #include "Thyra_BlockedLinearOpBase.hpp"
55 #include "Thyra_ProductVectorSpaceBase.hpp"
56 #include "Thyra_TpetraLinearOp.hpp"
57 
58 #include "Teko_TpetraThyraConverter.hpp"
59 #include "Teuchos_Ptr.hpp"
60 
61 
62 namespace Teko {
63 namespace TpetraHelpers {
64 
65 
66 using namespace Teuchos;
67 using namespace Thyra;
68 
69 DefaultMappingStrategy::DefaultMappingStrategy(const RCP<const Thyra::LinearOpBase<double> > & thyraOp,const Teuchos::Comm<Thyra::Ordinal> & comm)
70 {
71  RCP<Teuchos::Comm<Thyra::Ordinal> > newComm = comm.duplicate();
72 
73  // extract vector spaces from linear operator
74  domainSpace_ = thyraOp->domain();
75  rangeSpace_ = thyraOp->range();
76 
77  domainMap_ = Teko::TpetraHelpers::thyraVSToTpetraMap(*domainSpace_,newComm);
78  rangeMap_ = Teko::TpetraHelpers::thyraVSToTpetraMap(*rangeSpace_,newComm);
79 }
80 
81 void DefaultMappingStrategy::copyTpetraIntoThyra(const Tpetra::MultiVector<ST,LO,GO,NT>& x, const Ptr<Thyra::MultiVectorBase<ST> > & thyraVec) const
82 {
83  Teko::TpetraHelpers::blockTpetraToThyra(x,thyraVec);
84 }
85 
86 void DefaultMappingStrategy::copyThyraIntoTpetra(const RCP<const Thyra::MultiVectorBase<ST> > & thyraVec, Tpetra::MultiVector<ST,LO,GO,NT>& v) const
87 {
88  Teko::TpetraHelpers::blockThyraToTpetra(thyraVec,v);
89 }
90 
91 TpetraOperatorWrapper::TpetraOperatorWrapper()
92 {
93  useTranspose_ = false;
94  mapStrategy_ = Teuchos::null;
95  thyraOp_ = Teuchos::null;
96  comm_ = Teuchos::null;
97  label_ = Teuchos::null;
98 }
99 
100 TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<ST> > & thyraOp)
101 {
102  SetOperator(thyraOp);
103 }
104 
105 TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<ST> > & thyraOp,const RCP<const MappingStrategy> & mapStrategy)
106  : mapStrategy_(mapStrategy)
107 {
108  SetOperator(thyraOp);
109 }
110 
111 TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const MappingStrategy> & mapStrategy)
112  : mapStrategy_(mapStrategy)
113 {
114  useTranspose_ = false;
115  thyraOp_ = Teuchos::null;
116  comm_ = Teuchos::null;
117  label_ = Teuchos::null;
118 }
119 
120 void TpetraOperatorWrapper::SetOperator(const RCP<const LinearOpBase<ST> > & thyraOp,bool buildMap)
121 {
122  useTranspose_ = false;
123  thyraOp_ = thyraOp;
124  comm_ = getThyraComm(*thyraOp);
125  label_ = thyraOp_->description();
126  if(mapStrategy_==Teuchos::null && buildMap)
127  mapStrategy_ = Teuchos::rcp(new DefaultMappingStrategy(thyraOp,*comm_));
128 }
129 
130 double TpetraOperatorWrapper::NormInf() const
131 {
132  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
133  "TpetraOperatorWrapper::NormInf not implemated");
134  return 1.0;
135 }
136 
137 void TpetraOperatorWrapper::apply(const Tpetra::MultiVector<ST,LO,GO,NT>& X, Tpetra::MultiVector<ST,LO,GO,NT>& Y,Teuchos::ETransp mode,ST alpha, ST beta) const
138 {
139  if (!useTranspose_)
140  {
141  // allocate space for each vector
142  RCP<Thyra::MultiVectorBase<ST> > tX;
143  RCP<Thyra::MultiVectorBase<ST> > tY;
144 
145  tX = Thyra::createMembers(thyraOp_->domain(),X.getNumVectors());
146  tY = Thyra::createMembers(thyraOp_->range(),X.getNumVectors());
147 
148  Thyra::assign(tX.ptr(),0.0);
149  Thyra::assign(tY.ptr(),0.0);
150 
151  // copy epetra X into thyra X
152  mapStrategy_->copyTpetraIntoThyra(X, tX.ptr());
153  mapStrategy_->copyTpetraIntoThyra(Y, tY.ptr()); // if this matrix isn't block square, this probably won't work!
154 
155  // perform matrix vector multiplication
156  thyraOp_->apply(Thyra::NOTRANS,*tX,tY.ptr(),alpha,beta);
157 
158  // copy thyra Y into epetra Y
159  mapStrategy_->copyThyraIntoTpetra(tY, Y);
160  }
161  else
162  {
163  TEUCHOS_ASSERT(false);
164  }
165 }
166 
167 
168 void TpetraOperatorWrapper::applyInverse(const Tpetra::MultiVector<ST,LO,GO,NT>& X,
169  Tpetra::MultiVector<ST,LO,GO,NT>& Y, Teuchos::ETransp mode, ST alpha, ST beta) const
170 {
171  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
172  "TpetraOperatorWrapper::applyInverse not implemented");
173 }
174 
175 
176 RCP<const Teuchos::Comm<Thyra::Ordinal> >
177 TpetraOperatorWrapper::getThyraComm(const Thyra::LinearOpBase<ST>& inOp) const
178 {
179  RCP<const VectorSpaceBase<ST> > vs = inOp.domain();
180 
181  RCP<const SpmdVectorSpaceBase<ST> > spmd;
182  RCP<const VectorSpaceBase<ST> > current = vs;
183  while(current!=Teuchos::null) {
184  // try to cast to a product vector space first
185  RCP<const ProductVectorSpaceBase<ST> > prod
186  = rcp_dynamic_cast<const ProductVectorSpaceBase<ST> >(current);
187 
188  // figure out what type it is
189  if(prod==Teuchos::null) {
190  // hopfully this is a SPMD vector space
191  spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<ST> >(current);
192 
193  break;
194  }
195  else // get first convenient vector space
196  current = prod->getBlock(0);
197  }
198 
199  TEUCHOS_TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error,
200  "TpetraOperatorWrapper requires std::vector space "
201  "blocks to be SPMD std::vector spaces");
202 
203  return spmd->getComm();
204 /*
205  const Thyra::ConstLinearOperator<double> thyraOp = rcpFromRef(inOp);
206 
207  RCP<Epetra_Comm> rtn;
208  // VectorSpace<double> vs = thyraOp.domain().getBlock(0);
209  RCP<const VectorSpaceBase<double> > vs = thyraOp.domain().getBlock(0).constPtr();
210 
211  // search for an SpmdVectorSpaceBase object
212  RCP<const SpmdVectorSpaceBase<double> > spmd;
213  RCP<const VectorSpaceBase<double> > current = vs;
214  while(current!=Teuchos::null) {
215  // try to cast to a product vector space first
216  RCP<const ProductVectorSpaceBase<double> > prod
217  = rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(current);
218 
219  // figure out what type it is
220  if(prod==Teuchos::null) {
221  // hopfully this is a SPMD vector space
222  spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(current);
223 
224  break;
225  }
226  else {
227  // get first convenient vector space
228  current = prod->getBlock(0);
229  }
230  }
231 
232  TEUCHOS_TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error,
233  "TpetraOperatorWrapper requires std::vector space "
234  "blocks to be SPMD std::vector spaces");
235 
236  const SerialComm<Thyra::Ordinal>* serialComm
237  = dynamic_cast<const SerialComm<Thyra::Ordinal>*>(spmd->getComm().get());
238 
239 #ifdef HAVE_MPI
240  const MpiComm<Thyra::Ordinal>* mpiComm
241  = dynamic_cast<const MpiComm<Thyra::Ordinal>*>(spmd->getComm().get());
242 
243  TEUCHOS_TEST_FOR_EXCEPTION(mpiComm==0 && serialComm==0, std::runtime_error,
244  "SPMD std::vector space has a communicator that is "
245  "neither a serial comm nor an MPI comm");
246 
247  if (mpiComm != 0)
248  {
249  rtn = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
250  }
251  else
252  {
253  rtn = rcp(new Epetra_SerialComm());
254  }
255 #else
256  TEUCHOS_TEST_FOR_EXCEPTION(serialComm==0, std::runtime_error,
257  "SPMD std::vector space has a communicator that is "
258  "neither a serial comm nor an MPI comm");
259  rtn = rcp(new Epetra_SerialComm());
260 
261 #endif
262 
263  TEUCHOS_TEST_FOR_EXCEPTION(rtn.get()==0, std::runtime_error, "null communicator created");
264  return rtn;
265 */
266 }
267 
269 {
270  const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
271  = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());
272 
273  return blkOp->productRange()->numBlocks();
274 }
275 
277 {
278  const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
279  = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());
280 
281  return blkOp->productDomain()->numBlocks();
282 }
283 
284 Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > TpetraOperatorWrapper::GetBlock(int i,int j) const
285 {
286  const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
287  = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());
288 
289  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blkOp->getBlock(i,j));
290 
291  return tOp->getConstTpetraOperator();
292 }
293 
294 } // namespace TpetraHelpers
295 } // namespace Teko
const RCP< const Thyra::LinearOpBase< ST > > getThyraOp() const
Return the thyra operator associated with this wrapper.
Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > GetBlock(int i, int j) const
Grab the i,j block.
virtual void copyTpetraIntoThyra(const Tpetra::MultiVector< ST, LO, GO, NT > &tpetraX, const Teuchos::Ptr< Thyra::MultiVectorBase< ST > > &thyraX) const
Copy an Epetra_MultiVector into a Thyra::MultiVectorBase.
virtual void copyThyraIntoTpetra(const RCP< const Thyra::MultiVectorBase< ST > > &thyraX, Tpetra::MultiVector< ST, LO, GO, NT > &tpetraX) const
Copy an Thyra::MultiVectorBase into a Epetra_MultiVector.
virtual int GetBlockRowCount()
Get the number of block rows in this operator.
virtual int GetBlockColCount()
Get the number of block columns in this operator.