NAME¶
pminres -- conjugate gradient algorithm.
SYNOPSIS¶
template <class Matrix, class Vector, class Preconditioner, class Real>
int pminres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
int &max_iter, Real &tol, odiststream *p_derr=0);
EXAMPLE¶
The simplest call to 'pminres' has the folling form:
size_t max_iter = 100;
double tol = 1e-7;
int status = pminres(a, x, b, EYE, max_iter, tol, &derr);
DESCRIPTION¶
pminres solves the symmetric positive definite linear system Ax=b using
the Conjugate Gradient 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
NOTE¶
pminres follows the algorithm described in "Solution of sparse
indefinite systems of linear equations", C. C. Paige and M. A. Saunders,
SIAM J. Numer. Anal., 12(4), 1975. For more, see
http://www.stanford.edu/group/SOL/software.html and also the PhD
"Iterative methods for singular linear equations and least-squares
problems", S.-C. T. Choi, Stanford University, 2006,
http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf at
page 60. The present implementation style is inspired from
IML++ 1.2
iterative method library,
http://math.nist.gov/iml++.
IMPLEMENTATION¶
template <class Matrix, class Vector, class Preconditioner, class Real, class Size>
int pminres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
Size &max_iter, Real &tol, odiststream *p_derr = 0, std::string label = "minres")
{
Vector b = M.solve(Mb);
Real norm_b = sqrt(fabs(dot(Mb,b)));
if (norm_b == Real(0.)) norm_b = 1;
Vector Mr = Mb - A*x;
Vector z = M.solve(Mr);
Real beta2 = dot(Mr, z);
Real norm_r = sqrt(fabs(beta2));
if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
if (p_derr) (*p_derr) << "[" << label << "] 0 " << norm_r/norm_b << std::endl;
if (beta2 < 0 || norm_r <= tol*norm_b) {
tol = norm_r/norm_b;
max_iter = 0;
dis_warning_macro ("beta2 = " << beta2 << " < 0: stop");
return 0;
}
Real beta = sqrt(beta2);
Real eta = beta;
Vector Mv = Mr/beta;
Vector u = z/beta;
Real c_old = 1.;
Real s_old = 0.;
Real c = 1.;
Real s = 0.;
Vector u_old (x.ownership(), 0.);
Vector Mv_old (x.ownership(), 0.);
Vector w (x.ownership(), 0.);
Vector w_old (x.ownership(), 0.);
Vector w_old2 (x.ownership(), 0.);
for (Size n = 1; n <= max_iter; n++) {
// Lanczos
Mr = A*u;
z = M.solve(Mr);
Real alpha = dot(Mr, u);
Mr = Mr - alpha*Mv - beta*Mv_old;
z = z - alpha*u - beta*u_old;
beta2 = dot(Mr, z);
if (beta2 < 0) {
dis_warning_macro ("pminres: machine precision problem");
tol = norm_r/norm_b;
max_iter = n;
return 2;
}
Real beta_old = beta;
beta = sqrt(beta2);
// QR factorisation
Real c_old2 = c_old;
Real s_old2 = s_old;
c_old = c;
s_old = s;
Real rho0 = c_old*alpha - c_old2*s_old*beta_old;
Real rho2 = s_old*alpha + c_old2*c_old*beta_old;
Real rho1 = sqrt(sqr(rho0) + sqr(beta));
Real rho3 = s_old2 * beta_old;
// Givens rotation
c = rho0 / rho1;
s = beta / rho1;
// update
w_old2 = w_old;
w_old = w;
w = (u - rho2*w_old - rho3*w_old2)/rho1;
x += c*eta*w;
eta = -s*eta;
Mv_old = Mv;
u_old = u;
Mv = Mr/beta;
u = z/beta;
// check residue
norm_r *= s;
if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << norm_r/norm_b << std::endl;
if (norm_r <= tol*norm_b) {
tol = norm_r/norm_b;
max_iter = n;
return 0;
}
}
tol = norm_r/norm_b;
return 1;
}