function [s, z] = seplam_initial(A, B, prtlevel) % search along lines parallel to the real axis passing through all 2n eigenvalues % of A and B, looking for ALL real positive s and z on these real lines such that % sigma_min(A - zI) = sigma_min(B - zI) = s. Return the smallest such s and % the corresponding z. This ensures that we don't get stuck with a component of % A's eps-pseudospectrum strictly contained inside a component of B's eps-pseudospectrum, % for some eps, thus defeating "realimagcheck". % could make this more efficient % can save factor of 2 by looking along lines joining eigenvalues % can avoid repeat checks along same line, e.g. if eigenvalues are real eA = eig(A); eB = eig(B); eigtol = 1e-8; for j=1:length(A) [s,z] = linecheck(A, B, eA(j), eigtol); [val,indx] = min(s); sminA(j) = val; zminA(j) = z(indx); if prtlevel > 1 fprintf('\nfound initial upper bound %e on line %d for A', val, j) end end for j=1:length(B) [s,z] = linecheck(A, B, eB(j), eigtol); [val,indx] = min(s); sminB(j) = val; zminB(j) = z(indx); if prtlevel > 1 fprintf('\nfound initial upper bound %e on line %d for B', val, j) end end [valA,indxA] = min(sminA); [valB,indxB] = min(sminB); if valA < valB s = valA; z = zminA(indxA); else s = valB; z = zminB(indxB); end % possible that s is inf and z is nan, if no s and z were found on % any of these lines