% cd c:\working@PC\tex\Beijing2009\MATLAB
%
% This is an example of the Locally Optimal CG method for computing the smallest
% eigenvalue of a definite Hermitian matrix pencil A-lambda B with B positive definite.
%
% RCL 7/21/2009

addpath Yuji
global A ff hh;
load Yuji/data875/data875

[m,n]=size(A);
if m~= n
   disp('A is NOT square');
end

x=randn(n,1);
Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx); x0=(1/nrmx)*x; Bx=(1/nrmx)*Bx;
Ax=Axx(x0); rho0=real(x0'*Ax); 
r0=Ax-rho0*Bx; rBr=real(r0'*(B*r0)); nrmr=sqrt(rBr); r0n=(1/nrmr)*r0;

Ar0=Axx(r0n);
rhos=[rho0 real(r0n'*Ar0)];
lambs = rho0; resd=norm(r0);  % history
[x1, t]=onestep(Ax,Bx,x0,r0n,rhos);

Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx); x1=(1/nrmx)*x1; Bx=(1/nrmx)*Bx;
Ax=Axx(x1); rho1=real(x1'*Ax); lambs = [ lambs; rho1];
r=Ax-rho1*Bx; 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]=mgsB(B,[x0 x1 r]);
    AV=Axx(V); AA=V'*(AV); AA=0.5*(AA+AA'); 
    [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*B*x1; nrmr=norm(r); r=(1/nrmr)*r; resd=[resd; nrmr];
end

err=lambs-eigAB(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=(eigAB(2)-eigAB(1))/(eigAB(n)-eigAB(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
