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
|
def cg(A, b, x, imax=50, eps=1e-6):
"""
A, b, x : matrix, rhs, solution
imax : maximum allowed iterations
eps : tolerance for convergence
"""
# allocate work vectors
r = b.duplicate()
d = b.duplicate()
q = b.duplicate()
# initialization
i = 0
A.mult(x, r)
r.aypx(-1, b)
r.copy(d)
delta_0 = r.dot(r)
delta = delta_0
# enter iteration loop
while i < imax and \
delta > delta_0 * eps**2:
A.mult(d, q)
alpha = delta / d.dot(q)
x.axpy(+alpha, d)
r.axpy(-alpha, q)
delta_old = delta
delta = r.dot(r)
beta = delta / delta_old
d.aypx(beta, r)
i = i + 1
return i, delta**0.5
|
