1import numpy as np
2import numpy.linalg as npla
3
4def mgc(A, b, x0,
5 maxiter=100, atol=1.49e-8, rtol=1.49e-8):
6
7 n = b.size
8 x = x0.copy()
9 r = b - A@x
10 d = r.copy()
11 nr = npla.norm(r)
12 print(f'0: {x}, nr = {nr:.1e}')
13 info = -1
14 for k in range(maxiter):
15 rdr = np.dot(r, r)
16 Ad = A@d
17 alpha = rdr/np.dot(d,Ad)
18 x = x + alpha*d
19
20 r = r - alpha*Ad
21 beta = np.dot(r,r)/rdr
22 d = r + beta*d
23
24 nr = npla.norm(r)
25 print(f'{k+1}: {x}, nr = {nr:.1e}')
26
27 if (nr <= max(atol, rtol*npla.norm(b))):
28 info = 0
29 break
30
31 return x, info
32
33
34A = np.array([[2., -1., 0., 0.],
35 [-1., 2., -1., 0.],
36 [0., -1., 2., -1.],
37 [0., 0., -1., 2.]])
38
39b = np.array([-3., 2., 2., -3.])
40
41
42x0 = np.zeros_like(b)
43
44x, info = mgc(A, b, x0)