%------- Steepest Descent ---------------------
%
% Define matrix and solution
%
ndim = 100;
R    = randn(ndim);  % Random matrix
%
npot = 0.2;            % Steers the eigenvalue distribution
A    = (R'*R)^npot;    % A matrix, positive definite
lamb = eig(A);         % Eigenvalues
kappa= max(lamb)/min(lamb) % The condition number
b    = rand(ndim,1);   % The b-vector
xsol = -A\b;           % The solution
%
Norm2 = sqrt(xsol'*xsol);   % For normalization
NormA = sqrt(xsol'*A*xsol); % For normalization  
%
% Initialization of the SD iteration
%
x = zeros(size(b)); g = b; 
Niterations = 100;
for loop = 1 : Niterations
   Ag   = A*g;
   alfa = (g'*g)./(g'*Ag);
   x    = x - alfa*g;
   g    = g - alfa*Ag; %  g = b+A*x;
   err2(loop) = sqrt((x-xsol)'*(x-xsol))/Norm2;
   errA(loop) = sqrt((x-xsol)'*A*(x-xsol))/NormA;
   errE(loop) = ((kappa-1)/(kappa+1)).^loop;
end;
%
% Display results
subplot(5,2,5); 
plot(lamb,zeros(size(lamb)),'.'); 
axis([0 ceil(max(lamb)) -.5 .5]);
Tittel = ['Eigenvalues, ndim=' num2str(ndim)]; 
title(Tittel);
%
subplot(1,2,2); 
semilogy(1:Niterations,err2,'b',1:Niterations,errA,'r',...
         1:Niterations,errE,'k--' );
axis([0 ndim 1.0e-16 1])
legend( '2-norm' , 'A-norm', 'Error estimate' );
xlabel('Iteration number');
ylabel('||x_{Niter}-x*|| / ||x*||')
Tittel = ['npot= ' num2str(npot) '  \kappa=',num2str(kappa)];
title(Tittel);
%
%------------- Conjugate Gradient Standard Example -------------------
%
% Standard CG  
clf; 
clear;
%
ndim = 200; % dimension
nsim = 100; % number of iterations.
%
R = randn(ndim);
npot = 0.2; 
A = (R'*R)^npot;      % A-matrix
xsol = rand(ndim,1);  % x*
b    = A*xsol;        % RHS
kappa= max(eig(A))/min(eig(A)); % Condition number
%
errBound= (sqrt(kappa)-1)/(sqrt(kappa)+1);
Norm2 = sqrt(xsol'*xsol); 
NormA = sqrt(xsol'*A*xsol);
%
x = zeros(size(b)); 
g = A*x-b ; 
p = -g; 
% Iteration
for loop = 1:nsim
  Ap = A*p; % Only one matrix-vector product!
  alfa = -(p'*g)./(p'*Ap);
  x = x + alfa*p;
  g = g + alfa*Ap; % g = A*x-b;
  beta = (g'*Ap)./(p'*Ap);
  p = -g + beta*p;
  err2(loop) = sqrt((x-xsol)'*(x-xsol))/Norm2;
  errA(loop) = sqrt((x-xsol)'*A*(x-xsol))/NormA;
  errB(loop) = 2*errBound^loop;
end; 
%
semilogy(1:nsim, err2,1:nsim,errA,'g',1:nsim,errB,'r',...
        'Linewidth', 2);
set(gca,'FontSize' ,14)    
axis([0 nsim  1e-10 10])
legend( '2-norm' , 'A-norm', 'Err. bound' );
xlabel('Iteration number'); ylabel('Error')
Tittel = ['Size: ' num2str(ndim) ',   \kappa: ',...
                   num2str(kappa,'%8.2e\n')];
title(Tittel);

%
%------------- Conjugate Gradient Extended Matrix -------------------
%
clf; 
clear;
%
ndim = 200; 
nsim = 100;
%
R    = randn(ndim);
npot = .1; 
A0 = (R'*R)^npot;
L=randn(ndim,12);
mu =100; % Large eigenvalues
A  = A0 + mu*L*L'; 
kappa= max(eig(A))/min(eig(A));
errBound= 2*(sqrt(kappa)-1)/(sqrt(kappa)+1);
%
xsol = rand(ndim,1);
b    = A*xsol;%
%
Norm2 = sqrt(xsol'*xsol); 
NormA = sqrt(xsol'*A*xsol);
%
x = zeros(size(b)); 
g = A*x-b ; 
p = -g; 
for loop = 1:nsim
  Ap = A*p; % Only one matrix-vector product!
  alfa = -(p'*g)./(p'*Ap);
  x = x + alfa*p;
  g = g + alfa*Ap; % g = A*x-b;
  beta = (g'*Ap)./ (p'*Ap);
  p = -g + beta*p;
  err2(loop) = sqrt((x-xsol)'*(x-xsol))/Norm2;
  errA(loop) = sqrt((x-xsol)'*A*(x-xsol))/NormA;
  errB(loop) = 2*(errBound^loop);
end; 
%
semilogy(1:nsim, err2,1:nsim,errA,'g',1:nsim,errB,'r');
axis([0 nsim  1e-16 10])
legend( '2-norm' , 'A-norm','Error Bound');
xlabel('Iteration number'); ylabel('Error')
Tittel = ['Size: ' num2str(ndim) ',   \kappa: ',num2str(kappa)];
title(Tittel);