%
% This is an example of the CG method for computing the smallest 
% eigenvalue of a Hermitian matrix
%
% RCL 6/28/2009

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

whichbeta = 1; maxitn=75;

x=ones(n,1); x0=(1/norm(x))*x;
Ax=A*x0; rho0=x0'*Ax; lambs = []; resd=[];
r0=Ax-rho0*x0; p0=r0; p0n=(1/norm(p0))*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];
    if whichbeta == 1,
       beta=r1'*r1/(nrmr0*nrmr0); 
    else
       tmp=nrmx1*r1;
       beta=r1'*(tmp-r0)/(nrmr0*nrmr0); 
    end
    p1=r1+beta*p0;
    x0=x1; 
    rho0=rho1; 
    r0=r1; 
    p0=p1; 
    p0n=(1/norm(p1))*p1;
end
resd=[resd; norm(r0)];

err=lambs-eigA(1);
figure(1);
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('Conjugate Gradient'); %grid
hold off

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


