function BalData = hsvd(D,Options)
%HSVD  Computes the Hankel singular values.
%
%   Second input OPTIONS must be of class mor.BalancedTruncationOptions.

%	 Author(s): P. Gahinet, S. Gumussoy
%   Copyright 1986-2023 The MathWorks, Inc.
Ts = abs(D.Ts);
if hasInternalDelay(D)
   error(message('Control:transformation:BALROM14'))
end
[ny,nu] = size(D.d);
RELERR = strcmp(Options.Goal,'relative');
LIMITED = ~(isempty(Options.TimeIntervals) && isempty(Options.FreqIntervals));

% Eliminate structural cancellations and cancellations at infinity
nx0 = size(D.a,1);
if LIMITED
   % Note: Time- or frequency-limited requires explicit model
   [isProper,D] = isproper(sminreal(D),true);
   if ~isProper
      error(message('Control:transformation:BALROM12'))
   end
else
   % Eliminate cancellations but keep E matrix
   [~,D] = isproper(sminreal(D),false);
end
nsnm = nx0-size(D.a,1);  % number of cancelled states
   
% Stable/unstable partition
% Note: DS is always proper (poles at infinity go to DNS)
[Ds,Dns,~,sepInfo] = stabsep(D,Options);

% Compute Cholesky factors of reachability and observability Gramians
% RE: 1) Uses Hammarling's algorithm that solves directly for the
%        Cholesky factor of the Lyapunov equation solution (see LYAPCHOL)
%     2) Skip scaling (rigorous scaling already performed in STABSEP)
a = Ds.a;  b = Ds.b;  c = Ds.c;  d = Ds.d;  e = Ds.e;
if RELERR
   % Relative error control
   if ny>nu
      % Implicitly transpose (BST requires G wide)
      if Ts==0
         Ro = lyapchol(a',c',e','noscale');
      else
         Ro = dlyapchol(a',c',e','noscale');
      end
      [Rr,alpha] = gramBST(a',c',b',d',e',Ts,Ro,Options.Regularization);
   else
      % G square or wide
      if Ts==0
         Rr = lyapchol(a,b,e,'noscale');
      else
         Rr = dlyapchol(a,b,e,'noscale');
      end
      [Ro,alpha] = gramBST(a,b,c,d,e,Ts,Rr,Options.Regularization);
   end
else
   % Absolute error control
   ns = size(a,1);
   if isempty(Options.InputWeight)
      if Ts==0
         Rr = lyapchol(a,b,e,'noscale');
      else
         Rr = dlyapchol(a,b,e,'noscale');
      end
      if LIMITED
         % use time/frequency limited approach for balanced realization
         Rr = gramLimited(a,b,Ts,Options,Rr);
      end
   else
      % Account for input weight (Antoulas, p. 235)
      [ai,bi,ci,di] = ssdata(Options.InputWeight);
      ni = size(ai,1);
      aa = [a b*ci;zeros(ni,ns) ai];   bb = [b*di;bi];
      if isempty(e)
         ee = [];
      else 
         ee = blkdiag(e,eye(ni));
      end
      if Ts==0
         Rr = lyapchol(aa,bb,ee,'noscale');
      else
         Rr = dlyapchol(aa,bb,ee,'noscale');
      end
      Rr = Rr(1:ns,1:ns);
   end
   if isempty(Options.OutputWeight)
      if Ts==0
         Ro = lyapchol(a',c',e','noscale');
      else
         Ro = dlyapchol(a',c',e','noscale');
      end
      if LIMITED
         Ro = gramLimited(a',c',Ts,Options,Ro);
      end
   else
      % Account for output weight
      [ao,bo,co,do] = ssdata(Options.OutputWeight);
      no = size(ao,1);
      aa = [a zeros(ns,no);bo*c ao];   cc = [do*c co];
      if isempty(e)
         ee = [];
      else 
         ee = blkdiag(e,eye(no));
      end
      if Ts==0
         Ro = lyapchol(aa',cc',ee','noscale');
      else
         Ro = dlyapchol(aa',cc',ee','noscale');
      end
      Ro = Ro(1:ns,1:ns);
   end
   alpha = 0;
end

% Compute SVD of the "product of the Cholesky factors"
% NOTE: Numerically, the product SVD algorithm of Heath et al. (reference [12]
% of [1]) is superior to forming the product ro*rr' directly and then
% computing the SVD.  In other words, the following code should be used:
%      [u,g,v] = prodsvd(ro,rr')
if isempty(e)
   [u,s,v] = svd(Ro*Rr');
   tau = norm(abs(Ro)*abs(Rr)','fro');
else
   [u,s,v] = svd(Ro*e*Rr');
   tau = norm(abs(Ro)*abs(e)*abs(Rr)','fro');
end
gs = diag(s);
% Estimate threshold ZEROTOL below which the computed HSV should be
% considered zero (note threshold decided wrt to original gramians not
% limited ones)
if RELERR
   gs = min(gs,1);
   ZeroTol = 100*eps * tau;
else
   ZeroTol = 100*eps * (norm(d,1) + tau);
end

% Add unstable HSV
nns = order(Dns);
g = [inf(nns,1) ; gs ; zeros(nsnm,1)];

% Error bound (ERRBND are the error bounds for orders 0:NSV-1)
if RELERR
   % Relative error bound
   aux = (1+g)./(1-g);
   aux(isinf(g)) = inf;
   errbnd = flipud(abs(cumprod(flipud(aux)))-1);
else
   % Absolute error bound
   errbnd = 2*flipud(cumsum(flipud(g)));
end
   
% Return balancing data
BalData = mor.util.BalancedTruncationCache(sepInfo);
BalData.Sigma = g;
BalData.Error = errbnd;
if isempty(gs)
   BalData.Energy = g;
else
   BalData.Energy = (g/gs(1)).^2;
end
BalData.Regularization = alpha;
BalData.ZeroTol = ZeroTol;
BalData.Lr = Rr';
BalData.Lo = Ro';
BalData.u = u;
BalData.v = v;
