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

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

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

whichbeta = 1; maxitn=75;

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); lambs = []; resd=[];
r0=Ax-rho0*Bx; p0=r0; pBp=real(p0'*(B*p0)); nrmp=sqrt(pBp); p0n=(1/nrmp)*p0;
for i=1:maxitn,
    Ap=Axx(p0n); rhos=[rho0 real(p0n'*Ap)];
    lambs = [lambs ; rho0];
    [x1, t]=onestep(Ax,Bx,x0,p0n,rhos);
    Bx=B*x1; xBx=real(x1'*Bx); nrmx1=sqrt(xBx); x1=(1/nrmx1)*x1; Bx=(1/nrmx1)*Bx; 
    Ax=Axx(x1); rho1=real(x1'*Ax); r1=Ax-rho1*Bx;
    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; 
    pBp=real(p0'*(B*p0)); nrmp=sqrt(pBp); p0n=(1/nrmp)*p0;
end
resd=[resd; norm(r0)];

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('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
