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)