Teko Version of the Day
Loading...
Searching...
No Matches
Teko_TpetraBlockPreconditioner.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_TpetraBlockPreconditioner.hpp"
48#include "Teko_Preconditioner.hpp"
49
50// Thyra includes
51#include "Thyra_DefaultLinearOpSource.hpp"
52#include "Thyra_TpetraLinearOp.hpp"
53
54// Teuchos includes
55#include "Teuchos_Time.hpp"
56
57// Teko includes
58#include "Teko_TpetraBasicMappingStrategy.hpp"
59#include "Teko_Utilities.hpp"
60
61namespace Teko {
62namespace TpetraHelpers {
63
64using Teuchos::RCP;
65using Teuchos::rcp;
66using Teuchos::rcpFromRef;
67using Teuchos::rcp_dynamic_cast;
68
75TpetraBlockPreconditioner::TpetraBlockPreconditioner(const Teuchos::RCP<const PreconditionerFactory> & bfp)
76 : preconFactory_(bfp), firstBuildComplete_(false)
77{
78}
79
81{
82 if((not clearOld) && preconObj_!=Teuchos::null)
83 return;
84 preconObj_ = preconFactory_->createPrec();
85}
86
99void TpetraBlockPreconditioner::buildPreconditioner(const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & A,bool clear)
100{
101 Teko_DEBUG_SCOPE("TBP::buildPreconditioner",10);
102
103 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
104 RCP<const Thyra::LinearOpBase<ST> > thyraA = extractLinearOp(A);
105
106 // set the mapping strategy
107 // SetMapStrategy(rcp(new InverseMappingStrategy(eow->getMapStrategy())));
108 SetMapStrategy(rcp(new InverseMappingStrategy(extractMappingStrategy(A))));
109
110 // build preconObj_
111 initPreconditioner(clear);
112
113 // actually build the preconditioner
114 RCP<const Thyra::LinearOpSourceBase<ST> > lOpSrc = Thyra::defaultLinearOpSource(thyraA);
115 preconFactory_->initializePrec(lOpSrc,&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
116
117 // extract preconditioner operator
118 RCP<const Thyra::LinearOpBase<ST> > preconditioner = preconObj_->getUnspecifiedPrecOp();
119
120 SetOperator(preconditioner,false);
121
122 firstBuildComplete_ = true;
123
124 TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
125 TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
126 TEUCHOS_ASSERT(getMapStrategy()!=Teuchos::null);
127 TEUCHOS_ASSERT(firstBuildComplete_==true);
128}
129
143void TpetraBlockPreconditioner::buildPreconditioner(const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & A,const Tpetra::MultiVector<ST,LO,GO,NT> & tpetra_mv,bool clear)
144{
145 Teko_DEBUG_SCOPE("TBP::buildPreconditioner - with solution",10);
146
147 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
148 RCP<const Thyra::LinearOpBase<ST> > thyraA = extractLinearOp(A);
149
150 // set the mapping strategy
151 SetMapStrategy(rcp(new InverseMappingStrategy(extractMappingStrategy(A))));
152
153 TEUCHOS_ASSERT(getMapStrategy()!=Teuchos::null);
154
155 // build the thyra version of the source multivector
156 RCP<Thyra::MultiVectorBase<ST> > thyra_mv = Thyra::createMembers(thyraA->range(),tpetra_mv.getNumVectors());
157 getMapStrategy()->copyTpetraIntoThyra(tpetra_mv,thyra_mv.ptr());
158
159 // build preconObj_
160 initPreconditioner(clear);
161
162 // actually build the preconditioner
163 preconFactory_->initializePrec(Thyra::defaultLinearOpSource(thyraA),thyra_mv,&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
164 RCP<const Thyra::LinearOpBase<ST> > preconditioner = preconObj_->getUnspecifiedPrecOp();
165
166 SetOperator(preconditioner,false);
167
168 firstBuildComplete_ = true;
169
170 TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
171 TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
172 TEUCHOS_ASSERT(getMapStrategy()!=Teuchos::null);
173 TEUCHOS_ASSERT(firstBuildComplete_==true);
174}
175
189void TpetraBlockPreconditioner::rebuildPreconditioner(const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & A)
190{
191 Teko_DEBUG_SCOPE("TBP::rebuildPreconditioner",10);
192
193 // if the preconditioner hasn't been built yet, rebuild from scratch
194 if(not firstBuildComplete_) {
195 buildPreconditioner(A,false);
196 return;
197 }
198 Teko_DEBUG_EXPR(Teuchos::Time timer(""));
199
200 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
201 Teko_DEBUG_EXPR(timer.start(true));
202 RCP<const Thyra::LinearOpBase<ST> > thyraA = extractLinearOp(A);
203 Teko_DEBUG_EXPR(timer.stop());
204 Teko_DEBUG_MSG("TBP::rebuild get thyraop time = " << timer.totalElapsedTime(),2);
205
206 // reinitialize the preconditioner
207 Teko_DEBUG_EXPR(timer.start(true));
208 preconFactory_->initializePrec(Thyra::defaultLinearOpSource(thyraA),&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
209 RCP<const Thyra::LinearOpBase<ST> > preconditioner = preconObj_->getUnspecifiedPrecOp();
210 Teko_DEBUG_EXPR(timer.stop());
211 Teko_DEBUG_MSG("TBP::rebuild initialize prec time = " << timer.totalElapsedTime(),2);
212
213 Teko_DEBUG_EXPR(timer.start(true));
214 SetOperator(preconditioner,false);
215 Teko_DEBUG_EXPR(timer.stop());
216 Teko_DEBUG_MSG("TBP::rebuild set operator time = " << timer.totalElapsedTime(),2);
217
218 TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
219 TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
220 TEUCHOS_ASSERT(firstBuildComplete_==true);
221}
222
236void TpetraBlockPreconditioner::rebuildPreconditioner(const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & A,const Tpetra::MultiVector<ST,LO,GO,NT> & tpetra_mv)
237{
238 Teko_DEBUG_SCOPE("TBP::rebuildPreconditioner - with solution",10);
239
240 // if the preconditioner hasn't been built yet, rebuild from scratch
241 if(not firstBuildComplete_) {
242 buildPreconditioner(A,tpetra_mv,false);
243 return;
244 }
245 Teko_DEBUG_EXPR(Teuchos::Time timer(""));
246
247 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
248 Teko_DEBUG_EXPR(timer.start(true));
249 RCP<const Thyra::LinearOpBase<ST> > thyraA = extractLinearOp(A);
250 Teko_DEBUG_EXPR(timer.stop());
251 Teko_DEBUG_MSG("TBP::rebuild get thyraop time = " << timer.totalElapsedTime(),2);
252
253 // build the thyra version of the source multivector
254 Teko_DEBUG_EXPR(timer.start(true));
255 RCP<Thyra::MultiVectorBase<ST> > thyra_mv = Thyra::createMembers(thyraA->range(),tpetra_mv.getNumVectors());
256 getMapStrategy()->copyTpetraIntoThyra(tpetra_mv,thyra_mv.ptr());
257 Teko_DEBUG_EXPR(timer.stop());
258 Teko_DEBUG_MSG("TBP::rebuild vector copy time = " << timer.totalElapsedTime(),2);
259
260 // reinitialize the preconditioner
261 Teko_DEBUG_EXPR(timer.start(true));
262 preconFactory_->initializePrec(Thyra::defaultLinearOpSource(thyraA),thyra_mv,&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
263 RCP<const Thyra::LinearOpBase<ST> > preconditioner = preconObj_->getUnspecifiedPrecOp();
264 Teko_DEBUG_EXPR(timer.stop());
265 Teko_DEBUG_MSG("TBP::rebuild initialize prec time = " << timer.totalElapsedTime(),2);
266
267 Teko_DEBUG_EXPR(timer.start(true));
268 SetOperator(preconditioner,false);
269 Teko_DEBUG_EXPR(timer.stop());
270 Teko_DEBUG_MSG("TBP::rebuild set operator time = " << timer.totalElapsedTime(),2);
271
272 TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
273 TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
274 TEUCHOS_ASSERT(firstBuildComplete_==true);
275}
276
286Teuchos::RCP<PreconditionerState> TpetraBlockPreconditioner::getPreconditionerState()
287{
288 Teuchos::RCP<Preconditioner> bp = rcp_dynamic_cast<Preconditioner>(preconObj_);
289
290 if(bp!=Teuchos::null)
291 return bp->getStateObject();
292
293 return Teuchos::null;
294}
295
305Teuchos::RCP<const PreconditionerState> TpetraBlockPreconditioner::getPreconditionerState() const
306{
307 Teuchos::RCP<const Preconditioner> bp = rcp_dynamic_cast<const Preconditioner>(preconObj_);
308
309 if(bp!=Teuchos::null)
310 return bp->getStateObject();
311
312 return Teuchos::null;
313}
314
315Teuchos::RCP<const Thyra::LinearOpBase<ST> > TpetraBlockPreconditioner::extractLinearOp(const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & A) const
316{
317 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
318 const RCP<const TpetraOperatorWrapper> & tow = rcp_dynamic_cast<const TpetraOperatorWrapper>(A);
319
320 // if it is an EpetraOperatorWrapper, then get the Thyra operator
321 if(tow!=Teuchos::null)
322 return tow->getThyraOp();
323
324 // otherwise wrap it up as a thyra operator
325 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(A->getDomainMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(A->getRangeMap()),A);
326}
327
328Teuchos::RCP<const MappingStrategy> TpetraBlockPreconditioner::extractMappingStrategy(const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & A) const
329{
330 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
331 const RCP<const TpetraOperatorWrapper> & tow = rcp_dynamic_cast<const TpetraOperatorWrapper>(A);
332
333 // if it is an EpetraOperatorWrapper, then get the Thyra operator
334 if(tow!=Teuchos::null)
335 return tow->getMapStrategy();
336
337 // otherwise wrap it up as a thyra operator
338 RCP<const Tpetra::Map<LO,GO,NT> > range = A->getRangeMap();
339 RCP<const Tpetra::Map<LO,GO,NT> > domain = A->getDomainMap();
340 return rcp(new BasicMappingStrategy(range,domain,*Thyra::convertTpetraToThyraComm(range->getComm())));
341}
342
343} // end namespace Epetra
344} // end namespace Teko
Flip a mapping strategy object around to give the "inverse" mapping strategy.
virtual Teuchos::RCP< PreconditionerState > getPreconditionerState()
virtual void initPreconditioner(bool clearOld=false)
Build the underlying data structure for the preconditioner.
virtual void buildPreconditioner(const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &A, bool clear=true)
Build this preconditioner from an Epetra_Operator passed in to this object.
virtual void rebuildPreconditioner(const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &A)
Rebuild this preconditioner from an Epetra_Operator passed in this to object.
const RCP< const MappingStrategy > getMapStrategy() const
Get the mapping strategy for this wrapper (translate between Thyra and Epetra)
const RCP< const Thyra::LinearOpBase< ST > > getThyraOp() const
Return the thyra operator associated with this wrapper.