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 = "DD";
00079 Teuchos::RCP<
00080   Teuchos::StringToIntegralParameterEntryValidator<EMLProblemType>
00081   >
00082 BaseMethodDefaults_validator;
00083   
00084 const std::string MLSettings_name = "ML Settings";
00085 
00086 
00087 } // namespace
00088 
00089 
00090 namespace Thyra {
00091 
00092 
00093 using Teuchos::RCP;
00094 using Teuchos::ParameterList;
00095 
00096 
00097 // Constructors/initializers/accessors
00098 
00099   
00100 MLPreconditionerFactory::MLPreconditionerFactory()
00101   :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
00102 {}
00103 
00104 
00105 // Overridden from PreconditionerFactoryBase
00106 
00107 
00108 bool MLPreconditionerFactory::isCompatible(
00109   const LinearOpSourceBase<double> &fwdOpSrc
00110   ) const
00111 {
00112   using Teuchos::outArg;
00113   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00114   EOpTransp epetraFwdOpTransp;
00115   EApplyEpetraOpAs epetraFwdOpApplyAs;
00116   EAdjointEpetraOp epetraFwdOpAdjointSupport;
00117   double epetraFwdOpScalar;
00118   Teuchos::RCP<const LinearOpBase<double> >
00119     fwdOp = fwdOpSrc.getOp();
00120   epetraFwdOpViewExtractor_->getEpetraOpView(
00121     fwdOp,
00122     outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
00123     outArg(epetraFwdOpApplyAs),
00124     outArg(epetraFwdOpAdjointSupport),
00125     outArg(epetraFwdOpScalar)
00126     );
00127   if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
00128     return false;
00129   return true;
00130 }
00131 
00132 
00133 bool MLPreconditionerFactory::applySupportsConj(EConj conj) const
00134 {
00135   return true;
00136 }
00137 
00138 
00139 bool MLPreconditionerFactory::applyTransposeSupportsConj(EConj conj) const
00140 {
00141   return false; // See comment below
00142 }
00143 
00144 
00145 RCP<PreconditionerBase<double> >
00146 MLPreconditionerFactory::createPrec() const
00147 {
00148   return Teuchos::rcp(new DefaultPreconditioner<double>());
00149 }
00150 
00151 
00152 void MLPreconditionerFactory::initializePrec(
00153   const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
00154   PreconditionerBase<double> *prec,
00155   const ESupportSolveUse supportSolveUse
00156   ) const
00157 {
00158   using Teuchos::outArg;
00159   using Teuchos::OSTab;
00160   using Teuchos::dyn_cast;
00161   using Teuchos::RCP;
00162   using Teuchos::null;
00163   using Teuchos::rcp;
00164   using Teuchos::rcp_dynamic_cast;
00165   using Teuchos::rcp_const_cast;
00166   using Teuchos::set_extra_data;
00167   using Teuchos::get_optional_extra_data;
00168   using Teuchos::implicit_cast;
00169   Teuchos::Time totalTimer(""), timer("");
00170   totalTimer.start(true);
00171   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00172   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00173   Teuchos::OSTab tab(out);
00174   if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
00175     *out << "\nEntering Thyra::MLPreconditionerFactory::initializePrec(...) ...\n";
00176 
00177   Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
00178 #ifdef _DEBUG
00179   TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00180   TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
00181 #endif
00182   //
00183   // Unwrap and get the forward Epetra_Operator object
00184   //
00185   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00186   EOpTransp epetraFwdOpTransp;
00187   EApplyEpetraOpAs epetraFwdOpApplyAs;
00188   EAdjointEpetraOp epetraFwdOpAdjointSupport;
00189   double epetraFwdOpScalar;
00190   epetraFwdOpViewExtractor_->getEpetraOpView(
00191     fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
00192     outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
00193                                              );
00194   // Validate what we get is what we need
00195   RCP<const Epetra_RowMatrix>
00196     epetraFwdRowMat = rcp_dynamic_cast<const Epetra_RowMatrix>(epetraFwdOp,true);
00197   TEUCHOS_TEST_FOR_EXCEPTION(
00198     epetraFwdOpApplyAs != EPETRA_OP_APPLY_APPLY, std::logic_error
00199     ,"Error, incorrect apply mode for an Epetra_RowMatrix"
00200     );
00201 
00202   //
00203   // Get the concrete preconditioner object
00204   //
00205   DefaultPreconditioner<double>
00206     *defaultPrec = &Teuchos::dyn_cast<DefaultPreconditioner<double> >(*prec);
00207   //
00208   // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
00209   //
00210   RCP<EpetraLinearOp>
00211     epetra_precOp = rcp_dynamic_cast<EpetraLinearOp>(defaultPrec->getNonconstUnspecifiedPrecOp(),true);
00212   //
00213   // Get the embedded ML_Epetra::MultiLevelPreconditioner object if it exists
00214   //
00215   Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> ml_precOp;
00216   if(epetra_precOp.get())
00217     ml_precOp = rcp_dynamic_cast<ML_Epetra::MultiLevelPreconditioner>(epetra_precOp->epetra_op(),true);
00218   //
00219   // Get the attached forward operator if it exists and make sure that it matches
00220   //
00221   if(ml_precOp!=Teuchos::null) {
00222     // Get the forward operator and make sure that it matches what is
00223     // already being used!
00224     const Epetra_RowMatrix & rm = ml_precOp->RowMatrix();
00225    
00226     TEUCHOS_TEST_FOR_EXCEPTION(
00227        &rm!=&*epetraFwdRowMat, std::logic_error
00228        ,"ML requires Epetra_RowMatrix to be the same for each initialization of the preconditioner"
00229        );
00230   }
00231   //
00232   // Perform initialization if needed
00233   //
00234   const bool startingOver = (ml_precOp.get() == NULL);
00235   if(startingOver) 
00236   {
00237     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00238       *out << "\nCreating the initial ML_Epetra::MultiLevelPreconditioner object...\n";
00239     timer.start(true);
00240     // Create the initial preconditioner: DO NOT compute it yet
00241     ml_precOp = rcp(
00242       new ML_Epetra::MultiLevelPreconditioner(
00243         *epetraFwdRowMat, paramList_->sublist(MLSettings_name),false
00244         )
00245       );
00246     
00247     timer.stop();
00248     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00249       OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
00250     // RAB: Above, I am just passing a string to ML::Create(...) in order
00251     // get this code written.  However, in the future, it would be good to
00252     // copy the contents of what is in ML::Create(...) into a local
00253     // function and then use switch(...) to create the initial
00254     // ML_Epetra::MultiLevelPreconditioner object.  This would result in better validation
00255     // and faster code.
00256     // Set parameters if the list exists
00257     if(paramList_.get())
00258       TEUCHOS_TEST_FOR_EXCEPT(
00259         0!=ml_precOp->SetParameterList(paramList_->sublist(MLSettings_name))
00260         );
00261     // Initialize the structure for the preconditioner
00262     //        TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->Initialize());
00263   }
00264   //
00265   // Attach the epetraFwdOp to the ml_precOp to guarantee that it will not go away
00266   //
00267   set_extra_data(epetraFwdOp, "IFPF::epetraFwdOp", Teuchos::inOutArg(ml_precOp),
00268     Teuchos::POST_DESTROY, false);
00269   //
00270   // Update the factorization
00271   //
00272   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00273     *out << "\nComputing the preconditioner ...\n";
00274   timer.start(true);
00275   TEUCHOS_TEST_FOR_EXCEPT(0!=ml_precOp->ComputePreconditioner());
00276   timer.stop();
00277   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00278     OSTab(out).o() <<"=> Setup time = "<<timer.totalElapsedTime()<<" sec\n";
00279   //
00280   // Compute the conditioner number estimate if asked
00281   //
00282 
00283   // ToDo: Implement
00284 
00285   //
00286   // Attach fwdOp to the ml_precOp
00287   //
00288   set_extra_data(fwdOp, "IFPF::fwdOp", Teuchos::inOutArg(ml_precOp),
00289     Teuchos::POST_DESTROY, false);
00290   //
00291   // Initialize the output EpetraLinearOp
00292   //
00293   if(startingOver) {
00294     epetra_precOp = rcp(new EpetraLinearOp);
00295   }
00296   epetra_precOp->initialize(
00297     ml_precOp
00298     ,epetraFwdOpTransp
00299     ,EPETRA_OP_APPLY_APPLY_INVERSE
00300     ,EPETRA_OP_ADJOINT_UNSUPPORTED  // ToDo: Look into adjoints again.
00301     );
00302   //
00303   // Initialize the preconditioner
00304   //
00305   defaultPrec->initializeUnspecified(
00306     Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp)
00307     );
00308   totalTimer.stop();
00309   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00310     *out << "\nTotal time in MLPreconditionerFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
00311   if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
00312     *out << "\nLeaving Thyra::MLPreconditionerFactory::initializePrec(...) ...\n";
00313 }
00314 
00315 
00316 void MLPreconditionerFactory::uninitializePrec(
00317   PreconditionerBase<double> *prec,
00318   Teuchos::RCP<const LinearOpSourceBase<double> > *fwdOp,
00319   ESupportSolveUse *supportSolveUse
00320   ) const
00321 {
00322   TEUCHOS_TEST_FOR_EXCEPT(true);
00323 }
00324 
00325 
00326 // Overridden from ParameterListAcceptor
00327 
00328 
00329 void MLPreconditionerFactory::setParameterList(
00330   Teuchos::RCP<ParameterList> const& paramList
00331   )
00332 {
00333   TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
00334   paramList->validateParameters(*this->getValidParameters(),0);
00335   paramList_ = paramList;
00336   const EMLProblemType
00337     defaultType = BaseMethodDefaults_validator->getIntegralValue(
00338       *paramList_,BaseMethodDefaults_name,BaseMethodDefaults_default
00339       );
00340   if( ML_PROBTYPE_NONE != defaultType ) {
00341     const std::string
00342       defaultTypeStr = BaseMethodDefaults_valueNames[defaultType];
00343     Teuchos::ParameterList defaultParams;
00344     TEUCHOS_TEST_FOR_EXCEPTION(
00345       0!=ML_Epetra::SetDefaults(defaultTypeStr,defaultParams)
00346       ,Teuchos::Exceptions::InvalidParameterValue
00347       ,"Error, the ML problem type \"" << defaultTypeStr << "\' is not recongnised by ML!"
00348       );
00349     // Note, the only way the above exception message could be generated is if
00350     // a default problem type was removed from ML_Epetra::SetDefaults(...).
00351     // When a new problem type is added to this function, it must be added to
00352     // our enum EMLProblemType along with associated objects ...  In other
00353     // words, this adapter must be maintained as ML is maintained.  An
00354     // alternative design would be to just pass in whatever string to this
00355     // function.  This would improve maintainability but it would not generate
00356     // very good error messages when a bad string was passed in.  Currenly,
00357     // the error message attached to the exception will contain the list of
00358     // valid problem types.
00359     paramList_->sublist(MLSettings_name).setParametersNotAlreadySet(
00360       defaultParams);
00361   }
00362 #ifdef TEUCHOS_DEBUG
00363   paramList->validateParameters(*this->getValidParameters(),0);
00364 #endif
00365 }
00366 
00367 
00368 RCP<ParameterList>
00369 MLPreconditionerFactory::getNonconstParameterList()
00370 {
00371   return paramList_;
00372 }
00373 
00374 
00375 RCP<ParameterList>
00376 MLPreconditionerFactory::unsetParameterList()
00377 {
00378   Teuchos::RCP<ParameterList> _paramList = paramList_;
00379   paramList_ = Teuchos::null;
00380   return _paramList;
00381 }
00382 
00383 
00384 RCP<const ParameterList>
00385 MLPreconditionerFactory::getParameterList() const
00386 {
00387   return paramList_;
00388 }
00389 
00390 
00391 RCP<const ParameterList>
00392 MLPreconditionerFactory::getValidParameters() const
00393 {
00394 
00395   using Teuchos::rcp;
00396   using Teuchos::tuple;
00397   using Teuchos::implicit_cast;
00398   using Teuchos::rcp_implicit_cast;
00399   typedef Teuchos::ParameterEntryValidator PEV;
00400 
00401   static RCP<const ParameterList> validPL;
00402 
00403   if(is_null(validPL)) {
00404 
00405     RCP<ParameterList>
00406       pl = rcp(new ParameterList());
00407 
00408     BaseMethodDefaults_validator = rcp(
00409       new Teuchos::StringToIntegralParameterEntryValidator<EMLProblemType>(
00410         BaseMethodDefaults_valueNames,
00411         tuple<std::string>(
00412           "Do not set any default parameters",
00413           "Set default parameters for a smoothed aggregation method",
00414           "Set default parameters for a domain decomposition method",
00415           "Set default parameters for a domain decomposition method special to ML",
00416           "Set default parameters for a Maxwell-type of linear operator"
00417           ),
00418         tuple<EMLProblemType>(
00419           ML_PROBTYPE_NONE,
00420           ML_PROBTYPE_SMOOTHED_AGGREGATION,
00421           ML_PROBTYPE_DOMAIN_DECOMPOSITION,
00422           ML_PROBTYPE_DOMAIN_DECOMPOSITION_ML,
00423           ML_PROBTYPE_MAXWELL
00424           ),
00425         BaseMethodDefaults_name
00426         )
00427       );
00428 
00429     pl->set(BaseMethodDefaults_name,BaseMethodDefaults_default,
00430       "Select the default method type which also sets parameter defaults\n"
00431       "in the sublist \"" + MLSettings_name + "\"!",
00432       rcp_implicit_cast<const PEV>(BaseMethodDefaults_validator)
00433       );
00434 
00435 /* 2007/07/02: rabartl:  The statement below should be the correct way to
00436  * get the list of valid parameters but it seems to be causing problems so
00437  * I am commenting it out for now.
00438  */
00439 /*
00440     pl->sublist(
00441       MLSettings_name, false,
00442       "Parameters directly accepted by ML_Epetra interface."
00443       ).setParameters(*rcp(ML_Epetra::GetValidMLPParameters()));
00444 */
00445     
00446     {
00447       ParameterList &mlSettingsPL = pl->sublist(
00448         MLSettings_name, false,
00449         "Sampling of the parameters directly accepted by ML\n"
00450         "This list of parameters is generated by combining all of\n"
00451         "the parameters set for all of the default problem types supported\n"
00452         "by ML.  Therefore, do not assume these parameters are at values that\n"
00453         "are reasonable to ML.  This list is just to give a sense of some of\n"
00454         "the parameters that ML accepts.  Consult ML documentation on how to\n"
00455         "set these parameters.  Also, you can print the parameter list after\n"
00456         "it is used and see what defaults where set for each default problem\n"
00457         "type.  Warning! the parameters in this sublist are currently *not*\n"
00458         "being validated by ML!"
00459         );
00460       //std::cout << "\nMLSettings doc before = " << pl->getEntryPtr(MLSettings_name)->docString() << "\n";
00461       { // Set of valid parameters (but perhaps not accetable values)
00462         for (
00463           int i = 0;
00464           i < implicit_cast<int>(BaseMethodDefaults_valueNames.size());
00465           ++i
00466           )
00467         {
00468           ParameterList defaultParams;
00469           const std::string defaultTypeStr = BaseMethodDefaults_valueNames[i];
00470           if (defaultTypeStr != BaseMethodDefaults_valueNames_none) {
00471             TEUCHOS_TEST_FOR_EXCEPTION(
00472               0!=ML_Epetra::SetDefaults(defaultTypeStr,defaultParams)
00473               ,Teuchos::Exceptions::InvalidParameterValue
00474               ,"Error, the ML problem type \"" << defaultTypeStr
00475               << "\' is not recongnised by ML!"
00476               );
00477             mlSettingsPL.setParameters(defaultParams);
00478           }
00479         }
00480       }
00481     }
00482 
00483     validPL = pl;
00484 
00485   }
00486 
00487   return validPL;
00488 
00489 }
00490 
00491 
00492 // Public functions overridden from Teuchos::Describable
00493 
00494 
00495 std::string MLPreconditionerFactory::description() const
00496 {
00497   std::ostringstream oss;
00498   oss << "Thyra::MLPreconditionerFactory";
00499   return oss.str();
00500 }
00501 
00502 
00503 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines