% cd c:\working@PC\tex\Beijing2009\MATLAB
%
% This is an example of the Locally Optimal 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

x=randn(n,1); nrmx=norm(x); x0=(1/nrmx)*x;
Ax=A*x0; rho0=x0'*Ax; 
r0=Ax-rho0*x0; nrmr0=norm(r0); r0=(1/nrmr0)*r0; 
rhos=[rho0 r0'*(A*r0)];
lambs = rho0; resd=nrmr0;  % history
[x1, t]=onestep0(Ax,x0,r0,rhos); x1=(1/norm(x1))*x1;
Ax=A*x1; rho1=x1'*Ax; lambs = [ lambs; rho1];
r=Ax-rho1*x1; nrmr=norm(r); r=(1/nrmr)*r; resd=[resd; nrmr];
    
maxitn=70; 
for i=1:maxitn,
    % Project the problem onto span{x0 x1 r}
    [V, R]=mgs([x0 x1 r]);
    AV=A*V; AA=V'*(AV); 
    [Q,D]=schur(AA); D=diag(D);
    [dmin,imin]=min(D); lambs = [ lambs; dmin];
    x0=x1; x1=V*Q(:,imin);
    r=AV*Q(:,imin)-dmin*x1; nrmr=norm(r); r=(1/nrmr)*r; resd=[resd; nrmr];
end

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('Locally Optimal CG'); grid
hold off

figure(2)
m=length(err);
plot(err(2:m)./err(1:m-1)); hold on
xi=(eigA(2)-eigA(1))/(eigA(n)-eigA(1));
rate=( (1-xi)/(1+xi) )^2;
plot([1 maxitn],[rate rate],'r--');
title('Locally Optimal CG: (\rho_{i+1}-\lambda_1)/(\rho_i-\lambda_1)');
grid
hold off
