import numpy as np

A = np.array([[1.12,0.13,0.14,0.15],
              [0.16,1.17,0.18,0.19],
              [0.20,0.21,1.22,0.23],
              [0.24,0.25,0.26,1.27]])

b = np.array([0.05,-0.05,0.15,-0.15])

def cg(A, b, tol):
    it=0; x = 0;
    r = np.copy(b); r_prev = np.copy(b)
    rho = np.dot(r,r)
    p = np.copy(r)
    while np.sqrt(rho) > tol * np.sqrt(np.dot(b,b)):
        it += 1
        if it == 1:
            p = r
        else:
            beta = np.dot(r,r)/np.dot(r_prev,r_prev)
            p = r + beta*p
        w = np.dot(A, p)
        alpha = np.dot(r,r)/np.dot(p, w)
        x = x + alpha*p
        r_prev = r
        r = r - alpha*w
        rho = np.dot(r,r)
    return x, it

for i in range(2,8):
    x, it = cg(A,b,10**-i)
    print("eps = ", 10**-i, ", x = ", x, ", iteration = ", it)