Teko Version of the Day
Loading...
Searching...
No Matches
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
62namespace Teko {
63namespace TpetraHelpers {
64
65
66using namespace Teuchos;
67using namespace Thyra;
68
69DefaultMappingStrategy::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
81void 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
86void 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
91TpetraOperatorWrapper::TpetraOperatorWrapper()
92{
93 useTranspose_ = false;
94 mapStrategy_ = Teuchos::null;
95 thyraOp_ = Teuchos::null;
96 comm_ = Teuchos::null;
97 label_ = Teuchos::null;
98}
99
100TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<ST> > & thyraOp)
101{
102 SetOperator(thyraOp);
103}
104
105TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<ST> > & thyraOp,const RCP<const MappingStrategy> & mapStrategy)
106 : mapStrategy_(mapStrategy)
107{
108 SetOperator(thyraOp);
109}
110
111TpetraOperatorWrapper::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
120void 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
130double TpetraOperatorWrapper::NormInf() const
131{
132 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
133 "TpetraOperatorWrapper::NormInf not implemated");
134 return 1.0;
135}
136
137void 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
168void 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
176RCP<const Teuchos::Comm<Thyra::Ordinal> >
177TpetraOperatorWrapper::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
284Teuchos::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
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.
RCP< const Thyra::VectorSpaceBase< ST > > rangeSpace_
Range space object.
RCP< const Tpetra::Map< LO, GO, NT > > domainMap_
Pointer to the constructed domain map.
RCP< const Thyra::VectorSpaceBase< ST > > domainSpace_
Domain space object.
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.
RCP< const Tpetra::Map< LO, GO, NT > > rangeMap_
Pointer to the constructed range map.
virtual int GetBlockRowCount()
Get the number of block rows in this operator.
Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > GetBlock(int i, int j) const
Grab the i,j block.
virtual int GetBlockColCount()
Get the number of block columns in this operator.
const RCP< const Thyra::LinearOpBase< ST > > getThyraOp() const
Return the thyra operator associated with this wrapper.