function [a,b,c,t] = compbal(a,b,c)
%COMPBAL  Balancing for SIMO state-space realizations in 
%    companion form.
%
%    [A1,B1,C1,T] = COMPBAL(A,B,C)  returns the balanced 
%    realization:
%            -1                   -1
%      A1 = T  * A * T  ,   B1 = T  * B  ,   C1 = C * T 
%         
%    where T is a diagonal matrix.

%    Authors: P. Gahinet, C. Moler, and A. Laub, 4-5-96
%    Copyright 1986-2000 The MathWorks, Inc. 
%    $Revision: 1.10 $  $Date: 2000/06/14 17:14:06 $

nx = size(a,1);
ny = size(c,1);
if nx==0,
   return
end

% Replace C by the norms of its columns in the SIMO case
cnorms = max(abs(c),[],1);

% Balance [a b;cnorms 0]
% Watch for state permutation in BALANCE
[t,abcd] = balance([a b;cnorms 0]);
[dt,p] = max(t,[],2);
abcd = abcd(p,p);  % undo possible permutation
s = dt(nx+1);      % equalizes |B| and |C|
dt = dt(1:nx);     % balances the states

% Extract transformed matrices
a = abcd(1:nx,1:nx);
b = abcd(1:nx,nx+1);
% Perform s\c*t in o(nx^2) flops
c = (c.*repmat(dt',size(c,1),1))/s;

% Form T
t = s * diag(dt);

