Stratimikos Version of the Day
Thyra_AmesosLinearOpWithSolve.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 #include "Thyra_AmesosLinearOpWithSolve.hpp"
00032 #include "Thyra_EpetraThyraWrappers.hpp"
00033 #include "Thyra_MultiVectorStdOps.hpp"
00034 #include "Epetra_MultiVector.h"
00035 #include "Teuchos_TimeMonitor.hpp"
00036 
00037 
00038 namespace Thyra {
00039 
00040 
00041 // Constructors/initializers/accessors
00042 
00043 
00044 AmesosLinearOpWithSolve::AmesosLinearOpWithSolve()
00045 {}
00046 
00047 
00048 AmesosLinearOpWithSolve::AmesosLinearOpWithSolve(
00049   const Teuchos::RCP<const LinearOpBase<double> > &fwdOp,
00050   const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
00051   const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
00052   const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
00053   const EOpTransp amesosSolverTransp,
00054   const double amesosSolverScalar
00055   )
00056 {
00057   this->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,
00058     amesosSolverTransp,amesosSolverScalar);
00059 }
00060 
00061 
00062 void AmesosLinearOpWithSolve::initialize(
00063   const Teuchos::RCP<const LinearOpBase<double> > &fwdOp,
00064   const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
00065   const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
00066   const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
00067   const EOpTransp amesosSolverTransp,
00068   const double amesosSolverScalar
00069   )
00070 {
00071 #ifdef TEUCHOS_DEBUG
00072   TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00073   TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00074   TEUCHOS_TEST_FOR_EXCEPT(epetraLP.get()==NULL);
00075   TEUCHOS_TEST_FOR_EXCEPT(amesosSolver.get()==NULL);
00076   TEUCHOS_TEST_FOR_EXCEPT(epetraLP->GetLHS()!=NULL);
00077   TEUCHOS_TEST_FOR_EXCEPT(epetraLP->GetRHS()!=NULL);
00078 #endif
00079   fwdOp_ = fwdOp;
00080   fwdOpSrc_ = fwdOpSrc;
00081   epetraLP_ = epetraLP;
00082   amesosSolver_ = amesosSolver;
00083   amesosSolverTransp_ = amesosSolverTransp;
00084   amesosSolverScalar_ = amesosSolverScalar;
00085   const std::string fwdOpLabel = fwdOp_->getObjectLabel();
00086   if(fwdOpLabel.length())
00087     this->setObjectLabel( "lows("+fwdOpLabel+")" );
00088 }
00089 
00090 
00091 Teuchos::RCP<const LinearOpSourceBase<double> >
00092 AmesosLinearOpWithSolve::extract_fwdOpSrc()
00093 {
00094   Teuchos::RCP<const LinearOpSourceBase<double> >
00095     _fwdOpSrc = fwdOpSrc_;
00096   fwdOpSrc_ = Teuchos::null;
00097   return _fwdOpSrc;
00098 }
00099 
00100 
00101 void AmesosLinearOpWithSolve::uninitialize(
00102   Teuchos::RCP<const LinearOpBase<double> > *fwdOp,
00103   Teuchos::RCP<const LinearOpSourceBase<double> > *fwdOpSrc,
00104   Teuchos::RCP<Epetra_LinearProblem> *epetraLP,
00105   Teuchos::RCP<Amesos_BaseSolver> *amesosSolver,
00106   EOpTransp *amesosSolverTransp,
00107   double *amesosSolverScalar
00108   )
00109 {
00110 
00111   if(fwdOp) *fwdOp = fwdOp_;
00112   if(fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
00113   if(epetraLP) *epetraLP = epetraLP_;
00114   if(amesosSolver) *amesosSolver = amesosSolver_;
00115   if(amesosSolverTransp) *amesosSolverTransp = amesosSolverTransp_;
00116   if(amesosSolverScalar) *amesosSolverScalar = amesosSolverScalar_;
00117 
00118   fwdOp_ = Teuchos::null;
00119   fwdOpSrc_ = Teuchos::null;
00120   epetraLP_ = Teuchos::null;
00121   amesosSolver_ = Teuchos::null;
00122   amesosSolverTransp_ = NOTRANS;
00123   amesosSolverScalar_ = 0.0;
00124 
00125 }
00126 
00127 
00128 // Overridden from LinearOpBase
00129 
00130 
00131 Teuchos::RCP< const VectorSpaceBase<double> >
00132 AmesosLinearOpWithSolve::range() const
00133 {
00134   return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
00135 }
00136 
00137 
00138 Teuchos::RCP< const VectorSpaceBase<double> >
00139 AmesosLinearOpWithSolve::domain() const
00140 {
00141   return  ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
00142 }
00143 
00144 
00145 Teuchos::RCP<const LinearOpBase<double> >
00146 AmesosLinearOpWithSolve::clone() const
00147 {
00148   return Teuchos::null; // Not supported yet but could be
00149 }
00150 
00151 
00152 // Overridden from Teuchos::Describable
00153 
00154 
00155 std::string AmesosLinearOpWithSolve::description() const
00156 {
00157   std::ostringstream oss;
00158   oss << Teuchos::Describable::description();
00159   if(!is_null(amesosSolver_)) {
00160     oss
00161       << "{fwdOp="<<fwdOp_->description()
00162       << ",amesosSolver="<<typeName(*amesosSolver_)<<"}";
00163   }
00164   return oss.str();
00165 }
00166 
00167 
00168 void AmesosLinearOpWithSolve::describe(
00169   Teuchos::FancyOStream &out,
00170   const Teuchos::EVerbosityLevel verbLevel
00171   ) const
00172 {
00173   using Teuchos::OSTab;
00174   using Teuchos::typeName;
00175   using Teuchos::describe;
00176   switch(verbLevel) {
00177     case Teuchos::VERB_DEFAULT:
00178     case Teuchos::VERB_LOW:
00179       out << this->description() << std::endl;
00180       break;
00181     case Teuchos::VERB_MEDIUM:
00182     case Teuchos::VERB_HIGH:
00183     case Teuchos::VERB_EXTREME:
00184     {
00185       out
00186         << Teuchos::Describable::description() << "{"
00187         << "rangeDim=" << this->range()->dim()
00188         << ",domainDim="<< this->domain()->dim() << "}\n";
00189       OSTab tab(out);
00190       if(!is_null(fwdOp_)) {
00191         out << "fwdOp = " << describe(*fwdOp_,verbLevel);
00192       }
00193       if(!is_null(amesosSolver_)) {
00194         out << "amesosSolver=" << typeName(*amesosSolver_) << "\n";
00195       }
00196       break;
00197     }
00198     default:
00199       TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
00200   }
00201 }
00202 
00203 
00204 // protected
00205 
00206 
00207 // Overridden from LinearOpBase
00208 
00209 
00210 bool AmesosLinearOpWithSolve::opSupportedImpl(EOpTransp M_trans) const
00211 {
00212   return ::Thyra::opSupported(*fwdOp_,M_trans);
00213 }
00214 
00215 
00216 void AmesosLinearOpWithSolve::applyImpl(
00217   const EOpTransp M_trans,
00218   const MultiVectorBase<double> &X,
00219   const Ptr<MultiVectorBase<double> > &Y,
00220   const double alpha,
00221   const double beta
00222   ) const
00223 {
00224   Thyra::apply( *fwdOp_, M_trans, X, Y, alpha, beta );
00225 }
00226 
00227 
00228 // Overridden from LinearOpWithSolveBase
00229 
00230 
00231 bool AmesosLinearOpWithSolve::solveSupportsImpl(EOpTransp M_trans) const
00232 {
00233   if (Thyra::real_trans(M_trans) == Thyra::NOTRANS) {
00234     // Assume every amesos solver supports a basic forward solve!
00235     return true;
00236   }
00237   // Query the amesos solver to see if it supports the transpose operation.
00238   // NOTE: Amesos_BaseSolver makes you change the state of the object to
00239   // determine if the object supports an adjoint solver.  This is a bad design
00240   // but I have no control over that.  This is why you see this hacked
00241   // oldUseTranspose varible and logic.  NOTE: This function meets the basic
00242   // guarantee but if setUseTransplse(...) throws, then the state of
00243   // UseTranspose() may be different.
00244   const bool oldUseTranspose = amesosSolver_->UseTranspose();
00245   const bool supportsAdjoint = (amesosSolver_->SetUseTranspose(true) == 0);
00246   amesosSolver_->SetUseTranspose(oldUseTranspose);
00247   return supportsAdjoint;
00248 }
00249 
00250 
00251 bool AmesosLinearOpWithSolve::solveSupportsSolveMeasureTypeImpl(
00252   EOpTransp M_trans, const SolveMeasureType& solveMeasureType
00253   ) const
00254 {
00255   return true; // I am a direct solver so I should be able to do it all!
00256 }
00257 
00258 
00259 SolveStatus<double>
00260 AmesosLinearOpWithSolve::solveImpl(
00261   const EOpTransp M_trans,
00262   const MultiVectorBase<double> &B,
00263   const Ptr<MultiVectorBase<double> > &X,
00264   const Ptr<const SolveCriteria<double> > solveCriteria
00265   ) const
00266 {
00267   using Teuchos::rcpFromPtr;
00268   using Teuchos::rcpFromRef;
00269   using Teuchos::OSTab;
00270 
00271   Teuchos::Time totalTimer("");
00272   totalTimer.start(true);
00273 
00274   THYRA_FUNC_TIME_MONITOR("Stratimikos: AmesosLOWS");
00275 
00276   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00277   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00278   OSTab tab = this->getOSTab();
00279   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00280     *out << "\nSolving block system using Amesos solver "
00281          << typeName(*amesosSolver_) << " ...\n\n";
00282 
00283   //
00284   // Get the op(...) range and domain maps
00285   //
00286   const EOpTransp amesosOpTransp = real_trans(trans_trans(amesosSolverTransp_,M_trans));
00287   const Epetra_Operator *amesosOp = epetraLP_->GetOperator();
00288   const Epetra_Map
00289     &opRangeMap  = ( amesosOpTransp == NOTRANS
00290       ? amesosOp->OperatorRangeMap()  : amesosOp->OperatorDomainMap() ),
00291     &opDomainMap = ( amesosOpTransp == NOTRANS
00292       ? amesosOp->OperatorDomainMap() : amesosOp->OperatorRangeMap()  );
00293 
00294   //
00295   // Get Epetra_MultiVector views of B and X
00296   //
00297   Teuchos::RCP<const Epetra_MultiVector>
00298     epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
00299   Teuchos::RCP<Epetra_MultiVector>
00300     epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
00301 
00302   //
00303   // Set B and X in the linear problem
00304   //
00305   epetraLP_->SetLHS(&*epetra_X);
00306   epetraLP_->SetRHS(const_cast<Epetra_MultiVector*>(&*epetra_B));
00307   // Above should be okay but cross your fingers!
00308 
00309   //
00310   // Solve the linear system
00311   //
00312   const bool oldUseTranspose = amesosSolver_->UseTranspose();
00313   amesosSolver_->SetUseTranspose(amesosOpTransp==TRANS);
00314   const int err = amesosSolver_->Solve();
00315   TEUCHOS_TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
00316     "Error, the function Solve() on the amesos solver of type\n"
00317     "\'"<<typeName(*amesosSolver_)<<"\' failed with error code "<<err<<"!"
00318     );
00319   amesosSolver_->SetUseTranspose(oldUseTranspose);
00320 
00321   //
00322   // Unset B and X
00323   //
00324   epetraLP_->SetLHS(NULL);
00325   epetraLP_->SetRHS(NULL);
00326   epetra_X = Teuchos::null;
00327   epetra_B = Teuchos::null;
00328 
00329   //
00330   // Scale X if needed
00331   //
00332   if(amesosSolverScalar_!=1.0)
00333     Thyra::scale(1.0/amesosSolverScalar_, X);
00334 
00335   //
00336   // Set the solve status if requested
00337   //
00338   SolveStatus<double> solveStatus;
00339   solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00340   solveStatus.achievedTol = SolveStatus<double>::unknownTolerance();
00341   solveStatus.message =
00342     std::string("Solver ")+typeName(*amesosSolver_)+std::string(" converged!");
00343   
00344   //
00345   // Report the overall time
00346   //
00347   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00348     *out
00349       << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
00350 
00351   return solveStatus;
00352 
00353 }
00354 
00355 
00356 } // end namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends