Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Thyra_BelosLinearOpWithSolve_def.hpp
Go to the documentation of this file.
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_GeneralSolveCriteriaBelosStatusTest.hpp"
00007 #include "Thyra_LinearOpWithSolveHelpers.hpp"
00008 #include "Teuchos_DebugDefaultAsserts.hpp"
00009 #include "Teuchos_Assert.hpp"
00010 #include "Teuchos_TimeMonitor.hpp"
00011 
00012 namespace {
00013   // Set the Belos solver's parameter list to scale its residual norms
00014   // in the specified way.
00015   //
00016   // We break this out in a separate function because the parameters
00017   // to set depend on which parameters the Belos solver supports.  Not
00018   // all Belos solvers support both the "Implicit Residual Scaling"
00019   // and "Explicit Residual Scaling" parameters, so we have to check
00020   // the solver's list of valid parameters for the existence of these.
00021   //
00022   // Scaling options: Belos lets you decide whether the solver will
00023   // scale residual norms by the (left-)preconditioned initial
00024   // residual norms (residualScalingType = "Norm of Initial
00025   // Residual"), or by the unpreconditioned initial residual norms
00026   // (residualScalingType = "Norm of RHS").  Usually you want to scale
00027   // by the unpreconditioned initial residual norms.  This is because
00028   // preconditioning is just an optimization, and you really want to
00029   // make ||B - AX|| small, rather than ||M B - M (A X)||.  If you're
00030   // measuring ||B - AX|| and scaling by the initial residual, you
00031   // should use the unpreconditioned initial residual to match it.
00032   //
00033   // Note, however, that the implicit residual test computes
00034   // left-preconditioned residuals, if a left preconditioner was
00035   // provided.  That's OK because when Belos solvers (at least the
00036   // GMRES variants) are given a left preconditioner, they first check
00037   // the implicit residuals.  If those converge, they then check the
00038   // explicit residuals.  The explicit residual test does _not_ apply
00039   // the left preconditioner when computing the residual.  The
00040   // implicit residual test is just an optimization so that Belos
00041   // doesn't have to compute explicit residuals B - A*X at every
00042   // iteration.  This is why we use the same scaling factor for both
00043   // the implicit and explicit residuals.
00044   //
00045   // Arguments:
00046   //
00047   // solverParams [in/out] Parameters for the current solve.
00048   //
00049   // solverValidParams [in] Valid parameter list for the Belos solver.
00050   //   Result of calling the solver's getValidParameters() method.
00051   //
00052   // residualScalingType [in] String describing how the solver should
00053   //   scale residuals.  Valid values include "Norm of RHS" and "Norm
00054   //   of Initial Residual" (these are the only two options this file
00055   //   currently uses, though Belos offers other options).
00056   void
00057   setResidualScalingType (const Teuchos::RCP<Teuchos::ParameterList>& solverParams,
00058         const Teuchos::RCP<const Teuchos::ParameterList>& solverValidParams,
00059         const std::string& residualScalingType)
00060   {
00061     // Many Belos solvers (especially the GMRES variants) define both
00062     // "Implicit Residual Scaling" and "Explicit Residual Scaling"
00063     // options.
00064     //
00065     // "Implicit" means "the left-preconditioned approximate
00066     // a.k.a. 'recursive' residual as computed by the Krylov method."
00067     // 
00068     // "Explicit" means ||B - A*X||, the unpreconditioned, "exact"
00069     // residual.
00070     //
00071     // Belos' GMRES implementations chain these two tests in sequence.
00072     // Implicit comes first, and explicit is not evaluated unless
00073     // implicit passes.  In some cases (e.g., no left preconditioner),
00074     // GMRES _only_ uses the implicit tests.  This means that only
00075     // setting "Explicit Residual Scaling" won't change the solver's
00076     // behavior.  Stratimikos tends to prefer using a right
00077     // preconditioner, in which case setting only the "Explicit
00078     // Residual Scaling" argument has no effect.  Furthermore, if
00079     // "Explicit Residual Scaling" is set to something other than the
00080     // default (initial residual norm), without "Implicit Residual
00081     // Scaling" getting the same setting, then the implicit residual
00082     // test will be using a radically different scaling factor than
00083     // the user wanted.
00084     // 
00085     // Not all Belos solvers support both options.  We check the
00086     // solver's valid parameter list first before attempting to set
00087     // the option.
00088     if (solverValidParams->isParameter ("Implicit Residual Scaling")) {
00089       solverParams->set ("Implicit Residual Scaling", residualScalingType);
00090     }
00091     if (solverValidParams->isParameter ("Explicit Residual Scaling")) {
00092       solverParams->set ("Explicit Residual Scaling", residualScalingType);
00093     }
00094   }
00095 
00096 } // namespace (anonymous)
00097 
00098 
00099 namespace Thyra {
00100 
00101 
00102 // Constructors/initializers/accessors
00103 
00104 
00105 template<class Scalar>
00106 BelosLinearOpWithSolve<Scalar>::BelosLinearOpWithSolve()
00107   :convergenceTestFrequency_(-1),
00108   isExternalPrec_(false),
00109   supportSolveUse_(SUPPORT_SOLVE_UNSPECIFIED),
00110   defaultTol_(-1.0)
00111 {}
00112 
00113 
00114 template<class Scalar>
00115 void BelosLinearOpWithSolve<Scalar>::initialize(
00116   const RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > &lp,
00117   const RCP<Teuchos::ParameterList> &solverPL,
00118   const RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > &iterativeSolver,
00119   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00120   const RCP<const PreconditionerBase<Scalar> > &prec,
00121   const bool isExternalPrec,
00122   const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
00123   const ESupportSolveUse &supportSolveUse,
00124   const int convergenceTestFrequency
00125   )
00126 {
00127   this->setLinePrefix("BELOS/T");
00128   // ToDo: Validate input
00129   lp_ = lp;
00130   solverPL_ = solverPL;
00131   iterativeSolver_ = iterativeSolver;
00132   fwdOpSrc_ = fwdOpSrc;
00133   prec_ = prec;
00134   isExternalPrec_ = isExternalPrec;
00135   approxFwdOpSrc_ = approxFwdOpSrc;
00136   supportSolveUse_ = supportSolveUse;
00137   convergenceTestFrequency_ = convergenceTestFrequency;
00138   // Check if "Convergence Tolerance" is in the solver parameter list.  If
00139   // not, use the default from the solver.
00140   if ( !is_null(solverPL_) ) {
00141     if (solverPL_->isParameter("Convergence Tolerance")) {
00142       defaultTol_ = solverPL_->get<double>("Convergence Tolerance");
00143     }
00144   }
00145   else {
00146     RCP<const Teuchos::ParameterList> defaultPL =
00147       iterativeSolver->getValidParameters();
00148     defaultTol_ = defaultPL->get<double>("Convergence Tolerance");
00149   }
00150 }
00151 
00152 
00153 template<class Scalar>
00154 RCP<const LinearOpSourceBase<Scalar> >
00155 BelosLinearOpWithSolve<Scalar>::extract_fwdOpSrc()
00156 {
00157   RCP<const LinearOpSourceBase<Scalar> >
00158     _fwdOpSrc = fwdOpSrc_;
00159   fwdOpSrc_ = Teuchos::null;
00160   return _fwdOpSrc;
00161 }
00162 
00163 
00164 template<class Scalar>
00165 RCP<const PreconditionerBase<Scalar> >
00166 BelosLinearOpWithSolve<Scalar>::extract_prec()
00167 {
00168   RCP<const PreconditionerBase<Scalar> >
00169     _prec = prec_;
00170   prec_ = Teuchos::null;
00171   return _prec;
00172 }
00173 
00174 
00175 template<class Scalar>
00176 bool BelosLinearOpWithSolve<Scalar>::isExternalPrec() const
00177 {
00178   return isExternalPrec_;
00179 }
00180 
00181 
00182 template<class Scalar>
00183 RCP<const LinearOpSourceBase<Scalar> >
00184 BelosLinearOpWithSolve<Scalar>::extract_approxFwdOpSrc()
00185 {
00186   RCP<const LinearOpSourceBase<Scalar> >
00187     _approxFwdOpSrc = approxFwdOpSrc_;
00188   approxFwdOpSrc_ = Teuchos::null;
00189   return _approxFwdOpSrc;
00190 }
00191 
00192 
00193 template<class Scalar>
00194 ESupportSolveUse BelosLinearOpWithSolve<Scalar>::supportSolveUse() const
00195 {
00196   return supportSolveUse_;
00197 }
00198 
00199 
00200 template<class Scalar>
00201 void BelosLinearOpWithSolve<Scalar>::uninitialize(
00202   RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > *lp,
00203   RCP<Teuchos::ParameterList> *solverPL,
00204   RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > *iterativeSolver,
00205   RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
00206   RCP<const PreconditionerBase<Scalar> > *prec,
00207   bool *isExternalPrec,
00208   RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
00209   ESupportSolveUse *supportSolveUse
00210   )
00211 {
00212   if (lp) *lp = lp_;
00213   if (solverPL) *solverPL = solverPL_;
00214   if (iterativeSolver) *iterativeSolver = iterativeSolver_;
00215   if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
00216   if (prec) *prec = prec_;
00217   if (isExternalPrec) *isExternalPrec = isExternalPrec_;
00218   if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
00219   if (supportSolveUse) *supportSolveUse = supportSolveUse_;
00220 
00221   lp_ = Teuchos::null;
00222   solverPL_ = Teuchos::null;
00223   iterativeSolver_ = Teuchos::null;
00224   fwdOpSrc_ = Teuchos::null;
00225   prec_ = Teuchos::null;
00226   isExternalPrec_ = false;
00227   approxFwdOpSrc_ = Teuchos::null;
00228   supportSolveUse_ = SUPPORT_SOLVE_UNSPECIFIED;
00229 }
00230 
00231 
00232 // Overridden from LinearOpBase
00233 
00234 
00235 template<class Scalar>
00236 RCP< const VectorSpaceBase<Scalar> >
00237 BelosLinearOpWithSolve<Scalar>::range() const
00238 {
00239   if (!is_null(lp_))
00240     return lp_->getOperator()->range();
00241   return Teuchos::null;
00242 }
00243 
00244 
00245 template<class Scalar>
00246 RCP< const VectorSpaceBase<Scalar> >
00247 BelosLinearOpWithSolve<Scalar>::domain() const
00248 {
00249   if (!is_null(lp_))
00250     return lp_->getOperator()->domain();
00251   return Teuchos::null;
00252 }
00253 
00254 
00255 template<class Scalar>
00256 RCP<const LinearOpBase<Scalar> >
00257 BelosLinearOpWithSolve<Scalar>::clone() const
00258 {
00259   return Teuchos::null; // Not supported yet but could be
00260 }
00261 
00262 
00263 // Overridden from Teuchos::Describable
00264 
00265 
00266 template<class Scalar>
00267 std::string BelosLinearOpWithSolve<Scalar>::description() const
00268 {
00269   std::ostringstream oss;
00270   oss << Teuchos::Describable::description();
00271   if ( !is_null(lp_) && !is_null(lp_->getOperator()) ) {
00272     oss << "{";
00273     oss << "iterativeSolver=\'"<<iterativeSolver_->description()<<"\'";
00274     oss << ",fwdOp=\'"<<lp_->getOperator()->description()<<"\'";
00275     if (lp_->getLeftPrec().get())
00276       oss << ",leftPrecOp=\'"<<lp_->getLeftPrec()->description()<<"\'";
00277     if (lp_->getRightPrec().get())
00278       oss << ",rightPrecOp=\'"<<lp_->getRightPrec()->description()<<"\'";
00279     oss << "}";
00280   }
00281   // ToDo: Make Belos::SolverManager derive from Teuchos::Describable so
00282   // that we can get better information.
00283   return oss.str();
00284 }
00285 
00286 
00287 template<class Scalar>
00288 void BelosLinearOpWithSolve<Scalar>::describe(
00289   Teuchos::FancyOStream &out_arg,
00290   const Teuchos::EVerbosityLevel verbLevel
00291   ) const
00292 {
00293   typedef Teuchos::ScalarTraits<Scalar> ST;
00294   using Teuchos::FancyOStream;
00295   using Teuchos::OSTab;
00296   using Teuchos::describe;
00297   RCP<FancyOStream> out = rcp(&out_arg,false);
00298   OSTab tab(out);
00299   switch (verbLevel) {
00300     case Teuchos::VERB_DEFAULT:
00301     case Teuchos::VERB_LOW:
00302       *out << this->description() << std::endl;
00303       break;
00304     case Teuchos::VERB_MEDIUM:
00305     case Teuchos::VERB_HIGH:
00306     case Teuchos::VERB_EXTREME:
00307     {
00308       *out
00309         << Teuchos::Describable::description()<< "{"
00310         << "rangeDim=" << this->range()->dim()
00311         << ",domainDim=" << this->domain()->dim() << "}\n";
00312       if (lp_->getOperator().get()) {
00313         OSTab tab(out);
00314         *out
00315           << "iterativeSolver = "<<describe(*iterativeSolver_,verbLevel)
00316           << "fwdOp = " << describe(*lp_->getOperator(),verbLevel);
00317         if (lp_->getLeftPrec().get())
00318           *out << "leftPrecOp = "<<describe(*lp_->getLeftPrec(),verbLevel);
00319         if (lp_->getRightPrec().get())
00320           *out << "rightPrecOp = "<<describe(*lp_->getRightPrec(),verbLevel);
00321       }
00322       break;
00323     }
00324     default:
00325       TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
00326   }
00327 }
00328 
00329 
00330 // protected
00331 
00332 
00333 // Overridden from LinearOpBase
00334 
00335 
00336 template<class Scalar>
00337 bool BelosLinearOpWithSolve<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00338 {
00339   return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
00340 }
00341 
00342 
00343 template<class Scalar>
00344 void BelosLinearOpWithSolve<Scalar>::applyImpl(
00345   const EOpTransp M_trans,
00346   const MultiVectorBase<Scalar> &X,
00347   const Ptr<MultiVectorBase<Scalar> > &Y,
00348   const Scalar alpha,
00349   const Scalar beta
00350   ) const
00351 {
00352   ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
00353 }
00354 
00355 
00356 // Overridden from LinearOpWithSolveBase
00357 
00358 
00359 template<class Scalar>
00360 bool
00361 BelosLinearOpWithSolve<Scalar>::solveSupportsImpl(EOpTransp M_trans) const
00362 {
00363   return solveSupportsNewImpl(M_trans, Teuchos::null);
00364 }
00365 
00366 
00367 template<class Scalar>
00368 bool
00369 BelosLinearOpWithSolve<Scalar>::solveSupportsNewImpl(EOpTransp transp,
00370   const Ptr<const SolveCriteria<Scalar> > solveCriteria) const
00371 {
00372   // Only support forward solve right now!
00373   if (real_trans(transp)==NOTRANS) return true;
00374   return false; // ToDo: Support adjoint solves!
00375   // Otherwise, Thyra/Belos now supports every solve criteria type that exists
00376   // because of the class Thyra::GeneralSolveCriteriaBelosStatusTest!
00377   return true;
00378 /*
00379   if (real_trans(M_trans)==NOTRANS) {
00380     return (
00381       solveMeasureType.useDefault()
00382       ||
00383       solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
00384       ||
00385       solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL)
00386       );
00387   }
00388 */
00389 }
00390 
00391 
00392 template<class Scalar>
00393 bool
00394 BelosLinearOpWithSolve<Scalar>::solveSupportsSolveMeasureTypeImpl(
00395   EOpTransp M_trans, const SolveMeasureType& solveMeasureType) const
00396 {
00397   SolveCriteria<Scalar> solveCriteria(solveMeasureType, SolveCriteria<Scalar>::unspecifiedTolerance());
00398   return solveSupportsNewImpl(M_trans, Teuchos::constOptInArg(solveCriteria));
00399 }
00400 
00401 
00402 template<class Scalar>
00403 SolveStatus<Scalar>
00404 BelosLinearOpWithSolve<Scalar>::solveImpl(
00405   const EOpTransp M_trans,
00406   const MultiVectorBase<Scalar> &B,
00407   const Ptr<MultiVectorBase<Scalar> > &X,
00408   const Ptr<const SolveCriteria<Scalar> > solveCriteria
00409   ) const
00410 {
00411 
00412   THYRA_FUNC_TIME_MONITOR("Stratimikos: BelosLOWS");
00413 
00414   using Teuchos::rcp;
00415   using Teuchos::rcpFromRef;
00416   using Teuchos::rcpFromPtr;
00417   using Teuchos::FancyOStream;
00418   using Teuchos::OSTab;
00419   using Teuchos::ParameterList;
00420   using Teuchos::parameterList;
00421   using Teuchos::describe;
00422   typedef Teuchos::ScalarTraits<Scalar> ST;
00423   typedef typename ST::magnitudeType ScalarMag;
00424   Teuchos::Time totalTimer(""), timer("");
00425   totalTimer.start(true);
00426 
00427   assertSolveSupports(*this, M_trans, solveCriteria);
00428   // 2010/08/22: rabartl: Bug 4915 ToDo: Move the above into the NIV function
00429   // solve(...).
00430 
00431   const RCP<FancyOStream> out = this->getOStream();
00432   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00433   OSTab tab = this->getOSTab();
00434   if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW)) {
00435     *out << "\nStarting iterations with Belos:\n";
00436     OSTab tab2(out);
00437     *out << "Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
00438     *out << "Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
00439     *out << "With #Eqns="<<B.range()->dim()<<", #RHSs="<<B.domain()->dim()<<" ...\n";
00440   }
00441 
00442   //
00443   // Set RHS and LHS
00444   //
00445 
00446   bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
00447   TEUCHOS_TEST_FOR_EXCEPTION(
00448     ret == false, CatastrophicSolveFailure
00449     ,"Error, the Belos::LinearProblem could not be set for the current solve!"
00450     );
00451 
00452   //
00453   // Set the solution criteria
00454   //
00455 
00456   // Parameter list for the current solve.
00457   const RCP<ParameterList> tmpPL = Teuchos::parameterList();
00458 
00459   // The solver's valid parameter list.
00460   RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();
00461 
00462   SolveMeasureType solveMeasureType;
00463   RCP<GeneralSolveCriteriaBelosStatusTest<Scalar> > generalSolveCriteriaBelosStatusTest;
00464   if (nonnull(solveCriteria)) {
00465     solveMeasureType = solveCriteria->solveMeasureType;
00466     const ScalarMag requestedTol = solveCriteria->requestedTol;
00467     if (solveMeasureType.useDefault()) {
00468       tmpPL->set("Convergence Tolerance", defaultTol_);
00469     }
00470     else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
00471       if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
00472         tmpPL->set("Convergence Tolerance", requestedTol);
00473       }
00474       else {
00475         tmpPL->set("Convergence Tolerance", defaultTol_);
00476       }
00477       setResidualScalingType (tmpPL, validPL, "Norm of RHS");
00478     }
00479     else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
00480       if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
00481         tmpPL->set("Convergence Tolerance", requestedTol);
00482       }
00483       else {
00484         tmpPL->set("Convergence Tolerance", defaultTol_);
00485       }
00486       setResidualScalingType (tmpPL, validPL, "Norm of Initial Residual");
00487     }
00488     else {
00489       // Set the most generic (and inefficient) solve criteria
00490       generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
00491         *solveCriteria, convergenceTestFrequency_);
00492       // Set the verbosity level (one level down)
00493       generalSolveCriteriaBelosStatusTest->setOStream(out);
00494       generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
00495       // Set the default convergence tolerance to always converged to allow
00496       // the above status test to control things.
00497       tmpPL->set("Convergence Tolerance", 1.0);
00498     }
00499     // maximum iterations
00500     if (nonnull(solveCriteria->extraParameters)) {
00501       if (Teuchos::isParameterType<int>(*solveCriteria->extraParameters,"Maximum Iterations")) {
00502         tmpPL->set("Maximum Iterations", Teuchos::get<int>(*solveCriteria->extraParameters,"Maximum Iterations"));
00503       }
00504     }
00505   }
00506   else {
00507     // No solveCriteria was even passed in!
00508     tmpPL->set("Convergence Tolerance", defaultTol_);
00509   }
00510 
00511   //
00512   // Solve the linear system
00513   //
00514 
00515   Belos::ReturnType belosSolveStatus;
00516   {
00517     RCP<std::ostream>
00518       outUsed =
00519       ( static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW)
00520         ? out
00521         : rcp(new FancyOStream(rcp(new Teuchos::oblackholestream())))
00522         );
00523     Teuchos::OSTab tab(outUsed,1,"BELOS");
00524     tmpPL->set("Output Stream", outUsed);
00525     iterativeSolver_->setParameters(tmpPL);
00526     if (nonnull(generalSolveCriteriaBelosStatusTest)) {
00527       iterativeSolver_->setUserConvStatusTest(generalSolveCriteriaBelosStatusTest);
00528     }
00529     belosSolveStatus = iterativeSolver_->solve();
00530   }
00531 
00532   //
00533   // Report the solve status
00534   //
00535 
00536   totalTimer.stop();
00537 
00538   SolveStatus<Scalar> solveStatus;
00539 
00540   switch (belosSolveStatus) {
00541     case Belos::Unconverged: {
00542       solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
00543       // Set achievedTol even if the solver did not converge.  This is
00544       // helpful for things like nonlinear solvers, which might be
00545       // able to use a partially converged result, and which would
00546       // like to know the achieved convergence tolerance for use in
00547       // computing bounds.  It's also helpful for estimating whether a
00548       // small increase in the maximum iteration count might be
00549       // helpful next time.
00550       try {
00551   // Some solvers might not have implemented achievedTol(). 
00552   // The default implementation throws std::runtime_error.
00553   solveStatus.achievedTol = iterativeSolver_->achievedTol();
00554       } catch (std::runtime_error&) {
00555   // Do nothing; use the default value of achievedTol.
00556       }
00557       break;
00558     }
00559     case Belos::Converged: {
00560       solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00561       if (nonnull(generalSolveCriteriaBelosStatusTest)) {
00562   // The user set a custom status test.  This means that we
00563   // should ask the custom status test itself, rather than the
00564   // Belos solver, what the final achieved convergence tolerance
00565   // was.
00566         const ArrayView<const ScalarMag> achievedTol = 
00567           generalSolveCriteriaBelosStatusTest->achievedTol();
00568         solveStatus.achievedTol = ST::zero();
00569         for (Ordinal i = 0; i < achievedTol.size(); ++i) {
00570           solveStatus.achievedTol = std::max(solveStatus.achievedTol, achievedTol[i]);
00571         }
00572       }
00573       else {
00574   try {
00575     // Some solvers might not have implemented achievedTol(). 
00576     // The default implementation throws std::runtime_error.
00577     solveStatus.achievedTol = iterativeSolver_->achievedTol();
00578   } catch (std::runtime_error&) {
00579     // Use the default convergence tolerance.  This is a correct
00580     // upper bound, since we did actually converge.
00581     solveStatus.achievedTol = tmpPL->get("Convergence Tolerance", defaultTol_);
00582   }
00583       }
00584       break;
00585     }
00586     TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT();
00587   }
00588 
00589   std::ostringstream ossmessage;
00590   ossmessage
00591     << "The Belos solver of type \""<<iterativeSolver_->description()
00592     <<"\" returned a solve status of \""<< toString(solveStatus.solveStatus) << "\""
00593     << " in " << iterativeSolver_->getNumIters() << " iterations"
00594     << " with total CPU time of " << totalTimer.totalElapsedTime() << " sec" ;
00595   if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00596     *out << "\n" << ossmessage.str() << "\n";
00597 
00598   solveStatus.message = ossmessage.str();
00599 
00600   // Dump the getNumIters() and the achieved convergence tolerance
00601   // into solveStatus.extraParameters, as the "Belos/Iteration Count"
00602   // resp. "Belos/Achieved Tolerance" parameters.
00603   if (solveStatus.extraParameters.is_null()) {
00604     solveStatus.extraParameters = parameterList ();
00605   }
00606   solveStatus.extraParameters->set ("Belos/Iteration Count", 
00607             iterativeSolver_->getNumIters());\
00608   // package independent version of the same
00609   solveStatus.extraParameters->set ("Iteration Count", 
00610             iterativeSolver_->getNumIters());\
00611   // NOTE (mfh 13 Dec 2011) Though the most commonly used Belos
00612   // solvers do implement achievedTol(), some Belos solvers currently
00613   // do not.  In the latter case, if the solver did not converge, the
00614   // reported achievedTol() value may just be the default "invalid"
00615   // value -1, and if the solver did converge, the reported value will
00616   // just be the convergence tolerance (a correct upper bound).
00617   solveStatus.extraParameters->set ("Belos/Achieved Tolerance", 
00618             solveStatus.achievedTol);
00619 
00620 //  This information is in the previous line, which is printed anytime the verbosity
00621 //  is not set to Teuchos::VERB_NONE, so I'm commenting this out for now.
00622 //  if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00623 //    *out << "\nTotal solve time in Belos = "<<totalTimer.totalElapsedTime()<<" sec\n";
00624 
00625   return solveStatus;
00626 
00627 }
00628 
00629 
00630 } // end namespace Thyra
00631 
00632 
00633 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines