function [D,SingularFlag] = ssequiv(D,TL,eL,TR,eR)
% Equivalence transformation (TL*A*TR,TL*B,C*TR,TL*E*TR).

%   Copyright 2023 The MathWorks, Inc.
nL = size(TL,1);
nR = size(TR,1);
if nL==0 && nR==0
   SingularFlag = false;  return
end
a = D.a;  b = D.b;  c = D.c;  e = D.e;
x0 = D.x0;  dx0 = D.dx0;
nx = size(a,1);
if (nL>0 && nL~=nx) || (nR>0 && nR~=nx)
   error(message('Control:transformation:ssequiv1'))
end
[ny,nu] = size(D.d);
EXPLICIT = isempty(e);
if nL==0
   eL = 0;
end
if nR==0
   eR = 0;
end

if EXPLICIT && eL+eR==0 && isequal(TL,TR)
   % Special handling of T*A/T or T\A*T
   T = TL;
   if eL>0
      % T*A/T
      aux = matlab.internal.math.nowarn.mrdivide([a ; c],T);
      a = T*aux(1:nx,:);  b = T*b;  c = aux(nx+1:nx+ny,:);
      x0 = T*x0;  dx0 = T*dx0;
   else
      % T\A*T
      aux = matlab.internal.math.nowarn.mldivide(T,[a b x0 dx0]);
      a = aux(:,1:nx)*T;  b = aux(:,nx+1:nx+nu);  c = c*T;
      x0 = aux(:,nx+nu+1);  dx0 = aux(:,nx+nu+2);
   end
else
   % General case
   aux = [];
   if EXPLICIT
      e = eye(nx);
   end
   if eL>0
      a = TL*a;  b = TL*b;  e = TL*e;  dx0 = TL*dx0;
   end
   if eR>0
      a = a*TR;  c = c*TR;  e = e*TR;  x0 = TR\x0;
   end
   if eL<0
      aux = matlab.internal.math.nowarn.mldivide(TL,[a e b dx0]);
      a = aux(:,1:nx);  e = aux(:,nx+1:2*nx);  b = aux(:,2*nx+1:2*nx+nu);
      dx0 = aux(:,2*nx+nu+1);
   end
   if eR<0
      aux = matlab.internal.math.nowarn.mrdivide([a;e;c],TR);
      a = aux(1:nx,:);  e = aux(nx+1:2*nx,:);  c = aux(2*nx+1:2*nx+ny,:);
      x0 = TR*x0;
   end
end
SingularFlag = ~allfinite(aux);
D.a = a;  D.b = b;  D.c = c;  D.e = e;
D.dx0 = dx0;  D.x0 = x0;
D.StateName = strings(0,1);
D.StatePath = strings(0,1);
D.StateUnit = strings(0,1);
D.Scaled = false;
