% cd c:\working@PC\tex\Beijing2009\MATLAB
%
% This is an example of Jacobi-Divison method for computing the smallest
% eigenvalue of a symmetric. It runs, too, on nonsymmetric matrices.
%
% RCL 7/23/2009

n=1000;
A=diag(1:n)+.5*diag(ones(n-1,1),1)+.5*diag(ones(n-1,1),-1);
A(1,n)=0.5; A(n,1)=0.5;

Anrm=n+1;
kgmres=20; kJD=100; tol=1.0e-10;
%v0=zeros(n,1); v0(n)=1;
v0=ones(n,1);

V=v0/norm(v0); tau=n+1;
mu=V'*(A*V); r=A*V-mu*V; resd=norm(r)/Anrm; err=resd;
itn=0; z=V; kappa=tau; AA=mu;
while itn <= kJD && err >= tol,
    itn = itn + 1;
    %
    % Use GMRES; comment out to use the direct method
    v=gmresJD(A,z,kappa,r,kgmres);
    %
    % Solve correction equation directly; comment out to use GMRES
    %v=corrJD(A,z,kappa,r);
    %
    v=v-V*(V'*v); v=v/norm(v);
    Av=A*v;
    %AA=[V v]'*A*[V,v];
    AA=[ AA          V'*Av;
        (A'*v)'*V    v'*Av];
    V=[V v];
    [Q,D]=eig(AA); D=diag(D);
    [tmp,ii]=min(abs(D-tau)); mu=D(ii); z=V*Q(:,ii); z=z/norm(z);
    r=A*z-mu*z; err=norm(r)/Anrm; resd=[resd; err];
    kappa=mu;
end
