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