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

k = 5; % number of desired eigenvalues
X0=randn(n,k); [V, R1]=mgsB(B,X0); X0=V;
lambs=[]; resd=[]; % Y=zeros(n,k);

AX=Axx(X0); AA0=X0'*AX; AA0=0.5*(AA0+AA0');
[Q,D]=schur(AA0); D=diag(D); 
X0=X0*Q; R=AX*Q-B*X0*diag(D); 
lambs = D'; resd=sqrt(sum(R.*conj(R)));
resorth=norm(X0'*B*X0-eye(k),1);

[V, R1]=mgsB(B,[X0 R]); resorth=[resorth; norm(V'*B*V-eye(2*k),1)];
AV=Axx(V); AA=V'*(AV); AA=0.5*(AA+AA');
[Q,D]=schur(AA); D=diag(D); [tmp,iD]=sort(D);
X1=V*Q(:,iD(1:k)); R=AV*Q(:,iD(1:k))-B*X1*diag(D(iD(1:k)));
lambs = [lambs ; D(iD(1:k))']; 
resd = [resd ; sqrt(sum(R.*conj(R))) ];

maxitn=70;
for i=1:maxitn,
    [V, R1]=mgsB(B,[X1 X0 R]); resorth=[resorth; norm(V'*B*V-eye(3*k),1)];
    AV=Axx(V); AA=V'*(AV); AA=0.5*(AA+AA');
    [Q,D]=schur(AA); D=diag(D); [tmp,iD]=sort(D);
    X0 = X1;
    X1=V*Q(:,iD(1:k)); R=AV*Q(:,iD(1:k))-B*X1*diag(D(iD(1:k)));
    lambs = [lambs ; D(iD(1:k))']; 
    resd = [resd ; sqrt(sum(R.*conj(R))) ];
end

err=lambs-ones(maxitn+2,1)*eigAB(1:k)';
figure(1);
subplot(1,2,1); 
semilogy(abs(err)); title('eigenvalue errors'); grid
xlabel('i')
subplot(1,2,2); 
semilogy(resd,'-'); title('residual errors'); grid
hold off
