import numpy as np

A = np.array([[-2,3.01,0.12,-0.11],
              [2.92,-0.17,0.11,0.22],
              [0.66,0.52,3.17,2.11],
              [3.01,0.42,0.27,-0.15]])

b = np.array([0.2,0.1,0,-0.1])

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)