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