NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
All Classes Namespaces Files Functions Variables Typedefs Friends Macros Pages
Functions
cgs.h File Reference

Go to the source code of this file.

Functions

template<class Matrix , class Vector , class Preconditioner , class Real >
int CGS (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 CGS ( const Matrix &  A,
Vector &  x,
const Vector &  b,
const Preconditioner &  M,
int &  max_iter,
Real &  tol 
)

Definition at line 24 of file cgs.h.

26 {
27  Real resid;
28  Vector rho_1(1), rho_2(1), alpha(1), beta(1);
29  Vector p, phat, q, qhat, vhat, u, uhat;
30 
31  Real normb = norm(b);
32  Vector r = b - A*x;
33  Vector rtilde = r;
34 
35  if (normb == 0.0)
36  normb = 1;
37 
38  if ((resid = norm(r) / normb) <= tol) {
39  tol = resid;
40  max_iter = 0;
41  return 0;
42  }
43 
44  for (int i = 1; i <= max_iter; i++) {
45  rho_1(0) = dot(rtilde, r);
46  if (rho_1(0) == 0) {
47  tol = norm(r) / normb;
48  return 2;
49  }
50  if (i == 1) {
51  u = r;
52  p = u;
53  } else {
54  beta(0) = rho_1(0) / rho_2(0);
55  u = r + beta(0) * q;
56  p = u + beta(0) * (q + beta(0) * p);
57  }
58  phat = M.solve(p);
59  vhat = A*phat;
60  alpha(0) = rho_1(0) / dot(rtilde, vhat);
61  q = u - alpha(0) * vhat;
62  uhat = M.solve(u + q);
63  x += alpha(0) * uhat;
64  qhat = A * uhat;
65  r -= alpha(0) * qhat;
66  rho_2(0) = rho_1(0);
67  if ((resid = norm(r) / normb) < tol) {
68  tol = resid;
69  max_iter = i;
70  return 0;
71  }
72  }
73 
74  tol = resid;
75  return 1;
76 }
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: