NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
ir.h
Go to the documentation of this file.
1 //*****************************************************************
2 // Iterative template routine -- Preconditioned Richardson
3 //
4 // IR solves the unsymmetric linear system Ax = b using
5 // Iterative Refinement (preconditioned Richardson iteration).
6 //
7 // The return value indicates convergence within max_iter (input)
8 // iterations (0), or no convergence within max_iter iterations (1).
9 //
10 // Upon successful return, output arguments have the following values:
11 //
12 // x -- approximate solution to Ax = b
13 // max_iter -- the number of iterations performed before the
14 // tolerance was reached
15 // tol -- the residual after the final iteration
16 //
17 //*****************************************************************
18 
19 template < class Matrix, class Vector, class Preconditioner, class Real >
20 int
21 IR(const Matrix &A, Vector &x, const Vector &b,
22  const Preconditioner &M, int &max_iter, Real &tol)
23 {
24  Real resid;
25  Vector z;
26 
27  Real normb = norm(b);
28  Vector r = b - A*x;
29 
30  if (normb == 0.0)
31  normb = 1;
32 
33  if ((resid = norm(r) / normb) <= tol) {
34  tol = resid;
35  max_iter = 0;
36  return 0;
37  }
38 
39  for (int i = 1; i <= max_iter; i++) {
40  z = M.solve(r);
41  x += z;
42  r = b - A * x;
43 
44  if ((resid = norm(r) / normb) <= tol) {
45  tol = resid;
46  max_iter = i;
47  return 0;
48  }
49  }
50 
51  tol = resid;
52  return 1;
53 }
54 
55 
56 
int IR(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)
Definition: ir.h:21
A::value_type norm(const A &a)
Euklidean norm of array or expression.
Definition: util.hxx:54