function [A,B,C,E,StateMat] = smreal(A,B,C,E,StateMat)
%SMREAL  Compute a structurally minimal realization.
%
%   [A,B,C,E,NAMES] = SMREAL(A,B,C,E,NAMES) eliminates 
%   states that are not connected to any input or output.  
%   The resulting realization is structurally minimal, i.e., 
%   minimal when all nonzero entries of A,B,C,E are set to 1.
%   The argument NAMES can contains the state names or any 
%   array whose rows map to the states.

%   Author(s): P. Gahinet
%   Copyright 1986-2003 The MathWorks, Inc.
%   $Revision: 1.1.6.3 $  $Date: 2004/04/10 23:14:58 $

% Uses characterization of the controllable subspace as smallest
% subspace S such that Im(B)<AS+ES, dim(S)=dim(AS+ES).
% This subspace is generated by
%     S0 = B;  Sk+1 = Sk + ASk + ESk
% To compute a structurally minimal realization, this recursion 
% is implemented on the boolean algebra (+ means OR, * means AND) 
% with A = A(A~=0) and b = any(B,2):
%     s(0) = b;  s(k+1) = A s(k) | E s(k) | s(k)
% The limit (stationary value) of the boolean vector s(k) indicates
% which states to keep.
%
% RE: The propagation equation is implemented as
%       ds(k+1) = (A+E) ds(k) & ~s(k)  
%       s(k+1) = s(k) + ds(k+1)
%     (ds(k+1) = newly activated states at iteration k)

% Dimensionality
sizes = size(A);

% Derive structural connectivity information
Ab = any(A,3);
Bb = any(B,2);
Cb = any(C,1);
if ~isempty(E),
   Ab = Ab | any(E,3);
end

% Compute structurally minimal realization
% First identify structurally controllable states
xc = Bb;    % x(1)
dx = xc;    % dx(k) = x(k)-x(k-1)
while any(dx),
   Adx = any(Ab(:,dx),2);  % A dx(k) + E dx(k)
   dx = Adx & ~xc;         % dx(k+1)
   xc = xc | Adx;          % xc(k+1)
end

% Identify structurally observable states
xo = Cb;   
dx = xo;
while any(dx),
   Adx = any(Ab(dx,:),1);
   dx = Adx & ~xo;
   xo = xo | Adx;
end

% Minimal states are both controllable and observable
xco = xo & xc';

% Extract structurally minimal realization
if ~all(xco),
   A = A(xco,xco);
   B = B(xco,:);
   C = C(:,xco);
   if ~isempty(E),  
      E = E(xco,xco);
   end
   if nargin>4
       % Extract rows associated with remaining states
      StateMat = StateMat(xco,:);
   end
end