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

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