%function [d] = sblkbal(m,opt)
%	Balances a square matrix assuming only scalar blocks.
%	"opt" is an optional string specifying the method: 
%	'o' is for Osborne's method , 'p' is for a Perron vector method.
%	The methods give comparable results for most matrices, but 'p' exactly 
%	computes mu for positive matrices.  The default is 'p' for size(m) < 10
%	and 'o' otherwise, since 'p' is faster for small matrices but has a
%	growth rate of n^3 compared with less than n^2 for 'o'.  This is partly
%	due to the matlab implementation, which greatly favors 'p'.
function [d] = sblkbal(m,opt)
	[nr nc] = size(m);
	if nargin==2 & ( any(opt=='p') | any(opt=='o') )
		elseif nr < 10
			opt = 'p';
		   else
			opt = 'o';
		
	end %if nargin

		% First, add small numbers to the abs(m) to avoid degeneracies
		% due to zeros rows or columns.  Then get an initial balancing
		% of the matrix, which is now real and positive.

[db,sm] = balance(abs(m)+10*eps*ones(m)); 
dsqr=diag(db)'; dsqr = ones(1,size(dsqr))./(dsqr.*dsqr);

if any(opt=='p') 
	[vec,ev] = eig(sm);		% This is a Perron eigenvector approach
	[ev,indx] = max(diag(ev));	% x and y are right and left eigenvecs
	x = vec(:,indx);		  	%	for largest eig (ev)
	y = zeros(x'); y(indx)=1;	% x,y and ev are all positive
	y = y/vec;			
	dsqr= abs(y./x').*dsqr; 	% balance r and l eigenvectors
  else			
	j = 0;				% Osborne's method
	smoffd = sm - diag(diag(sm));   		%zero the diagonal
	sqr = smoffd.*smoffd;				
	percd = 1; oldcost = sum(sum(sqr)');
	while j < 30 & percd > 1e-4
		j = j + 1;
		rowsum = sum(sqr') ; colsum = sum(sqr);
		newdsqr = sqrt(colsum./rowsum);  	%balance row and
		dsqr = newdsqr.*dsqr;			%	column norms
		newdupdate = ones(nr,1)*newdsqr;
		sqr = newdupdate'.*sqr./newdupdate;
		cost = sum(sum(sqr)');			%cost=frobenius norm
		percd = (oldcost-cost)/oldcost;		%percd= percentage
		oldcost = cost;				%        drop in cost
	end % while j 
end  %if any(opt=='p')
 
d = sqrt(dsqr);
%
% Copyright MUSYN INC 1991,  All Rights Reserved
