Stratimikos Version of the Day
Loading...
Searching...
No Matches
Thyra_IfpackPreconditionerFactory.cpp
1// @HEADER
2// ***********************************************************************
3//
4// Stratimikos: Thyra-based strategies for linear solvers
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Thyra_IfpackPreconditionerFactory.hpp"
43#include "Thyra_EpetraOperatorViewExtractorStd.hpp"
44#include "Thyra_EpetraLinearOp.hpp"
45#include "Thyra_DefaultPreconditioner.hpp"
46#include "Ifpack_ValidParameters.h"
47#include "Ifpack_Preconditioner.h"
48#include "Ifpack.h"
49#include "Epetra_RowMatrix.h"
50#include "Teuchos_TimeMonitor.hpp"
51#include "Teuchos_dyn_cast.hpp"
52#include "Teuchos_implicit_cast.hpp"
53#include "Teuchos_StandardParameterEntryValidators.hpp"
54#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
55#include "Teuchos_ValidatorXMLConverterDB.hpp"
56#include "Teuchos_StaticSetupMacro.hpp"
57
58
59namespace {
60
61Teuchos::RCP<Teuchos::Time> overallTimer, creationTimer, factorizationTimer;
62
63const std::string Ifpack_name = "Ifpack";
64
65const std::string IfpackSettings_name = "Ifpack Settings";
66
67const std::string PrecType_name = "Prec Type";
68Teuchos::RCP<Teuchos::StringToIntegralParameterEntryValidator<Ifpack::EPrecType> >
69precTypeValidator; // Will be setup below!
70const Ifpack::EPrecType PrecType_default = Ifpack::ILU;
71const std::string PrecTypeName_default = Ifpack::precTypeNames[PrecType_default];
72
73const std::string Overlap_name = "Overlap";
74const int Overlap_default = 0;
75
76
77TEUCHOS_STATIC_SETUP()
78{
79 TEUCHOS_ADD_STRINGTOINTEGRALVALIDATOR_CONVERTER(Ifpack::EPrecType);
80}
81
82
83} // namespace
84
85namespace Thyra {
86
87// Constructors/initializers/accessors
88
90 :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
91 ,precType_(PrecType_default)
92 ,overlap_(Overlap_default)
93{
94 initializeTimers();
95 getValidParameters(); // Make sure validators get created!
96}
97
98// Overridden from PreconditionerFactoryBase
99
101 const LinearOpSourceBase<double> &fwdOpSrc
102 ) const
103{
104 using Teuchos::outArg;
105 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
106 EOpTransp epetraFwdOpTransp;
107 EApplyEpetraOpAs epetraFwdOpApplyAs;
108 EAdjointEpetraOp epetraFwdOpAdjointSupport;
109 double epetraFwdOpScalar;
110 epetraFwdOpViewExtractor_->getEpetraOpView(
111 fwdOpSrc.getOp(),
112 outArg(epetraFwdOp), outArg(epetraFwdOpTransp),
113 outArg(epetraFwdOpApplyAs), outArg(epetraFwdOpAdjointSupport),
114 outArg(epetraFwdOpScalar)
115 );
116 if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
117 return false;
118 return true;
119}
120
122{
123 return true;
124}
125
127{
128 return false; // See comment below
129}
130
131Teuchos::RCP<PreconditionerBase<double> >
133{
134 return Teuchos::rcp(new DefaultPreconditioner<double>());
135}
136
138 const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc
139 ,PreconditionerBase<double> *prec
140 ,const ESupportSolveUse /* supportSolveUse */
141 ) const
142{
143 using Teuchos::outArg;
144 using Teuchos::OSTab;
145 using Teuchos::dyn_cast;
146 using Teuchos::RCP;
147 using Teuchos::null;
148 using Teuchos::rcp;
149 using Teuchos::rcp_dynamic_cast;
150 using Teuchos::rcp_const_cast;
151 using Teuchos::set_extra_data;
152 using Teuchos::get_optional_extra_data;
153 using Teuchos::implicit_cast;
154 Teuchos::Time totalTimer(""), timer("");
155 totalTimer.start(true);
156#ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
157 Teuchos::TimeMonitor overallTimeMonitor(*overallTimer);
158#endif
159 const Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
160 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
161 Teuchos::OSTab tab(out);
162 if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
163 *out << "\nEntering Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
164#ifdef TEUCHOS_DEBUG
165 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
166 TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
167#endif
168 Teuchos::RCP<const LinearOpBase<double> >
169 fwdOp = fwdOpSrc->getOp();
170#ifdef TEUCHOS_DEBUG
171 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
172#endif
173 //
174 // Unwrap and get the forward Epetra_Operator object
175 //
176 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
177 EOpTransp epetraFwdOpTransp;
178 EApplyEpetraOpAs epetraFwdOpApplyAs;
179 EAdjointEpetraOp epetraFwdOpAdjointSupport;
180 double epetraFwdOpScalar;
181 epetraFwdOpViewExtractor_->getEpetraOpView(
182 fwdOp,
183 outArg(epetraFwdOp), outArg(epetraFwdOpTransp),
184 outArg(epetraFwdOpApplyAs), outArg(epetraFwdOpAdjointSupport),
185 outArg(epetraFwdOpScalar)
186 );
187 // Validate what we get is what we need
188 RCP<const Epetra_RowMatrix>
189 epetraFwdRowMat = rcp_dynamic_cast<const Epetra_RowMatrix>(epetraFwdOp,true);
190 TEUCHOS_TEST_FOR_EXCEPTION(
191 epetraFwdOpApplyAs != EPETRA_OP_APPLY_APPLY, std::logic_error
192 ,"Error, incorrect apply mode for an Epetra_RowMatrix"
193 );
194 //
195 // Get the concrete preconditioner object
196 //
197 DefaultPreconditioner<double>
198 *defaultPrec = &Teuchos::dyn_cast<DefaultPreconditioner<double> >(*prec);
199 //
200 // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
201 //
202 RCP<EpetraLinearOp>
203 epetra_precOp = rcp_dynamic_cast<EpetraLinearOp>(defaultPrec->getNonconstUnspecifiedPrecOp(),true);
204 //
205 // Get the embedded Ifpack_Preconditioner object if it exists
206 //
207 Teuchos::RCP<Ifpack_Preconditioner>
208 ifpack_precOp;
209 if(epetra_precOp.get())
210 ifpack_precOp = rcp_dynamic_cast<Ifpack_Preconditioner>(epetra_precOp->epetra_op(),true);
211 //
212 // Get the attached forward operator if it exists and make sure that it matches
213 //
214 if(ifpack_precOp.get()) {
215 // ToDo: Get the forward operator and make sure that it matches what is
216 // already being used!
217 }
218 //
219 // Permform initialization if needed
220 //
221 //const bool startingOver = (ifpack_precOp.get() == NULL);
222 const bool startingOver = true;
223 // ToDo: Comment back in the above original version of startingOver to allow
224 // for resuse. Rob H. just pointed out to me that a memory leak is being
225 // created when you just call Ifpack_ILU::Compute() over and over again.
226 // Rob H. said that he will check in a fix the the development branch when
227 // he can.
228 if(startingOver) {
229 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
230 *out << "\nCreating the initial Ifpack_Preconditioner object of type \'"<<Ifpack::toString(precType_)<<"\' ...\n";
231 timer.start(true);
232#ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
233 Teuchos::TimeMonitor creationTimeMonitor(*creationTimer);
234#endif
235 // Create the initial preconditioner
236 ifpack_precOp = rcp(
238 precType_
239 ,const_cast<Epetra_RowMatrix*>(&*epetraFwdRowMat)
240 ,overlap_
241 )
242 );
243 timer.stop();
244 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
245 OSTab(out).o() <<"=> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
246 // Set parameters if the list exists
247 if(paramList_.get()) {
248 Teuchos::ParameterList
249 &ifpackSettingsPL = paramList_->sublist(IfpackSettings_name);
250 // Above will create new sublist if it does not exist!
251 TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->SetParameters(ifpackSettingsPL));
252 // Above, I have not idea how any error messages for a mistake will be
253 // reported back to the user!
254 }
255 // Initialize the structure for the preconditioner
256 TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->Initialize());
257 }
258 //
259 // Attach the epetraFwdOp to the ifpack_precOp to guarantee that it will not go away
260 //
261 set_extra_data(epetraFwdOp, "IFPF::epetraFwdOp", Teuchos::inOutArg(ifpack_precOp),
262 Teuchos::POST_DESTROY, false);
263 //
264 // Update the factorization
265 //
266 {
267 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
268 *out << "\nComputing the preconditioner ...\n";
269#ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
270 Teuchos::TimeMonitor factorizationTimeMonitor(*factorizationTimer);
271#endif
272 timer.start(true);
273 TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->Compute());
274 timer.stop();
275 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
276 OSTab(out).o() <<"=> Setup time = "<<timer.totalElapsedTime()<<" sec\n";
277 }
278 //
279 // Compute the conditioner number estimate if asked
280 //
281
282 // ToDo: Implement
283
284 //
285 // Attach fwdOp to the ifpack_precOp
286 //
287 set_extra_data(fwdOpSrc, "IFPF::fwdOpSrc", Teuchos::inOutArg(ifpack_precOp),
288 Teuchos::POST_DESTROY, false);
289 //
290 // Initialize the output EpetraLinearOp
291 //
292 if(startingOver) {
293 epetra_precOp = rcp(new EpetraLinearOp);
294 }
295 epetra_precOp->initialize(
296 ifpack_precOp
297 ,epetraFwdOpTransp
298 ,EPETRA_OP_APPLY_APPLY_INVERSE
299 ,EPETRA_OP_ADJOINT_SUPPORTED
300 );
301 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_MEDIUM)) {
302 *out << "\nDescription of created preconditioner:\n";
303 OSTab tab2(out);
304 ifpack_precOp->Print(*out);
305 }
306
307 //
308 // Initialize the preconditioner
309 //
310 defaultPrec->initializeUnspecified(
311 Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp)
312 );
313 totalTimer.stop();
314 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
315 *out << "\nTotal time in IfpackPreconditionerFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
316 if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
317 *out << "\nLeaving Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
318}
319
321 PreconditionerBase<double> * /* prec */
322 ,Teuchos::RCP<const LinearOpSourceBase<double> > * /* fwdOpSrc */
323 ,ESupportSolveUse * /* supportSolveUse */
324 ) const
325{
326 TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
327}
328
329// Overridden from ParameterListAcceptor
330
331void IfpackPreconditionerFactory::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
332{
333 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
334 paramList->validateParameters(*this->getValidParameters(),2);
335 // Note: The above validation will validate right down into the the sublist
336 // named IfpackSettings_name!
337 paramList_ = paramList;
338 overlap_ = paramList_->get(Overlap_name,Overlap_default);
339 std::ostringstream oss;
340 oss << "(sub)list \""<<paramList->name()<<"\"parameter \"Prec Type\"";
341 precType_ =
342 ( paramList_.get()
343 ? precTypeValidator->getIntegralValue(*paramList_,PrecType_name,PrecTypeName_default)
344 : PrecType_default
345 );
346 Teuchos::readVerboseObjectSublist(&*paramList_,this);
347#ifdef TEUCHOS_DEBUG
348 // Validate my use of the parameters!
349 paramList->validateParameters(*this->getValidParameters(),1);
350#endif
351}
352
353Teuchos::RCP<Teuchos::ParameterList>
358
359Teuchos::RCP<Teuchos::ParameterList>
361{
362 Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
363 paramList_ = Teuchos::null;
364 return _paramList;
365}
366
367Teuchos::RCP<const Teuchos::ParameterList>
369{
370 return paramList_;
371}
372
373Teuchos::RCP<const Teuchos::ParameterList>
375{
376 using Teuchos::rcp;
377 using Teuchos::rcp_implicit_cast;
378 typedef Teuchos::ParameterEntryValidator PEV;
379 static Teuchos::RCP<Teuchos::ParameterList> validParamList;
380 if(validParamList.get()==NULL) {
381 validParamList = Teuchos::rcp(new Teuchos::ParameterList(Ifpack_name));
382 {
383 // Create the validator for the preconditioner type!
384 Teuchos::Array<std::string>
385 precTypeNames;
386 precTypeNames.insert(
387 precTypeNames.begin(),
390 );
391 Teuchos::Array<Ifpack::EPrecType>
392 precTypeValues;
393 precTypeValues.insert(
394 precTypeValues.begin(),
397 );
398 precTypeValidator = rcp(
399 new Teuchos::StringToIntegralParameterEntryValidator<Ifpack::EPrecType>(
400 precTypeNames,precTypeValues,PrecType_name
401 )
402 );
403 }
404 validParamList->set(
405 PrecType_name, PrecTypeName_default,
406 "Type of Ifpack preconditioner to use.",
407 rcp_implicit_cast<const PEV>(precTypeValidator)
408 );
409 validParamList->set(
410 Overlap_name, Overlap_default,
411 "Number of rows/columns overlapped between subdomains in different"
412 "\nprocesses in the additive Schwarz-type domain-decomposition preconditioners."
413 );
414 validParamList->set(
415 IfpackSettings_name, Ifpack_GetValidParameters(),
416 "Preconditioner settings that are passed onto the Ifpack preconditioners themselves."
417 );
418 // Note that in the above setParameterList(...) function that we actually
419 // validate down into the first level of this sublist. Really the
420 // Ifpack_Preconditioner objects themselves should do validation but we do
421 // it ourselves taking the return from the Ifpack_GetValidParameters()
422 // function as gospel!
423 Teuchos::setupVerboseObjectSublist(&*validParamList);
424 }
425 return validParamList;
426}
427
428// Public functions overridden from Teuchos::Describable
429
431{
432 std::ostringstream oss;
433 oss << "Thyra::IfpackPreconditionerFactory{";
434 oss << "precType=\"" << ::Ifpack::toString(precType_) << "\"";
435 oss << ",overlap=" << overlap_;
436 oss << "}";
437 return oss.str();
438}
439
440// private
441
442void IfpackPreconditionerFactory::initializeTimers()
443{
444 if(!overallTimer.get()) {
445#ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
446 overallTimer = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF");
447 creationTimer = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF:Creation");
448 factorizationTimer = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF:Factorization");
449#endif
450 }
451}
452
453} // namespace Thyra
static const int numPrecTypes
static const char * precTypeNames[numPrecTypes]
static const char * toString(const EPrecType precType)
Ifpack_Preconditioner * Create(const std::string PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
static const EPrecType precTypeValues[numPrecTypes]
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
void uninitializePrec(PreconditionerBase< double > *prec, Teuchos::RCP< const LinearOpSourceBase< double > > *fwdOpSrc, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< PreconditionerBase< double > > createPrec() const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
bool isCompatible(const LinearOpSourceBase< double > &fwdOpSrc) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< double > > &fwdOpSrc, PreconditionerBase< double > *prec, const ESupportSolveUse supportSolveUse) const

Generated for Stratimikos by doxygen 1.10.0