Thyra_MLPreconditionerFactory.cpp

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