Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Thyra_MLPreconditionerFactory.cpp
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Thyra_MLPreconditionerFactory.hpp"
00031 
00032 #include "Thyra_EpetraOperatorViewExtractorStd.hpp"
00033 #include "Thyra_EpetraLinearOp.hpp"
00034 #include "Thyra_DefaultPreconditioner.hpp"
00035 #include "ml_MultiLevelPreconditioner.h"
00036 #include "ml_MultiLevelOperator.h"
00037 #include "ml_ValidateParameters.h"
00038 #include "Epetra_RowMatrix.h"
00039 #include "Teuchos_StandardParameterEntryValidators.hpp"
00040 #include "Teuchos_dyn_cast.hpp"
00041 #include "Teuchos_TimeMonitor.hpp"
00042 #include "Teuchos_implicit_cast.hpp"
00043 #include "Teuchos_ValidatorXMLConverterDB.hpp"
00044 #include "Teuchos_StaticSetupMacro.hpp"
00045 #include "Teuchos_iostream_helpers.hpp"
00046 
00047 
00048 namespace {
00049 
00050 
00051 enum EMLProblemType {
00052   ML_PROBTYPE_NONE,
00053   ML_PROBTYPE_SMOOTHED_AGGREGATION, 
00054   ML_PROBTYPE_DOMAIN_DECOMPOSITION,
00055   ML_PROBTYPE_DOMAIN_DECOMPOSITION_ML,
00056   ML_PROBTYPE_MAXWELL
00057 };
00058 const std::string BaseMethodDefaults_valueNames_none = "none";
00059 const Teuchos::Array<std::string> BaseMethodDefaults_valueNames
00060 = Teuchos::tuple<std::string>(
00061   BaseMethodDefaults_valueNames_none,
00062   "SA", 
00063   "DD",
00064   "DD-ML",
00065   "maxwell"
00066   );
00067 
00068 
00069 TEUCHOS_ENUM_INPUT_STREAM_OPERATOR(EMLProblemType)
00070 
00071 
00072 TEUCHOS_STATIC_SETUP()
00073 {
00074   TEUCHOS_ADD_STRINGTOINTEGRALVALIDATOR_CONVERTER(EMLProblemType);
00075 }
00076 
00077 const std::string BaseMethodDefaults_name = "Base Method Defaults";
00078 const std::string BaseMethodDefaults_default = "SA";
00079 Teuchos::RCP<
00080   Teuchos::StringToIntegralParameterEntryValidator<EMLProblemType>
00081   >
00082 BaseMethodDefaults_validator;
00083 
00084 const std::string ReuseFineLevelSmoother_name = "Reuse Fine Level Smoother";
00085 const bool ReuseFineLevelSmoother_default = false;
00086   
00087 const std::string MLSettings_name = "ML Settings";
00088 
00089 
00090 } // namespace
00091 
00092 
00093 namespace Thyra {
00094 
00095 
00096 using Teuchos::RCP;
00097 using Teuchos::ParameterList;
00098 
00099 
00100 // Constructors/initializers/accessors
00101 
00102   
00103 MLPreconditionerFactory::MLPreconditionerFactory()
00104   :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
00105 {}
00106 
00107 
00108 // Overridden from PreconditionerFactoryBase
00109 
00110 
00111 bool MLPreconditionerFactory::isCompatible(
00112   const LinearOpSourceBase<double> &fwdOpSrc
00113   ) const
00114 {
00115   using Teuchos::outArg;
00116   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00117   EOpTransp epetraFwdOpTransp;
00118   EApplyEpetraOpAs epetraFwdOpApplyAs;
00119   EAdjointEpetraOp epetraFwdOpAdjointSupport;
00120   double epetraFwdOpScalar;
00121   Teuchos::RCP<const LinearOpBase<double> >
00122     fwdOp = fwdOpSrc.getOp();
00123   epetraFwdOpViewExtractor_->getEpetraOpView(
00124     fwdOp,
00125     outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
00126     outArg(epetraFwdOpApplyAs),
00127     outArg(epetraFwdOpAdjointSupport),
00128     outArg(epetraFwdOpScalar)
00129     );
00130   if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
00131     return false;
00132   return true;
00133 }
00134 
00135 
00136 bool MLPreconditionerFactory::applySupportsConj(EConj conj) const
00137 {
00138   return true;
00139 }
00140 
00141 
00142 bool MLPreconditionerFactory::applyTransposeSupportsConj(EConj conj) const
00143 {
00144   return false; // See comment below
00145 }
00146 
00147 
00148 RCP<PreconditionerBase<double> >
00149 MLPreconditionerFactory::createPrec() const
00150 {
00151   return Teuchos::rcp(new DefaultPreconditioner<double>());
00152 }
00153 
00154 
00155 void MLPreconditionerFactory::initializePrec(
00156   const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
00157   PreconditionerBase<double> *prec,
00158   const ESupportSolveUse supportSolveUse
00159   ) const
00160 {
00161   using Teuchos::outArg;
00162   using Teuchos::OSTab;
00163   using Teuchos::dyn_cast;
00164   using Teuchos::RCP;
00165   using Teuchos::null;
00166   using Teuchos::rcp;
00167   using Teuchos::rcp_dynamic_cast;
00168   using Teuchos::rcp_const_cast;
00169   using Teuchos::set_extra_data;
00170   using Teuchos::get_optional_extra_data;
00171   using Teuchos::implicit_cast;
00172   Teuchos::Time totalTimer(""), timer("");
00173   totalTimer.start(true);
00174   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00175   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00176   Teuchos::OSTab tab(out);
00177   if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
00178     *out << "\nEntering Thyra::MLPreconditionerFactory::initializePrec(...) ...\n";
00179 
00180   Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
00181 #ifdef _DEBUG
00182   TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00183   TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
00184 #endif
00185   //
00186   // Unwrap and get the forward Epetra_Operator object
00187   //
00188   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00189   EOpTransp epetraFwdOpTransp;
00190   EApplyEpetraOpAs epetraFwdOpApplyAs;
00191   EAdjointEpetraOp epetraFwdOpAdjointSupport;
00192   double epetraFwdOpScalar;
00193   epetraFwdOpViewExtractor_->getEpetraOpView(
00194     fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
00195     outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
00196                                              );
00197   // Validate what we get is what we need
00198   RCP<const Epetra_RowMatrix>
00199     epetraFwdRowMat = rcp_dynamic_cast<const Epetra_RowMatrix>(epetraFwdOp,true);
00200   TEUCHOS_TEST_FOR_EXCEPTION(
00201     epetraFwdOpApplyAs != EPETRA_OP_APPLY_APPLY, std::logic_error
00202     ,"Error, incorrect apply mode for an Epetra_RowMatrix"
00203     );
00204 
00205   //
00206   // Get the concrete preconditioner object
00207   //
00208   DefaultPreconditioner<double>
00209     *defaultPrec = &Teuchos::dyn_cast<DefaultPreconditioner<double> >(*prec);
00210   //
00211   // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
00212   //
00213   RCP<EpetraLinearOp>
00214     epetra_precOp = rcp_dynamic_cast<EpetraLinearOp>(defaultPrec->getNonconstUnspecifiedPrecOp(),true);
00215   //
00216   // Get the embedded ML_Epetra::MultiLevelPreconditioner object if it exists
00217   //
00218   Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> ml_precOp;
00219   if(epetra_precOp.get())
00220     ml_precOp = rcp_dynamic_cast<ML_Epetra::MultiLevelPreconditioner>(epetra_precOp->epetra_op(),true);
00221   //
00222   // Get the attached forward operator if it exists and make sure that it matches
00223   //
00224   if(ml_precOp!=Teuchos::null) {
00225     // Get the forward operator and make sure that it matches what is
00226     // already being used!
00227     const Epetra_RowMatrix & rm = ml_precOp->RowMatrix();
00228    
00229     TEUCHOS_TEST_FOR_EXCEPTION(
00230        &rm!=&*epetraFwdRowMat, std::logic_error
00231        ,"ML requires Epetra_RowMatrix to be the same for each initialization of the preconditioner"
00232        );
00233   }
00234   //
00235   // Perform initialization if needed
00236   //
00237   const bool startingOver = (ml_precOp.get() == NULL);
00238   if(startingOver) 
00239   {
00240     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00241       *out << "\nCreating the initial ML_Epetra::MultiLevelPreconditioner object...\n";
00242     timer.start(true);
00243     // Create the initial preconditioner: DO NOT compute it yet
00244     ml_precOp = rcp(
00245       new ML_Epetra::MultiLevelPreconditioner(
00246         *epetraFwdRowMat, paramList_->sublist(MLSettings_name),false
00247         )
00248       );
00249     
00250     timer.stop();
00251     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00252       OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
00253     // RAB: Above, I am just passing a string to ML::Create(...) in order
00254     // get this code written.  However, in the future, it would be good to
00255     // copy the contents of what is in ML::Create(...) into a local
00256     // function and then use switch(...) to create the initial
00257     // ML_Epetra::MultiLevelPreconditioner object.  This would result in better validation
00258     // and faster code.
00259     // Set parameters if the list exists
00260     if(paramList_.get())
00261       TEUCHOS_TEST_FOR_EXCEPT(
00262         0!=ml_precOp->SetParameterList(paramList_->sublist(MLSettings_name))
00263         );
00264     // Initialize the structure for the preconditioner
00265     //        TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->Initialize());
00266   }
00267   //
00268   // Attach the epetraFwdOp to the ml_precOp to guarantee that it will not go away
00269   //
00270   set_extra_data(epetraFwdOp, "IFPF::epetraFwdOp", Teuchos::inOutArg(ml_precOp),
00271     Teuchos::POST_DESTROY, false);
00272   //
00273   // Update the factorization
00274   //
00275   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00276     *out << "\nComputing the preconditioner ...\n";
00277   timer.start(true);
00278   if (startingOver) {
00279     TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->ComputePreconditioner());
00280   }
00281   else {
00282     TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->ReComputePreconditioner(paramList_->get<bool>(ReuseFineLevelSmoother_name)));
00283   }
00284   timer.stop();
00285   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00286     OSTab(out).o() <<"=> Setup time = "<<timer.totalElapsedTime()<<" sec\n";
00287   //
00288   // Compute the conditioner number estimate if asked
00289   //
00290 
00291   // ToDo: Implement
00292 
00293   //
00294   // Attach fwdOp to the ml_precOp
00295   //
00296   set_extra_data(fwdOp, "IFPF::fwdOp", Teuchos::inOutArg(ml_precOp),
00297     Teuchos::POST_DESTROY, false);
00298   //
00299   // Initialize the output EpetraLinearOp
00300   //
00301   if(startingOver) {
00302     epetra_precOp = rcp(new EpetraLinearOp);
00303   }
00304   epetra_precOp->initialize(
00305     ml_precOp
00306     ,epetraFwdOpTransp
00307     ,EPETRA_OP_APPLY_APPLY_INVERSE
00308     ,EPETRA_OP_ADJOINT_UNSUPPORTED  // ToDo: Look into adjoints again.
00309     );
00310   //
00311   // Initialize the preconditioner
00312   //
00313   defaultPrec->initializeUnspecified(
00314     Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp)
00315     );
00316   totalTimer.stop();
00317   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00318     *out << "\nTotal time in MLPreconditionerFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
00319   if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
00320     *out << "\nLeaving Thyra::MLPreconditionerFactory::initializePrec(...) ...\n";
00321 }
00322 
00323 
00324 void MLPreconditionerFactory::uninitializePrec(
00325   PreconditionerBase<double> *prec,
00326   Teuchos::RCP<const LinearOpSourceBase<double> > *fwdOp,
00327   ESupportSolveUse *supportSolveUse
00328   ) const
00329 {
00330   TEUCHOS_TEST_FOR_EXCEPT(true);
00331 }
00332 
00333 
00334 // Overridden from ParameterListAcceptor
00335 
00336 
00337 void MLPreconditionerFactory::setParameterList(
00338   Teuchos::RCP<ParameterList> const& paramList
00339   )
00340 {
00341   TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
00342   paramList->validateParameters(*this->getValidParameters(),0);
00343   paramList_ = paramList;
00344 
00345   // set default for reuse of fine level smoother
00346   if(!paramList_->isType<bool>(ReuseFineLevelSmoother_name))
00347     paramList_->set<bool>(ReuseFineLevelSmoother_name,ReuseFineLevelSmoother_default);
00348 
00349   const EMLProblemType
00350     defaultType = BaseMethodDefaults_validator->getIntegralValue(
00351       *paramList_,BaseMethodDefaults_name,BaseMethodDefaults_default
00352       );
00353   if( ML_PROBTYPE_NONE != defaultType ) {
00354     const std::string
00355       defaultTypeStr = BaseMethodDefaults_valueNames[defaultType];
00356     Teuchos::ParameterList defaultParams;
00357     TEUCHOS_TEST_FOR_EXCEPTION(
00358       0!=ML_Epetra::SetDefaults(defaultTypeStr,defaultParams)
00359       ,Teuchos::Exceptions::InvalidParameterValue
00360       ,"Error, the ML problem type \"" << defaultTypeStr << "\' is not recognized by ML!"
00361       );
00362     // Note, the only way the above exception message could be generated is if
00363     // a default problem type was removed from ML_Epetra::SetDefaults(...).
00364     // When a new problem type is added to this function, it must be added to
00365     // our enum EMLProblemType along with associated objects ...  In other
00366     // words, this adapter must be maintained as ML is maintained.  An
00367     // alternative design would be to just pass in whatever string to this
00368     // function.  This would improve maintainability but it would not generate
00369     // very good error messages when a bad string was passed in.  Currently,
00370     // the error message attached to the exception will contain the list of
00371     // valid problem types.
00372     paramList_->sublist(MLSettings_name).setParametersNotAlreadySet(
00373       defaultParams);
00374   }
00375 #ifdef TEUCHOS_DEBUG
00376   paramList->validateParameters(*this->getValidParameters(),0);
00377 #endif
00378 }
00379 
00380 
00381 RCP<ParameterList>
00382 MLPreconditionerFactory::getNonconstParameterList()
00383 {
00384   return paramList_;
00385 }
00386 
00387 
00388 RCP<ParameterList>
00389 MLPreconditionerFactory::unsetParameterList()
00390 {
00391   Teuchos::RCP<ParameterList> _paramList = paramList_;
00392   paramList_ = Teuchos::null;
00393   return _paramList;
00394 }
00395 
00396 
00397 RCP<const ParameterList>
00398 MLPreconditionerFactory::getParameterList() const
00399 {
00400   return paramList_;
00401 }
00402 
00403 
00404 RCP<const ParameterList>
00405 MLPreconditionerFactory::getValidParameters() const
00406 {
00407 
00408   using Teuchos::rcp;
00409   using Teuchos::tuple;
00410   using Teuchos::implicit_cast;
00411   using Teuchos::rcp_implicit_cast;
00412   typedef Teuchos::ParameterEntryValidator PEV;
00413 
00414   static RCP<const ParameterList> validPL;
00415 
00416   if(is_null(validPL)) {
00417 
00418     RCP<ParameterList>
00419       pl = rcp(new ParameterList());
00420 
00421     BaseMethodDefaults_validator = rcp(
00422       new Teuchos::StringToIntegralParameterEntryValidator<EMLProblemType>(
00423         BaseMethodDefaults_valueNames,
00424         tuple<std::string>(
00425           "Do not set any default parameters",
00426           "Set default parameters for a smoothed aggregation method",
00427           "Set default parameters for a domain decomposition method",
00428           "Set default parameters for a domain decomposition method special to ML",
00429           "Set default parameters for a Maxwell-type of linear operator"
00430           ),
00431         tuple<EMLProblemType>(
00432           ML_PROBTYPE_NONE,
00433           ML_PROBTYPE_SMOOTHED_AGGREGATION,
00434           ML_PROBTYPE_DOMAIN_DECOMPOSITION,
00435           ML_PROBTYPE_DOMAIN_DECOMPOSITION_ML,
00436           ML_PROBTYPE_MAXWELL
00437           ),
00438         BaseMethodDefaults_name
00439         )
00440       );
00441 
00442     pl->set(BaseMethodDefaults_name,BaseMethodDefaults_default,
00443       "Select the default method type which also sets parameter defaults\n"
00444       "in the sublist \"" + MLSettings_name + "\"!",
00445       rcp_implicit_cast<const PEV>(BaseMethodDefaults_validator)
00446       );
00447 
00448     pl->set(ReuseFineLevelSmoother_name,ReuseFineLevelSmoother_default,
00449       "Enables/disables the reuse of the fine level smoother.");
00450 
00451 /* 2007/07/02: rabartl:  The statement below should be the correct way to
00452  * get the list of valid parameters but it seems to be causing problems so
00453  * I am commenting it out for now.
00454  */
00455 /*
00456     pl->sublist(
00457       MLSettings_name, false,
00458       "Parameters directly accepted by ML_Epetra interface."
00459       ).setParameters(*rcp(ML_Epetra::GetValidMLPParameters()));
00460 */
00461     
00462     {
00463       ParameterList &mlSettingsPL = pl->sublist(
00464         MLSettings_name, false,
00465         "Sampling of the parameters directly accepted by ML\n"
00466         "This list of parameters is generated by combining all of\n"
00467         "the parameters set for all of the default problem types supported\n"
00468         "by ML.  Therefore, do not assume these parameters are at values that\n"
00469         "are reasonable to ML.  This list is just to give a sense of some of\n"
00470         "the parameters that ML accepts.  Consult ML documentation on how to\n"
00471         "set these parameters.  Also, you can print the parameter list after\n"
00472         "it is used and see what defaults where set for each default problem\n"
00473         "type.  Warning! the parameters in this sublist are currently *not*\n"
00474         "being validated by ML!"
00475         );
00476       //std::cout << "\nMLSettings doc before = " << pl->getEntryPtr(MLSettings_name)->docString() << "\n";
00477       { // Set of valid parameters (but perhaps not acceptable values)
00478         for (
00479           int i = 0;
00480           i < implicit_cast<int>(BaseMethodDefaults_valueNames.size());
00481           ++i
00482           )
00483         {
00484           ParameterList defaultParams;
00485           const std::string defaultTypeStr = BaseMethodDefaults_valueNames[i];
00486           if (defaultTypeStr != BaseMethodDefaults_valueNames_none) {
00487             TEUCHOS_TEST_FOR_EXCEPTION(
00488               0!=ML_Epetra::SetDefaults(defaultTypeStr,defaultParams)
00489               ,Teuchos::Exceptions::InvalidParameterValue
00490               ,"Error, the ML problem type \"" << defaultTypeStr
00491               << "\' is not recognized by ML!"
00492               );
00493             mlSettingsPL.setParameters(defaultParams);
00494           }
00495         }
00496       }
00497     }
00498 
00499     validPL = pl;
00500 
00501   }
00502 
00503   return validPL;
00504 
00505 }
00506 
00507 
00508 // Public functions overridden from Teuchos::Describable
00509 
00510 
00511 std::string MLPreconditionerFactory::description() const
00512 {
00513   std::ostringstream oss;
00514   oss << "Thyra::MLPreconditionerFactory";
00515   return oss.str();
00516 }
00517 
00518 
00519 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines