Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
//*****************************************************************
// Iterative template routine -- BiCGSTAB
//
// BiCGSTAB solves the unsymmetric linear system Ax = b
// using the Preconditioned BiConjugate Gradient Stabilized method
//
// BiCGSTAB follows the algorithm described on p. 27 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
//
// Slightly modified version of IML++ template. See ReadMe file.
//
// Jesper Andersson
//
#ifndef bicgstab_h
#define bicgstab_h
namespace MISCMATHS {
template < class Matrix, class Vector, class Preconditioner, class Real >
int
BiCGSTAB(const Matrix &A, Vector &x, const Vector &b,
const Preconditioner &M, int &max_iter, Real &tol)
{
Real resid;
Vector rho_1(1), rho_2(1), alpha(1), beta(1), omega(1);
Vector p, phat, s, shat, t, v;
Real normb = b.NormFrobenius();
Vector r = b - A * x;
Vector rtilde = r;
if (normb == 0.0)
normb = 1;
if ((resid = r.NormFrobenius() / normb) <= tol) {
tol = resid;
max_iter = 0;
return 0;
}
for (int i = 1; i <= max_iter; i++) {
rho_1(1) = DotProduct(rtilde, r);
if (rho_1(1) == 0) {
tol = r.NormFrobenius() / normb;
return 2;
}
if (i == 1)
p = r;
else {
beta(1) = (rho_1(1)/rho_2(1)) * (alpha(1)/omega(1));
p = r + beta(1) * (p - omega(1) * v);
}
phat = M.solve(p);
v = A * phat;
alpha(1) = rho_1(1) / DotProduct(rtilde, v);
s = r - alpha(1) * v;
if ((resid = s.NormFrobenius()/normb) < tol) {
x += alpha(1) * phat;
tol = resid;
return 0;
}
shat = M.solve(s);
t = A * shat;
omega = DotProduct(t,s) / DotProduct(t,t);
x += alpha(1) * phat + omega(1) * shat;
r = s - omega(1) * t;
rho_2(1) = rho_1(1);
if ((resid = r.NormFrobenius() / normb) < tol) {
tol = resid;
max_iter = i;
return 0;
}
if (omega(1) == 0) {
tol = r.NormFrobenius() / normb;
return 3;
}
}
tol = resid;
return 1;
}
} // End namespace MISCMATHS
#endif // End #ifndef bicgstab_h