% CGtest3.m is a script for comparing {pcg, minres, symmlq}
% on UF sparse matrix number 1360 (PARSEC/Si2, n=769).
% See http://www.cise.ufl.edu/research/sparse/matrices/ (Tim Davis).
% The matrix A is an example from quantum chemistry.
% It is symmetric indefinite with lambda_min = -0.3844.
%
% 30 Mar 2010: MS&E 318 / CME 338 Homework 2.
%----------------------------------------------------------

load Si2.mat
A     = Problem.A;
[m,n] = size(A);
x     = 1./(1:n)';

%----------------------------------------------------------
fprintf('\nPOSITIVE DEFINITE SYSTEM')
fprintf('\n                flag  iter   relres    error\n')
P     = A + 0.4*speye(n);
b     = P*x;
tol   = 1e-12;
maxit = 600;

[xL,flagL,relresL,iterL,resvecL] = symmlq(P,b,tol,maxit);
[xC,flagC,relresC,iterC,resvecC] = pcg   (P,b,tol,maxit);
[xM,flagM,relresM,iterM,resvecM] = minres(P,b,tol,maxit);

errL  = norm(xL-x,inf);
errC  = norm(xC-x,inf);
errM  = norm(xM-x,inf);
fprintf(' SYMMLQ   Px = b%4g %5g %8.1e %8.1e  r\n', flagL,iterL,relresL,errL)
fprintf(' CG       Px = b%4g %5g %8.1e %8.1e  b\n', flagC,iterC,relresC,errC)
fprintf(' MINRES   Px = b%4g %5g %8.1e %8.1e  g\n', flagM,iterM,relresM,errM)

figure(1)
hold off;  plot(log10(resvecL),'ro')
hold  on;  plot(log10(resvecC),'b-')
hold  on;  plot(log10(resvecM),'g.')


%----------------------------------------------------------
fprintf('\nINDEFINITE SYSTEM')
fprintf('\n                flag  iter   relres    error\n')
Afun  = @(x) A*x;       b  = A*x;      % Treat  A  as a function
Afun2 = @(x) A*(A*x);   b2 = A*b;      % Treat A*A as a function

[xL,flagL,relresL,iterL,resvecL] = symmlq(Afun ,b, tol,maxit);
[xC,flagC,relresC,iterC,resvecC] = pcg   (Afun ,b ,tol,maxit);
[xM,flagM,relresM,iterM,resvecM] = minres(Afun ,b ,tol,maxit);
[xN,flagN,relresN,iterN,resvecN] = pcg   (Afun2,b2,tol,maxit);
[xS,flagS,relresS,iterS,resvecS] = lsqr  (A    ,b ,tol,maxit);

errL  = norm(xL-x,inf);
errC  = norm(xC-x,inf);
errM  = norm(xM-x,inf);
errN  = norm(xN-x,inf);
errS  = norm(xS-x,inf);
fprintf(' SYMMLQ   Ax = b%4g %5g %8.1e %8.1e  r\n', flagL,iterL,relresL,errL)
fprintf(' CG       Ax = b%4g %5g %8.1e %8.1e  b\n', flagC,iterC,relresC,errC)
fprintf(' MINRES   Ax = b%4g %5g %8.1e %8.1e  g\n', flagM,iterM,relresM,errM)
fprintf(' CG     A^2x =Ab%4g %5g %8.1e %8.1e  k\n', flagN,iterN,relresN,errN)
fprintf(' LSQR     Ax = b%4g %5g %8.1e %8.1e  m\n', flagS,iterS,relresS,errS)

figure(2)
hold off;  plot(log10(resvecL),'ro')
hold  on;  plot(log10(resvecC),'b-')
hold  on;  plot(log10(resvecM),'g.')
hold  on;  plot(log10(resvecN),'k.')
hold  on;  plot(log10(resvecS),'m-')

%return

%----------------------------------------------------------
% Plot the eigenvalues of A.
%----------------------------------------------------------
lambda = eig(full(A));
figure(3)
hold off;    plot(lambda,'b.')
xlabel('Eigenvalue number')
ylabel('\lambda(A)')
title('Eigenvalues of A')

% Try to show if the eigenvalues are clustered.
figure(4);
hold off;    plot(lambda,10*ones(n,1),'b.')
hold on;

y1    = -5;
yn    = 44.9;
step  = 2;

for y = y1:step:yn
  y2  = y1 + step;
  nlambda = length( find(lambda>y1 & lambda<=y2) );
  bar(y1+0.5*step,nlambda)
  y1  = y2;
end

xlabel('\lambda(A)')
ylabel('No. of \lambda(A)')
title('Distribution of eigenvalues of A')