% cd c:\working@PC\tex\Beijing2009\MATLAB
%
% This is an example of the Locally Optimal Block Pre-conditioned CG method for computing the smallest
% few eigenvalues of a Hermitian matrix
%
% RCL 6/30/2009

load SiNa
load eigSiNa
A=Problem.A;
[m,n]=size(A);
if m~= n
   disp('A is NOT square');
end

% construct a preconditioner
bds=bdeigest(A); P=A-bds(3)*speye(n);
RR=cholinc(P,'0'); clear P;

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

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

KR=RR\( (RR')\R );
[V, R1]=mgs([X0 KR]); resorth=[resorth; norm(V'*V-eye(2*k),1)];
AV=A*V; AA=V'*(AV);
[Q,D]=schur(AA); D=diag(D); [tmp,iD]=sort(D);
X1=V*Q(:,iD(1:k)); R=A*X1-X1*diag(D(iD(1:k)));
lambs = [lambs ; D(iD(1:k))'];
resd = [resd ; sqrt(sum(R.*conj(R))) ];

maxitn=35;
for i=1:maxitn,
    KR=RR\( (RR')\R );
    [V, R1]=mgs([X1 X0 KR]); resorth=[resorth; norm(V'*V-eye(3*k),1)];
    AV=A*V; AA=V'*(AV);
    [Q,D]=schur(AA); D=diag(D); [tmp,iD]=sort(D);
    X0 = X1;
    X1=V*Q(:,iD(1:k)); R=A*X1-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)*eigA(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
