Classes | Public Types | Public Member Functions | Static Public Attributes | List of all members
LinearSolver::Serial::Gmres< VarsType > Struct Template Reference

A serial GMRES iterative solver for nonsymmetric linear systems of equations. More...

#include <Gmres.hpp>

Public Types

using options = tmpl::list< ConvergenceCriteria, Verbosity, Restart >
 

Public Member Functions

 Gmres (Convergence::Criteria convergence_criteria, ::Verbosity verbosity, size_t restart=0) noexcept
 
 Gmres (const Gmres &)=default
 
Gmresoperator= (const Gmres &)=default
 
 Gmres (Gmres &&)=default
 
Gmresoperator= (Gmres &&)=default
 
void initialize () noexcept
 
const Convergence::Criteriaconvergence_criteria () const noexcept
 
void pup (PUP::er &p) noexcept
 
template<typename LinearOperator , typename SourceType , typename Preconditioner = IdentityPreconditioner<VarsType>, typename IterationCallback = NoIterationCallback>
std::pair< Convergence::HasConverged, VarsType > operator() (LinearOperator &&linear_operator, const SourceType &source, const VarsType &initial_guess, Preconditioner &&preconditioner=IdentityPreconditioner< VarsType >{}, IterationCallback &&iteration_callback=NoIterationCallback{}) const noexcept
 Iteratively solve the problem \(Ax=b\) for \(x\) where \(A\) is the linear_operator and \(b\) is the source, starting \(x\) at initial_guess. More...
 

Static Public Attributes

static constexpr OptionString help
 

Detailed Description

template<typename VarsType>
struct LinearSolver::Serial::Gmres< VarsType >

A serial GMRES iterative solver for nonsymmetric linear systems of equations.

This is an iterative algorithm to solve general linear equations \(Ax=b\) where \(A\) is a linear operator. See [55], chapter 6.5 for a description of the GMRES algorithm and Algorithm 9.6 for this implementation. It is matrix-free, which means the operator \(A\) needs not be provided explicity as a matrix but only the operator action \(A(x)\) must be provided for an argument \(x\).

The GMRES algorithm does not require the operator \(A\) to be symmetric or positive-definite. Note that other algorithms such as conjugate gradients may be more efficient for symmetric positive-definite operators.

Convergence:
Given a set of \(N_A\) equations (e.g. through an \(N_A\times N_A\) matrix) the GMRES algorithm will converge to numerical precision in at most \(N_A\) iterations. However, depending on the properties of the linear operator, an approximate solution can ideally be obtained in only a few iterations. See [55], section 6.11.4 for details on the convergence of the GMRES algorithm.
Restarting:
This implementation of the GMRES algorithm supports restarting, as detailed in [55], section 6.5.5. Since the GMRES algorithm iteratively builds up an orthogonal basis of the solution space the cost of each iteration increases linearly with the number of iterations. Therefore it is sometimes helpful to restart the algorithm every \(N_\mathrm{restart}\) iterations, discarding the set of basis vectors and starting again from the current solution estimate. This strategy can improve the performance of the solver, but note that the solver can stagnate for non-positive-definite operators and is not guaranteed to converge within \(N_A\) iterations anymore. Set the restart argument of the constructor to \(N_\mathrm{restart}\) to activate restarting, or set it to zero to deactivate restarting (default behaviour).
Preconditioning:
This implementation of the GMRES algorithm also supports preconditioning. You can provide a linear operator \(P\) that approximates the inverse of the operator \(A\) to accelerate the convergence of the linear solve. The algorithm is right-preconditioned, which allows the preconditioner to change in every iteration ("flexible" variant). See [55], sections 9.3.2 and 9.4.1 for details. This implementation follows Algorithm 9.6 in [55].
Improvements:
Further improvements can potentially be implemented for this algorithm, see e.g. [3].

Example

