NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
cg.h
Go to the documentation of this file.
1 //*****************************************************************
2 // Iterative template routine -- CG
3 //
4 // CG solves the symmetric positive definite linear
5 // system Ax=b using the Conjugate Gradient method.
6 //
7 // CG follows the algorithm described on p. 15 in the
8 // SIAM Templates book.
9 //
10 // The return value indicates convergence within max_iter (input)
11 // iterations (0), or no convergence within max_iter iterations (1).
12 //
13 // Upon successful return, output arguments have the following values:
14 //
15 // x -- approximate solution to Ax = b
16 // max_iter -- the number of iterations performed before the
17 // tolerance was reached
18 // tol -- the residual after the final iteration
19 //
20 //*****************************************************************
21 
22 template < class Matrix, class Vector, class Preconditioner, class Real >
23 int
24 CG(const Matrix &A, Vector &x, const Vector &b,
25  const Preconditioner &M, int &max_iter, Real &tol)
26 {
27  Real resid;
28  Vector p, z, q;
29  Vector alpha(1), beta(1), rho(1), rho_1(1);
30 
31  Real normb = norm(b);
32  Vector r = b - A*x;
33 
34  if (normb == 0.0)
35  normb = 1;
36 
37  if ((resid = norm(r) / normb) <= tol) {
38  tol = resid;
39  max_iter = 0;
40  return 0;
41  }
42 
43  for (int i = 1; i <= max_iter; i++) {
44  z = M.solve(r);
45  rho(0) = dot(r, z);
46 
47  if (i == 1)
48  p = z;
49  else {
50  beta(0) = rho(0) / rho_1(0);
51  p = z + beta(0) * p;
52  }
53 
54  q = A*p;
55  alpha(0) = rho(0) / dot(p, q);
56 
57  x += alpha(0) * p;
58  r -= alpha(0) * q;
59 
60  if ((resid = norm(r) / normb) <= tol) {
61  tol = resid;
62  max_iter = i;
63  return 0;
64  }
65 
66  rho_1(0) = rho(0);
67  }
68 
69  tol = resid;
70  return 1;
71 }
72 
int CG(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)
Definition: cg.h:24
A::value_type dot(const A &a, const B &b)
Dot product of array or expression.
Definition: util.ixx:91
A::value_type norm(const A &a)
Euklidean norm of array or expression.
Definition: util.hxx:54