% cd c:\working@PC\tex\Beijing2009\MATLAB
%
% This is an example of the Steep Decent method for computing the smallest 
% eigenvalue of a definite Hermitian matrix pencil A-lambda B with B positive definite.
%
% RCL 7/7/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); x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;
lambs=[]; maxitn=200; resd=[];
for i=1:maxitn,
    Ax=Axx(x); rho=real(x'*Ax); lambs = [ lambs; rho];
    r=Ax-rho*Bx; rBr=real(r'*(B*r)); nrmr=sqrt(rBr); r=(1/nrmr)*r; resd=[resd; nrmr];
    Ar=Axx(r); rhos=[rho real(r'*Ar)];
    [y, t]=onestep(Ax,Bx,x,r,rhos);
    x=y; 
    Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx); x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;
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 B x_i||');
xlabel('i')
title('Steep Decent'); 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('SD: (\rho_{i+1}-\lambda_1)/(\rho_i-\lambda_1)');
grid
hold off


