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