Thyra_BelosLinearOpWithSolve_def.hpp

00001 
00002 #ifndef THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
00003 #define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
00004 
00005 #include "Thyra_BelosLinearOpWithSolve_decl.hpp"
00006 #include "Thyra_LinearOpWithSolveHelpers.hpp"
00007 #include "Teuchos_Assert.hpp"
00008 #include "Teuchos_TimeMonitor.hpp"
00009 
00010 
00011 namespace Thyra {
00012 
00013 
00014 // Constructors/initializers/accessors
00015 
00016 
00017 template<class Scalar>
00018 BelosLinearOpWithSolve<Scalar>::BelosLinearOpWithSolve()
00019   :isExternalPrec_(false),
00020   supportSolveUse_(SUPPORT_SOLVE_UNSPECIFIED),
00021   defaultTol_(-1.0)
00022 {}
00023 
00024 
00025 template<class Scalar>
00026 BelosLinearOpWithSolve<Scalar>::BelosLinearOpWithSolve(
00027   const RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > &lp,
00028   const RCP<Teuchos::ParameterList> &solverPL,
00029   const RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > &iterativeSolver,
00030   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00031   const RCP<const PreconditionerBase<Scalar> > &prec,
00032   const bool isExternalPrec,
00033   const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
00034   const ESupportSolveUse &supportSolveUse
00035   )
00036 {
00037   initialize(
00038     lp,solverPL,iterativeSolver,
00039     fwdOpSrc,prec,isExternalPrec,approxFwdOpSrc,supportSolveUse
00040     );
00041 }
00042 
00043 
00044 template<class Scalar>
00045 void BelosLinearOpWithSolve<Scalar>::initialize(
00046   const RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > &lp,
00047   const RCP<Teuchos::ParameterList> &solverPL,
00048   const RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > &iterativeSolver,
00049   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00050   const RCP<const PreconditionerBase<Scalar> > &prec,
00051   const bool isExternalPrec,
00052   const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
00053   const ESupportSolveUse &supportSolveUse
00054   )
00055 {
00056   this->setLinePrefix("BELOS/T");
00057   // ToDo: Validate input
00058   lp_ = lp;
00059   solverPL_ = solverPL;
00060   iterativeSolver_ = iterativeSolver;
00061   fwdOpSrc_ = fwdOpSrc;
00062   prec_ = prec;
00063   isExternalPrec_ = isExternalPrec;
00064   approxFwdOpSrc_ = approxFwdOpSrc;
00065   supportSolveUse_ = supportSolveUse;
00066   // Check if "Convergence Tolerance" is in the solver parameter list.  If
00067   // not, use the default from the solver.
00068   if ( !is_null(solverPL_) ) {
00069     if (solverPL_->isParameter("Convergence Tolerance")) {
00070       defaultTol_ = solverPL_->get<double>("Convergence Tolerance");
00071     }
00072   }
00073   else {
00074     RCP<const Teuchos::ParameterList> defaultPL =
00075       iterativeSolver->getValidParameters();
00076     defaultTol_ = defaultPL->get<double>("Convergence Tolerance");
00077   }
00078 }
00079 
00080 
00081 template<class Scalar>
00082 RCP<const LinearOpSourceBase<Scalar> >
00083 BelosLinearOpWithSolve<Scalar>::extract_fwdOpSrc()
00084 {
00085   RCP<const LinearOpSourceBase<Scalar> >
00086     _fwdOpSrc = fwdOpSrc_;
00087   fwdOpSrc_ = Teuchos::null;
00088   return _fwdOpSrc;
00089 }
00090 
00091 
00092 template<class Scalar>
00093 RCP<const PreconditionerBase<Scalar> >
00094 BelosLinearOpWithSolve<Scalar>::extract_prec()
00095 {
00096   RCP<const PreconditionerBase<Scalar> >
00097     _prec = prec_;
00098   prec_ = Teuchos::null;
00099   return _prec;
00100 }
00101 
00102 
00103 template<class Scalar>
00104 bool BelosLinearOpWithSolve<Scalar>::isExternalPrec() const
00105 {
00106   return isExternalPrec_;
00107 }
00108 
00109 
00110 template<class Scalar>
00111 RCP<const LinearOpSourceBase<Scalar> >
00112 BelosLinearOpWithSolve<Scalar>::extract_approxFwdOpSrc()
00113 {
00114   RCP<const LinearOpSourceBase<Scalar> >
00115     _approxFwdOpSrc = approxFwdOpSrc_;
00116   approxFwdOpSrc_ = Teuchos::null;
00117   return _approxFwdOpSrc;
00118 }
00119 
00120 
00121 template<class Scalar>
00122 ESupportSolveUse BelosLinearOpWithSolve<Scalar>::supportSolveUse() const
00123 {
00124   return supportSolveUse_;
00125 }
00126 
00127 
00128 template<class Scalar>
00129 void BelosLinearOpWithSolve<Scalar>::uninitialize(
00130   RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > *lp,
00131   RCP<Teuchos::ParameterList> *solverPL,
00132   RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > *iterativeSolver,
00133   RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
00134   RCP<const PreconditionerBase<Scalar> > *prec,
00135   bool *isExternalPrec,
00136   RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
00137   ESupportSolveUse *supportSolveUse
00138   )
00139 {
00140   if (lp) *lp = lp_;
00141   if (solverPL) *solverPL = solverPL_;
00142   if (iterativeSolver) *iterativeSolver = iterativeSolver_;
00143   if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
00144   if (prec) *prec = prec_;
00145   if (isExternalPrec) *isExternalPrec = isExternalPrec_;
00146   if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
00147   if (supportSolveUse) *supportSolveUse = supportSolveUse_;
00148 
00149   lp_ = Teuchos::null;
00150   solverPL_ = Teuchos::null;
00151   iterativeSolver_ = Teuchos::null;
00152   fwdOpSrc_ = Teuchos::null;
00153   prec_ = Teuchos::null;
00154   isExternalPrec_ = false;
00155   approxFwdOpSrc_ = Teuchos::null;
00156   supportSolveUse_ = SUPPORT_SOLVE_UNSPECIFIED;
00157 }
00158 
00159 
00160 // Overridden from LinearOpBase
00161 
00162 
00163 template<class Scalar>
00164 RCP< const VectorSpaceBase<Scalar> >
00165 BelosLinearOpWithSolve<Scalar>::range() const
00166 {
00167   if (!is_null(lp_))
00168     return lp_->getOperator()->range();
00169   return Teuchos::null;
00170 }
00171 
00172 
00173 template<class Scalar>
00174 RCP< const VectorSpaceBase<Scalar> >
00175 BelosLinearOpWithSolve<Scalar>::domain() const
00176 {
00177   if (!is_null(lp_))
00178     return lp_->getOperator()->domain();
00179   return Teuchos::null;
00180 }
00181 
00182 
00183 template<class Scalar>
00184 RCP<const LinearOpBase<Scalar> >
00185 BelosLinearOpWithSolve<Scalar>::clone() const
00186 {
00187   return Teuchos::null; // Not supported yet but could be
00188 }
00189 
00190 
00191 // Overridden from Teuchos::Describable
00192 
00193 
00194 template<class Scalar>
00195 std::string BelosLinearOpWithSolve<Scalar>::description() const
00196 {
00197   std::ostringstream oss;
00198   oss << Teuchos::Describable::description();
00199   if ( !is_null(lp_) && !is_null(lp_->getOperator()) ) {
00200     oss << "{";
00201     oss << "iterativeSolver=\'"<<iterativeSolver_->description()<<"\'";
00202     oss << ",fwdOp=\'"<<lp_->getOperator()->description()<<"\'";
00203     if (lp_->getLeftPrec().get())
00204       oss << ",leftPrecOp=\'"<<lp_->getLeftPrec()->description()<<"\'";
00205     if (lp_->getRightPrec().get())
00206       oss << ",rightPrecOp=\'"<<lp_->getRightPrec()->description()<<"\'";
00207     oss << "}";
00208   }
00209   // ToDo: Make Belos::SolverManager derive from Teuchos::Describable so
00210   // that we can get better information.
00211   return oss.str();
00212 }
00213 
00214 
00215 template<class Scalar>
00216 void BelosLinearOpWithSolve<Scalar>::describe(
00217   Teuchos::FancyOStream &out_arg,
00218   const Teuchos::EVerbosityLevel verbLevel
00219   ) const
00220 {
00221   typedef Teuchos::ScalarTraits<Scalar> ST;
00222   using Teuchos::FancyOStream;
00223   using Teuchos::OSTab;
00224   using Teuchos::describe;
00225   RCP<FancyOStream> out = rcp(&out_arg,false);
00226   OSTab tab(out);
00227   switch (verbLevel) {
00228     case Teuchos::VERB_DEFAULT:
00229     case Teuchos::VERB_LOW:
00230       *out << this->description() << std::endl;
00231       break;
00232     case Teuchos::VERB_MEDIUM:
00233     case Teuchos::VERB_HIGH:
00234     case Teuchos::VERB_EXTREME:
00235     {
00236       *out
00237         << Teuchos::Describable::description()<< "{"
00238         << "rangeDim=" << this->range()->dim()
00239         << ",domainDim=" << this->domain()->dim() << "}\n";
00240       if (lp_->getOperator().get()) {
00241         OSTab tab(out);
00242         *out
00243           << "iterativeSolver = "<<describe(*iterativeSolver_,verbLevel)
00244           << "fwdOp = " << describe(*lp_->getOperator(),verbLevel);
00245         if (lp_->getLeftPrec().get())
00246           *out << "leftPrecOp = "<<describe(*lp_->getLeftPrec(),verbLevel);
00247         if (lp_->getRightPrec().get())
00248           *out << "rightPrecOp = "<<describe(*lp_->getRightPrec(),verbLevel);
00249       }
00250       break;
00251     }
00252     default:
00253       TEST_FOR_EXCEPT(true); // Should never get here!
00254   }
00255 }
00256 
00257 
00258 // protected
00259 
00260 
00261 // Overridden from LinearOpBase
00262 
00263 
00264 template<class Scalar>
00265 bool BelosLinearOpWithSolve<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00266 {
00267   return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
00268 }
00269 
00270 
00271 template<class Scalar>
00272 void BelosLinearOpWithSolve<Scalar>::applyImpl(
00273   const EOpTransp M_trans,
00274   const MultiVectorBase<Scalar> &X,
00275   const Ptr<MultiVectorBase<Scalar> > &Y,
00276   const Scalar alpha,
00277   const Scalar beta
00278   ) const
00279 {
00280   ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
00281 }
00282 
00283 
00284 // Overridden from LinearOpWithSolveBase
00285 
00286 
00287 template<class Scalar>
00288 bool
00289 BelosLinearOpWithSolve<Scalar>::solveSupportsImpl(EOpTransp M_trans) const
00290 {
00291   if (real_trans(M_trans)==NOTRANS) return true;
00292   return false; // ToDo: Support adjoint solves!
00293 }
00294 
00295 
00296 template<class Scalar>
00297 bool
00298 BelosLinearOpWithSolve<Scalar>::solveSupportsSolveMeasureTypeImpl(
00299   EOpTransp M_trans, const SolveMeasureType& solveMeasureType) const
00300 {
00301   if (real_trans(M_trans)==NOTRANS) {
00302     return (
00303       solveMeasureType.useDefault()
00304       ||
00305       solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
00306       ||
00307       solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL)
00308       );
00309   }
00310   // TRANS
00311   return false; // ToDo: Support adjoint solves!
00312 }
00313 
00314 
00315 template<class Scalar>
00316 SolveStatus<Scalar>
00317 BelosLinearOpWithSolve<Scalar>::solveImpl(
00318   const EOpTransp M_trans,
00319   const MultiVectorBase<Scalar> &B,
00320   const Ptr<MultiVectorBase<Scalar> > &X,
00321   const Ptr<const SolveCriteria<Scalar> > solveCriteria
00322   ) const
00323 {
00324 
00325   TEUCHOS_FUNC_TIME_MONITOR("BelosLOWS");
00326 
00327   using Teuchos::rcp;
00328   using Teuchos::rcpFromRef;
00329   using Teuchos::rcpFromPtr;
00330   using Teuchos::FancyOStream;
00331   using Teuchos::OSTab;
00332   using Teuchos::describe;
00333   typedef Teuchos::ScalarTraits<Scalar> ST;
00334   typedef typename ST::magnitudeType ScalarMag;
00335   Teuchos::Time totalTimer(""), timer("");
00336   totalTimer.start(true);
00337   
00338   TEUCHOS_ASSERT(this->solveSupportsImpl(M_trans));
00339 
00340   const int numRhs = B.domain()->dim();
00341   const int numEquations = B.range()->dim();
00342 
00343   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00344   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00345   OSTab tab = this->getOSTab();
00346   if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE)) {
00347     *out << "\nStarting iterations with Belos:\n";
00348     OSTab tab2(out);
00349     *out << "Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
00350     *out << "Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
00351     *out << "With #Eqns="<<numEquations<<", #RHSs="<<numRhs<<" ...\n";
00352   }
00353 
00354   //
00355   // Set RHS and LHS
00356   //
00357 
00358   bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
00359   TEST_FOR_EXCEPTION(
00360     ret == false, CatastrophicSolveFailure
00361     ,"Error, the Belos::LinearProblem could not be set for the current solve!"
00362     );
00363 
00364   //
00365   // Set the solution criteria
00366   //
00367 
00368   const RCP<Teuchos::ParameterList> tmpPL = Teuchos::parameterList();
00369 
00370   SolveMeasureType solveMeasureType;
00371   if (nonnull(solveCriteria)) {
00372     solveMeasureType = solveCriteria->solveMeasureType;
00373     const ScalarMag requestedTol = solveCriteria->requestedTol;
00374     assertSupportsSolveMeasureType(*this,M_trans,solveMeasureType);
00375     if ( solveMeasureType.useDefault() ) {
00376       tmpPL->set("Convergence Tolerance", defaultTol_);
00377     }
00378     else if ( solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS) ) {
00379       if ( requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance() )
00380         tmpPL->set("Convergence Tolerance", requestedTol);
00381       else
00382         tmpPL->set("Convergence Tolerance", defaultTol_);
00383       tmpPL->set("Explicit Residual Scaling", "Norm of RHS");
00384     }
00385     else if ( solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL) ) {
00386       if ( requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance() )
00387         tmpPL->set("Convergence Tolerance", requestedTol);
00388       else
00389         tmpPL->set("Convergence Tolerance", defaultTol_);
00390       tmpPL->set("Explicit Residual Scaling", "Norm of Initial Residual");
00391     }
00392     else {
00393       TEST_FOR_EXCEPT(true); // Should never get there.
00394     }
00395   }
00396   else {
00397     tmpPL->set("Convergence Tolerance", defaultTol_);
00398   }
00399 
00400   //
00401   // Reset the blocksize if we adding more vectors than half the number of equations,
00402   // orthogonalization will fail on the first iteration!
00403   //
00404 
00405   RCP<const Teuchos::ParameterList> solverParams = iterativeSolver_->getCurrentParameters();
00406   const int currBlockSize = Teuchos::getParameter<int>(*solverParams, "Block Size");
00407   bool isNumBlocks = false;
00408   int currNumBlocks = 0;
00409   if (Teuchos::isParameterType<int>(*solverParams, "Num Blocks")) {
00410     currNumBlocks = Teuchos::getParameter<int>(*solverParams, "Num Blocks");
00411     isNumBlocks = true;
00412   }
00413   const int newBlockSize = TEUCHOS_MIN(currBlockSize,numEquations/2);
00414   if (nonnull(out)
00415     && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE)
00416     && newBlockSize != currBlockSize)
00417   {
00418     *out << "\nAdjusted block size = " << newBlockSize << "\n";
00419   }
00420   //
00421   tmpPL->set("Block Size",newBlockSize);
00422 
00423   //
00424   // Set the number of Krylov blocks if we are using a GMRES solver, or a solver
00425   // that recognizes "Num Blocks". Otherwise the solver will throw an error!
00426   //
00427 
00428   if (isNumBlocks) {
00429     const int Krylov_length = (currNumBlocks*currBlockSize)/newBlockSize;
00430     tmpPL->set("Num Blocks",Krylov_length);
00431   
00432     if (newBlockSize != currBlockSize) {
00433       if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00434         *out
00435           << "\nAdjusted max number of Krylov basis blocks = " << Krylov_length << "\n";
00436     }
00437   }
00438 
00439   //
00440   // Solve the linear system
00441   //
00442 
00443   Belos::ReturnType belosSolveStatus;
00444   {
00445     RCP<std::ostream>
00446       outUsed =
00447       ( static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE)
00448         ? out
00449         : rcp(new FancyOStream(rcp(new Teuchos::oblackholestream())))
00450         );
00451     Teuchos::OSTab tab(outUsed,1,"BELOS");
00452     tmpPL->set("Output Stream", outUsed);
00453     iterativeSolver_->setParameters( tmpPL );
00454     belosSolveStatus = iterativeSolver_->solve();
00455   }
00456 
00457   //
00458   // Report the solve status
00459   //
00460 
00461   totalTimer.stop();
00462 
00463   SolveStatus<Scalar> solveStatus;
00464 
00465   switch (belosSolveStatus) {
00466     case Belos::Unconverged:
00467       solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
00468       break;
00469     case Belos::Converged:
00470       solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00471       solveStatus.achievedTol = tmpPL->get("Convergence Tolerance", defaultTol_);
00472       break;
00473     default:
00474       TEST_FOR_EXCEPT(true); // Should never get here!
00475   }
00476 
00477   std::ostringstream ossmessage;
00478   ossmessage
00479     << "The Belos solver of type \""<<iterativeSolver_->description()
00480     <<"\" returned a solve status of \""<< toString(solveStatus.solveStatus)
00481     /* "\" in " << iterations << " iterations and achieved an approximate max tolerance of "
00482        << belosAchievedTol << 
00483     */
00484     << " with total CPU time of " << totalTimer.totalElapsedTime() << " sec" ;
00485   if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00486     *out << "\n" << ossmessage.str() << "\n";
00487 
00488   solveStatus.message = ossmessage.str();
00489 
00490   if (out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00491     *out << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
00492 
00493   return solveStatus;
00494 
00495 }
00496 
00497 
00498 } // end namespace Thyra
00499 
00500 
00501 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
 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