function [x,relres,resvec]=mycg(A,b,tol)

% simple minded Conjugate Gradient algorithm
%  x=mycg(A,b)  solves Ax=b using CG where A is and N-by-N
%  symmetric positive definite matrix and b is a vector of
%  length N. The iteration uses a "while" loop that might not
%  terminate (but try the example below!).
%
%  [x,relres,resvec]=mycg(A,b,tol) specifies relative tolerance (tol)
%  and returns x, the final relative residual (relres) and the relres
%  history (resvec)
%
%  Fabian Waleffe MIT, May 5, 1995 demo for 18.330
% "functionized" for UW Math 801, Oct 12, 2000.
% 
%  interesting example: (note the obscene tolerance!!)
%
%   N=20; A=randn(N); A=A'*A; b=randn(N,1);
%   [x,relres,resvec]=mycg(A,b,1.e-200); 
%   semilogy(resvec,'-*')
%
%   (Hint: now compare this "iterated residual" to the actual residual
%          norm(b-A*x)/norm(b)

% Initialize
  
if nargin<3 
   tol =1.e-6;
end 

  x=zeros(size(b)); bTb = b'*b;   
  r= b-A*x; rTr = r'*r;     
  p=r; 

  relres=  sqrt(rTr/bTb);
  resvec=[relres];

%% Conjugate gradient algorithm

 while  relres > tol
    
    Ap= A*p;  
    ak= rTr/(p'*Ap);
    x = x + ak*p;       roTro=rTr;
    r = r - ak*Ap;      rTr= r'*r;
%
    bk= rTr/roTro;  
    p = r + bk*p; 

    relres=  sqrt(rTr/bTb);
    resvec=[resvec;relres];
  end



