Provided by: librheolef-dev_6.7-6_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(Real(1.),norm_b)) {
tol = residue/norm_b;
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));
residue = abs(s(i+1));
if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << residue/norm_b << std::endl;
if (residue < tol*max(Real(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(Real(1.),norm_b)) {
tol = residue/norm_b;
max_iter = j;
delete [] v;
return 0;
}
}
tol = residue/norm_b;
delete [] v;
return 1;
}
rheolef-6.7 rheolef-6.7 pgmres(4rheolef)