Thyra_AmesosLinearOpWithSolveFactory.cpp

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_CXX
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 #include "Teuchos_TypeNameTraits.hpp"
00041 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00042 
00043 #ifdef HAVE_AMESOS_KLU
00044 #include "Amesos_Klu.h"
00045 #endif
00046 #ifdef HAVE_AMESOS_PASTIX
00047 #include "Amesos_Pastix.h"
00048 #endif
00049 #ifdef HAVE_AMESOS_LAPACK
00050 #include "Amesos_Lapack.h"
00051 #endif
00052 #ifdef HAVE_AMESOS_MUMPS
00053 #include "Amesos_Mumps.h"
00054 #endif
00055 #ifdef HAVE_AMESOS_SCALAPACK
00056 #include "Amesos_Scalapack.h"
00057 #endif
00058 #ifdef HAVE_AMESOS_UMFPACK
00059 #include "Amesos_Umfpack.h"
00060 #endif
00061 #ifdef HAVE_AMESOS_SUPERLUDIST
00062 #include "Amesos_Superludist.h"
00063 #endif
00064 #ifdef HAVE_AMESOS_SUPERLU
00065 #include "Amesos_Superlu.h"
00066 #endif
00067 #ifdef HAVE_AMESOS_DSCPACK
00068 #include "Amesos_Dscpack.h"
00069 #endif
00070 #ifdef HAVE_AMESOS_PARDISO
00071 #include "Amesos_Pardiso.h"
00072 #endif
00073 #ifdef HAVE_AMESOS_TAUCS
00074 #include "Amesos_Taucs.h"
00075 #endif
00076 #ifdef HAVE_AMESOS_PARAKLETE
00077 #include "Amesos_Paraklete.h"
00078 #endif
00079 
00080 namespace {
00081 
00082 const std::string epetraFwdOp_str = "epetraFwdOp";
00083 
00084 } // namespace
00085 
00086 namespace Thyra {
00087 
00088 
00089 // Parameter names for Paramter List
00090 
00091 const std::string AmesosLinearOpWithSolveFactory::SolverType_name = "Solver Type";
00092 
00093 const std::string AmesosLinearOpWithSolveFactory::RefactorizationPolicy_name = "Refactorization Policy";
00094 
00095 const std::string AmesosLinearOpWithSolveFactory::ThrowOnPreconditionerInput_name = "Throw on Preconditioner Input";
00096 
00097 const std::string AmesosLinearOpWithSolveFactory::Amesos_Settings_name = "Amesos Settings";
00098 
00099 // Constructors/initializers/accessors
00100 
00101 AmesosLinearOpWithSolveFactory::~AmesosLinearOpWithSolveFactory()
00102 {
00103 #ifdef TEUCHOS_DEBUG
00104   if(paramList_.get())
00105     paramList_->validateParameters(
00106       *this->getValidParameters(),0  // Only validate this level for now!
00107       );
00108 #endif
00109 }
00110 
00111 AmesosLinearOpWithSolveFactory::AmesosLinearOpWithSolveFactory(
00112   const Amesos::ESolverType                            solverType
00113   ,const Amesos::ERefactorizationPolicy                refactorizationPolicy
00114   ,const bool                                          throwOnPrecInput
00115     )
00116   :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
00117   ,solverType_(solverType)
00118   ,refactorizationPolicy_(refactorizationPolicy)
00119   ,throwOnPrecInput_(throwOnPrecInput)
00120 {}
00121 
00122 // Overridden from LinearOpWithSolveFactoryBase
00123 
00124 bool AmesosLinearOpWithSolveFactory::isCompatible(
00125   const LinearOpSourceBase<double> &fwdOpSrc
00126   ) const
00127 {
00128   RCP<const LinearOpBase<double> >
00129     fwdOp = fwdOpSrc.getOp();
00130   RCP<const Epetra_Operator> epetraFwdOp;
00131   EOpTransp                                     epetraFwdOpTransp;
00132   EApplyEpetraOpAs                            epetraFwdOpApplyAs;
00133   EAdjointEpetraOp                            epetraFwdOpAdjointSupport;
00134   double                                      epetraFwdOpScalar;
00135   epetraFwdOpViewExtractor_->getEpetraOpView(
00136     fwdOp
00137     ,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
00138     );
00139   if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
00140     return false;
00141   return true;
00142 }
00143 
00144 RCP<LinearOpWithSolveBase<double> >
00145 AmesosLinearOpWithSolveFactory::createOp() const
00146 {
00147   return Teuchos::rcp(new AmesosLinearOpWithSolve());
00148 }
00149 
00150 void AmesosLinearOpWithSolveFactory::initializeOp(
00151   const RCP<const LinearOpSourceBase<double> >    &fwdOpSrc
00152   ,LinearOpWithSolveBase<double>                                   *Op
00153   ,const ESupportSolveUse                                          supportSolveUse
00154   ) const
00155 {
00156   TEUCHOS_FUNC_TIME_MONITOR("AmesosLOWSF");
00157 #ifdef TEUCHOS_DEBUG
00158   TEST_FOR_EXCEPT(Op==NULL);
00159 #endif
00160   const RCP<const LinearOpBase<double> > 
00161     fwdOp = fwdOpSrc->getOp();
00162   //
00163   // Unwrap and get the forward Epetra_Operator object
00164   //
00165   RCP<const Epetra_Operator> epetraFwdOp;
00166   EOpTransp                                     epetraFwdOpTransp;
00167   EApplyEpetraOpAs                            epetraFwdOpApplyAs;
00168   EAdjointEpetraOp                            epetraFwdOpAdjointSupport;
00169   double                                      epetraFwdOpScalar;
00170   epetraFwdOpViewExtractor_->getEpetraOpView(
00171     fwdOp,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
00172     );
00173   // Get the AmesosLinearOpWithSolve object
00174   AmesosLinearOpWithSolve
00175     *amesosOp = &Teuchos::dyn_cast<AmesosLinearOpWithSolve>(*Op);
00176   //
00177   // Determine if we must start over or not
00178   //
00179   bool startOver = ( amesosOp->get_amesosSolver()==Teuchos::null );
00180   if(!startOver) {
00181     startOver =
00182       (
00183         epetraFwdOpTransp != amesosOp->get_amesosSolverTransp() ||
00184         epetraFwdOp.get() != amesosOp->get_epetraLP()->GetOperator()
00185         // We must start over if the matrix object changes.  This is a
00186         // weakness of Amesos but there is nothing I can do about this right
00187         // now!
00188         );
00189   }
00190   //
00191   // Update the amesos solver
00192   //
00193   if(startOver) {
00194     //
00195     // This LOWS object has not be initialized yet or is not compatible with the existing
00196     // 
00197     // so this is where we setup everything from the ground up.
00198     //
00199     // Create the linear problem and set the operator with memory of RCP to Epetra_Operator view!
00200     RCP<Epetra_LinearProblem>
00201       epetraLP = Teuchos::rcp(new Epetra_LinearProblem());
00202     epetraLP->SetOperator(const_cast<Epetra_Operator*>(&*epetraFwdOp));
00203     Teuchos::set_extra_data< RCP<const Epetra_Operator> >( epetraFwdOp, epetraFwdOp_str, &epetraLP );
00204     // Create the concrete solver
00205     RCP<Amesos_BaseSolver>
00206       amesosSolver;
00207     {
00208       TEUCHOS_FUNC_TIME_MONITOR("AmesosLOWSF:InitConstruct");
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_PASTIX
00219         case Thyra::Amesos::PASTIX :
00220           amesosSolver = Teuchos::rcp(new Amesos_Pastix(*epetraLP));
00221           break;
00222 #endif
00223 #ifdef HAVE_AMESOS_MUMPS
00224         case Thyra::Amesos::MUMPS :
00225           amesosSolver = Teuchos::rcp(new Amesos_Mumps(*epetraLP));
00226           break;
00227 #endif
00228 #ifdef HAVE_AMESOS_SCALAPACK
00229         case Thyra::Amesos::SCALAPACK :
00230           amesosSolver = Teuchos::rcp(new Amesos_Scalapack(*epetraLP));
00231           break;
00232 #endif
00233 #ifdef HAVE_AMESOS_UMFPACK
00234         case Thyra::Amesos::UMFPACK :
00235           amesosSolver = Teuchos::rcp(new Amesos_Umfpack(*epetraLP));
00236           break;
00237 #endif
00238 #ifdef HAVE_AMESOS_SUPERLUDIST
00239         case Thyra::Amesos::SUPERLUDIST :
00240           amesosSolver = Teuchos::rcp(new Amesos_Superludist(*epetraLP));
00241           break;
00242 #endif
00243 #ifdef HAVE_AMESOS_SUPERLU
00244         case Thyra::Amesos::SUPERLU :
00245           amesosSolver = Teuchos::rcp(new Amesos_Superlu(*epetraLP));
00246           break;
00247 #endif
00248 #ifdef HAVE_AMESOS_DSCPACK
00249         case Thyra::Amesos::DSCPACK :
00250           amesosSolver = Teuchos::rcp(new Amesos_Dscpack(*epetraLP));
00251           break;
00252 #endif
00253 #ifdef HAVE_AMESOS_PARDISO
00254         case Thyra::Amesos::PARDISO :
00255           amesosSolver = Teuchos::rcp(new Amesos_Pardiso(*epetraLP));
00256           break;
00257 #endif
00258 #ifdef HAVE_AMESOS_TAUCS
00259         case Thyra::Amesos::TAUCS :
00260           amesosSolver = Teuchos::rcp(new Amesos_Taucs(*epetraLP));
00261           break;
00262 #endif
00263 #ifdef HAVE_AMESOS_PARAKLETE
00264         case Thyra::Amesos::PARAKLETE :
00265           amesosSolver = Teuchos::rcp(new Amesos_Paraklete(*epetraLP));
00266           break;
00267 #endif
00268         default:
00269           TEST_FOR_EXCEPTION(
00270             true, std::logic_error
00271             ,"Error, the solver type ID = " << solverType_ << " is invalid!"
00272             );
00273       }
00274     }
00275     // Set the parameters
00276     if(paramList_.get()) amesosSolver->SetParameters(paramList_->sublist("Amesos Settings"));
00277     // Do the initial factorization
00278     {
00279       TEUCHOS_FUNC_TIME_MONITOR("AmesosLOWSF:Symbolic");
00280       const int err = amesosSolver->SymbolicFactorization();
00281       TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
00282         "Error, SymbolicFactorization() on amesos solver of type \'"<<Teuchos::typeName(*amesosSolver)<<"\'\n"
00283         "returned error code "<<err<<"!" );
00284     }
00285     {
00286       TEUCHOS_FUNC_TIME_MONITOR("AmesosLOWSF:Factor");
00287       const int err = amesosSolver->NumericFactorization();
00288       TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
00289         "Error, NumericFactorization() on amesos solver of type \'"<<Teuchos::typeName(*amesosSolver)<<"\'\n"
00290         "returned error code "<<err<<"!" );
00291     }
00292     // Initialize the LOWS object and we are done!
00293     amesosOp->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,epetraFwdOpTransp,epetraFwdOpScalar);
00294   }
00295   else {
00296     //
00297     // This LOWS object has already be initialized once so we must just reset
00298     // the matrix and refactor it.
00299     //
00300     // Get non-const pointers to the linear problem and the amesos solver.
00301     // These const-casts are just fine since the amesosOp in non-const.
00302     RCP<Epetra_LinearProblem>
00303       epetraLP = Teuchos::rcp_const_cast<Epetra_LinearProblem>(amesosOp->get_epetraLP());
00304     RCP<Amesos_BaseSolver>
00305       amesosSolver = amesosOp->get_amesosSolver();
00306     // Reset the forward operator with memory of RCP to Epetra_Operator view!
00307     epetraLP->SetOperator(const_cast<Epetra_Operator*>(&*epetraFwdOp));
00308     Teuchos::get_nonconst_extra_data<RCP<const Epetra_Operator> >(epetraLP,epetraFwdOp_str) = epetraFwdOp;
00309     // Reset the parameters
00310     if(paramList_.get()) amesosSolver->SetParameters(paramList_->sublist(Amesos_Settings_name));
00311     // Repivot if asked
00312     if(refactorizationPolicy_==Amesos::REPIVOT_ON_REFACTORIZATION) {
00313       TEUCHOS_FUNC_TIME_MONITOR("AmesosLOWSF:Symbolic");
00314       const int err = amesosSolver->SymbolicFactorization();
00315       TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
00316         "Error, SymbolicFactorization() on amesos solver of type \'"<<Teuchos::typeName(*amesosSolver)<<"\'\n"
00317         "returned error code "<<err<<"!" );
00318     }
00319     {
00320       TEUCHOS_FUNC_TIME_MONITOR("AmesosLOWSF::Factor");
00321       const int err = amesosSolver->NumericFactorization();
00322       TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
00323         "Error, NumericFactorization() on amesos solver of type \'"<<Teuchos::typeName(*amesosSolver)<<"\'\n"
00324         "returned error code "<<err<<"!" );
00325     }
00326     /* ToDo: Put this back in once PrintStatus accepts an std::ostream!
00327     OsTab tab(out);
00328     amesosSolver->PrintStatus()
00329     */
00330     // Reinitialize the LOWS object and we are done! (we must do this to get the
00331     // possibly new transpose and scaling factors back in)
00332     amesosOp->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,epetraFwdOpTransp,epetraFwdOpScalar);
00333   }
00334   amesosOp->setOStream(this->getOStream());
00335   amesosOp->setVerbLevel(this->getVerbLevel());
00336 }
00337 
00338 bool AmesosLinearOpWithSolveFactory::supportsPreconditionerInputType(const EPreconditionerInputType precOpType) const
00339 {
00340   return false;
00341 }
00342 
00343 void AmesosLinearOpWithSolveFactory::initializePreconditionedOp(
00344   const RCP<const LinearOpSourceBase<double> >       &fwdOpSrc
00345   ,const RCP<const PreconditionerBase<double> >      &prec
00346   ,LinearOpWithSolveBase<double>                                      *Op
00347   ,const ESupportSolveUse                                             supportSolveUse
00348   ) const
00349 {
00350   TEST_FOR_EXCEPTION(
00351     this->throwOnPrecInput_, std::logic_error
00352     ,"Error, the concrete implementation described as \'"<<this->description()<<"\' does not support precondtioners "
00353     "and has been configured to throw this exception when the  initializePreconditionedOp(...) function is called!"
00354     );
00355   this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the precondtioner!
00356 }
00357 
00358 void AmesosLinearOpWithSolveFactory::initializePreconditionedOp(
00359   const RCP<const LinearOpSourceBase<double> >       &fwdOpSrc
00360   ,const RCP<const LinearOpSourceBase<double> >      &approxFwdOpSrc
00361   ,LinearOpWithSolveBase<double>                                      *Op
00362   ,const ESupportSolveUse                                             supportSolveUse
00363   ) const
00364 {
00365   TEST_FOR_EXCEPTION(
00366     this->throwOnPrecInput_, std::logic_error
00367     ,"Error, the concrete implementation described as \'"<<this->description()<<"\' does not support precondtioners "
00368     "and has been configured to throw this exception when the  initializePreconditionedOp(...) function is called!"
00369     );
00370   this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the precondtioner!
00371 }
00372 
00373 void AmesosLinearOpWithSolveFactory::uninitializeOp(
00374   LinearOpWithSolveBase<double>                               *Op
00375   ,RCP<const LinearOpSourceBase<double> >    *fwdOpSrc
00376   ,RCP<const PreconditionerBase<double> >    *prec
00377   ,RCP<const LinearOpSourceBase<double> >    *approxFwdOpSrc
00378   ,ESupportSolveUse                                           *supportSolveUse
00379   ) const
00380 {
00381 #ifdef TEUCHOS_DEBUG
00382   TEST_FOR_EXCEPT(Op==NULL);
00383 #endif
00384   AmesosLinearOpWithSolve
00385     *amesosOp = &Teuchos::dyn_cast<AmesosLinearOpWithSolve>(*Op);
00386   RCP<const LinearOpSourceBase<double> >
00387     _fwdOpSrc = amesosOp->extract_fwdOpSrc(); // Will be null if uninitialized!
00388   if(_fwdOpSrc.get()) {
00389     // Erase the Epetra_Operator view of the forward operator!
00390     RCP<Epetra_LinearProblem> epetraLP = amesosOp->get_epetraLP();
00391     Teuchos::get_nonconst_extra_data< RCP<const Epetra_Operator> >(
00392       epetraLP,epetraFwdOp_str
00393       )
00394       = Teuchos::null;
00395     // Note, we did not erase the address of the operator in
00396     // epetraLP->GetOperator() since it seems that the amesos solvers do not
00397     // recheck the value of GetProblem()->GetOperator() so you had better not
00398     // rest this!
00399   }
00400   if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc; // It is fine if the client does not want this object back!
00401   if(prec) *prec = Teuchos::null; // We never keep a preconditioner!
00402   if(approxFwdOpSrc) *approxFwdOpSrc = Teuchos::null; // We never keep an approximate fwd operator!
00403 }
00404 
00405 // Overridden from ParameterListAcceptor
00406 
00407 void AmesosLinearOpWithSolveFactory::setParameterList(
00408   RCP<Teuchos::ParameterList> const& paramList
00409   )
00410 {
00411   TEST_FOR_EXCEPT(paramList.get()==NULL);
00412   paramList->validateParameters(*this->getValidParameters(),0); // Only validate this level for now!
00413   paramList_ = paramList;
00414   solverType_ =
00415     Amesos::solverTypeNameToEnumMap.get<Amesos::ESolverType>(
00416       paramList_->get(
00417         SolverType_name
00418         ,Amesos::toString(solverType_)
00419         )
00420       ,paramList_->name()+"->"+SolverType_name
00421       );
00422   refactorizationPolicy_ = 
00423     Amesos::refactorizationPolicyNameToEnumMap.get<Amesos::ERefactorizationPolicy>(
00424       paramList_->get(
00425         RefactorizationPolicy_name
00426         ,Amesos::toString(refactorizationPolicy_)
00427         )
00428       ,paramList_->name()+"->"+RefactorizationPolicy_name
00429       );
00430   throwOnPrecInput_ = paramList_->get(ThrowOnPreconditionerInput_name,throwOnPrecInput_);
00431   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00432 }
00433 
00434 RCP<Teuchos::ParameterList>
00435 AmesosLinearOpWithSolveFactory::getNonconstParameterList()
00436 {
00437   return paramList_;
00438 }
00439 
00440 RCP<Teuchos::ParameterList>
00441 AmesosLinearOpWithSolveFactory::unsetParameterList()
00442 {
00443   RCP<Teuchos::ParameterList> _paramList = paramList_;
00444   paramList_ = Teuchos::null;
00445   return _paramList;
00446 }
00447 
00448 RCP<const Teuchos::ParameterList>
00449 AmesosLinearOpWithSolveFactory::getParameterList() const
00450 {
00451   return paramList_;
00452 }
00453 
00454 RCP<const Teuchos::ParameterList>
00455 AmesosLinearOpWithSolveFactory::getValidParameters() const
00456 {
00457   return generateAndGetValidParameters();
00458 }
00459 
00460 // Public functions overridden from Teuchos::Describable
00461 
00462 std::string AmesosLinearOpWithSolveFactory::description() const
00463 {
00464   std::ostringstream oss;
00465   oss << "Thyra::AmesosLinearOpWithSolveFactory{";
00466   oss << "solverType=" << toString(solverType_);
00467   oss << "}";
00468   return oss.str();
00469 }
00470 
00471 // private
00472 
00473 RCP<const Teuchos::ParameterList>
00474 AmesosLinearOpWithSolveFactory::generateAndGetValidParameters()
00475 {
00476   static RCP<Teuchos::ParameterList> validParamList;
00477   if(validParamList.get()==NULL) {
00478     validParamList = Teuchos::rcp(new Teuchos::ParameterList("Amesos"));
00479     validParamList->set(
00480       SolverType_name
00481 #ifdef HAVE_AMESOS_KLU
00482       ,Amesos::toString(Amesos::KLU)
00483 #else
00484       ,Amesos::toString(Amesos::LAPACK)
00485 #endif
00486       );
00487     validParamList->set(RefactorizationPolicy_name,Amesos::toString(Amesos::REPIVOT_ON_REFACTORIZATION));
00488     validParamList->set(ThrowOnPreconditionerInput_name,bool(true));
00489     validParamList->sublist(Amesos_Settings_name).setParameters(::Amesos::GetValidParameters());
00490     Teuchos::setupVerboseObjectSublist(&*validParamList);
00491   }
00492   return validParamList;
00493 }
00494 
00495 } // namespace Thyra
00496 
00497 #endif // SUN_CXX

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