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