Stratimikos Version of the Day
Thyra_IfpackPreconditionerFactory.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_IfpackPreconditionerFactory.hpp"
00031 #include "Thyra_EpetraOperatorViewExtractorStd.hpp"
00032 #include "Thyra_EpetraLinearOp.hpp"
00033 #include "Thyra_DefaultPreconditioner.hpp"
00034 #include "Ifpack_ValidParameters.h"
00035 #include "Ifpack_Preconditioner.h"
00036 #include "Ifpack.h"
00037 #include "Epetra_RowMatrix.h"
00038 #include "Teuchos_TimeMonitor.hpp"
00039 #include "Teuchos_dyn_cast.hpp"
00040 #include "Teuchos_implicit_cast.hpp"
00041 #include "Teuchos_StandardParameterEntryValidators.hpp"
00042 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00043 #include "Teuchos_ValidatorXMLConverterDB.hpp"
00044 #include "Teuchos_StaticSetupMacro.hpp"
00045 
00046 
00047 namespace {
00048 
00049 Teuchos::RCP<Teuchos::Time> overallTimer, creationTimer, factorizationTimer;
00050 
00051 const std::string Ifpack_name = "Ifpack";
00052 
00053 const std::string IfpackSettings_name = "Ifpack Settings";
00054 
00055 const std::string PrecType_name = "Prec Type";
00056 Teuchos::RCP<Teuchos::StringToIntegralParameterEntryValidator<Ifpack::EPrecType> >
00057 precTypeValidator; // Will be setup below!
00058 const Ifpack::EPrecType PrecType_default = Ifpack::ILU;
00059 const std::string PrecTypeName_default = Ifpack::precTypeNames[PrecType_default];
00060 
00061 const std::string Overlap_name = "Overlap";
00062 const int Overlap_default = 0;
00063 
00064 
00065 TEUCHOS_STATIC_SETUP()
00066 {
00067   TEUCHOS_ADD_STRINGTOINTEGRALVALIDATOR_CONVERTER(Ifpack::EPrecType);
00068 }
00069 
00070 
00071 } // namespace
00072 
00073 namespace Thyra {
00074 
00075 // Constructors/initializers/accessors
00076 
00077 IfpackPreconditionerFactory::IfpackPreconditionerFactory()
00078   :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
00079   ,precType_(PrecType_default)
00080   ,overlap_(Overlap_default)
00081 {
00082   initializeTimers();
00083   getValidParameters(); // Make sure validators get created!
00084 }
00085 
00086 // Overridden from PreconditionerFactoryBase
00087 
00088 bool IfpackPreconditionerFactory::isCompatible(
00089   const LinearOpSourceBase<double> &fwdOpSrc
00090   ) const
00091 {
00092   using Teuchos::outArg;
00093   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00094   EOpTransp epetraFwdOpTransp;
00095   EApplyEpetraOpAs epetraFwdOpApplyAs;
00096   EAdjointEpetraOp epetraFwdOpAdjointSupport;
00097   double epetraFwdOpScalar;
00098   epetraFwdOpViewExtractor_->getEpetraOpView(
00099     fwdOpSrc.getOp(), 
00100     outArg(epetraFwdOp), outArg(epetraFwdOpTransp),
00101     outArg(epetraFwdOpApplyAs), outArg(epetraFwdOpAdjointSupport),
00102     outArg(epetraFwdOpScalar)
00103     );
00104   if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
00105     return false;
00106   return true;
00107 }
00108 
00109 bool IfpackPreconditionerFactory::applySupportsConj(EConj conj) const
00110 {
00111   return true;
00112 }
00113 
00114 bool IfpackPreconditionerFactory::applyTransposeSupportsConj(EConj conj) const
00115 {
00116   return false; // See comment below
00117 }
00118 
00119 Teuchos::RCP<PreconditionerBase<double> >
00120 IfpackPreconditionerFactory::createPrec() const
00121 {
00122   return Teuchos::rcp(new DefaultPreconditioner<double>());
00123 }
00124 
00125 void IfpackPreconditionerFactory::initializePrec(
00126   const Teuchos::RCP<const LinearOpSourceBase<double> >    &fwdOpSrc
00127   ,PreconditionerBase<double>                                      *prec
00128   ,const ESupportSolveUse                                           supportSolveUse
00129   ) const
00130 {
00131   using Teuchos::outArg;
00132   using Teuchos::OSTab;
00133   using Teuchos::dyn_cast;
00134   using Teuchos::RCP;
00135   using Teuchos::null;
00136   using Teuchos::rcp;
00137   using Teuchos::rcp_dynamic_cast;
00138   using Teuchos::rcp_const_cast;
00139   using Teuchos::set_extra_data;
00140   using Teuchos::get_optional_extra_data;
00141   using Teuchos::implicit_cast;
00142   Teuchos::Time totalTimer(""), timer("");
00143   totalTimer.start(true);
00144 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
00145   Teuchos::TimeMonitor overallTimeMonitor(*overallTimer);
00146 #endif
00147   const Teuchos::RCP<Teuchos::FancyOStream> out       = this->getOStream();
00148   const Teuchos::EVerbosityLevel                    verbLevel = this->getVerbLevel();
00149   Teuchos::OSTab tab(out);
00150   if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
00151     *out << "\nEntering Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
00152 #ifdef TEUCHOS_DEBUG
00153   TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00154   TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
00155 #endif
00156   Teuchos::RCP<const LinearOpBase<double> >
00157     fwdOp = fwdOpSrc->getOp();
00158 #ifdef TEUCHOS_DEBUG
00159   TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00160 #endif
00161   //
00162   // Unwrap and get the forward Epetra_Operator object
00163   //
00164   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00165   EOpTransp epetraFwdOpTransp;
00166   EApplyEpetraOpAs epetraFwdOpApplyAs;
00167   EAdjointEpetraOp epetraFwdOpAdjointSupport;
00168   double epetraFwdOpScalar;
00169   epetraFwdOpViewExtractor_->getEpetraOpView(
00170     fwdOp,
00171     outArg(epetraFwdOp), outArg(epetraFwdOpTransp),
00172     outArg(epetraFwdOpApplyAs), outArg(epetraFwdOpAdjointSupport),
00173     outArg(epetraFwdOpScalar)
00174     );
00175   // Validate what we get is what we need
00176   RCP<const Epetra_RowMatrix>
00177     epetraFwdRowMat = rcp_dynamic_cast<const Epetra_RowMatrix>(epetraFwdOp,true);
00178   TEUCHOS_TEST_FOR_EXCEPTION(
00179     epetraFwdOpApplyAs != EPETRA_OP_APPLY_APPLY, std::logic_error
00180     ,"Error, incorrect apply mode for an Epetra_RowMatrix"
00181     );
00182   //
00183   // Get the concrete preconditioner object
00184   //
00185   DefaultPreconditioner<double>
00186     *defaultPrec = &Teuchos::dyn_cast<DefaultPreconditioner<double> >(*prec);
00187   //
00188   // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
00189   //
00190   RCP<EpetraLinearOp>
00191     epetra_precOp = rcp_dynamic_cast<EpetraLinearOp>(defaultPrec->getNonconstUnspecifiedPrecOp(),true);
00192   //
00193   // Get the embedded Ifpack_Preconditioner object if it exists
00194   //
00195   Teuchos::RCP<Ifpack_Preconditioner>
00196     ifpack_precOp;
00197   if(epetra_precOp.get())
00198     ifpack_precOp = rcp_dynamic_cast<Ifpack_Preconditioner>(epetra_precOp->epetra_op(),true);
00199   //
00200   // Get the attached forward operator if it exists and make sure that it matches
00201   //
00202   if(ifpack_precOp.get()) {
00203     // ToDo: Get the forward operator and make sure that it matches what is
00204     // already being used!
00205   }
00206   //
00207   // Permform initialization if needed
00208   //
00209   //const bool startingOver = (ifpack_precOp.get() == NULL);
00210   const bool startingOver = true;
00211   // ToDo: Comment back in the above original version of startingOver to allow
00212   // for resuse.  Rob H. just pointed out to me that a memory leak is being
00213   // created when you just call Ifpack_ILU::Compute() over and over again.
00214   // Rob H. said that he will check in a fix the the development branch when
00215   // he can.
00216   if(startingOver) {
00217     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00218       *out << "\nCreating the initial Ifpack_Preconditioner object of type \'"<<Ifpack::toString(precType_)<<"\' ...\n";
00219     timer.start(true);
00220 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
00221     Teuchos::TimeMonitor creationTimeMonitor(*creationTimer);
00222 #endif
00223     // Create the initial preconditioner
00224     ifpack_precOp = rcp(
00225       ::Ifpack::Create(
00226         precType_
00227         ,const_cast<Epetra_RowMatrix*>(&*epetraFwdRowMat)
00228         ,overlap_
00229         )
00230       );
00231     timer.stop();
00232     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00233       OSTab(out).o() <<"=> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
00234     // Set parameters if the list exists
00235     if(paramList_.get()) {
00236       Teuchos::ParameterList
00237         &ifpackSettingsPL = paramList_->sublist(IfpackSettings_name);
00238       // Above will create new sublist if it does not exist!
00239       TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->SetParameters(ifpackSettingsPL));
00240       // Above, I have not idea how any error messages for a mistake will be
00241       // reported back to the user!
00242     }
00243     // Initialize the structure for the preconditioner
00244     TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->Initialize());
00245   }
00246   //
00247   // Attach the epetraFwdOp to the ifpack_precOp to guarantee that it will not go away
00248   //
00249   set_extra_data(epetraFwdOp, "IFPF::epetraFwdOp", Teuchos::inOutArg(ifpack_precOp),
00250     Teuchos::POST_DESTROY, false);
00251   //
00252   // Update the factorization
00253   //
00254   {
00255     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00256       *out << "\nComputing the preconditioner ...\n";
00257 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
00258     Teuchos::TimeMonitor factorizationTimeMonitor(*factorizationTimer);
00259 #endif
00260     timer.start(true);
00261     TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->Compute());
00262     timer.stop();
00263     if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00264       OSTab(out).o() <<"=> Setup time = "<<timer.totalElapsedTime()<<" sec\n";
00265   }
00266   //
00267   // Compute the conditioner number estimate if asked
00268   //
00269 
00270   // ToDo: Implement
00271 
00272   //
00273   // Attach fwdOp to the ifpack_precOp
00274   //
00275   set_extra_data(fwdOpSrc, "IFPF::fwdOpSrc", Teuchos::inOutArg(ifpack_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     ifpack_precOp
00285     ,epetraFwdOpTransp
00286     ,EPETRA_OP_APPLY_APPLY_INVERSE
00287     ,EPETRA_OP_ADJOINT_SUPPORTED
00288     );
00289   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_MEDIUM)) {
00290     *out << "\nDescription of created preconditioner:\n";
00291     OSTab tab2(out);
00292     ifpack_precOp->Print(*out);
00293   }
00294 
00295   //
00296   // Initialize the preconditioner
00297   //
00298   defaultPrec->initializeUnspecified(
00299     Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp)
00300     );
00301   totalTimer.stop();
00302   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00303     *out << "\nTotal time in IfpackPreconditionerFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
00304   if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
00305     *out << "\nLeaving Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
00306 }
00307 
00308 void IfpackPreconditionerFactory::uninitializePrec(
00309   PreconditionerBase<double>                                *prec
00310   ,Teuchos::RCP<const LinearOpSourceBase<double> >  *fwdOpSrc
00311   ,ESupportSolveUse                                         *supportSolveUse
00312   ) const
00313 {
00314   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
00315 }
00316 
00317 // Overridden from ParameterListAcceptor
00318 
00319 void IfpackPreconditionerFactory::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
00320 {
00321   TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
00322   paramList->validateParameters(*this->getValidParameters(),2);
00323   // Note: The above validation will validate right down into the the sublist
00324   // named IfpackSettings_name!
00325   paramList_ = paramList;
00326   overlap_ = paramList_->get(Overlap_name,Overlap_default);
00327   std::ostringstream oss;
00328   oss << "(sub)list \""<<paramList->name()<<"\"parameter \"Prec Type\"";
00329   precType_ =
00330     ( paramList_.get()
00331       ? precTypeValidator->getIntegralValue(*paramList_,PrecType_name,PrecTypeName_default)
00332       : PrecType_default
00333       );
00334   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00335 #ifdef TEUCHOS_DEBUG
00336   // Validate my use of the parameters!
00337   paramList->validateParameters(*this->getValidParameters(),1);
00338 #endif
00339 }
00340 
00341 Teuchos::RCP<Teuchos::ParameterList>
00342 IfpackPreconditionerFactory::getNonconstParameterList()
00343 {
00344   return paramList_;
00345 }
00346 
00347 Teuchos::RCP<Teuchos::ParameterList>
00348 IfpackPreconditionerFactory::unsetParameterList()
00349 {
00350   Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
00351   paramList_ = Teuchos::null;
00352   return _paramList;
00353 }
00354 
00355 Teuchos::RCP<const Teuchos::ParameterList>
00356 IfpackPreconditionerFactory::getParameterList() const
00357 {
00358   return paramList_;
00359 }
00360 
00361 Teuchos::RCP<const Teuchos::ParameterList>
00362 IfpackPreconditionerFactory::getValidParameters() const
00363 {
00364   using Teuchos::rcp;
00365   using Teuchos::rcp_implicit_cast;
00366   typedef Teuchos::ParameterEntryValidator PEV;
00367   static Teuchos::RCP<Teuchos::ParameterList> validParamList;
00368   if(validParamList.get()==NULL) {
00369     validParamList = Teuchos::rcp(new Teuchos::ParameterList(Ifpack_name));
00370     {
00371       // Create the validator for the preconditioner type!
00372       Teuchos::Array<std::string>
00373         precTypeNames;
00374       precTypeNames.insert(
00375         precTypeNames.begin(),
00376         &Ifpack::precTypeNames[0],
00377         &Ifpack::precTypeNames[0] + Ifpack::numPrecTypes
00378         );
00379       Teuchos::Array<Ifpack::EPrecType>
00380         precTypeValues;
00381       precTypeValues.insert(
00382         precTypeValues.begin(),
00383         &Ifpack::precTypeValues[0],
00384         &Ifpack::precTypeValues[0] + Ifpack::numPrecTypes
00385         );
00386       precTypeValidator = rcp(
00387         new Teuchos::StringToIntegralParameterEntryValidator<Ifpack::EPrecType>(
00388           precTypeNames,precTypeValues,PrecType_name
00389           )
00390         );
00391     }
00392     validParamList->set(
00393       PrecType_name, PrecTypeName_default,
00394       "Type of Ifpack preconditioner to use.",
00395       rcp_implicit_cast<const PEV>(precTypeValidator)
00396       );
00397     validParamList->set(
00398       Overlap_name, Overlap_default,
00399       "Number of rows/columns overlapped between subdomains in different"
00400       "\nprocesses in the additive Schwarz-type domain-decomposition preconditioners."
00401       );
00402     validParamList->sublist(
00403       IfpackSettings_name, false,
00404       "Preconditioner settings that are passed onto the Ifpack preconditioners themselves."
00405       ).setParameters(Ifpack_GetValidParameters());
00406     // Note that in the above setParameterList(...) function that we actually
00407     // validate down into the first level of this sublist.  Really the
00408     // Ifpack_Preconditioner objects themselves should do validation but we do
00409     // it ourselves taking the return from the Ifpack_GetValidParameters()
00410     // function as gospel!
00411     Teuchos::setupVerboseObjectSublist(&*validParamList);
00412   }
00413   return validParamList;
00414 }
00415 
00416 // Public functions overridden from Teuchos::Describable
00417 
00418 std::string IfpackPreconditionerFactory::description() const
00419 {
00420   std::ostringstream oss;
00421   oss << "Thyra::IfpackPreconditionerFactory{";
00422   oss << "precType=\"" << ::Ifpack::toString(precType_) << "\"";
00423   oss << ",overlap=" << overlap_;
00424   oss << "}";
00425   return oss.str();
00426 }
00427 
00428 // private
00429 
00430 void IfpackPreconditionerFactory::initializeTimers()
00431 {
00432   if(!overallTimer.get()) {
00433 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
00434     overallTimer       = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF");
00435     creationTimer      = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF:Creation");
00436     factorizationTimer = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF:Factorization");
00437 #endif
00438   }
00439 }
00440 
00441 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends