function [upper, lower, z] = seplambda(A, B, abstol, prtlevel)
% call:  [upper, lower, z] = seplambda(A, B, abstol, prtlevel)
% bisection method to compute seplambda(A,B): the largest epsilon such that the 
% epsilon-pseudospectra of A and B do not overlap, or equivalently
% the global minimum of max{sigmin(A - zI), sigmin(B - zI) over all
% z in the complex plane
% input: 
%   A: square complex matrix
%   B: square complex matrix, not necessarily same dimension as A
%   abstol: compute upper and lower bounds within abstol of each other
%   prtlevel: print level (0 none, 1 print upper and lower bounds as they
%     are determined by bisection, 2 also print info from initial phase)
% output: 
%   upper: strict upper bound on seplambda(A,B)
%   lower: lower bound on seplambda(A,B)
%   z: complex number certifying upper, i.e.,
%       with upper = sigmin(A-zI) = sigmin(B-zI) within a tolerance

% caveat: lower bound may be invalid if an error is made
%   deciding whether an eigenvalue is real or imaginary 
%   (tolerances are used in this decision). To get any guarantees would need to 
%   use specialized eigenvalue software instead of "eig".
%   Upper bound is guaranteed as it is certified by z (within a tolerance).
%   However, since upper bound is obtained by bisection, not optimization,
%   it may not be a good upper bound.

% method:
%  (1) initially: check along lines through eigenvalues to ensure epsilon is small
%       enough that a component of A's epsilon-pseudospectrum can't be strictly
%       inside a component of B's
%  (2) bisection: at each step, do real/imaginary eigenvalue check to find whether
%     eps-pseudospectra of A and B intersect for a given eps
%  Written by Michael Overton, overton@cs.nyu.edu, Jan 2004, for the SVG meeting.
%  Last revision Jan 2005.
%  Reference: An Algorithm to Compute Sep_Lambda, Ming Gu and Michael L. Overton
%  Submitted to SIAM J. Matrix Anal Appl, Jan 2005.
%  http://www.cs.nyu.edu/overton/papers/pdffiles/seplambda.pdf

%  M-files called: seplam_initial, seplam_bisect, linecheck, realimagcheck
if nargin < 4
   prtlevel = 1;
end
if nargin < 3
   abstol = 1e-6;   % absolute tolerance for bisection
end
if size(A,1) ~= size(A,2) | size(B,1) ~= size(B,2)
    error('A and B must be square')
end
maxnorm = max([norm(A,inf); norm(B,inf)]);
if abstol < 1e-16*maxnorm
    abstol = 1e-16*maxnorm;
    if prtlevel > 0
        fprintf('\abstol too small, resetting it to %e', abstol)
    end
end
% initial phase: checking along lines
[f,z] = seplam_initial(A, B, prtlevel);  
if f < inf
   upper = f;
else
   upper = max(min(svd(A)), min(svd(B)));
end
if prtlevel > 0
    fprintf('\nstarting bisection with upper bound equal to %e', upper)
end
lower = 0;
% bisection
[lower, upper, z] = seplam_bisect(A, B, lower, upper, z, abstol, prtlevel);
if prtlevel > 0
    fprintf('\nfinished: lower = %e, upper = %e\n', lower, upper)
end
