Provided by: librheolef-dev_6.6-1build2_amd64 #### NAME

```       pgmres -- generalized minimum residual method

```

#### SYNOPSIS

```          template <class Matrix, class Vector, class Preconditioner,
class SmallMatrix, class SmallVector, class Real, class Size>
int pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
SmallMatrix &H, const SmallVector& dummy,
Size m, Size &max_iter, Real &tol);

```

#### EXAMPLE

```       The simplest call to pgmres has the folling form:

double tol = 1e-7;
size_t max_iter = 100;
size_t m = 6;
boost::numeric::ublas::matrix<double> H(m+1,m+1);
vec<double,sequential> dummy;
int status = pgmres (a, x, b, ic0(a), H, dummy, m, max_iter, tol);

```

#### DESCRIPTION

```       pgmres  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 pgmres iterations, m specifies the number of iterations for each
restart.

pgmres  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 pgmres iterations. Within pgmres, 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,  pgmres  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.

pgmres is an iterative template routine.

pgmres  follows  the  algorithm described on p. 20 in 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.

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

```

#### IMPLEMENTATION

```       template <class SmallMatrix, class Vector, class SmallVector, class Size>
void
Update(Vector &x, Size k, SmallMatrix &h, SmallVector &s, Vector 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 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;
}
template <class Matrix, class Vector, class Preconditioner,
class SmallMatrix, class SmallVector, class Real, class Size>
int
pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
SmallMatrix &H, const SmallVector&, const Size &m, Size &max_iter, Real &tol,
odiststream* p_derr = 0, std::string label = "pgmres")
{
Vector w;
SmallVector s(m+1), cs(m+1), sn(m+1);
Size i;
Size j = 1;
Size k;
Real residue;
Real norm_b = norm(M.solve(b));
Vector r = M.solve(b - A * x);
Real beta = norm(r);
if (p_derr) (*p_derr) << "[" << label << "] # norm_b=" << norm_b << std::endl;
if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
if (norm_b == Real(0)) {
norm_b = 1;
}
residue = norm(r);
if (residue  <= tol*max(1.,norm_b)) {
tol = residue/norm_b;
max_iter = 0;
return 0;
}
Vector *v = new Vector[m+1];
while (j <= max_iter) {
v = 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));
residue = abs(s(i+1));
if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << residue/norm_b << std::endl;
if (residue  < tol*max(1.,norm_b)) {
Update(x, i, H, s, v);
tol = residue/norm_b;
max_iter = j;
delete [] v;
return 0;
}
}
Update(x, m - 1, H, s, v);
r = M.solve(b - A * x);
beta = norm(r);
residue = beta;
if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << residue/norm_b << std::endl;
if (residue < tol*max(1.,norm_b)) {
tol = residue/norm_b;
max_iter = j;
delete [] v;
return 0;
}
}
tol = residue/norm_b;
delete [] v;
return 1;
}
```