function cgdemo1(A)
%CGDEMO1  simple demonstration of SD and CG
%
% A simple demonstration comparing the convergence of
% steepest descents and conjugate gradients
%
% The argument, if given, should be a sparse SPD matrix (e.g.,
% gallery('poisson',50) ).  If the first argument is omitted,
% a previously computed finite element stiffness matrix is used.
% Douglas N. Arnold, 2009-10-01

if nargin == 0
  load StiffnessMatrix.mat
elseif nargin ~= 1
  fprintf('\nUsage: pcgdemo1 or pcgdemo1(A)\n');
end

n=length(A);
b=ones(n,1);
clf

% number of iterations
niter = 249;

%%% Steepest descent
% initial guess
x=zeros(n,1);
% store norms of residuals for steepest descents
ressd=zeros(niter+1,1);

r=b-A*x;
r2 = r'*r;
ressd(1)=sqrt(r2);
for i = 1:niter
  p = A*r;
  lambda = r2/(r'*p);
  x = x + lambda*r;
  r = r - lambda*p;
  r2 = r'*r;
  ressd(i+1)=sqrt(r2);
%  fprintf('iteration %3d   x(1)  %20.18f   rel l2 resid  %5.3g%%\n',i,x(1),ressd(i+1)/ressd(1)*100)
end

%%% Conjugate gradients
% initial guess
x=zeros(n,1);
% store norms of residuals for conjugate gradients
rescg=zeros(niter+1,1);

r = b-A*x;
s = r;
r2 = r'*r;
rescg(1) = sqrt(r2);
for i = 1:niter
  t = A*s;
  s2 = s'*t;
  lambda = r2/s2;
  x = x + lambda*s;
  r = r-lambda*t;
  r2old = r2;
  r2 = r'*r;
  s=r+(r2/r2old)*s;
  rescg(i+1)=sqrt(r2);
%  fprintf('iteration %3d   x(1)  %20.18f   rel l2 resid  %5.3g%%\n',i,x(1),rescg(i+1)/rescg(1)*100)
end

% plot norms of residuals: linear convergence means a straight line
p=semilogy(ressd,'r');
set(p,'color',[.75,.25,.75],'linewidth',2,'marker','o')
hold on
p=semilogy(rescg);
set(p,'color',[.75,.25,.25],'linewidth',2,'marker','o')
grid on
set(gca,'fontsize',20)
xlabel('iterations')
ylabel('norm of residual')
legend('SD','CG','Location','SouthWest')

% compute observed rates of convergence by linear fit of log of norm of residuals
lfit=polyfit(1:niter+1,log(ressd)',1);
ratesd = exp(lfit(1));
lfit=polyfit(1:niter+1,log(rescg)',1);
ratecg = exp(lfit(1));
fprintf('Observed rates of linear convergence\n')
fprintf('SD:         %6.4f\nCG:         %6.5f\n',...
  ratesd,ratecg)
