function [a,b,c,e,t] = abcbalance(a,b,c,e,condt,NoPerm)
%ABCBALANCE  Balancing of state-space models.
%
%   [A,B,C,E,T] = ABCBALANCE(A,B,C,E,CONDT) uses BALANCE to
%   compute a diagonal similarity transformation T such that 
%   [T*A/T,T*B;C/T,0] has approximately equal row and column 
%   norms.  CONDT specifies an upper bound on the condition 
%   number of T.
%
%   LOW-LEVEL UTILITY, CALLED BY SSBAL.

%   Authors: P. Gahinet and C. Moler
%   Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.6 $  $Date: 2000/06/14 17:14:06 $

% RE: Expects A,B,C,E to be matrices or ND arrays

NoPerm = (nargin>5);  % disables state permutation when = 1
   
% Get dimensions
asizes = size(a);
nx = asizes(1);
ne = size(e,1);
nu = size(b,2);
ny = size(c,1);
nsys = prod(asizes(3:end));

% Quick exit when no state
if nx==0 | nsys==0,
   t = eye(nx);   
   return
end

% Form 2D matrix M = [|A|+|E| |B|;|C| 0] to be balanced
mae = max(abs(a(:,:,:)),[],3);
if ne,
   mae = mae + max(abs(e(:,:,:)),[],3);
end
mb = max([max(abs(b(:,:,:)),[],3) , zeros(nx,1)],[],2);
mc = max([max(abs(c(:,:,:)),[],3) ; zeros(1,nx)],[],1);
% To activate balancing when M=[A B;C 0] is triangular, set
% zero entries of B and C to
%   * C(j)=TINY*|C| if C~=0, |B| otherwise 
%   * B(j)=TINY*|B| if B~=0, |C| otherwise
% RE: This perturbation will be visible only when C.ej=0 and 
%     A.ej=s.ej, i.e., x(j) is unobservable or uncontrollable
bmax = max(mb); 
cmax = max(mc);
tiny = max(1e-100,(1/condt)^2);
mb(~mb) = tiny * bmax + (bmax==0)*cmax;
mc(~mc) = tiny * cmax + (cmax==0)*bmax;


% Perform the balancing with BALANCE
[t,Mb] = balance([mae mb;mc 0]);
dt = max(t(1:nx,:),[],2);     % balances the states
s = max(t(nx+1,:));           % equalizes |B| and |C|

% BALANCE(A) may permute the rows/cols of A to enforce triangularity
% (desirable for balancing prior to Hessenberg reduction). Acquire
% this permutation when NOPERM=0 
if NoPerm
   p = ':';  % No permutation
else
   [tae,junk] = balance(mae);
   [junk,p] = max(tae,[],2);
end

% If cond(T) exceeds CONDT, rescale diag(T) to match the bound on CONDT 
if isfinite(condt) & max(dt)>10*condt*min(dt),
   dt = log2(dt);
   scalf = log2(condt)/(max(dt)-min(dt));
   dt = pow2(round(scalf*dt));
end

% Compute balanced realization
dti = 1./dt;
dtt = dt.';
for k=1:nsys,
   a(p,p,k) = (dti(:,ones(1,nx)) .* a(:,:,k)) .* dtt(ones(1,nx),:);
   b(p,:,k) = (dti(:,ones(1,nu)) .* b(:,:,k)) * s;
   c(:,p,k) = (c(:,:,k) .* dtt(ones(1,ny),:)) / s;
   if ne,
      e(p,p,k) = (dti(:,ones(1,nx)) .* e(:,:,k)) .* dtt(ones(1,nx),:);
   end
end

% Return inv(T) to be consistent with SS2SS
t = zeros(nx);
t(p,:) = s * diag(dti);
