Loading...
Searching...
No Matches
<tt>gmres</tt>

generalized minimum residual algorithm

# Synopsis

template <class Matrix, class Vector, class Preconditioner,
class SmallMatrix, class SmallVector>
int gmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
SmallMatrix &H, const SmallVector& V, const solver_option& sopt = solver_option())
rheolef::gmres
int gmres(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, SmallMatrix &H, const SmallVector &V, const solver_option &sopt=solver_option())
Definition gmres.h:171
A
Definition diffusion_isotropic.h:25
M
Expr1::memory_type M
Definition vec_expr_v2.h:416

# Example

        solver_option sopt;
        sopt.tol = 1e-7;
        sopt.max_iter = 100;
        size_t m = sopt.krylov_dimension = 6;
        Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> H(m+1,m+1);
        Eigen::Matrix<T,Eigen::Dynamic,1>              V(m);
        int status = gmres (A, x, b, ilut(a), H, V, sopt);

# Description

This function solves the unsymmetric linear system `A*x=b` with the generalized minimum residual algorithm. The `gmres` function follows the algorithm described on p. 20 in

        R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato,
        J. Dongarra, V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst,
        Templates for the solution of linear systems: building blocks for iterative methods,
        SIAM, 1994.

The fourth argument of `gmres` is a preconditionner: here, the `ilut` one is used, for simplicity.

Next, `H` specifies a matrix to hold the coefficients of the upper Hessenberg matrix constructed by the `gmres` iterations, `m` specifies the number of iterations for each restart. We have used here the `eigen` dense matrix and vector types for the `H` and `V` vectors, with sizes related to the Krylov space dimension `m`. Finally, the `solver_option` variable `sopt` transmits the stopping criterion with `sopt.tol` and `sopt.max_iter`.

On return, the `sopt.residue` and `sopt.n_iter` indicate the reached residue and the number of iterations effectively performed. The return `status` is zero when the prescribed tolerance `tol` has been obtained, and non-zero otherwise. Also, the `x` variable contains the approximate solution. See also the `solver_option` for more controls upon the stopping criterion.

# Remarks

`gmres` requires two matrices as input: `A` and `H`. The `A` matrix, which will typically be sparse, corresponds to the matrix involved in the linear system A*x=b. Conversely, the `H` matrix, which will typically be dense, corresponds to the upper Hessenberg matrix that is constructed during the `gmres` iterations. Within `gmres`, `H` is used in a different way than `A`, so its class must supply different functionality. That is, `A` is only accessed though its matrix-vector and transpose-matrix-vector multiplication functions. On the other hand, `gmres` solves a dense upper triangular linear system of equations on `H`. Therefore, the class to which `H` belongs must provide `H(i,j)` accessors.

It is important to remember that we use the convention that indices are 0-based. That is `H(0,0)` is the first component of the matrix. Also, the type of the matrix must be compatible with the type of single vector entry. That is, operations such as H(i,j)*x(j) must be able to be carried out.

# Implementation

This documentation has been generated from file linalg/lib/gmres.h

The present template implementation is inspired from the `IML++ 1.2` iterative method library, http://math.nist.gov/iml++

template <class Matrix, class Vector, class Preconditioner,
class SmallMatrix, class SmallVector>
int gmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
SmallMatrix &H, const SmallVector& V, const solver_option& sopt = solver_option())
{
typedef typename Vector::size_type Size;
typedef typename Vector::float_type Real;
std::string label = (sopt.label != "" ? sopt.label : "gmres");
Size m = sopt.krylov_dimension;
Vector w;
SmallVector s(m+1), cs(m+1), sn(m+1);
Real residue;
Real norm_b = norm(M.solve(b));
Vector r = M.solve(b - A * x);
Real beta = norm(r);
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] # norm_b=" << norm_b << std::endl
<< "[" << label << "] #iteration residue" << std::endl;
if (sopt.absolute_stopping || norm_b == Real(0)) norm_b = 1;
sopt.n_iter = 0;
sopt.residue = norm(r)/norm_b;
if (sopt.residue <= sopt.tol) return 0;
std::vector<Vector> v (m+1);
for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; ) {
v[0] = r/beta;
for (Size i = 0; i < m+1; i++) s(i) = 0; // std::numeric_limits<Float>::max();
s(0) = beta;
for (Size i = 0; i < m && sopt.n_iter <= sopt.max_iter; i++, sopt.n_iter++) {
w = M.solve(A * v[i]);
for (Size k = 0; k <= i; k++) {
H(k, i) = dot(w, v[k]);
w -= H(k, i) * v[k];
}
H(i+1, i) = norm(w);
v[i+1] = w/H(i+1,i);
for (Size k = 0; k < i; k++) {
details::apply_plane_rotation (H(k,i), H(k+1,i), cs(k), sn(k));
}
details::generate_plane_rotation (H(i,i), H(i+1,i), cs(i), sn(i));
details::apply_plane_rotation (H(i,i), H(i+1,i), cs(i), sn(i));
details::apply_plane_rotation (s(i), s(i+1), cs(i), sn(i));
sopt.residue = abs(s(i+1))/norm_b;
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
if (sopt.residue <= sopt.tol) {
details::update (x, i, H, s, v);
return 0;
}
}
details::update (x, m - 1, H, s, v);
r = M.solve(b - A * x);
beta = norm(r);
sopt.residue = beta/norm_b;
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
if (sopt.residue < sopt.tol) return 0;
}
return 1;
}
rheolef::details::update
void update(Vector &x, Size k, const SmallMatrix &h, const SmallVector &s, Vector2 &v)
Definition gmres.h:132
rheolef::details::apply_plane_rotation
void apply_plane_rotation(Real &dx, Real &dy, const Real &cs, const Real &sn)
Definition gmres.h:160
rheolef::details::generate_plane_rotation
void generate_plane_rotation(const Real &dx, const Real &dy, Real &cs, Real &sn)
Definition gmres.h:145
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition vec_expr_v2.h:415
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
rk::beta
Float beta[][pmax+1]
Definition runge_kutta_semiimplicit.icc:60
residue
field residue(Float p, const field &uh)
Definition p_laplacian_post.cc:35
rheolef::Vector::size_type
DATA::size_type size_type
Definition Vector.h:188
template <class SmallMatrix, class SmallVector, class Vector, class Vector2, class Size>
void update (Vector& x, Size k, const SmallMatrix& h, const SmallVector& s, Vector2& v) {
SmallVector y = s;
// back solve:
for (int i = k; i >= 0; i--) {
y(i) /= h(i,i);
for (int j = i - 1; j >= 0; j--)
y(j) -= h(j,i) * y(i);
}
for (Size j = 0; j <= k; j++) {
x += v[j] * y(j);
}
}
template<class Real>
void generate_plane_rotation (const Real& dx, const Real& dy, Real& cs, Real& sn) {
if (dy == Real(0)) {
cs = 1.0;
sn = 0.0;
} else if (abs(dy) > abs(dx)) {
Real temp = dx / dy;
sn = 1.0 / sqrt( 1.0 + temp*temp );
cs = temp * sn;
} else {
Real temp = dy / dx;
cs = 1.0 / sqrt( 1.0 + temp*temp );
sn = temp * cs;
}
}
template<class Real>
void apply_plane_rotation (Real& dx, Real& dy, const Real& cs, const Real& sn) {
Real temp = cs * dx + sn * dy;
dy = -sn * dx + cs * dy;
dx = temp;
}