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