GNU.WIKI: The GNU/Linux Knowledge Base

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 = &derr;
std::string label = "pgmres";
if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
Vector w;
SmallVector s(m+1), cs(m+1), sn(m+1);
Size i;
Size j = 1;
Size k;
Real resid;
Real normb = norm(M.solve(b));
Vector r = M.solve(b - A * x);
Real beta = norm(r);
if (normb == Real(0)) {
normb = 1;
}
resid = norm(r) / normb;
if (resid  <= 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));
resid = abs(s(i+1)) / normb;
if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
if (resid  < 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);
resid = beta / normb;
if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
if (resid < tol) {
tol = resid;
max_iter = j;
delete [] v;
return 0;
}
}
tol = resid;
delete [] v;
return 1;
}```

All copyrights belong to their respective owners. Other content (c) 2014-2018, GNU.WIKI. Please report site errors to webmaster@gnu.wiki.
Page load time: 0.128 seconds. Last modified: November 04 2018 12:49:43.