function [upper, lower, z, localmins] = dunc_fast(A, B, tol)
% call:  [upper, lower, z, localmins] = dunc_fast(A, B, tol)
% distance to uncontrollability via BFGS, followed by a call to 
% Gu's algorithm to check whether the global minimum was found,
% within a tolerance, and if not, restart BFGS at a lower point. 
% The d to u is 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: relative tolerance (default: 1e-8*sigma_max([A B]))
% output parameters
%   upper: best approximation to (and strict upper bound on) true value
%   lower: lower bound on true value, hopefully with lower = (1-tol)*upper
%   z: satisfies sigma_min([A - zI  B]) = upper
%   localmins: number of BFGS runs

% 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.eud (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;
done = 0;
localmins = 0;
if nargin < 3
   tol = 1e-8*max(sv);
end
while ~done
   localmins = localmins + 1;
   bfgspars.A = A;
   bfgspars.B = B;
   bfgspars.nvar = 2;  % 1 complex variable is 2 real variables
   bfgspars.fgname = 'sigmin';
   bfgsoptions.x0 = [real(z); imag(z)];
   bfgsoptions.gradtol = tol;% not purpose of tol, but must put something
                             % that demands better accuracy for smaller tol
   [x, f, g] = bfgs(bfgspars, bfgsoptions);   % find local minimizer
   if f >= upper             % no improvement - seems to be best can do
      done = 1;              % in this case f should be very small
      disp('give up as bfgs did not find lower point - upper should be small')
   else
      z = x(1) + i*x(2);
      upper = f;
   end
   delta = f*(1 - 0.5*tol);      % try to certify local min is global
   eta = f*tol;   % then eta/2 is 0.5*f*tol, so "delta_2" is f(1-tol)
   [bigger, lambda] = gu_test(delta, eta, A, B);
   if bigger  % distance to uncontrollability <= f*(1 - 0.5*tol)
      % disp('Gu: can do better')
      z = lambda;
      upper = f*(1 - 0.5*tol);   % relative test
   else       % distance to uncontrollability >= f*(1 - tol)
      done = 1;
      lower = f*(1 - tol);
   end
   if localmins > 100 
      disp('dunc_fast: infinite loop? Too many local mins.')
      keyboard
   end
end