INFO("Solve a symmetric 2x2 matrix");
DenseMatrix<double> matrix(2, 2);
matrix(0, 0) = 4.;
matrix(0, 1) = 1.;
matrix(1, 0) = 1.;
matrix(1, 1) = 3.;
const DenseVector<double> source{1., 2.};
const DenseVector<double> initial_guess{2., 1.};
const DenseVector<double> expected_solution{0.0909090909090909,
0.6363636363636364};
const Convergence::Criteria convergence_criteria{2, 1.e-14, 0.};
const Gmres<DenseVector<double>> gmres{convergence_criteria,
::Verbosity::Verbose};
const auto linear_operator =
[&matrix](const DenseVector<double>& arg) noexcept {
return matrix * arg;
};
std::vector<double> recorded_residuals;
const auto result = gmres(
linear_operator, source, initial_guess,
IdentityPreconditioner<DenseVector<double>>{},
[&recorded_residuals](const Convergence::HasConverged& has_converged) {
recorded_residuals.push_back(has_converged.residual_magnitude());
});
const auto& has_converged = result.first;
const auto& solution = result.second;
REQUIRE(has_converged);
CHECK(has_converged.reason() == Convergence::Reason::AbsoluteResidual);
CHECK(has_converged.num_iterations() == 2);
CHECK(has_converged.residual_magnitude() <= 1.e-14);
CHECK(has_converged.initial_residual_magnitude() ==
approx(8.54400374531753));
CHECK_ITERABLE_APPROX(solution, expected_solution);
// The residuals should decrease monotonically
CHECK(recorded_residuals[0] == has_converged.initial_residual_magnitude());
for (size_t i = 1; i < has_converged.num_iterations(); ++i) {
CHECK(recorded_residuals[i] <= recorded_residuals[i - 1]);
}

Member Function Documentation

◆ operator()()

template<typename VarsType >
template<typename LinearOperator , typename SourceType , typename Preconditioner , typename IterationCallback >
std::pair< Convergence::HasConverged, VarsType > LinearSolver::Serial::Gmres< VarsType >::operator() ( LinearOperator &&  linear_operator,
const SourceType &  source,
const VarsType &  initial_guess,
Preconditioner &&  preconditioner = IdentityPreconditioner<VarsType>{},
IterationCallback &&  iteration_callback = NoIterationCallback{} 
) const
noexcept

Iteratively solve the problem \(Ax=b\) for \(x\) where \(A\) is the linear_operator and \(b\) is the source, starting \(x\) at initial_guess.

Optionally provide a preconditioner (see class documentation).

Returns: An instance of Convergence::HasConverged that provides information on the convergence status of the completed solve, and the approximate solution \(x\).

Member Data Documentation

◆ help

template<typename VarsType >
constexpr OptionString LinearSolver::Serial::Gmres< VarsType >::help
staticconstexpr
Initial value:
=
"A serial GMRES iterative solver for nonsymmetric linear systems of\n"
"equations Ax=b. It will converge to numerical precision in at most N_A\n"
"iterations, where N_A is the number of equations represented by the\n"
"linear operator A, but will ideally converge to a reasonable\n"
"approximation of the solution x in only a few iterations.\n"
"\n"
"Restarting: It is sometimes helpful to restart the algorithm every\n"
"N_restart iterations to speed it up. Note that it can stagnate for\n"
"non-positive-definite matrices and is not guaranteed to converge\n"
"within N_A iterations anymore when restarting is activated.\n"
"Activate restarting by setting the 'Restart' option to N_restart, or\n"
"deactivate restarting by setting it to zero (default)."

The documentation for this struct was generated from the following file:
std::vector< double >
Convergence::Criteria
Criteria that determine an iterative algorithm has converged.
Definition: Criteria.hpp:35
DenseVector< double >
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:59
CHECK_ITERABLE_APPROX
#define CHECK_ITERABLE_APPROX(a, b)
A wrapper around Catch's CHECK macro that checks approximate equality of entries in iterable containe...
Definition: TestingFramework.hpp:139
DenseMatrix< double >