function [a,b,c,e,s,p] = aebalance(a,b,c,e,varargin)
%AEBALANCE  Balances (A,E) pair.
%
%   [A,B,C,E,S,P] = AEBALANCE(A,B,C,E) balances the (A,E) pair
%   and returns the balanced system (T\A*T,T\E*T,T\B,C*T) together 
%   with the scaling S and permutation P such that T(:,P) = diag(S).  
%  
%   [A,B,C,E,S,P] = AEBALANCE(A,B,C,E,Option1,Option2,...) specifies 
%   additional options as strings:
%     'noperm'    Prevents state permutation during balancing
%     'perm'      Enables state permutation (default)
%     'fullbal'   Performs full balancing of A using BALANCE
%                 (default)
%     'safebal'   Performs row/col scaling only when it decreases 
%                 ||A||
%
%   LOW-LEVEL UTILITY.

%   Authors: P. Gahinet
%   Copyright 1986-2002 The MathWorks, Inc. 
%   $Revision: 1.1.6.2 $  $Date: 2004/12/10 19:25:53 $

% RE: Expects A,B,C,E to be matrices 2D arrays

% Get dimensions
nx = size(a,1);
ne = size(e,1);

% Quick exit when no state
if nx==0,
   s = ones(0,1); 
   p = zeros(0,1);
   return
end

% Form 2D matrix M = [|A|+|E| |B|;|C| 0] to be balanced
if ne==0,
   mae = abs(a);
else
   anorm = norm(a-diag(diag(a)),1);
   enorm = norm(e-diag(diag(e)),1);
   if anorm>0 && enorm>0
      mae = enorm*abs(a)+anorm*abs(e);
   else
      mae = abs(a)+abs(e);
   end
end

% Perform full balancing
[s,p] = xbalance(mae,varargin{:}); 

% Apply scaling and permutation to (A,E)
is = 1./s;
a(p,p) = lrscale(a,is,s);
if ne,
   e(p,p) = lrscale(e,is,s);
end

% Scale B,C if nonempty
if ~isempty(b) && ~isempty(c)
   mb = max(abs(b),[],2);
   mc = max(abs(c),[],1);
   if any(mb) && any(mc)
      % Equalize ||B|| and ||C||
      bnorm = max(lrscale(mb,is,[]));
      cnorm = max(lrscale(mc,[],s));
      sbc = pow2(round(log2(cnorm/bnorm)/2));
      s = s / sbc;
      is = is * sbc;
   end
   b(p,:) = lrscale(b,is,[]);
   c(:,p) = lrscale(c,[],s);
end