
function [X,info] = CG(A,B,X,TOL)

% function [X,info] = CG(A,B[,X,TOL])
% solution of AX=B with conjugate gradients
% X [=0]: start vector
% TOL [=eps]: relative error reduction
% info.error: error norm(AX-B)
% info.iterations: number of iterations

if nargin < 3; 
   X = 0*B;
end;
if nargin < 4;
  TOL = eps;
end; 

G = A*X - B;
U = -G; 
gamma = G'*G;
TOL = gamma*TOL^2;

for k=1:length(X);
   V = A*U;
   alpha = gamma/(U'*V);
   X = X + alpha*U;
   G = G + alpha*V;
   beta = gamma;
   gamma = G'*G;
   if gamma < TOL; 
      break;
   end;
   beta = gamma/beta;
   U = beta*U - G;
end;

info.error = sqrt(gamma);
info.iterations = k;

