cg

def cg(A, b, x0=None, maxiter=None, tol=1e-5):
    x = x0 or np.zeros_like(b)
    k = 0
    kmax = maxiter or len(b)

    r = b - A @ x
    p = r

    d_o = np.dot(r, r)
    d_n = d_o
    while k < kmax and d_n > tol ** 2 * d_o:
        a = float(d_n / (p.T @ (A @ p)))
        x += a * p
        r -= a * (A @ p)
        d_o = d_n
        d_n = np.dot(r, r)
        beta = float(d_n / d_o)
        p = r + beta * p
        k += 1
    return x