function R = dlyapchol(A, B, E, ~)
%DLYAPCHOL  Square-root solver for discrete-time Lyapunov equations.
%
%   R = DLYAPCHOL(A,B) computes a Cholesky factorization X = R'*R of 
%   the solution X to the Lyapunov matrix equation:
%
%       A*X*A'- X + B*B' = 0
%
%   All eigenvalues of A must lie in the open unit disk for R to exist.
%
%   R = DLYAPCHOL(A,B,E) computes a Cholesky factorization X = R'*R of
%   X solving the generalized Lyapunov equation:
%
%       A*X*A' - E*X*E' + B*B' = 0
%
%   All generalized eigenvalues of (A,E) must lie in the open unit disk 
%   for R to exist.
%
%   See also DLYAP, LYAPCHOL.

%	 Author(s): V. Sima, P. Gahinet
%   Copyright 1986-2023 The MathWorks, Inc.

%  DLYAPCHOL is based on the SLICOT routines SB03OD and SG03BD.
arguments 
   A(:,:) double {mustBeFinite}
   B(:,:) double {mustBeFinite}
   E(:,:) double {mustBeFinite} = [];
   ~
end
DoScaling = (nargin<4);

% Validate data 
[nx,nc] = size(A);
if nc~=nx
   error(message('Control:foundation:LyapChol1','dlyapchol'))
elseif ~isempty(E) && any(size(E)~=nx)
   error(message('Control:foundation:LyapChol2','dlyapchol'))
elseif size(B,1)~=nx
   error(message('Control:foundation:LyapChol3','dlyapchol'))
end

% Convert to full
A = full(A);  B = full(B);  E = full(E);

% Balance to minimize spectrum distorsions in Schur/QZ factorizations
if DoScaling
   [A,~,~,E,s] = aebalance(A,[],[],E,'noperm','fullbal');  % S\A*S, S\E*S
   B = s.\B;  % B=S\B
end

% Solve equation
try
   % Call SLICOT routines SB03OD/Z (E=[]) or SG03BD/Z (descriptor)
   % See LYAPCHOL for details
   if isempty(E)
      % Tranpose when A is not upper Hessenberg or its transpose also is.
      TRANS = ~isbanded(A,1,nx) || isbanded(A,nx,1);
   else
      TRANS = ~(isbanded(A,1,nx) && isbanded(E,1,nx)) || ...
         (isbanded(A,nx,1) && isbanded(E,nx,1));
   end
   if TRANS
      R = hlyapslv('DO',A',E',B');
   else
      R = hlyapslv('DC',A,E,B);
   end
catch ME
   switch ME.identifier
      case 'Control:foundation:hlyapslv4'
         error(message('Control:foundation:LyapChol5'))
      otherwise
         throw(ME)
   end
end

% Undo scaling
if DoScaling
   R = R.*s'; % R->R*S
end
