Thyra_IfpackPreconditionerFactory.cpp

Go to the documentation of this file.
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 
00041 namespace {
00042 
00043 Teuchos::RefCountPtr<Teuchos::Time> overallTimer, creationTimer, factorizationTimer;
00044 
00045 const std::string Ifpack_name = "Ifpack";
00046 const std::string IfpackSettings_name = "Ifpack Settings";
00047 const std::string PrecType_name = "Prec Type";
00048 const std::string Overlap_name = "Overlap";
00049 
00050 } // namespace
00051 
00052 namespace Thyra {
00053 
00054 // Constructors/initializers/accessors
00055 
00056 IfpackPreconditionerFactory::IfpackPreconditionerFactory()
00057   :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
00058    ,precType_(Ifpack::ILU)
00059    ,overlap_(0)
00060 {
00061   initializeTimers();
00062 }
00063 
00064 // Overridden from PreconditionerFactoryBase
00065 
00066 bool IfpackPreconditionerFactory::isCompatible(
00067   const LinearOpSourceBase<double> &fwdOpSrc
00068   ) const
00069 {
00070   Teuchos::RefCountPtr<const Epetra_Operator> epetraFwdOp;
00071   ETransp                                     epetraFwdOpTransp;
00072   EApplyEpetraOpAs                            epetraFwdOpApplyAs;
00073   EAdjointEpetraOp                            epetraFwdOpAdjointSupport;
00074   double                                      epetraFwdOpScalar;
00075   epetraFwdOpViewExtractor_->getEpetraOpView(
00076     fwdOpSrc.getOp()
00077     ,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
00078     );
00079   if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
00080     return false;
00081   return true;
00082 }
00083 
00084 bool IfpackPreconditionerFactory::applySupportsConj(EConj conj) const
00085 {
00086   return true;
00087 }
00088 
00089 bool IfpackPreconditionerFactory::applyTransposeSupportsConj(EConj conj) const
00090 {
00091   return false; // See comment below
00092 }
00093 
00094 Teuchos::RefCountPtr<PreconditionerBase<double> >
00095 IfpackPreconditionerFactory::createPrec() const
00096 {
00097   return Teuchos::rcp(new DefaultPreconditioner<double>());
00098 }
00099 
00100 void IfpackPreconditionerFactory::initializePrec(
00101   const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >    &fwdOpSrc
00102   ,PreconditionerBase<double>                                      *prec
00103   ,const ESupportSolveUse                                           supportSolveUse
00104   ) const
00105 {
00106   using Teuchos::OSTab;
00107   using Teuchos::dyn_cast;
00108   using Teuchos::RefCountPtr;
00109   using Teuchos::null;
00110   using Teuchos::rcp;
00111   using Teuchos::rcp_dynamic_cast;
00112   using Teuchos::rcp_const_cast;
00113   using Teuchos::set_extra_data;
00114   using Teuchos::get_optional_extra_data;
00115   Teuchos::Time totalTimer(""), timer("");
00116   totalTimer.start(true);
00117   Teuchos::TimeMonitor overallTimeMonitor(*overallTimer);
00118   const Teuchos::RefCountPtr<Teuchos::FancyOStream> out       = this->getOStream();
00119   const Teuchos::EVerbosityLevel                    verbLevel = this->getVerbLevel();
00120   Teuchos::OSTab tab(out);
00121   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00122     *out << "\nEntering Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
00123 #ifdef TEUCHOS_DEBUG
00124   TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00125   TEST_FOR_EXCEPT(prec==NULL);
00126 #endif
00127   Teuchos::RefCountPtr<const LinearOpBase<double> >
00128     fwdOp = fwdOpSrc->getOp();
00129 #ifdef TEUCHOS_DEBUG
00130   TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00131 #endif
00132   //
00133   // Unwrap and get the forward Epetra_Operator object
00134   //
00135   Teuchos::RefCountPtr<const Epetra_Operator> epetraFwdOp;
00136   ETransp                                     epetraFwdOpTransp;
00137   EApplyEpetraOpAs                            epetraFwdOpApplyAs;
00138   EAdjointEpetraOp                            epetraFwdOpAdjointSupport;
00139   double                                      epetraFwdOpScalar;
00140   epetraFwdOpViewExtractor_->getEpetraOpView(
00141     fwdOp,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
00142     );
00143   // Validate what we get is what we need
00144   RefCountPtr<const Epetra_RowMatrix>
00145     epetraFwdRowMat = rcp_dynamic_cast<const Epetra_RowMatrix>(epetraFwdOp,true);
00146   TEST_FOR_EXCEPTION(
00147     epetraFwdOpApplyAs != EPETRA_OP_APPLY_APPLY, std::logic_error
00148     ,"Error, incorrect apply mode for an Epetra_RowMatrix"
00149     );
00150   //
00151   // Get the concrete precondtioner object
00152   //
00153   DefaultPreconditioner<double>
00154     *defaultPrec = &Teuchos::dyn_cast<DefaultPreconditioner<double> >(*prec);
00155   //
00156   // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
00157   //
00158   RefCountPtr<EpetraLinearOp>
00159     epetra_precOp = rcp_dynamic_cast<EpetraLinearOp>(defaultPrec->getNonconstUnspecifiedPrecOp(),true);
00160   //
00161   // Get the embedded Ifpack_Preconditioner object if it exists
00162   //
00163   Teuchos::RefCountPtr<Ifpack_Preconditioner>
00164     ifpack_precOp;
00165   if(epetra_precOp.get())
00166     ifpack_precOp = rcp_dynamic_cast<Ifpack_Preconditioner>(epetra_precOp->epetra_op(),true);
00167   //
00168   // Get the attached forward operator if it exists and make sure that it matches
00169   //
00170   if(ifpack_precOp.get()) {
00171     // ToDo: Get the forward operator and make sure that it matches what is
00172     // already being used!
00173   }
00174   //
00175   // Permform initialization if needed
00176   //
00177   //const bool startingOver = (ifpack_precOp.get() == NULL);
00178   const bool startingOver = true;
00179   // ToDo: Comment back in the above original version of startingOver to allow
00180   // for resuse.  Rob H. just pointed out to me that a memory leak is being
00181   // created when you just call Ifpack_ILU::Compute() over and over again.
00182   // Rob H. said that he will check in a fix the the development branch when
00183   // he can.
00184   if(startingOver) {
00185     if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00186       *out << "\nCreating the initial Ifpack_Preconditioner object of type \'"<<toString(precType_)<<"\' ...\n";
00187     timer.start(true);
00188     Teuchos::TimeMonitor creationTimeMonitor(*creationTimer);
00189     // Create the initial preconditioner
00190     ::Ifpack ifpackFcty; // Should be a static function!
00191     const std::string &precTypeName = toString(precType_);
00192     ifpack_precOp = rcp(
00193       ifpackFcty.Create(
00194         precTypeName
00195         ,const_cast<Epetra_RowMatrix*>(&*epetraFwdRowMat)
00196         ,overlap_
00197         )
00198       );
00199     timer.stop();
00200     if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00201       *OSTab(out).getOStream() <<"\n=> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
00202     // RAB: Above, I am just passing a string to Ifpack::Create(...) in order
00203     // get this code written.  However, in the future, it would be good to
00204     // copy the contents of what is in Ifpack::Create(...) into a local
00205     // function and then use switch(...) to create the initial
00206     // Ifpack_Preconditioner object.  This would result in better validation
00207     // and faster code.
00208     TEST_FOR_EXCEPTION(
00209       ifpack_precOp.get()==NULL, std::logic_error
00210       ,"Error, Ifpack::Create(precTypeName,...) returned NULL for precType = \""<<precTypeName<<"\"!"
00211       );
00212     // Set parameters if the list exists
00213     if(paramList_.get())
00214       TEST_FOR_EXCEPT(0!=ifpack_precOp->SetParameters(paramList_->sublist(IfpackSettings_name))); // This will create new sublist if it does not exist!
00215     // Initailize the structure for the preconditioner
00216     TEST_FOR_EXCEPT(0!=ifpack_precOp->Initialize());
00217   }
00218   //
00219   // Attach the epetraFwdOp to the ifpack_precOp to guarantee that it will not go away
00220   //
00221   set_extra_data(epetraFwdOp,"IFPF::epetraFwdOp",&ifpack_precOp,Teuchos::POST_DESTROY,false);
00222   //
00223   // Update the factorization
00224   //
00225   if(1){
00226     if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00227       *out << "\nComputing the factorization of the preconditioner ...\n";
00228     Teuchos::TimeMonitor factorizationTimeMonitor(*factorizationTimer);
00229     timer.start(true);
00230     TEST_FOR_EXCEPT(0!=ifpack_precOp->Compute());
00231     timer.stop();
00232     if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00233       *OSTab(out).getOStream() <<"\n=> Factorization time = "<<timer.totalElapsedTime()<<" sec\n";
00234   }
00235   //
00236   // Compute the conditioner number estimate if asked
00237   //
00238 
00239   // ToDo: Implement
00240 
00241   //
00242   // Attach fwdOp to the ifpack_precOp
00243   //
00244   set_extra_data(fwdOpSrc,"IFPF::fwdOpSrc",&ifpack_precOp,Teuchos::POST_DESTROY,false);
00245   //
00246   // Initialize the output EpetraLinearOp
00247   //
00248   if(startingOver) {
00249     epetra_precOp = rcp(new EpetraLinearOp);
00250   }
00251   epetra_precOp->initialize(
00252     ifpack_precOp
00253     ,epetraFwdOpTransp
00254     ,EPETRA_OP_APPLY_APPLY_INVERSE
00255     ,EPETRA_OP_ADJOINT_UNSUPPORTED  // ToDo: Look into adjoints again.
00256     );
00257   //
00258   // Initialize the preconditioner
00259   //
00260   defaultPrec->initializeUnspecified(
00261     Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp)
00262     );
00263   totalTimer.stop();
00264   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00265     *out
00266       << "\nTotal time = "<<totalTimer.totalElapsedTime()<<" sec\n"
00267       << "\nLeaving Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
00268 }
00269 
00270 void IfpackPreconditionerFactory::uninitializePrec(
00271   PreconditionerBase<double>                                *prec
00272   ,Teuchos::RefCountPtr<const LinearOpSourceBase<double> >  *fwdOpSrc
00273   ,ESupportSolveUse                                         *supportSolveUse
00274   ) const
00275 {
00276   TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
00277 }
00278 
00279 // Overridden from ParameterListAcceptor
00280 
00281 void IfpackPreconditionerFactory::setParameterList(Teuchos::RefCountPtr<Teuchos::ParameterList> const& paramList)
00282 {
00283   TEST_FOR_EXCEPT(paramList.get()==NULL);
00284   paramList->validateParameters(*this->getValidParameters(),1);
00285   paramList_ = paramList;
00286   overlap_ = paramList_->get(Overlap_name,overlap_);
00287   std::ostringstream oss;
00288   oss << "(sub)list \""<<paramList->name()<<"\"parameter \"Prec Type\"";
00289   precType_ = Ifpack::precTypeNameToEnum(
00290     paramList_->get(PrecType_name,toString(precType_))
00291     ,oss.str()
00292     );
00293 }
00294 
00295 Teuchos::RefCountPtr<Teuchos::ParameterList>
00296 IfpackPreconditionerFactory::getParameterList()
00297 {
00298   return paramList_;
00299 }
00300 
00301 Teuchos::RefCountPtr<Teuchos::ParameterList>
00302 IfpackPreconditionerFactory::unsetParameterList()
00303 {
00304   Teuchos::RefCountPtr<Teuchos::ParameterList> _paramList = paramList_;
00305   paramList_ = Teuchos::null;
00306   return _paramList;
00307 }
00308 
00309 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00310 IfpackPreconditionerFactory::getParameterList() const
00311 {
00312   return paramList_;
00313 }
00314 
00315 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00316 IfpackPreconditionerFactory::getValidParameters() const
00317 {
00318   return generateAndGetValidParameters();
00319 }
00320 
00321 // Public functions overridden from Teuchos::Describable
00322 
00323 std::string IfpackPreconditionerFactory::description() const
00324 {
00325   std::ostringstream oss;
00326   oss << "Thyra::IfpackPreconditionerFactory{";
00327   oss << "precType=\"" << toString(precType_) << "\"";
00328   oss << ",overlap=" << overlap_;
00329   oss << "}";
00330   return oss.str();
00331 }
00332 
00333 // private
00334 
00335 void IfpackPreconditionerFactory::initializeTimers()
00336 {
00337   if(!overallTimer.get()) {
00338     overallTimer       = Teuchos::TimeMonitor::getNewTimer("IfpackPF");
00339     creationTimer      = Teuchos::TimeMonitor::getNewTimer("IfpackPF:Creation");
00340     factorizationTimer = Teuchos::TimeMonitor::getNewTimer("IfpackPF:Factorization");
00341   }
00342 }
00343 
00344 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00345 IfpackPreconditionerFactory::generateAndGetValidParameters()
00346 {
00347   static Teuchos::RefCountPtr<Teuchos::ParameterList> validParamList;
00348   if(validParamList.get()==NULL) {
00349     validParamList = Teuchos::rcp(new Teuchos::ParameterList(Ifpack_name));
00350     validParamList->set(PrecType_name,"ILU");
00351     validParamList->set(Overlap_name,0);
00352     validParamList->sublist(IfpackSettings_name).setParameters(Ifpack_GetValidParameters());
00353   }
00354   return validParamList;
00355 }
00356 
00357 } // namespace Thyra

Generated on Thu Sep 18 12:37:22 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1