function z = realimagcheck(A, B, eps, all, eigtol)
% Seplambda test: look for a z such that
% sigma_min(A - zI) = sigma_min(B - zI) = eps
% If there isn't one, return nan.
% all: 1 if all z's are wanted (may return many duplicates), 0 if one is enough
% eigtol: tolerance used for real, imaginary parts zero (could cause trouble)
if nargin < 4
all = 0;
end
if nargin < 5
eigtol = 1e-8; % very arbitrary
end
nA = length(A);
nB = length(B);
IA = eye(nA);
IB = eye(nB);
OA = zeros(nA);
OB = zeros(nB);
% we are looking for z = alpha + i*beta for which
% [ 0 A - (alpha + i*beta)I
% A' - (alpha - i beta)I 0 ] has eigenvalue eps, i.e.
% [ -eps A - (alpha + i*beta)I
% A' - (alpha - i beta)I -eps ] has determinant 0, i.e.
% AA - alpha JJA has eigenvalue i*beta, where
AA = [ A -eps*IA % this is Hamiltonian
eps*IA -A' ];
IIA = eye(2*nA);
IIB = eye(2*nB);
JJA = [ IA OA
OA -IA ];
% and likewise for B, i.e. BB - alpha JJB has eigenvalue i*beta, where
BB = [ B -eps*IB % this is Hamiltonian
eps*IB -B' ];
JJB = [ IB OB
OB -IB ];
% Thus we are looking for alpha for which AA and BB have the SAME eigenvalue
% i*beta, i.e. det(K - alpha L) = 0, which is a generalized eigenvalue
% problem with half the eigenvalues being infinite (because of cancellation
% in L). This pencil is known to be regular (not singular) as long as
% A and B have no common eigenvalue. See Section 3 of Gu/Overton.
% The finite eigenvalues have skew Hamiltonian symmetry around the
% real axis (i.e., Hamiltonian symmetry if multiplied by i). Should worry
% about rounding; ideally would use a special skew-Hamiltonian code.
%
% Previously we had the following, for which eig(K,L) are exactly the
% same, but then the notation is not consistent with the paper
%%% K = kron(AA,IIB) - kron(IIA,conj(BB'));
%%% L = kron(JJA,IIB) - kron(IIA,JJB');
K = kron(IIB,AA) - kron(conj(BB'),IIA); % transpose, NOT conjugate transpose!!!
L = kron(IIB,JJA) - kron(JJB',IIA);
lambda = eig(K,L); % can't use "eigs" as K, L are indefinite
% Now we have to check ALL the finite REAL eigenvalues alpha, and for each one,
% find out if the common eigenvalues of AA - alpha JJ and BB - alpha JJ
% include an IMAGINARY number i*beta, and furthermore, if so, whether the
% singular values of A - z I and B - z I for z = alpha + i*beta that
% are, by construction, equal to eps are actually in both cases the MINIMUM
% singular value - the others are of no interest.
real_eigs = find(abs(imag(lambda)) <= eigtol & abs(real(lambda)) < inf);
z = nan; % default if don't find a candidate z
found = 0;
for j=1:length(real_eigs);
alpha = real(lambda(real_eigs(j))); % remove imaginary rounding
AAshift = [ A - alpha*IA -eps*IA
eps*IA -A' + alpha*IA ];
BBshift = [ B - alpha*IB -eps*IB
eps*IB -B' + alpha*IB ];
lamA = eig(AAshift);
lamB = eig(BBshift);
% disp('computed lamA, lamB')
imag_eigsA = find(abs(real(lamA)) <= eigtol);
imag_eigsB = find(abs(real(lamB)) <= eigtol);
if ~isempty(imag_eigsA) & ~isempty(imag_eigsB)
imag_valsA = imag(lamA(imag_eigsA)); % remove real rounding
imag_valsB = imag(lamB(imag_eigsB)); % remove real rounding
diff = imag_valsA*ones(1, length(imag_valsB)) - ...
ones(length(imag_valsA), 1)*imag_valsB';
[indA,indB] = find(abs(diff) < eigtol);
if ~isempty(indA) % go through all matching pairs
for iA = indA' % index vector must be row vector
for iB = indB'
beta = (imag_valsA(iA) + imag_valsB(iB))/2; % same to rounding
zcandidate = alpha + beta*i;
% found a common eigenvalue
svalsA = svd(A - zcandidate*IA); % singular value check
svalsB = svd(B - zcandidate*IB);
if max([abs(min(svalsA) - eps) abs(min(svalsB) - eps)]) <= eigtol
found = found + 1;
z(found) = zcandidate;
if ~all
return; % one z is enough
end
else
% disp('min singular value not close to eps')
end
end
end
% disp('after double loop')
end
end
end