NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
qmr.h
Go to the documentation of this file.
1 //*****************************************************************
2 // Iterative template routine -- QMR
3 //
4 // QMR.h solves the unsymmetric linear system Ax = b using the
5 // Quasi-Minimal Residual method following the algorithm as described
6 // on p. 24 in the SIAM Templates book.
7 //
8 // -------------------------------------------------------------
9 // return value indicates
10 // ------------ ---------------------
11 // 0 convergence within max_iter iterations
12 // 1 no convergence after max_iter iterations
13 // breakdown in:
14 // 2 rho
15 // 3 beta
16 // 4 gamma
17 // 5 delta
18 // 6 ep
19 // 7 xi
20 // -------------------------------------------------------------
21 //
22 // Upon successful return, output arguments have the following values:
23 //
24 // x -- approximate solution to Ax=b
25 // max_iter -- the number of iterations performed before the
26 // tolerance was reached
27 // tol -- the residual after the final iteration
28 //
29 //*****************************************************************
30 
31 
32 #include <math.h>
33 
34 template < class Matrix, class Vector, class Preconditioner1,
35  class Preconditioner2, class Real >
36 int
37 QMR(const Matrix &A, Vector &x, const Vector &b, const Preconditioner1 &M1,
38  const Preconditioner2 &M2, int &max_iter, Real &tol)
39 {
40  Real resid;
41 
42  Vector rho(1), rho_1(1), xi(1), gamma(1), gamma_1(1), theta(1), theta_1(1);
43  Vector eta(1), delta(1), ep(1), beta(1);
44 
45  Vector r, v_tld, y, w_tld, z;
46  Vector v, w, y_tld, z_tld;
47  Vector p, q, p_tld, d, s;
48 
49  Real normb = norm(b);
50 
51  r = b - A * x;
52 
53  if (normb == 0.0)
54  normb = 1;
55 
56  if ((resid = norm(r) / normb) <= tol) {
57  tol = resid;
58  max_iter = 0;
59  return 0;
60  }
61 
62  v_tld = r;
63  y = M1.solve(v_tld);
64  rho(0) = norm(y);
65 
66  w_tld = r;
67  z = M2.trans_solve(w_tld);
68  xi(0) = norm(z);
69 
70  gamma(0) = 1.0;
71  eta(0) = -1.0;
72  theta(0) = 0.0;
73 
74  for (int i = 1; i <= max_iter; i++) {
75 
76  if (rho(0) == 0.0)
77  return 2; // return on breakdown
78 
79  if (xi(0) == 0.0)
80  return 7; // return on breakdown
81 
82  v = (1. / rho(0)) * v_tld;
83  y = (1. / rho(0)) * y;
84 
85  w = (1. / xi(0)) * w_tld;
86  z = (1. / xi(0)) * z;
87 
88  delta(0) = dot(z, y);
89  if (delta(0) == 0.0)
90  return 5; // return on breakdown
91 
92  y_tld = M2.solve(y); // apply preconditioners
93  z_tld = M1.trans_solve(z);
94 
95  if (i > 1) {
96  p = y_tld - (xi(0) * delta(0) / ep(0)) * p;
97  q = z_tld - (rho(0) * delta(0) / ep(0)) * q;
98  } else {
99  p = y_tld;
100  q = z_tld;
101  }
102 
103  p_tld = A * p;
104  ep(0) = dot(q, p_tld);
105  if (ep(0) == 0.0)
106  return 6; // return on breakdown
107 
108  beta(0) = ep(0) / delta(0);
109  if (beta(0) == 0.0)
110  return 3; // return on breakdown
111 
112  v_tld = p_tld - beta(0) * v;
113  y = M1.solve(v_tld);
114 
115  rho_1(0) = rho(0);
116  rho(0) = norm(y);
117  w_tld = A.trans_mult(q) - beta(0) * w;
118  z = M2.trans_solve(w_tld);
119 
120  xi(0) = norm(z);
121 
122  gamma_1(0) = gamma(0);
123  theta_1(0) = theta(0);
124 
125  theta(0) = rho(0) / (gamma_1(0) * beta(0));
126  gamma(0) = 1.0 / sqrt(1.0 + theta(0) * theta(0));
127 
128  if (gamma(0) == 0.0)
129  return 4; // return on breakdown
130 
131  eta(0) = -eta(0) * rho_1(0) * gamma(0) * gamma(0) /
132  (beta(0) * gamma_1(0) * gamma_1(0));
133 
134  if (i > 1) {
135  d = eta(0) * p + (theta_1(0) * theta_1(0) * gamma(0) * gamma(0)) * d;
136  s = eta(0) * p_tld + (theta_1(0) * theta_1(0) * gamma(0) * gamma(0)) * s;
137  } else {
138  d = eta(0) * p;
139  s = eta(0) * p_tld;
140  }
141 
142  x += d; // update approximation vector
143  r -= s; // compute residual
144 
145  if ((resid = norm(r) / normb) <= tol) {
146  tol = resid;
147  max_iter = i;
148  return 0;
149  }
150  }
151 
152  tol = resid;
153  return 1; // no convergence
154 }
int QMR(const Matrix &A, Vector &x, const Vector &b, const Preconditioner1 &M1, const Preconditioner2 &M2, int &max_iter, Real &tol)
Definition: qmr.h:37
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