function t = getTmax(R,S,H)
% Computes the largest t for which 
%     H'*(R+t*S)*H < 0
% where R,S symmetric and S>=0. Assumes H=I when omitted.

%   Copyright 1986-2017 The MathWorks, Inc.

% Absorb H and set zero tolerances. Tolerances set based on expected 
% roundoff level in H'*R*H and H'*S*H (needed for insensitivity to 
% orthogonal transformation (R,S,H) -> (U'*R*U,U'*S*U,U'*H))
rtol = 1e2*eps;
tolR = rtol * norm(R,1);
tolS = rtol * norm(S,1);
if nargin>2
   nH2 = norm(H,1)^2;
   tolR = tolR * nH2;
   tolS = tolS * nH2;
   R = H'*R*H;
   S = H'*S*H;
end
n = size(R,1);
if ~(allfinite(R) && allfinite(S))
   t = NaN;  return
elseif n==0
   t = Inf;  return
end

% Compute null space U2 of S
[U,T] = schur(S+S');
[sv,is] = sort(abs(diag(T)),'descend');
m = sum(sv>tolS); % rank of S
if m==n
   % S full rank
   t = min(real(eig(-R,S)));
else
   % Reduce to [R11-t*S11 R12;R12' R22]<0
   U1 = U(:,is(1:m));
   U2 = U(:,is(m+1:n));
   RU2 = R*U2;
   R22 = U2'*RU2;
   
   % Check R22<0
   if any(eig(R22+R22')>=-tolR)
      t = -Inf;
   elseif m==0
      % R<0, S=0
      t = Inf;
   else
      % Row compress RU2 = [R12;R22] -> [0;*]
      [Q,~,~] = qr(RU2,'vector');
      Q = Q(:,n-m+1:n);
      t = min(real(eig(-Q'*R*U1,Q'*S*U1)));
   end
end