Amesos Package Browser (Single Doxygen Collection) Development
Thyra_AmesosLinearOpWithSolveFactory.cpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //                Amesos: Direct Sparse Solver Package
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 */
00030 
00031 #ifndef __sun
00032 
00033 #include "Thyra_AmesosLinearOpWithSolveFactory.hpp"
00034 
00035 #include "Thyra_AmesosLinearOpWithSolve.hpp"
00036 #include "Thyra_EpetraOperatorViewExtractorStd.hpp"
00037 #include "Amesos.h"
00038 #include "Teuchos_dyn_cast.hpp"
00039 #include "Teuchos_TimeMonitor.hpp"
00040 
00041 #ifdef HAVE_AMESOS_KLU
00042 #include "Amesos_Klu.h"
00043 #endif
00044 #ifdef HAVE_AMESOS_LAPACK
00045 #include "Amesos_Lapack.h"
00046 #endif
00047 #ifdef HAVE_AMESOS_MUMPS
00048 #include "Amesos_Mumps.h"
00049 #endif
00050 #ifdef HAVE_AMESOS_SCALAPACK
00051 #include "Amesos_Scalapack.h"
00052 #endif
00053 #ifdef HAVE_AMESOS_UMFPACK
00054 #include "Amesos_Umfpack.h"
00055 #endif
00056 #ifdef HAVE_AMESOS_SUPERLUDIST
00057 #include "Amesos_Superludist.h"
00058 #endif
00059 #ifdef HAVE_AMESOS_SUPERLU
00060 #include "Amesos_Superlu.h"
00061 #endif
00062 #ifdef HAVE_AMESOS_DSCPACK
00063 #include "Amesos_Dscpack.h"
00064 #endif
00065 #ifdef HAVE_AMESOS_PARDISO
00066 #include "Amesos_Pardiso.h"
00067 #endif
00068 #ifdef HAVE_AMESOS_TAUCS
00069 #include "Amesos_Taucs.h"
00070 #endif
00071 #ifdef HAVE_AMESOS_PARAKLETE
00072 #include "Amesos_Paraklete.h"
00073 #endif
00074 
00075 namespace {
00076 
00077 Teuchos::RCP<Teuchos::Time> overallTimer, constructTimer, symbolicTimer, factorTimer;
00078 
00079 const std::string epetraFwdOp_str = "epetraFwdOp";
00080 
00081 } // namespace
00082 
00083 namespace Thyra {
00084 
00085 
00086 // Parameter names for Paramter List
00087 
00088 const std::string AmesosLinearOpWithSolveFactory::SolverType_name = "Solver Type";
00089 
00090 const std::string AmesosLinearOpWithSolveFactory::RefactorizationPolicy_name = "Refactorization Policy";
00091 
00092 const std::string AmesosLinearOpWithSolveFactory::ThrowOnPreconditionerInput_name = "Throw on Preconditioner Input";
00093 
00094 const std::string AmesosLinearOpWithSolveFactory::Amesos_Settings_name = "Amesos Settings";
00095 
00096 // Constructors/initializers/accessors
00097 
00098 AmesosLinearOpWithSolveFactory::~AmesosLinearOpWithSolveFactory()
00099 {
00100 #ifdef TEUCHOS_DEBUG
00101   if(paramList_.get())
00102     paramList_->validateParameters(
00103       *this->getValidParameters(),0  // Only validate this level for now!
00104       );
00105 #endif
00106 }
00107 
00108 AmesosLinearOpWithSolveFactory::AmesosLinearOpWithSolveFactory(
00109   const Amesos::ESolverType                            solverType
00110   ,const Amesos::ERefactorizationPolicy                refactorizationPolicy
00111   ,const bool                                          throwOnPrecInput
00112     )
00113   :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
00114   ,solverType_(solverType)
00115   ,refactorizationPolicy_(refactorizationPolicy)
00116   ,throwOnPrecInput_(throwOnPrecInput)
00117 {
00118   initializeTimers();
00119 }
00120 
00121 // Overridden from LinearOpWithSolveFactoryBase
00122 
00123 bool AmesosLinearOpWithSolveFactory::isCompatible(
00124   const LinearOpSourceBase<double> &fwdOpSrc
00125   ) const
00126 {
00127   Teuchos::RCP<const LinearOpBase<double> >
00128     fwdOp = fwdOpSrc.getOp();
00129   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00130   ETransp                                     epetraFwdOpTransp;
00131   EApplyEpetraOpAs                            epetraFwdOpApplyAs;
00132   EAdjointEpetraOp                            epetraFwdOpAdjointSupport;
00133   double                                      epetraFwdOpScalar;
00134   epetraFwdOpViewExtractor_->getEpetraOpView(
00135     fwdOp
00136     ,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
00137     );
00138   if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
00139     return false;
00140   return true;
00141 }
00142 
00143 Teuchos::RCP<LinearOpWithSolveBase<double> >
00144 AmesosLinearOpWithSolveFactory::createOp() const
00145 {
00146   return Teuchos::rcp(new AmesosLinearOpWithSolve());
00147 }
00148 
00149 void AmesosLinearOpWithSolveFactory::initializeOp(
00150   const Teuchos::RCP<const LinearOpSourceBase<double> >    &fwdOpSrc
00151   ,LinearOpWithSolveBase<double>                                   *Op
00152   ,const ESupportSolveUse                                          supportSolveUse
00153   ) const
00154 {
00155   Teuchos::TimeMonitor overallTimeMonitor(*overallTimer);
00156 #ifdef TEUCHOS_DEBUG
00157   TEST_FOR_EXCEPT(Op==NULL);
00158 #endif
00159   const Teuchos::RCP<const LinearOpBase<double> > 
00160     fwdOp = fwdOpSrc->getOp();
00161   //
00162   // Unwrap and get the forward Epetra_Operator object
00163   //
00164   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00165   ETransp                                     epetraFwdOpTransp;
00166   EApplyEpetraOpAs                            epetraFwdOpApplyAs;
00167   EAdjointEpetraOp                            epetraFwdOpAdjointSupport;
00168   double                                      epetraFwdOpScalar;
00169   epetraFwdOpViewExtractor_->getEpetraOpView(
00170     fwdOp,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
00171     );
00172   // Get the AmesosLinearOpWithSolve object
00173   AmesosLinearOpWithSolve
00174     *amesosOp = &Teuchos::dyn_cast<AmesosLinearOpWithSolve>(*Op);
00175   //
00176   // Determine if we must start over or not
00177   //
00178   bool startOver = ( amesosOp->get_amesosSolver()==Teuchos::null );
00179   if(!startOver) {
00180     startOver =
00181       (
00182         epetraFwdOpTransp != amesosOp->get_amesosSolverTransp() ||
00183         epetraFwdOp.get() != amesosOp->get_epetraLP()->GetOperator()
00184         // We must start over if the matrix object changes.  This is a
00185         // weakness of Amesos but there is nothing I can do about this right
00186         // now!
00187         );
00188   }
00189   //
00190   // Update the amesos solver
00191   //
00192   if(startOver) {
00193     //
00194     // This LOWS object has not be initialized yet or is not compatible with the existing
00195     // 
00196     // so this is where we setup everything from the ground up.
00197     //
00198     // Create the linear problem and set the operator with memory of RCP to Epetra_Operator view!
00199     Teuchos::RCP<Epetra_LinearProblem>
00200       epetraLP = Teuchos::rcp(new Epetra_LinearProblem());
00201     epetraLP->SetOperator(const_cast<Epetra_Operator*>(&*epetraFwdOp));
00202     Teuchos::set_extra_data< Teuchos::RCP<const Epetra_Operator> >( epetraFwdOp,
00203       epetraFwdOp_str, Teuchos::inOutArg(epetraLP) );
00204     // Create the concrete solver
00205     Teuchos::RCP<Amesos_BaseSolver>
00206       amesosSolver;
00207     {
00208       Teuchos::TimeMonitor constructTimeMonitor(*constructTimer);
00209       switch(solverType_) {
00210         case Thyra::Amesos::LAPACK :
00211           amesosSolver = Teuchos::rcp(new Amesos_Lapack(*epetraLP));
00212           break;
00213 #ifdef HAVE_AMESOS_KLU
00214         case Thyra::Amesos::KLU :
00215           amesosSolver = Teuchos::rcp(new Amesos_Klu(*epetraLP));
00216           break;
00217 #endif
00218 #ifdef HAVE_AMESOS_MUMPS
00219         case Thyra::Amesos::MUMPS :
00220           amesosSolver = Teuchos::rcp(new Amesos_Mumps(*epetraLP));
00221           break;
00222 #endif
00223 #ifdef HAVE_AMESOS_SCALAPACK
00224         case Thyra::Amesos::SCALAPACK :
00225           amesosSolver = Teuchos::rcp(new Amesos_Scalapack(*epetraLP));
00226           break;
00227 #endif
00228 #ifdef HAVE_AMESOS_UMFPACK
00229         case Thyra::Amesos::UMFPACK :
00230           amesosSolver = Teuchos::rcp(new Amesos_Umfpack(*epetraLP));
00231           break;
00232 #endif
00233 #ifdef HAVE_AMESOS_SUPERLUDIST
00234         case Thyra::Amesos::SUPERLUDIST :
00235           amesosSolver = Teuchos::rcp(new Amesos_Superludist(*epetraLP));
00236           break;
00237 #endif
00238 #ifdef HAVE_AMESOS_SUPERLU
00239         case Thyra::Amesos::SUPERLU :
00240           amesosSolver = Teuchos::rcp(new Amesos_Superlu(*epetraLP));
00241           break;
00242 #endif
00243 #ifdef HAVE_AMESOS_DSCPACK
00244         case Thyra::Amesos::DSCPACK :
00245           amesosSolver = Teuchos::rcp(new Amesos_Dscpack(*epetraLP));
00246           break;
00247 #endif
00248 #ifdef HAVE_AMESOS_PARDISO
00249         case Thyra::Amesos::PARDISO :
00250           amesosSolver = Teuchos::rcp(new Amesos_Pardiso(*epetraLP));
00251           break;
00252 #endif
00253 #ifdef HAVE_AMESOS_TAUCS
00254         case Thyra::Amesos::TAUCS :
00255           amesosSolver = Teuchos::rcp(new Amesos_Taucs(*epetraLP));
00256           break;
00257 #endif
00258 #ifdef HAVE_AMESOS_PARAKLETE
00259         case Thyra::Amesos::PARAKLETE :
00260           amesosSolver = Teuchos::rcp(new Amesos_Paraklete(*epetraLP));
00261           break;
00262 #endif
00263         default:
00264           TEST_FOR_EXCEPTION(
00265             true, std::logic_error
00266             ,"Error, the solver type ID = " << solverType_ << " is invalid!"
00267             );
00268       }
00269     }
00270     // Set the parameters
00271     if(paramList_.get()) amesosSolver->setParameterList(sublist(paramList_,"Amesos Settings"));
00272     // Do the initial factorization
00273     {
00274       Teuchos::TimeMonitor symbolicTimeMonitor(*symbolicTimer);
00275       amesosSolver->SymbolicFactorization();
00276     }
00277     {
00278       Teuchos::TimeMonitor factorTimeMonitor(*factorTimer);
00279       amesosSolver->NumericFactorization();
00280     }
00281     // Initialize the LOWS object and we are done!
00282     amesosOp->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,epetraFwdOpTransp,epetraFwdOpScalar);
00283   }
00284   else {
00285     //
00286     // This LOWS object has already be initialized once so we must just reset
00287     // the matrix and refactor it.
00288     //
00289     // Get non-const pointers to the linear problem and the amesos solver.
00290     // These const-casts are just fine since the amesosOp in non-const.
00291     Teuchos::RCP<Epetra_LinearProblem>
00292       epetraLP = Teuchos::rcp_const_cast<Epetra_LinearProblem>(amesosOp->get_epetraLP());
00293     Teuchos::RCP<Amesos_BaseSolver>
00294       amesosSolver = amesosOp->get_amesosSolver();
00295     // Reset the forward operator with memory of RCP to Epetra_Operator view!
00296     epetraLP->SetOperator(const_cast<Epetra_Operator*>(&*epetraFwdOp));
00297     Teuchos::get_extra_data< Teuchos::RCP<const Epetra_Operator> >(epetraLP,epetraFwdOp_str) = epetraFwdOp;
00298     // Reset the parameters
00299     if(paramList_.get()) amesosSolver->setParameterList(sublist(paramList_,Amesos_Settings_name));
00300     // Repivot if asked
00301     if(refactorizationPolicy_==Amesos::REPIVOT_ON_REFACTORIZATION) {
00302       Teuchos::TimeMonitor symbolicTimeMonitor(*symbolicTimer);
00303       amesosSolver->SymbolicFactorization();
00304     }
00305     {
00306       Teuchos::TimeMonitor factorTimeMonitor(*factorTimer);
00307       amesosSolver->NumericFactorization();
00308     }
00309     // Reinitialize the LOWS object and we are done! (we must do this to get the
00310     // possibly new transpose and scaling factors back in)
00311     amesosOp->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,epetraFwdOpTransp,epetraFwdOpScalar);
00312   }
00313   amesosOp->setOStream(this->getOStream());
00314   amesosOp->setVerbLevel(this->getVerbLevel());
00315 }
00316 
00317 bool AmesosLinearOpWithSolveFactory::supportsPreconditionerInputType(const EPreconditionerInputType precOpType) const
00318 {
00319   return false;
00320 }
00321 
00322 void AmesosLinearOpWithSolveFactory::initializePreconditionedOp(
00323   const Teuchos::RCP<const LinearOpSourceBase<double> >       &fwdOpSrc
00324   ,const Teuchos::RCP<const PreconditionerBase<double> >      &prec
00325   ,LinearOpWithSolveBase<double>                                      *Op
00326   ,const ESupportSolveUse                                             supportSolveUse
00327   ) const
00328 {
00329   TEST_FOR_EXCEPTION(
00330     this->throwOnPrecInput_, std::logic_error
00331     ,"Error, the concrete implementation described as \'"<<this->description()<<"\' does not support precondtioners "
00332     "and has been configured to throw this exception when the  initializePreconditionedOp(...) function is called!"
00333     );
00334   this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the precondtioner!
00335 }
00336 
00337 void AmesosLinearOpWithSolveFactory::initializePreconditionedOp(
00338   const Teuchos::RCP<const LinearOpSourceBase<double> >       &fwdOpSrc
00339   ,const Teuchos::RCP<const LinearOpSourceBase<double> >      &approxFwdOpSrc
00340   ,LinearOpWithSolveBase<double>                                      *Op
00341   ,const ESupportSolveUse                                             supportSolveUse
00342   ) const
00343 {
00344   TEST_FOR_EXCEPTION(
00345     this->throwOnPrecInput_, std::logic_error
00346     ,"Error, the concrete implementation described as \'"<<this->description()<<"\' does not support precondtioners "
00347     "and has been configured to throw this exception when the  initializePreconditionedOp(...) function is called!"
00348     );
00349   this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the precondtioner!
00350 }
00351 
00352 void AmesosLinearOpWithSolveFactory::uninitializeOp(
00353   LinearOpWithSolveBase<double>                               *Op
00354   ,Teuchos::RCP<const LinearOpSourceBase<double> >    *fwdOpSrc
00355   ,Teuchos::RCP<const PreconditionerBase<double> >    *prec
00356   ,Teuchos::RCP<const LinearOpSourceBase<double> >    *approxFwdOpSrc
00357   ,ESupportSolveUse                                           *supportSolveUse
00358   ) const
00359 {
00360 #ifdef TEUCHOS_DEBUG
00361   TEST_FOR_EXCEPT(Op==NULL);
00362 #endif
00363   AmesosLinearOpWithSolve
00364     *amesosOp = &Teuchos::dyn_cast<AmesosLinearOpWithSolve>(*Op);
00365   Teuchos::RCP<const LinearOpSourceBase<double> >
00366     _fwdOpSrc = amesosOp->extract_fwdOpSrc(); // Will be null if uninitialized!
00367   if(_fwdOpSrc.get()) {
00368     // Erase the Epetra_Operator view of the forward operator!
00369     Teuchos::RCP<Epetra_LinearProblem> epetraLP = amesosOp->get_epetraLP();
00370     Teuchos::get_extra_data< Teuchos::RCP<const Epetra_Operator> >(
00371       epetraLP,epetraFwdOp_str
00372       )
00373       = Teuchos::null;
00374     // Note, we did not erase the address of the operator in
00375     // epetraLP->GetOperator() since it seems that the amesos solvers do not
00376     // recheck the value of GetProblem()->GetOperator() so you had better not
00377     // rest this!
00378   }
00379   if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc; // It is fine if the client does not want this object back!
00380   if(prec) *prec = Teuchos::null; // We never keep a preconditioner!
00381   if(approxFwdOpSrc) *approxFwdOpSrc = Teuchos::null; // We never keep an approximate fwd operator!
00382 }
00383 
00384 // Overridden from ParameterListAcceptor
00385 
00386 void AmesosLinearOpWithSolveFactory::setParameterList(
00387   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00388   )
00389 {
00390   TEST_FOR_EXCEPT(paramList.get()==NULL);
00391   paramList->validateParameters(*this->getValidParameters(),0); // Only validate this level for now!
00392   paramList_ = paramList;
00393   solverType_ =
00394     Amesos::solverTypeNameToEnumMap.get<Amesos::ESolverType>(
00395       paramList_->get(
00396         SolverType_name
00397         ,Amesos::toString(solverType_)
00398         )
00399       ,paramList_->name()+"->"+SolverType_name
00400       );
00401   refactorizationPolicy_ = 
00402     Amesos::refactorizationPolicyNameToEnumMap.get<Amesos::ERefactorizationPolicy>(
00403       paramList_->get(
00404         RefactorizationPolicy_name
00405         ,Amesos::toString(refactorizationPolicy_)
00406         )
00407       ,paramList_->name()+"->"+RefactorizationPolicy_name
00408       );
00409   throwOnPrecInput_ = paramList_->get(ThrowOnPreconditionerInput_name,throwOnPrecInput_);
00410 }
00411 
00412 Teuchos::RCP<Teuchos::ParameterList>
00413 AmesosLinearOpWithSolveFactory::getNonconstParameterList()
00414 {
00415   return paramList_;
00416 }
00417 
00418 Teuchos::RCP<Teuchos::ParameterList>
00419 AmesosLinearOpWithSolveFactory::unsetParameterList()
00420 {
00421   Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
00422   paramList_ = Teuchos::null;
00423   return _paramList;
00424 }
00425 
00426 Teuchos::RCP<const Teuchos::ParameterList>
00427 AmesosLinearOpWithSolveFactory::getParameterList() const
00428 {
00429   return paramList_;
00430 }
00431 
00432 Teuchos::RCP<const Teuchos::ParameterList>
00433 AmesosLinearOpWithSolveFactory::getValidParameters() const
00434 {
00435   return generateAndGetValidParameters();
00436 }
00437 
00438 // Public functions overridden from Teuchos::Describable
00439 
00440 std::string AmesosLinearOpWithSolveFactory::description() const
00441 {
00442   std::ostringstream oss;
00443   oss << "Thyra::AmesosLinearOpWithSolveFactory{";
00444   oss << "solverType=" << toString(solverType_);
00445   oss << "}";
00446   return oss.str();
00447 }
00448 
00449 // private
00450 
00451 
00452 void AmesosLinearOpWithSolveFactory::initializeTimers()
00453 {
00454   if(!overallTimer.get()) {
00455     overallTimer    = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF");
00456     constructTimer  = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF:InitConstruct");
00457     symbolicTimer   = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF:Symbolic");
00458     factorTimer     = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF:Factor");
00459   }
00460 }
00461 
00462 Teuchos::RCP<const Teuchos::ParameterList>
00463 AmesosLinearOpWithSolveFactory::generateAndGetValidParameters()
00464 {
00465   static Teuchos::RCP<Teuchos::ParameterList> validParamList;
00466   if(validParamList.get()==NULL) {
00467     validParamList = Teuchos::rcp(new Teuchos::ParameterList("Amesos"));
00468     validParamList->set(
00469       SolverType_name
00470 #ifdef HAVE_AMESOS_KLU
00471       ,Amesos::toString(Amesos::KLU)
00472 #else
00473       ,Amesos::toString(Amesos::LAPACK)
00474 #endif
00475       );
00476     validParamList->set(RefactorizationPolicy_name,Amesos::toString(Amesos::REPIVOT_ON_REFACTORIZATION));
00477     validParamList->set(ThrowOnPreconditionerInput_name,bool(true));
00478     validParamList->sublist(Amesos_Settings_name).setParameters(::Amesos::GetValidParameters());
00479   }
00480   return validParamList;
00481 }
00482 
00483 } // namespace Thyra
00484 
00485 #endif // __sun
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines