function [mu,x,k,resd]=eJD(A,U,tau,Anrm,v0,tol)
%
% Extended JD method, given U whose columns give a partial Schur decomposition
%
%      A U = U T
%
% input
%
%         A      (n-by-n)
%         U      (n-by-p)
%         tau    (scalar) near which an eigenvalue to be found
%         Anrm  scalar approx ||A||
%         v0     (n-by-1) initial approximation for the next Schur vector
%         tol    scalar   tolerance
%
% Output
%
%         mu     scalar   computed eigenvalue
%         x      (n-by-1) computed Schur vector
%         k      (integer) number of JD step taken
%         resd   (k-by-1)  residual history
%
%    RCL  7/27/2009

% parameters to be set
kgmres=20; kJD=100;

[n,p]=size(U); 

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