global A;

CGparams = struct('maxit',1000,'toler',1.0e-6);

n=20;
CGparams.maxit = min(1000, 10*n);
A=hilb(n); x0=zeros(n,1); b = ones(n,1);
[inform, x, residual_norm] = CG(@Aexplicit, x0, b, CGparams);

if inform.status==1
  fprintf(1,' n=%d, convergence in %d iterations\n', n, inform.iter);
  
  % plot the logs of the residuals
  figure(1);
  residual_norm=log10(residual_norm);
  plot(residual_norm,'LineWidth',1.1);
  top=max(residual_norm(:)); bottom=min(residual_norm(:));
  v=[0 length(residual_norm) bottom-.05*(top-bottom) top+.05*(top-bottom)];
  axis(v);
  xlabel('Iteration'); ylabel('log10 residual');
  title('CG on Hilbert matrix for n=20');
else
  fprintf(1,' n=%d, convergence not achieved in %d iterations\n',...
      n, CGparams.maxit);
end


n=100;
CGparams.maxit = min(1000, 5*n);
A=hilb(n); x0=zeros(n,1); b = ones(n,1);
[inform, x, residual_norm] = CG(@Aexplicit, x0, b, CGparams);

if inform.status==1
  fprintf(1,' n=%d, convergence in %d iterations\n', n, inform.iter);

  % plot the logs of the residuals
  figure(2);
  residual_norm=log10(residual_norm);
  plot(residual_norm,'LineWidth',1.1);
  top=max(residual_norm(:)); bottom=min(residual_norm(:));
  v=[0 length(residual_norm) bottom-.05*(top-bottom) top+.05*(top-bottom)];
  axis(v);
  xlabel('Iteration'); ylabel('log10 residual');
  title('CG on Hilbert matrix for n=100');
else
  fprintf(1,' n=%d, convergence not achieved in %d iterations\n',...
      n, CGparams.maxit);
end


n=100;
CGparams.maxit = min(1000, 5*n);
A=hilb(n) + 0.2*eye(n); x0=zeros(n,1); b = ones(n,1);
[inform, x, residual_norm] = CG(@Aexplicit, x0, b, CGparams);

if inform.status==1
  fprintf(1,' n=%d, convergence in %d iterations\n', n, inform.iter);

  % plot the logs of the residuals
  figure(3);
  residual_norm=log10(residual_norm);
  plot(residual_norm,'LineWidth',1.1);
  top=max(residual_norm(:)); bottom=min(residual_norm(:));
  v=[0 length(residual_norm) bottom-.05*(top-bottom) top+.05*(top-bottom)];
  axis(v);
  xlabel('Iteration'); ylabel('log10 residual');
  title('CG on Hilbert matrix plus .2*I for n=100');
else
  fprintf(1,' n=%d, convergence not achieved in %d iterations\n',...
      n, CGparams.maxit);
end