NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
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: