function [ x,r,i ] = conj_grad_mod(A,x,b,tol)
% x and b must be column vectors

if nargin<4
    tol = 1e-9;
end

n = length(b);
r = b - A*x;
s = r;

% r = A*x - b;
% s = -r;

i = 0;
rtr = r'*r;

while norm(r) > tol
%     i = i+1;
%     As = A*s;
%     alpha = (s'*r)/(s'*As);
%     x = x + alpha*s;
%     r = b - A*x;
%     beta = -(r'*As)/(s'*As);
%     s = r + beta*s;

%     i = i+1;
%     As = A*s;
%     alpha = rtr/(s'*As);
%     xn = x + alpha*s;
%     rn = r - alpha*As;
%     rntrn = rn'*rn;
%     beta = rntrn/rtr;
%     sn = rn + beta*s;
%     rtr = rntrn;
%     r = rn;
%     x = xn;
%     s = sn;

    i = i+1;
    As = A*s;
    % alpha = rtr/(s'*As);
    alpha = (s'*r)/(s'*As);    
    x = x + alpha*s;
    r = r - alpha*As;
    rntrn = r'*r;
    beta = rntrn/rtr;
    s = r + beta*s;
    rtr = rntrn;

end

end

