function [IC,B,xkeep] = sminrealLFT(IC,B)
% Eliminates structurally nonminimal states, delays, and blocks.

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

% Note: The rows and columns of IC.d are partitioned as [y;z;zd] and [u;w;wd] 
% where u,y are the external I/Os, w,z are the block channels, and wd,zd
% are the delay channels.

% States and delays
Delay = IC.Delay;
nfd = length(Delay.Internal);
nx = size(IC.a,1);

% Signal dimensions
[rs,cs] = size(IC.d);
[nw,nz] = iosize(B);
ny = rs-nz-nfd;  
nu = cs-nw-nfd;
nb = numel(B);

% Replace each block by a scalar in IC matrix
ICS = ([IC.a IC.b;IC.c IC.d]~=0); 
ir = nx+ny;   ic = nx+nu;
for ct=1:nb
   bs = iosize(B(ct));
   ICS(:,nx+nu+ct) = any(ICS(:,ic+1:ic+bs(1)),2);
   ICS(nx+ny+ct,:) = any(ICS(ir+1:ir+bs(2),:),1);
   ir = ir+bs(2);   ic = ic+bs(1);
end
nxyb = nx+ny+nb;   nxub = nx+nu+nb;
ICS(:,nxub+1:nxub+nfd) = ICS(:,ic+1:ic+nfd);
ICS(nxyb+1:nxyb+nfd,:) = ICS(ir+1:ir+nfd,:);
% D now has ny+nb+nfd rows and nu+nb+nfd columns
   
% Identify components to eliminate
if isempty(IC.e)
   E = [];
else
   E = blkdiag(IC.e,eye(nb+nfd));
end
ix = nx+1:nx+ny;  ixc = 1:nxyb+nfd;  ixc(ix) = [];
jx = nx+1:nx+nu;  jxc = 1:nxub+nfd;  jxc(jx) = [];
[~,~,~,~,xbdkeep,ebdkeep] = smreal(ICS(ixc,jxc),ICS(ixc,jx),ICS(ix,jxc),E);
xkeep = xbdkeep(1:nx);
bkeep = xbdkeep(nx+1:nx+nb);
ix = find(xkeep);
ie = find(ebdkeep(1:nx,:));
id = find(xbdkeep(nx+nb+1:nx+nb+nfd,:))';
   
% Perform simplification
[jkeep,ikeep] = getRowColSelection(B,bkeep);
iy = [1:ny , ny+ikeep];  iyz = [iy , ny+nz+id];
iu = [1:nu , nu+jkeep];  iuw = [iu , nu+nw+id];
IC.a = IC.a(ie,ix);
if ~isempty(IC.e)
   IC.e = IC.e(ie,ix);
end
IC.dx0 = IC.dx0(ie,:);
IC.x0 = IC.x0(ix,:);
IC.b = IC.b(ie,iuw);
IC.u0 = IC.u0(iuw,:);
IC.c = IC.c(iyz,ix);
IC.y0 = IC.y0(iyz,:);
IC.d = IC.d(iyz,iuw);
Delay.Output = Delay.Output(iy,:);
Delay.Input = Delay.Input(iu,:);
Delay.Internal = Delay.Internal(id,:);
IC.Delay = Delay;
if ~isempty(IC.StateName)
   IC.StateName = IC.StateName(ix,:);
end
if ~isempty(IC.StatePath)
   IC.StatePath = IC.StatePath(ix,:);
end
if ~isempty(IC.StateUnit)
   IC.StateUnit = IC.StateUnit(ix,:);
end
% Expand XKEEP to include block states and match model order
xkeep = [xkeep ; getStateSelection(B,bkeep)];
B = B(bkeep,:);
