NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
Functions
cg.h File Reference

Go to the source code of this file.

Functions

template<class Matrix , class Vector , class Preconditioner , class Real >
int CG (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)
 

Function Documentation

template<class Matrix , class Vector , class Preconditioner , class Real >
int CG ( const Matrix &  A,
Vector &  x,
const Vector &  b,
const Preconditioner &  M,
int &  max_iter,
Real &  tol 
)

Definition at line 24 of file cg.h.

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 }
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

+ Here is the call graph for this function: