function s7a2

mu = [0.3 0.1 0.03];	% parameter to control cond(A)
n = 10;			% problem size
b = -ones(n,1);		% -rhs

for i=1:length(mu)
  % construct random matrix A, spd
  while (1)
    a = rand(n); A = eye(n)+mu(i)*(a'+a);
    if (min(eig(A)) > 0) break; end	% A is pd
  end
  % solve with cg
  [x,R] = cg(A,b);
  % plot
  subplot(1,length(mu),i);
  plot(0:length(R)-1,log10(R),'.:'); 
  title(['\mu=',num2str(mu(i))]);
  xlabel('iteration k'); ylabel('||r_{k+1}||');
end

function [x,Rsave]=cg(A,b)

% conjugate gradients to solve Ax + b = 0
% in:	A	the (n x n)-Matrix
%	b	the negative rhs
% out:	x	solution
%	Rsave	norm(residuals,2) 

A,

tol = 1e-10; tol2=tol*tol;	% stop if norm(residual)<tol
Rsave = [];			% save norm of resdiduals 

% init
n=length(b); x = zeros(n,1); r = b; p = x; Rp = r'*r; Rsave = Rp;
%  loop
for k = 0:10*n
  if k==0 sigma = 0; else sigma = Rp/R; end
  R     = Rp;
  p     = sigma*p - r;
  z     = A*p;
  rho   = R/(p'*z);
  x     = x + rho*p;
  r     = r + rho*z;
  Rp    = r'*r;
  Rsave = [Rsave sqrt(Rp)];
  if (Rp < tol2) break; end
end

