Homepage > Man Pages > Category > File Formats
Homepage > Man Pages > Name > G

# gmres

## man page of gmres

### gmres: generalized minimum residual method

```NAME
gmres - generalized minimum residual method

SYNOPSIS
template <class Operator, class Vector, class Preconditioner,
class Matrix, class Real, class Int>
int gmres (const Operator &A, Vector &x, const Vector &b,
const Preconditioner &M, Matrix &H, Int m, Int &max_iter, Real &tol);

EXAMPLE
The simplest call to gmres has the folling form:

int m = 6;
dns H(m+1,m+1);
int status = gmres(a, x, b, EYE, H, m, 100, 1e-7);

DESCRIPTION
gmres solves the unsymmetric linear system Ax = b using the generalized
minimum residual method.

The  return  value  indicates  convergence  within   max_iter   (input)
iterations (0), or no convergence within max_iter iterations (1).  Upon
successful return, output arguments have the following values:

x      approximate solution to Ax = b

max_iter
the number of iterations  performed  before  the  tolerance  was
reached

tol    the  residual after the final iteration In addition, M specifies
a preconditioner, 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.

gmres  requires  two  matrices  as input, A and H.  The matrix A, which
will typically be a sparse matrix) corresponds to  the  matrix  in  the
linear  system  Ax=b.   The  matrix  H, which will be typically a dense
matrix,  corresponds  to  the  upper  Hessenberg  matrix  H   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)  operator
for element acess.

NOTE
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  H.  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.

gmres is an iterative template routine.

gmres  follows the algorithm described on p. 20 in @quotation Templates
for the Solution of  Linear  Systems:  Building  Blocks  for  Iterative
Methods,  2nd  Edition, R. Barrett, M. Berry, T. F. Chan, J. Demmel, J.
Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. Van der Vorst,
SIAM, 1994, ftp.netlib.org/templates/templates.ps.  @end quotation

The  present implementation is inspired from IML++ 1.2 iterative method
library, //math.nist.gov/iml++.

IMPLEMENTATION
template < class Matrix, class Vector, class Int >
void
Update(Vector &x, Int k, Matrix &h, Vector &s, Vector v[])
{
Vector y(s);

// Backsolve:
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 (Int j = 0; j <= k; j++)
x += v[j] * y(j);
}

#ifdef TO_CLEAN
template < class Real >
Real
abs(Real x)
{
return (x > Real(0) ? x : -x);
}
#endif // TO_CLEAN

template < class Operator, class Vector, class Preconditioner,
class Matrix, class Real, class Int >
int
gmres(const Operator &A, Vector &x, const Vector &b,
const Preconditioner &M, Matrix &H, const Int &m, Int &max_iter,
Real &tol)
{
Real resid;
Int i, j = 1, k;
Vector s(m+1), cs(m+1), sn(m+1), w;

Real normb = norm(M.solve(b));
Vector r = M.solve(b - A * x);
Real beta = norm(r);

if (normb == Real(0))
normb = 1;

if ((resid = norm(r) / normb) <= tol) {
tol = resid;
max_iter = 0;
return 0;
}

Vector *v = new Vector[m+1];

while (j <= max_iter) {
v[0] = r * (1.0 / beta);    // ??? r / beta
s = 0.0;
s(0) = beta;

for (i = 0; i < m && j <= max_iter; i++, j++) {
w = M.solve(A * v[i]);
for (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 * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)

for (k = 0; k < i; k++)
ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));

GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));

if ((resid = abs(s(i+1)) / normb) < tol) {
Update(x, i, H, s, v);
tol = resid;
max_iter = j;
delete [] v;
return 0;
}
}
Update(x, m - 1, H, s, v);
r = M.solve(b - A * x);
beta = norm(r);
if ((resid = beta / normb) < tol) {
tol = resid;
max_iter = j;
delete [] v;
return 0;
}
}

tol = resid;
delete [] v;
return 1;
}
template<class Real>
void GeneratePlaneRotation(Real &dx, 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 ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
Real temp  =  cs * dx + sn * dy;
dy = -sn * dx + cs * dy;
dx = temp;
}

GMRES(5)
```