function [G,SingularFlag] = imp2exp(G,yidx,uidx)
% Implicit to explicit conversion. Given the implicit model
%    0 = G11(s) Y + G12(s) U
%    Z = G21(s) Y + G22(s) U
% IMP2EXP computes the explicit model H(s) mapping U to [Y;Z].
% The index vectors YIDX,UIDX partition the columns of G into Y,U and
% the first NUMEL(YIDX) outputs of G are constrained to be zero.

%   Copyright 2003-2020 The MathWorks, Inc.
if any(G.Delay.Input(yidx)>0)
   error(message('Control:transformation:imp2exp6'))
end
[nOut,nIn] = iosize(G);
nX = size(G.a,1);
nDelay = numel(G.Delay.Internal);
nU = numel(uidx);
nY = numel(yidx);  % Note: nY + numel(uidx) = nIn
iU = [uidx nIn+1:nIn+nDelay];  % U input channels + internal delays
d = G.d;  
u0 = G.u0;

% Update delays
G.Delay.Input = G.Delay.Input(uidx);
G.Delay.Output = [zeros(nY,1) ; G.Delay.Output(nY+1:nOut)];

if nX==0 && nDelay==0
   % Static problem
   D = [zeros(nY,nIn-nY,'like',d) u0(yidx,:); d(nY+1:nOut,iU) G.y0(nY+1:nOut,:)] - ...
      [eye(nY,'like',d) ; d(nY+1:nOut,yidx)] * ...
      matlab.internal.math.nowarn.mldivide(d(1:nY,yidx),[d(1:nY,iU) G.y0(1:nY,:)]);
   if allfinite(D)
      SingularFlag = false;
      G = createGain(G,D(:,1:end-1));
      G.u0 = u0(iU,:);
      G.y0 = D(:,end);
   else
      SingularFlag = true;
      G = createGain(G,NaN(nOut,nU));
   end
else
   % Form implicit DAE with state vector [X ; Y]
   a = G.a;  b = G.b;  c = G.c;  e = G.e;
   A = [a b(:,yidx);c(1:nY,:) d(1:nY,yidx)];
   if isempty(e)
      E = blkdiag(eye(nX,'like',a),zeros(nY,'like',a));
   else
      E = blkdiag(e,zeros(nY,'like',a));
   end
   p = dmperm((A~=0) | (E~=0));
   if any(p==0)
      % Pencil A-s*E is structurally singular (G11 normal-rank-deficient)
      SingularFlag = true;
      G = createGain(G,NaN(nOut,nU));
   else
      % Build implicit model
      SingularFlag = false;
      G.a = A;   G.e = E;
      G.b = [b(:,iU) ; d(1:nY,iU)];
      G.c = [zeros(nY,nX,'like',a) eye(nY,'like',a) ; ...
         c(nY+1:nOut+nDelay,:) d(nY+1:nOut+nDelay,yidx)];
      G.d = [zeros(nY,nIn-nY+nDelay,'like',a) ; d(nY+1:nOut+nDelay,iU)];
      G.dx0 = [G.dx0 ; G.y0(1:nY,:)];
      G.x0 = [G.x0 ; u0(yidx,:)];
      G.u0 = u0(iU,:);
      G.y0 = [u0(yidx,:);G.y0(nY+1:nOut+nDelay,:)];
      G.Scaled = false;
   end
end
