%% simple minded Conjugate Gradient Demo
%% Fabian Waleffe MIT, May 5, 1995  for 18.330

%% define a symmetric positive definite matrix A and RHS b

%  A=[1,-1;-1,10];  b=[1,2]'; 

% Initialize
  
  x=[0,0]'; bTb = b'*b;               k=0; hist=[k,x'];
  r= b-A*x; rTr = r'*r;     
  y=r; 

%% Conjugate gradient algorithm

 while  rTr > (1.e-4*bTb) 
    
    Ay= A*y;  
    ak= rTr/(y'*Ay);
    x = x + ak*y;       roTro=rTr;
    r = r - ak*Ay;      rTr= r'*r;
%
    bk= rTr/roTro;  
    y = r + bk*y;           % conjugate gradient y instead of r
                             k= k+1; hist=[hist;k,x'];
  end

%% plot iterates

  plot(hist(:,2),hist(:,3)), axis image
