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