NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
bicg.h
Go to the documentation of this file.
1 //*****************************************************************
2 // Iterative template routine -- BiCG
3 //
4 // BiCG solves the unsymmetric linear system Ax = b
5 // using the Preconditioned BiConjugate Gradient method
6 //
7 // BiCG follows the algorithm described on p. 22 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 int
25 BiCG(const Matrix &A, Vector &x, const Vector &b,
26  const Preconditioner &M, int &max_iter, Real &tol)
27 {
28  Real resid;
29  Vector rho_1(1), rho_2(1), alpha(1), beta(1);
30  Vector z, ztilde, p, ptilde, q, qtilde;
31 
32  Real normb = norm(b);
33  Vector r = b - A * x;
34  Vector rtilde = r;
35 
36  if (normb == 0.0)
37  normb = 1;
38 
39  if ((resid = norm(r) / normb) <= tol) {
40  tol = resid;
41  max_iter = 0;
42  return 0;
43  }
44 
45  for (int i = 1; i <= max_iter; i++) {
46  z = M.solve(r);
47  ztilde = M.trans_solve(rtilde);
48  rho_1(0) = dot(z, rtilde);
49  if (rho_1(0) == 0) {
50  tol = norm(r) / normb;
51  max_iter = i;
52  return 2;
53  }
54  if (i == 1) {
55  p = z;
56  ptilde = ztilde;
57  } else {
58  beta(0) = rho_1(0) / rho_2(0);
59  p = z + beta(0) * p;
60  ptilde = ztilde + beta(0) * ptilde;
61  }
62  q = A * p;
63  qtilde = A.trans_mult(ptilde);
64  alpha(0) = rho_1(0) / dot(ptilde, q);
65  x += alpha(0) * p;
66  r -= alpha(0) * q;
67  rtilde -= alpha(0) * qtilde;
68 
69  rho_2(0) = rho_1(0);
70  if ((resid = norm(r) / normb) < tol) {
71  tol = resid;
72  max_iter = i;
73  return 0;
74  }
75  }
76 
77  tol = resid;
78  return 1;
79 }
80 
int BiCG(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)
Definition: bicg.h:25
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