Skip to content

Commit

Permalink
Implements the achievedTol interface for the GmresPolySolMgr
Browse files Browse the repository at this point in the history
A user request has been received for the Gmres polynomial preconditioned
solver to return the achieved tolerance from the outer iterative method.
This is best done through the achievedTol() method in the abstract solver
manager interface, that is utilized by most of the iterative methods.
This commit implements that function call and includes a test of that call
to ensure it works.
  • Loading branch information
hkthorn committed Apr 13, 2022
1 parent b09a0a0 commit 19c3699
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 4 deletions.
30 changes: 27 additions & 3 deletions packages/belos/src/BelosGmresPolySolMgr.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -152,9 +152,9 @@ template<class ScalarType, class MV, class OP>
class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP> {
private:

typedef MultiVecTraits<ScalarType,MV> MVT;
typedef Teuchos::ScalarTraits<ScalarType> STS;
typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;

typedef Teuchos::ScalarTraits<MagnitudeType> MTS;
typedef Belos::GmresPolyOp<ScalarType,MV,OP> gmres_poly_t;
typedef Belos::GmresPolyMv<ScalarType,MV> gmres_poly_mv_t;

Expand Down Expand Up @@ -217,6 +217,26 @@ class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP> {
*/
Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }

/*! \brief Tolerance achieved by the last \c solve() invocation.
This is the maximum over all right-hand sides' achieved
convergence tolerances, and is set whether or not the solve
actually managed to achieve the desired convergence tolerance.
\warning This solver manager may be use as either a polynomial
preconditioned iterative method or a polynomial preconditioner.
In the later case, where a static polynomial is being applied
through each call to solve(), there is not an outer solve that
can provide the achieved tolerance.
\warning This result may not be meaningful if there was a loss
of accuracy during the outer solve. You should first call \c
isLOADetected() to check for a loss of accuracy during the
last solve.
*/
MagnitudeType achievedTol() const override {
return achievedTol_;
}

/*! \brief Return the timers for this object.
*
* The timers are ordered as follows:
Expand Down Expand Up @@ -331,7 +351,7 @@ class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP> {
#endif

// Current solver values.
MagnitudeType polyTol_;
MagnitudeType polyTol_, achievedTol_;
int maxDegree_, numIters_;
int verbosity_;
bool hasOuterSolver_;
Expand Down Expand Up @@ -363,6 +383,7 @@ template<class ScalarType, class MV, class OP>
GmresPolySolMgr<ScalarType,MV,OP>::GmresPolySolMgr () :
outputStream_ (Teuchos::rcp(outputStream_default_,false)),
polyTol_ (DefaultSolverParameters::polyTol),
achievedTol_(MTS::zero()),
maxDegree_ (maxDegree_default_),
numIters_ (0),
verbosity_ (verbosity_default_),
Expand Down Expand Up @@ -689,6 +710,7 @@ ReturnType GmresPolySolMgr<ScalarType,MV,OP>::solve ()
ret = solver->solve();
numIters_ = solver->getNumIters();
loaDetected_ = solver->isLOADetected();
achievedTol_ = solver->achievedTol();

} // if (hasOuterSolver_ && maxDegree_)
else if (hasOuterSolver_) {
Expand All @@ -704,12 +726,14 @@ ReturnType GmresPolySolMgr<ScalarType,MV,OP>::solve ()
ret = solver->solve();
numIters_ = solver->getNumIters();
loaDetected_ = solver->isLOADetected();
achievedTol_ = solver->achievedTol();

}
else if (maxDegree_) {

// Apply the polynomial to the current linear system
poly_Op_->ApplyPoly( *problem_->getRHS(), *problem_->getLHS() );
achievedTol_ = MTS::one();

}

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -249,8 +249,9 @@ int main(int argc, char *argv[]) {
//
Belos::ReturnType ret = solver->solve();
//
// Compute actual residuals.
std::cout << "Belos::GmresPolySolMgr returned achievedTol: " << solver->achievedTol() << std::endl << std::endl;
//
// Compute actual residuals.
RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
OPT::Apply( *A, *soln, *temp );
MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
Expand Down

0 comments on commit 19c3699

Please sign in to comment.