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

n = length(b);
u = u0;
r = b - A*u;
z = feval(invCFunc,r); % nove
p = z;                 % nove
normk = norm(z);       % nove
initNorm = normk;

for k=1:n
    Ap = A*p;
    invCAp = feval(invCFunc,Ap); % nove
    pinvCAp = p'*invCAp;         % nove
    alpha = (z'*p)/(pinvCAp);    % nove
    u = u + alpha*p;
    r = b - A*u;
    z = feval(invCFunc,r);       % nove
    normk = norm(z);             % nove
    if normk/initNorm <= prec
        break;
    else
        beta = -(z'*invCAp)/pinvCAp; % nove
        p = z + beta*p;              % nove
    end

end