function [s,p,b] = xbalance(a,varargin)
%XBALANCE  Balances A matrix.
%
%   [S,P,B] = XBALANCE(A) returns the balanced matrix B = T\A*T 
%   together with the scaling S and permutation P such that 
%   T(:,P) = diag(S).  
%
%   [S,P,B] = XBALANCE(A,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.10.1 $  $Date: 2004/12/10 19:33:32 $

% Differences with BALANCE:
%  * Use two-step approach when PERM is required to prevent scaling 
%    from operating only on a submatrix (can lead to poor overall 
%    scaling in some cases) (g209996)
%  * Skip permutation when A or A' is upper Hessenberg (g241217)
XPerm = ~any(strcmp('noperm',varargin));       % state permutation enabled?
SafeBalance = any(strcmp('safebal',varargin)); % perform safe balancing?
n = size(a,1);

% Compute scaling
[s,p,b] = balance(a,'noperm');

% Safe scaling option
if SafeBalance
   % Only balance (row/col) pair when this reduces ||A||
   rhoA = max(abs(b(:)));  % target norm for A
   ma = abs(a);
   ma(1:n+1:n^2) = 0;

   % Balancing loop
   s = ones(n,1);
   cnt = 0;
   while cnt<n
      for j=1:n
         % Rescale j-th (row,col) pair to bring down max
         % element close to rhoA
         muj = max(ma(j,:));
         nuj = max(ma(:,j));
         theta = max(rhoA,sqrt(muj*nuj));
         theta2 = 2*theta;
         if muj>theta2
            % Bring row norm down
            sj = pow2(fix(log2(muj/theta)));
            cnt = 0;
            % Update scaling
            s(j) = s(j) * sj;
            ma(j,:) = ma(j,:)/sj;
            ma(:,j) = ma(:,j)*sj;
         elseif nuj>theta2
            % Bring column norm down
            sj = pow2(fix(log2(theta/nuj)));
            cnt = 0;
            % Update scaling
            s(j) = s(j) * sj;
            ma(j,:) = ma(j,:)/sj;
            ma(:,j) = ma(:,j)*sj;
         else
            % No progress (sj=1)
            cnt = cnt+1;
            if cnt>n
               break
            end
         end
      end
   end
   % Apply "safe" scaling
   b = lrscale(a,1./s,s);
end

% Compute permutation
if XPerm
   if norm(tril(a,-2),1)==0
      % Real upper Hessenberg: no permutation
      p = (1:n)';
   elseif norm(triu(a,2),1)==0
      % Real lower Hessenberg
      p = (n:-1:1)';
   else
      [junk,p,junk] = balance(a);
   end
   b(p,p) = b;
else
   p = (1:n)';  % No permutation
end