NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
cheby.h
Go to the documentation of this file.
1 //*****************************************************************
2 // Iterative template routine -- CHEBY
3 //
4 // CHEBY solves the symmetric positive definite linear
5 // system Ax = b using the Preconditioned Chebyshev Method
6 //
7 // CHEBY follows the algorithm described on p. 30 of 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 
23 template < class Matrix, class Vector, class Preconditioner, class Real,
24  class Type >
25 int
26 CHEBY(const Matrix &A, Vector &x, const Vector &b,
27  const Preconditioner &M, int &max_iter, Real &tol,
28  Type eigmin, Type eigmax)
29 {
30  Real resid;
31  Type alpha, beta, c, d;
32  Vector p, q, z;
33 
34  Real normb = norm(b);
35  Vector r = b - A * x;
36 
37  if (normb == 0.0)
38  normb = 1;
39 
40  if ((resid = norm(r) / normb) <= tol) {
41  tol = resid;
42  max_iter = 0;
43  return 0;
44  }
45 
46  c = (eigmax - eigmin) / 2.0;
47  d = (eigmax + eigmin) / 2.0;
48 
49  for (int i = 1; i <= max_iter; i++) {
50  z = M.solve(r); // apply preconditioner
51 
52  if (i == 1) {
53  p = z;
54  alpha = 2.0 / d;
55  } else {
56  beta = c * alpha / 2.0; // calculate new beta
57  beta = beta * beta;
58  alpha = 1.0 / (d - beta); // calculate new alpha
59  p = z + beta * p; // update search direction
60  }
61 
62  q = A * p;
63  x += alpha * p; // update approximation vector
64  r -= alpha * q; // compute residual
65 
66  if ((resid = norm(r) / normb) <= tol) {
67  tol = resid;
68  max_iter = i;
69  return 0; // convergence
70  }
71  }
72 
73  tol = resid;
74  return 1; // no convergence
75 }
A::value_type norm(const A &a)
Euklidean norm of array or expression.
Definition: util.hxx:54
int CHEBY(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol, Type eigmin, Type eigmax)
Definition: cheby.h:26