function [yes, z] = gu_test(delta, eta, A, B)
% Gu's test on rho(A,B), the distance to uncontrollability for (A, B)   
% ref: SIAM J Matrix Anal Appl 21, pp 989-1003 (2000)
% If the test returns yes = 1 we conclude that
%            delta >= rho(A,B), and we have certifying z such that
%                  sigma_min([A - zI   B]) = delta
% If the test returns yes = 0 we conclude that
%            rho(A,B) > delta - eta/2.  (and we have z = NaN)
% The bigger eta is, the more both inequalities could be true.
% If lower and upper bounds L and U are already known on rho(A,B), then
% by choosing delta = (L + 2U)/3 and eta = 2(U-L)/3, so that
% delta - eta/2 = (2L + U)/3.  Then, if the result is yes=1, the
% new bounds are [L, (L + 2U)/3] and if the result is yes=0, the
% new bounds are [(2L + U)/3, U].  Either way, we can reduce the interval
% length by a factor of 2/3.  This is implemented in dunc_trisect.m
% Compare with with gu_test1.m, where eta is not an input parameter but
% is chosen very small, depending on a tolerance, so as to achieve bisection
% up to a tolerance; this leads to numerical trouble.
% The lines containing NOT in the comments implement the actual published
% version of Gu's test, with sign errors.
% Written by Michael Overton, overton@cs.nyu.edu (last updated Nov 2003)

if(delta < 0)
   disp('called with delta negative: trivially returns yes=0')
   yes = 0;
   z = nan;
   return
elseif(delta < 1e-15)
   disp('delta is tiny, quitting Gu test')
   yes = 0;
   z = nan;
   return
end
if(eta <= 0)
   error('eta is not allowed to be negative or zero')
end
n = size(A,1);
m = size(B,2);
if n ~= size(A,2) | m > n
   error('A must be square and B should be tall and skinny')
end
Baug = [B               % n+m by m
        -delta*eye(m)];  
[Q,R] = qr(Baug);       % Q is n+m by n+m
Q12 = Q(1:n, m+1:m+n);  % Q12 is n by n (square)
Q22 = Q(n+1:n+m, m+1:m+n);  % Q22 is m by n (short and fat)
Zn2 = zeros(n^2);
In = eye(n);
scriptB = [Zn2 Zn2 -kron(Q12,In) Zn2            % 2n^2 by 4n^2
           Zn2 Zn2      Zn2      kron(In,Q12)];
hatB = B*Q22 - delta*Q12;   % n by n
In2 = eye(n^2);
Ashift = A - eta*In;     % eta is real
scriptA = ...            % 2n^2 by 4n^2
[ delta*In2 -kron(Q12,hatB)  -kron(Q12,A)-kron(Ashift'*Q12,In)   Zn2
 -delta*In2  kron(hatB,Q12)  Zn2   kron(Ashift,Q12)+kron(In,A'*Q12)];
%%NOT [ delta*In2 -kron(Q12,hatB)  kron(Q12,A)-kron(Ashift'*Q12,In)   Zn2
%%NOT -delta*In2  kron(hatB,Q12)  Zn2  kron(Ashift,Q12)-kron(In,A'*Q12)];
scriptH = ...            % 2n^2 by 4n^2
[ kron(In,A)-kron(Ashift,In) Zn2  -kron(hatB,In)  kron(In,hatB)
Zn2 kron(Ashift'*Q12,Q12)-kron(Q12,A'*Q12) delta*kron(Q12,In) -delta*kron(In,Q12)];
[scriptR, scriptQ] = rq(scriptH);  % scriptQ is 4n^2 by 4n^2
two_nsq = 2*n^2;
indx2 = two_nsq + 1 : 2*two_nsq;
scriptQ2 = scriptQ(indx2,:);   
scriptQ22 = scriptQ(indx2,indx2);  % 2n^2 by 2n^2
scriptA1 = scriptA*scriptQ2';
D = [-kron(Q12,In)  Zn2
     Zn2      kron(In,Q12)];
scriptB1 = D*scriptQ22';
lambda = eig(scriptA1, scriptB1);
real_eigs = find(abs(imag(lambda)) <= 1e-6);  % may need to worry about rounding
Zn = zeros(n);
for k=1:size(real_eigs)
   alpha = 0.5*lambda(real_eigs(k));   % note the factor of 1/2 in (3.13)
                                       % "real" removes imaginary rounding
   M = [A - alpha*In  hatB
           - delta*In  (A' - alpha*In)*Q12];
   N = [A - (alpha + eta)*In   hatB
           - delta*In  (A' - (alpha + eta)*In)*Q12];
   RHS = [In  Zn
         Zn  -Q12];
   % verify that e-vector satisfies (3.10)
   %% [evecs,lambda] = eig(scriptA1, scriptB1);
   %%lambda = diag(lambda);
   %%gu_verify(evecs(:,real_eigs(k)), scriptQ, A, B, Q22, Q12, alpha, ...
   %%   eta, In, Zn, M, N, RHS, scriptH, scriptA, scriptB, scriptA1, scriptB1)
   if isnaninf(M)| isnaninf(RHS)  % may happen if delta too small
      break
   end
   lam_M = eig(M, RHS);
   lam_N = eig(N, RHS);
   imag_eigs_M = find(abs(real(lam_M)) <= 1e-7);
   imag_eigs_N = find(abs(real(lam_N)) <= 1e-7);
   imag_vals_M = imag(lam_M(imag_eigs_M)); % "imag" removes real rounding
   imag_vals_N = imag(lam_N(imag_eigs_N));
   diff = imag_vals_M*ones(1,length(imag_vals_N)) - ... % matrix of differences
          ones(length(imag_vals_M),1)*imag_vals_N';
   [val,col_ind] = min(min(abs(diff)));      % looking for equality
   if val < 1e-7
      z = alpha + imag_vals_N(col_ind)*i;
      svals = svd([A - z*eye(n), B]);
      if (min(svals) - delta) > 1e-7*delta
         disp('warning: sigma_min not very close to delta?')
         difference = min(svals) - delta
      end
      yes = 1;
      return
   end
end
yes = 0;
z = nan;
return
