A serial GMRES iterative solver for nonsymmetric linear systems of equations.
More...
|
|
| Gmres (Convergence::Criteria convergence_criteria, ::Verbosity verbosity, std::optional< size_t > restart=std::nullopt, std::optional< typename Base::PreconditionerType > local_preconditioner=std::nullopt, const Options::Context &context={}) |
| |
|
| Gmres (Gmres &&)=default |
| |
|
Gmres & | operator= (Gmres &&)=default |
| |
|
| Gmres (const Gmres &rhs) |
| |
|
Gmres & | operator= (const Gmres &rhs) |
| |
|
std::unique_ptr< LinearSolver< LinearSolverRegistrars > > | get_clone () const override |
| |
|
void | initialize () |
| |
|
const Convergence::Criteria & | convergence_criteria () const |
| |
|
::Verbosity | verbosity () const |
| |
|
size_t | restart () const |
| |
|
void | pup (PUP::er &p) override |
| |
|
template<typename LinearOperator , typename SourceType , typename... OperatorArgs, typename IterationCallback = NoIterationCallback> |
| Convergence::HasConverged | solve (gsl::not_null< VarsType * > initial_guess_in_solution_out, const LinearOperator &linear_operator, const SourceType &source, const std::tuple< OperatorArgs... > &operator_args=std::tuple{}, const IterationCallback &iteration_callback=NoIterationCallback{}) const |
| |
|
void | reset () override |
| |
|
template<typename LinearOperator , typename SourceType , typename... OperatorArgs, typename IterationCallback > |
| Convergence::HasConverged | solve (const gsl::not_null< VarsType * > initial_guess_in_solution_out, const LinearOperator &linear_operator, const SourceType &source, const std::tuple< OperatorArgs... > &operator_args, const IterationCallback &iteration_callback) const |
| |
template<typename VarsType, typename Preconditioner = NoPreconditioner, typename LinearSolverRegistrars = tmpl::list<Registrars::Gmres<VarsType>>>
class LinearSolver::Serial::Gmres< VarsType, Preconditioner, LinearSolverRegistrars >
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 [160], 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 [160], 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 [160], 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 'None' to deactivate restarting.
- 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 [160], sections 9.3.2 and 9.4.1 for details. This implementation follows Algorithm 9.6 in [160].
- Improvements:
- Further improvements can potentially be implemented for this algorithm, see e.g. [7].
Example
INFO("Solve a symmetric 2x2 matrix");
blaze::DynamicMatrix<double> matrix{{4., 1.}, {1., 3.}};
const helpers::ApplyMatrix linear_operator{std::move(matrix)};
const blaze::DynamicVector<double> source{1., 2.};
blaze::DynamicVector<double> initial_guess_in_solution_out{2., 1.};
const blaze::DynamicVector<double> expected_solution{0.0909090909090909,
0.6363636363636364};
const Gmres<blaze::DynamicVector<double>> gmres{convergence_criteria,
::Verbosity::Verbose};
CHECK_FALSE(gmres.has_preconditioner());
const auto has_converged = gmres.solve(
make_not_null(&initial_guess_in_solution_out), linear_operator, source,
[&recorded_residuals](
recorded_residuals.push_back(
local_has_converged.residual_magnitude());
});
REQUIRE(has_converged);
CHECK(linear_operator.invocations == 3);
CHECK(has_converged.num_iterations() == 2);
CHECK(has_converged.residual_magnitude() <= 1.e-14);
CHECK(has_converged.initial_residual_magnitude() ==
approx(8.54400374531753));
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]);
}
#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:91
@ AbsoluteResidual
Residual converged below absolute tolerance.
Criteria that determine an iterative algorithm has converged.
Definition: Criteria.hpp:35
Signals convergence or termination of the algorithm.
Definition: HasConverged.hpp:71