% cd c:\working@PC\tex\Beijing2009\MATLAB
%
% This is an example of the Pre-conditioned CG method for computing the smallest
% eigenvalue of a Hermitian matrix
%
% RCL 7/3/2009

load SiNa
load eigSiNa
A=Problem.A;
[m,n]=size(A);
if m~= n
   disp('A is NOT square');
end

% construct a preconditioner
bds=bdeigest(A); P=A-bds(3)*speye(n);
RR=cholinc(P,'0'); clear P;

whichbeta = 1; maxitn=32;

x=randn(n,1); x0=(1/norm(x))*x;
Ax=A*x0; rho0=x0'*Ax; r0=Ax-rho0*x0; 
lambs = []; resd=[];
% P\r approx (RR'*RR)\r = RR\(RR'\r);
p0=RR\( (RR')\r0 ); p0n=(1/norm(p0))*p0;
r0Kr0=r0'*p0;
for i=1:maxitn,
    rhos=[rho0 p0n'*(A*p0n)];
    lambs = [lambs ; rho0];
    [x1, t]=onestep0(Ax,x0,p0n,rhos);
    nrmx1=norm(x1); x1=(1/nrmx1)*x1; 
    Ax=A*x1; rho1=x1'*Ax; r1=Ax-rho1*x1;
    nrmr0=norm(r0); resd=[resd; nrmr0];
    % P\r approx (RR'*RR)\r = RR\(RR'\r);
    pwk=RR\( (RR')\r1 ); r1Kr1=r1'*pwk;
    if whichbeta == 1,
       beta=r1Kr1/r0Kr0;
    else
       tmp=nrmx1*r1;
       beta=pwk'*(tmp-r0)/r0Kr0;
    end
    p1=pwk+beta*p0;
    x0=x1; rho0=rho1; r0=r1; p0=p1; p0n=(1/norm(p1))*p1; r0Kr0=r1Kr1;
end
resd=[resd; norm(r0)];

err=lambs-eigA(1);
figure(3);
plot1=semilogy(err,'r-'); hold on
plot2=semilogy(resd,'g--');
legend([plot1 plot2],'\rho_i-\lambda_1', '||A x_i-\rho_i x_i||');
xlabel('i')
title('Preconditioned Conjugate Gradient'); grid
hold off

figure(4)
m=length(err);
plot(err(2:m)./err(1:m-1));
title('CG: (\rho_{i+1}-\lambda_1)/(\rho_i-\lambda_1)');
grid


