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

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

%  LYAPCHOL 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','lyapchol'))
elseif ~isempty(E) && any(size(E)~=nx)
   error(message('Control:foundation:LyapChol2','lyapchol'))
elseif size(B,1)~=nx
   error(message('Control:foundation:LyapChol3','lyapchol'))
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)
   % Note: Solving A*X+X*A'+B*B'=0 is more expensive than solving 
   % A'*X+X*A+C'*C=0 in SLICOT (extra QR factorization of R). For 
   % maximum performance, transpose the data except when (A,E) is upper 
   % Hessenberg (g148601: When (A,E) is already in Schur form, transposing
   % causes SLICOT to recompute the eigenvalues via qz(A',E'), which is 
   % inefficient and may disagree with eig(A,E) on stability).
   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('CO',A',E',B');
   else
      R = hlyapslv('CC',A,E,B);
   end
catch ME
   switch ME.identifier
      case 'Control:foundation:hlyapslv4'
         error(message('Control:foundation:LyapChol4'))
      otherwise
         throw(ME)
   end
end

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

