function [upper, lower ,z] = dunc_trisect(A, B, tol)
% call:  [upper, lower, z] = dunc_trisect(A, B, tol)
% distance to uncontrollability via a trisection version of Gu's algorithm, 
% i.e., the minimum value, over all z in complex plane, of
%   sigma_min([A - zI   B])        
% input parameters
%   A: square matrix
%   B: rectangular matrix (n by p with n >= p)
%   tol: absolute tolerance (default: 1e-8*sigma_max([A B]))
% output parameters
%   upper: upper bound on true value
%   lower: lower bound on true value, hopefully with lower >= upper - tol
%   z: satsfies sigma_min([A - zI  B]) = upper

% ref: Ming Gu, SIAM J Matrix Anal Appl 21, pp 989-1003 (2000)
% and: Jim Burke, Adrian Lewis, Michael Overton, SIAM J Matrix Anal Appl, submitted
% Written by Michael Overton, overton@cs.nyu.edu (last updated Nov 2003)

n = size(A,1);
m = size(B,1);
if n ~= size(A,2) | m > n
   error('A must be square and B should be tall and skinny')
end
sv = svd([A B]);
upper = min(sv);
z = 0;
lower = 0;
done = 0;
count = 0;
if nargin < 3
   tol = 1e-8*max(sv);
end
while ~done
   count = count + 1;
   delta = (lower + 2*upper)/3;   % compare with "upper" in dunc_bisect.m
   eta = 2*(upper - lower)/3;     % compare with small tolerance in dunc_bisect.m
   [bigger, lambda] = gu_test(delta, eta, A, B);
   if bigger
      upper = (lower + 2*upper)/3;
      z = lambda;
   else
      % delta <= 2*distance to uncontrollability 
      lower = (2*lower + upper)/3;   % don't update z, which is NaN
   end
   if count > 100 | upper - lower <= tol;
      done = 1;
   end
end
