function [x] = conj_grad(A,b,x)
%  Conjugate gradient steps starting from input point x 
%  Returns: vector x, the point of convergence

acc  = 10^-13;   % accuracy goal
% initial values
r    = b - A*x; % residual
p    = r;
k    = 0;
% iteration till convergence
dot  = r'*r;
while dot > acc,
    pp      = A*p;           % A used only here
    alpha   = dot/(p'*pp);
    x       = x + alpha*p;
    
    r       = r - alpha*pp;
    dot2    = r'*r;
    beta    = dot2/dot;
    p       = r + beta*p;
    dot     = dot2; 
    k       = k + 1;
end 
disp(['converged after ' num2str(k) ' steps']);