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())
int gmres(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, SmallMatrix &H, const SmallVector &V, const solver_option &sopt=solver_option())
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::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 norm_b =
norm(
M.solve(b));
Vector r =
M.solve(b -
A * x);
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;
s(0) = beta;
for (Size i = 0; i < m && sopt.n_iter <= sopt.max_iter; i++, sopt.n_iter++) {
for (Size k = 0; k <= i; k++) {
w -= H(k, i) * v[k];
}
v[i+1] = w/H(i+1,i);
for (Size k = 0; k < i; k++) {
}
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) {
return 0;
}
}
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;
}
void update(Vector &x, Size k, const SmallMatrix &h, const SmallVector &s, Vector2 &v)
void apply_plane_rotation(Real &dx, Real &dy, const Real &cs, const Real &sn)
void generate_plane_rotation(const Real &dx, const Real &dy, Real &cs, Real &sn)
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
field residue(Float p, const field &uh)
DATA::size_type size_type
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;
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>
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>
Real temp = cs * dx + sn * dy;
dy = -sn * dx + cs * dy;
dx = temp;
}