function [s,z] = linecheck(A, B, z0, eigtol)
% look for ALL points z on the line parallel to the real axis passing through
% the point z0 such that sigma_min(A - zI) = sigma_min(B - zI)
% i.e. all real numbers t such that
% sigma_min(A - (z0 + t)I)  = sigma_min(B - (z0 + t)I) = s for some s.
% Return all such s and z.
% If there aren't any such points, return inf for s and nan for z

if nargin < 4
   eigtol = 1e-8;
end
nA = length(A); 
nB = length(B); 
IA = eye(nA);
IB = eye(nB);
OA = zeros(nA);
OB = zeros(nB);
% we are looking for z = z0 + t for which
% [      0                 A - (z0 + t)I
%    A' - (z0' + t)I              0            ]   
% and
% [      0                 B - (z0 + t)I
%    B' - (z0' + t)I              0            ]   
% have same eigenvalue eps, i.e. AA - t JJA and BB - t JJB have same eigenvalue eps,
% where
AA = [   OA        A - z0*IA
      A' - z0'*IA     OA     ];   
BB = [   OB        B - z0*IB
      B' - z0'*IB     OB     ];   
JJA = [   OA   IA
          IA    OA ];   
JJB = [   OB   IB
          IB    OB ];   
% i.e. det(K - tL) = 0 for some real t, where

IIA = eye(2*nA);
IIB = eye(2*nB);
K = kron(AA,IIB) - kron(IIA,conj(BB'));  % transpose, NOT conjugate transpose!!!
L = kron(JJA,IIB) - kron(IIA,JJB');
% Compute eigenvalues of pencil, but could potentially be in trouble if pencil is 
% singular, i.e. ALL REAL t are eigenvalues.  See Section 3 of Gu/Overton. There it
% is explained how to use an additional parameter theta to ensure that the
% pencil is regular, but in the presence of rounding this does not seem necessary.
%
% Half the eigenvalues are infinite.  The finite ones have skew Hamiltonian symmetry
% around the real axis.  Should worry about rounding - ideally would use
% special skew-Hamiltonian code.
lambda = eig(K,L);        % can't use "eigs" as K, L are indefinite
% Now we have to check ALL the finite REAL eigenvalues t, and for each one, 
% check whether the least singular value of A - z I and B - z I for z = z0 + t 
% are the same - they might not be because (a) the common eigenvalue might not
% be real and (b) even if it is, it might not correspond to the LEAST singular values.

real_eigs = find(abs(imag(lambda)) <= eigtol & abs(real(lambda)) < inf);
z = nan;       % default if don't find a candidate z
s = inf;       % because want to take min of several such s's in eiglines
found = 0;
for j=1:length(real_eigs);
   zcandidate = z0 + lambda(real_eigs(j));
   sA = min(svd(A - zcandidate*IA));
   sB = min(svd(B - zcandidate*IB));
   if abs(sA - sB) < eigtol
      found = found + 1;
      z(found) = zcandidate;
      s(found) = (sA + sB)/2;  % sA and sB may not be quite the same
   end
end
