function BalData = hsvd(D,Options)
%HSVD  Computes low-rank approximations of gramians and corresponding
%      Hankel singular values.
%
%   Second input OPTIONS must be of class mor.SparseBalancedTruncationOptions.

%	 Author(s): P. Gahinet
%   Copyright 2023 The MathWorks, Inc.

% Note: CT only at this point
Ts = abs(D.Ts);
if hasInternalDelay(D)
   error(message('Control:transformation:BALROM14'))
end
a = D.a;  b = D.b;  c = D.c;  d = D.d;  e = D.e;
nx = size(a,1);

% Apply stabilizing offset
alpha = Options.Offset;
if alpha>0
   if isempty(e)
      a = a - alpha * speye(nx);
   else
      a = a - alpha * e;
   end
end

% Compoute low-rank approximations
[Lr,Lo,FF,LyapRes] = LRADI(a,b,c,d,e,Ts,Options);

% Compute SVD of the "product of the Cholesky factors"
if isempty(e)
   [u,g,v] = svd(Lo'*Lr,"econ","vector");
else
   [u,g,v] = svd(Lo'*e*Lr,"econ","vector");
end
if numel(g)>nx
   % Rr or Ro can be tall for small problems
   g = g(1:nx,:);  u = u(:,1:nx);  v = v(:,1:nx);
end

% HSVs below ZeroTol are treated as zero (noisy).
% Note: In principle, ZeroTol should be greater than Options.RankTol, but  
% this is too conservative for models with integrators (large first HSV). 
% Relax at the cost of error bound violations when g(r)<<RankTol*g(1).
ZeroTol = 1e-12*g(1);

% Absolute error bound
errbnd = 2*flipud(cumsum(flipud(g)));
   
% Save balancing data
BalData = mor.util.BalancedTruncationCache();
BalData.Split = [0 numel(g)];
BalData.Sigma = g;
BalData.Error = errbnd;
if isempty(g)
   BalData.Energy = zeros(0,1);
else
   BalData.Energy = (g/g(1)).^2;
end
BalData.ZeroTol = ZeroTol;
BalData.Lr = Lr;
BalData.Lo = Lo;
BalData.u = u;
BalData.v = v;
BalData.Feedthrough = FF;     % 1x2 cell
BalData.Residuals = LyapRes;  % 1x2 vector
