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

Generated on Wed Jul 22 12:56:30 2009 for Stratimikos by  doxygen 1.5.8