function [x,k,res,negCurv] = cg(A,b,optTol,maxIter)
% Linear Conjugate Gradient, where products with A are given by a function

x = zeros(size(b));
r = -b;

rr = r'*r;
p = -r;
k = 0;

res = norm(r);
while res > optTol & k < maxIter
    % Compute Matrix-vector product
    Ap = A(p);
    pAp = p'*Ap;

    % Conjugate Gradient
    alpha = rr/(pAp);
    x = x + alpha*p;
    r = r + alpha*Ap;    
    rr_new = r'*r;
    beta = rr_new/rr;
    p = -r + beta*p;
    k = k + 1;
    rr = rr_new;
    res = norm(r);
end