Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Thyra_MLPreconditionerFactory.cpp
Go to the documentation of this file.
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
43
44#include "Thyra_EpetraOperatorViewExtractorStd.hpp"
45#include "Thyra_EpetraLinearOp.hpp"
46#include "Thyra_DefaultPreconditioner.hpp"
47#include "ml_MultiLevelPreconditioner.h"
48#include "ml_MultiLevelOperator.h"
49#include "ml_ValidateParameters.h"
50#include "ml_RefMaxwell.h"
51#include "Epetra_RowMatrix.h"
52#include "Teuchos_StandardParameterEntryValidators.hpp"
53#include "Teuchos_dyn_cast.hpp"
54#include "Teuchos_TimeMonitor.hpp"
55#include "Teuchos_implicit_cast.hpp"
56#include "Teuchos_ValidatorXMLConverterDB.hpp"
57#include "Teuchos_StaticSetupMacro.hpp"
58#include "Teuchos_iostream_helpers.hpp"
59
60
61namespace {
62
63
64enum EMLProblemType {
65 ML_PROBTYPE_NONE,
66 ML_PROBTYPE_SMOOTHED_AGGREGATION,
67 ML_PROBTYPE_NONSYMMETRIC_SMOOTHED_AGGREGATION,
68 ML_PROBTYPE_DOMAIN_DECOMPOSITION,
69 ML_PROBTYPE_DOMAIN_DECOMPOSITION_ML,
70 ML_PROBTYPE_MAXWELL,
71 ML_PROBTYPE_REFMAXWELL
72};
73const std::string BaseMethodDefaults_valueNames_none = "none";
74const Teuchos::Array<std::string> BaseMethodDefaults_valueNames
75= Teuchos::tuple<std::string>(
76 BaseMethodDefaults_valueNames_none,
77 "SA",
78 "NSSA",
79 "DD",
80 "DD-ML",
81 "maxwell",
82 "refmaxwell"
83 );
84
85
86TEUCHOS_ENUM_INPUT_STREAM_OPERATOR(EMLProblemType)
87
88
89TEUCHOS_STATIC_SETUP()
90{
91 TEUCHOS_ADD_STRINGTOINTEGRALVALIDATOR_CONVERTER(EMLProblemType);
92}
93
94const std::string BaseMethodDefaults_name = "Base Method Defaults";
95const std::string BaseMethodDefaults_default = "SA";
96Teuchos::RCP<
97 Teuchos::StringToIntegralParameterEntryValidator<EMLProblemType>
98 >
99BaseMethodDefaults_validator;
100
101const std::string ReuseFineLevelSmoother_name = "Reuse Fine Level Smoother";
102const bool ReuseFineLevelSmoother_default = false;
103
104const std::string MLSettings_name = "ML Settings";
105
106
107} // namespace
108
109
110namespace Thyra {
111
112
113using Teuchos::RCP;
114using Teuchos::ParameterList;
115
116
117// Constructors/initializers/accessors
118
119
121 :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
122{}
123
124
125// Overridden from PreconditionerFactoryBase
126
127
129 const LinearOpSourceBase<double> &fwdOpSrc
130 ) const
131{
132 using Teuchos::outArg;
133 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
134 EOpTransp epetraFwdOpTransp;
135 EApplyEpetraOpAs epetraFwdOpApplyAs;
136 EAdjointEpetraOp epetraFwdOpAdjointSupport;
137 double epetraFwdOpScalar;
138 Teuchos::RCP<const LinearOpBase<double> >
139 fwdOp = fwdOpSrc.getOp();
140 epetraFwdOpViewExtractor_->getEpetraOpView(
141 fwdOp,
142 outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
143 outArg(epetraFwdOpApplyAs),
144 outArg(epetraFwdOpAdjointSupport),
145 outArg(epetraFwdOpScalar)
146 );
147 if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
148 return false;
149 return true;
150}
151
152
154{
155 return true;
156}
157
158
160{
161 return false; // See comment below
162}
163
164
165RCP<PreconditionerBase<double> >
167{
168 return Teuchos::rcp(new DefaultPreconditioner<double>());
169}
170
171
173 const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
174 PreconditionerBase<double> *prec,
175 const ESupportSolveUse /* supportSolveUse */
176 ) const
177{
178 using Teuchos::outArg;
179 using Teuchos::OSTab;
180 using Teuchos::dyn_cast;
181 using Teuchos::RCP;
182 using Teuchos::null;
183 using Teuchos::rcp;
184 using Teuchos::rcp_dynamic_cast;
185 using Teuchos::rcp_const_cast;
186 using Teuchos::set_extra_data;
187 using Teuchos::get_optional_extra_data;
188 using Teuchos::implicit_cast;
189 Teuchos::Time totalTimer(""), timer("");
190 totalTimer.start(true);
191 const RCP<Teuchos::FancyOStream> out = this->getOStream();
192 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
193 Teuchos::OSTab tab(out);
194 if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
195 *out << "\nEntering Thyra::MLPreconditionerFactory::initializePrec(...) ...\n";
196
197 // Get the problem type
198 const EMLProblemType problemType = BaseMethodDefaults_validator->getIntegralValue(*paramList_,BaseMethodDefaults_name,BaseMethodDefaults_default);
199
200 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
201#ifdef _DEBUG
202 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
203 TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
204#endif
205 //
206 // Unwrap and get the forward Epetra_Operator object
207 //
208 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
209 EOpTransp epetraFwdOpTransp;
210 EApplyEpetraOpAs epetraFwdOpApplyAs;
211 EAdjointEpetraOp epetraFwdOpAdjointSupport;
212 double epetraFwdOpScalar;
213 epetraFwdOpViewExtractor_->getEpetraOpView(
214 fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
215 outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
216 );
217 // Validate what we get is what we need
218 RCP<const Epetra_RowMatrix>
219 epetraFwdRowMat = rcp_dynamic_cast<const Epetra_RowMatrix>(epetraFwdOp,true);
220 TEUCHOS_TEST_FOR_EXCEPTION(
221 epetraFwdOpApplyAs != EPETRA_OP_APPLY_APPLY, std::logic_error
222 ,"Error, incorrect apply mode for an Epetra_RowMatrix"
223 );
224 RCP<const Epetra_CrsMatrix> epetraFwdCrsMat = rcp_dynamic_cast<const Epetra_CrsMatrix>(epetraFwdRowMat);
225
226 //
227 // Get the concrete preconditioner object
228 //
229 DefaultPreconditioner<double>
230 *defaultPrec = &Teuchos::dyn_cast<DefaultPreconditioner<double> >(*prec);
231 //
232 // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
233 //
234 RCP<EpetraLinearOp>
235 epetra_precOp = rcp_dynamic_cast<EpetraLinearOp>(defaultPrec->getNonconstUnspecifiedPrecOp(),true);
236 //
237 // Get the embedded ML_Epetra Preconditioner object if it exists
238 //
239 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> ml_precOp;
240 Teuchos::RCP<ML_Epetra::RefMaxwellPreconditioner> rm_precOp;
241 if(epetra_precOp.get()) {
242 if(problemType == ML_PROBTYPE_REFMAXWELL)
243 rm_precOp = rcp_dynamic_cast<ML_Epetra::RefMaxwellPreconditioner>(epetra_precOp->epetra_op(),true);
244 else
245 ml_precOp = rcp_dynamic_cast<ML_Epetra::MultiLevelPreconditioner>(epetra_precOp->epetra_op(),true);
246 }
247 //
248 // Get the attached forward operator if it exists and make sure that it matches
249 //
250 if(ml_precOp!=Teuchos::null) {
251 // Get the forward operator and make sure that it matches what is
252 // already being used!
253 const Epetra_RowMatrix & rm = ml_precOp->RowMatrix();
254
255 TEUCHOS_TEST_FOR_EXCEPTION(
256 &rm!=&*epetraFwdRowMat, std::logic_error
257 ,"ML requires Epetra_RowMatrix to be the same for each initialization of the preconditioner"
258 );
259 }
260 // NOTE: No such check exists for RefMaxwell
261
262 //
263 // Perform initialization if needed
264 //
265 const bool startingOver = (ml_precOp.get() == NULL && rm_precOp.get() == NULL);
266 if(startingOver)
267 {
268 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
269 *out << "\nCreating the initial ML_Epetra::MultiLevelPreconditioner object...\n";
270 timer.start(true);
271 // Create the initial preconditioner: DO NOT compute it yet
272
273 if(problemType==ML_PROBTYPE_REFMAXWELL)
274 rm_precOp = rcp(new ML_Epetra::RefMaxwellPreconditioner(*epetraFwdCrsMat, paramList_->sublist(MLSettings_name),false));
275 else
276 ml_precOp = rcp(new ML_Epetra::MultiLevelPreconditioner(*epetraFwdRowMat, paramList_->sublist(MLSettings_name),false));
277
278
279 timer.stop();
280 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
281 OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
282 // RAB: Above, I am just passing a string to ML::Create(...) in order
283 // get this code written. However, in the future, it would be good to
284 // copy the contents of what is in ML::Create(...) into a local
285 // function and then use switch(...) to create the initial
286 // ML_Epetra::MultiLevelPreconditioner object. This would result in better validation
287 // and faster code.
288 // Set parameters if the list exists
289 if(paramList_.get()) {
290 if (problemType==ML_PROBTYPE_REFMAXWELL) {
291 TEUCHOS_TEST_FOR_EXCEPT(0!=rm_precOp->SetParameterList(paramList_->sublist(MLSettings_name)));
292 }
293 else {
294 TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->SetParameterList(paramList_->sublist(MLSettings_name)));
295 }
296 }
297 }
298 //
299 // Attach the epetraFwdOp to the ml_precOp to guarantee that it will not go away
300 //
301 if (problemType==ML_PROBTYPE_REFMAXWELL)
302 set_extra_data(epetraFwdOp, "IFPF::epetraFwdOp", Teuchos::inOutArg(rm_precOp),
303 Teuchos::POST_DESTROY, false);
304 else
305 set_extra_data(epetraFwdOp, "IFPF::epetraFwdOp", Teuchos::inOutArg(ml_precOp),
306 Teuchos::POST_DESTROY, false);
307 //
308 // Update the factorization
309 //
310 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
311 *out << "\nComputing the preconditioner ...\n";
312 timer.start(true);
313 if (problemType==ML_PROBTYPE_REFMAXWELL) {
314 if (startingOver) {
315 TEUCHOS_TEST_FOR_EXCEPT(0!=rm_precOp->ComputePreconditioner());
316 }
317 else {
318 TEUCHOS_TEST_FOR_EXCEPT(0!=rm_precOp->ReComputePreconditioner());
319 }
320 }
321 else {
322 if (startingOver) {
323 TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->ComputePreconditioner());
324 }
325 else {
326 TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->ReComputePreconditioner(paramList_->get<bool>(ReuseFineLevelSmoother_name)));
327 }
328 }
329 timer.stop();
330 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
331 OSTab(out).o() <<"=> Setup time = "<<timer.totalElapsedTime()<<" sec\n";
332 //
333 // Compute the conditioner number estimate if asked
334 //
335
336 // ToDo: Implement
337
338 //
339 // Attach fwdOp to the ml_precOp
340 //
341 if (problemType==ML_PROBTYPE_REFMAXWELL)
342 set_extra_data(fwdOp, "IFPF::fwdOp", Teuchos::inOutArg(rm_precOp),
343 Teuchos::POST_DESTROY, false);
344 else
345 set_extra_data(fwdOp, "IFPF::fwdOp", Teuchos::inOutArg(ml_precOp),
346 Teuchos::POST_DESTROY, false);
347 //
348 // Initialize the output EpetraLinearOp
349 //
350 if(startingOver) {
351 epetra_precOp = rcp(new EpetraLinearOp);
352 }
353 // ToDo: Look into adjoints again.
354 if (problemType==ML_PROBTYPE_REFMAXWELL)
355 epetra_precOp->initialize(rm_precOp,epetraFwdOpTransp,EPETRA_OP_APPLY_APPLY_INVERSE,EPETRA_OP_ADJOINT_UNSUPPORTED);
356 else
357 epetra_precOp->initialize(ml_precOp,epetraFwdOpTransp,EPETRA_OP_APPLY_APPLY_INVERSE,EPETRA_OP_ADJOINT_UNSUPPORTED);
358 //
359 // Initialize the preconditioner
360 //
361 defaultPrec->initializeUnspecified(
362 Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp)
363 );
364 totalTimer.stop();
365 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
366 *out << "\nTotal time in MLPreconditionerFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
367 if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
368 *out << "\nLeaving Thyra::MLPreconditionerFactory::initializePrec(...) ...\n";
369}
370
371
373 PreconditionerBase<double> * /* prec */,
374 Teuchos::RCP<const LinearOpSourceBase<double> > * /* fwdOp */,
375 ESupportSolveUse * /* supportSolveUse */
376 ) const
377{
378 TEUCHOS_TEST_FOR_EXCEPT(true);
379}
380
381
382// Overridden from ParameterListAcceptor
383
384
386 Teuchos::RCP<ParameterList> const& paramList
387 )
388{
389 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
390
391 // Do not recurse
392 paramList->validateParameters(*this->getValidParameters(),0);
393 paramList_ = paramList;
394
395 // set default for reuse of fine level smoother
396 if(!paramList_->isType<bool>(ReuseFineLevelSmoother_name))
397 paramList_->set<bool>(ReuseFineLevelSmoother_name,ReuseFineLevelSmoother_default);
398
399 const EMLProblemType
400 defaultType = BaseMethodDefaults_validator->getIntegralValue(
401 *paramList_,BaseMethodDefaults_name,BaseMethodDefaults_default
402 );
403 if( ML_PROBTYPE_NONE != defaultType ) {
404 const std::string defaultTypeStr = BaseMethodDefaults_valueNames[defaultType];
405
406 // ML will do validation on its own. We don't need to duplicate that here.
407 Teuchos::ParameterList defaultParams;
408 if(defaultType == ML_PROBTYPE_REFMAXWELL) {
409 ML_Epetra::SetDefaultsRefMaxwell(defaultParams);
410 }
411 else {
412 TEUCHOS_TEST_FOR_EXCEPTION(0!=ML_Epetra::SetDefaults(defaultTypeStr,defaultParams)
413 ,Teuchos::Exceptions::InvalidParameterValue
414 ,"Error, the ML problem type \"" << defaultTypeStr << "\' is not recognized by ML!"
415 );
416 }
417
418 // Note, the only way the above exception message could be generated is if
419 // a default problem type was removed from ML_Epetra::SetDefaults(...).
420 // When a new problem type is added to this function, it must be added to
421 // our enum EMLProblemType along with associated objects ... In other
422 // words, this adapter must be maintained as ML is maintained. An
423 // alternative design would be to just pass in whatever string to this
424 // function. This would improve maintainability but it would not generate
425 // very good error messages when a bad string was passed in. Currently,
426 // the error message attached to the exception will contain the list of
427 // valid problem types.
428 paramList_->sublist(MLSettings_name).setParametersNotAlreadySet(
429 defaultParams);
430 }
431
432}
433
434
435RCP<ParameterList>
437{
438 return paramList_;
439}
440
441
442RCP<ParameterList>
444{
445 Teuchos::RCP<ParameterList> _paramList = paramList_;
446 paramList_ = Teuchos::null;
447 return _paramList;
448}
449
450
451RCP<const ParameterList>
453{
454 return paramList_;
455}
456
457
458RCP<const ParameterList>
460{
461 // NOTE: We're only going to use this function to genrate valid *Stratimikos* parameters.
462 // Since ML's parameters can be validated internally, we'll handle those separarely.
463
464
465 using Teuchos::rcp;
466 using Teuchos::tuple;
467 using Teuchos::implicit_cast;
468 using Teuchos::rcp_implicit_cast;
469 typedef Teuchos::ParameterEntryValidator PEV;
470
471 static RCP<const ParameterList> validPL;
472
473 if(is_null(validPL)) {
474
475 RCP<ParameterList>
476 pl = rcp(new ParameterList());
477
478 BaseMethodDefaults_validator = rcp(
479 new Teuchos::StringToIntegralParameterEntryValidator<EMLProblemType>(
480 BaseMethodDefaults_valueNames,
481 tuple<std::string>(
482 "Do not set any default parameters",
483 "Set default parameters for a smoothed aggregation method",
484 "Set default parameters for a nonsymmetric smoothed aggregation method",
485 "Set default parameters for a domain decomposition method",
486 "Set default parameters for a domain decomposition method special to ML",
487 "Set default parameters for a Maxwell-type of preconditioner",
488 "Set default parameters for a RefMaxwell-type preconditioner"
489 ),
490 tuple<EMLProblemType>(
491 ML_PROBTYPE_NONE,
492 ML_PROBTYPE_SMOOTHED_AGGREGATION,
493 ML_PROBTYPE_NONSYMMETRIC_SMOOTHED_AGGREGATION,
494 ML_PROBTYPE_DOMAIN_DECOMPOSITION,
495 ML_PROBTYPE_DOMAIN_DECOMPOSITION_ML,
496 ML_PROBTYPE_MAXWELL,
497 ML_PROBTYPE_REFMAXWELL
498 ),
499 BaseMethodDefaults_name
500 )
501 );
502
503 pl->set(BaseMethodDefaults_name,BaseMethodDefaults_default,
504 "Select the default method type which also sets parameter defaults\n"
505 "in the sublist \"" + MLSettings_name + "\"!",
506 rcp_implicit_cast<const PEV>(BaseMethodDefaults_validator)
507 );
508
509 pl->set(ReuseFineLevelSmoother_name,ReuseFineLevelSmoother_default,
510 "Enables/disables the reuse of the fine level smoother.");
511
512 ParameterList mlpl;
513 pl->set(MLSettings_name,mlpl,
514 "Sampling of the parameters directly accepted by ML\n"
515 "This list of parameters is generated by combining all of\n"
516 "the parameters set for all of the default problem types supported\n"
517 "by ML. Therefore, do not assume these parameters are at values that\n"
518 "are reasonable to ML. This list is just to give a sense of some of\n"
519 "the parameters that ML accepts. Consult ML documentation on how to\n"
520 "set these parameters. Also, you can print the parameter list after\n"
521 "it is used and see what defaults where set for each default problem\n"
522 "type. Warning! the parameters in this sublist are currently *not*\n"
523 "being validated by ML!"
524 );
525
526 validPL = pl;
527
528 }
529
530 return validPL;
531
532}
533
534
535// Public functions overridden from Teuchos::Describable
536
537
539{
540 std::ostringstream oss;
541 oss << "Thyra::MLPreconditionerFactory";
542 return oss.str();
543}
544
545
546} // namespace Thyra
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
bool isCompatible(const LinearOpSourceBase< double > &fwdOp) const
void uninitializePrec(PreconditionerBase< double > *prec, Teuchos::RCP< const LinearOpSourceBase< double > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< double > > &fwdOp, PreconditionerBase< double > *prec, const ESupportSolveUse supportSolveUse) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
Teuchos::RCP< PreconditionerBase< double > > createPrec() const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Teuchos::RCP< Teuchos::ParameterList > paramList_