function [u,k] = conjugate_gradients(A,b,u0,prec)

n = length(b);
u = u0;
r = b - A*u;
p = r;
normk = norm(r);
initNorm = normk;

for k=1:n
        
    Ap = A*p;
    pAp = p'*Ap;
    alpha = (r'*p)/(pAp);
    u = u + alpha*p;
    r = b - A*u;
    normk = norm(r);
    
    if normk/initNorm <= prec
        break;
    else
        beta = -(r'*Ap)/pAp;
        p = r + beta*p;
    end

end