home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / SBLKBAL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.0 KB  |  55 lines
  1. %function [d] = sblkbal(m,opt)
  2. %    Balances a square matrix assuming only scalar blocks.
  3. %    "opt" is an optional string specifying the method: 
  4. %    'o' is for Osborne's method , 'p' is for a Perron vector method.
  5. %    The methods give comparable results for most matrices, but 'p' exactly 
  6. %    computes mu for positive matrices.  The default is 'p' for size(m) < 10
  7. %    and 'o' otherwise, since 'p' is faster for small matrices but has a
  8. %    growth rate of n^3 compared with less than n^2 for 'o'.  This is partly
  9. %    due to the matlab implementation, which greatly favors 'p'.
  10. function [d] = sblkbal(m,opt)
  11.     [nr nc] = size(m);
  12.     if nargin==2 & ( any(opt=='p') | any(opt=='o') )
  13.         elseif nr < 10
  14.             opt = 'p';
  15.            else
  16.             opt = 'o';
  17.         
  18.     end %if nargin
  19.  
  20.         % First, add small numbers to the abs(m) to avoid degeneracies
  21.         % due to zeros rows or columns.  Then get an initial balancing
  22.         % of the matrix, which is now real and positive.
  23.  
  24. [db,sm] = balance(abs(m)+10*eps*ones(m)); 
  25. dsqr=diag(db)'; dsqr = ones(1,size(dsqr))./(dsqr.*dsqr);
  26.  
  27. if any(opt=='p') 
  28.     [vec,ev] = eig(sm);        % This is a Perron eigenvector approach
  29.     [ev,indx] = max(diag(ev));    % x and y are right and left eigenvecs
  30.     x = vec(:,indx);              %    for largest eig (ev)
  31.     y = zeros(x'); y(indx)=1;    % x,y and ev are all positive
  32.     y = y/vec;            
  33.     dsqr= abs(y./x').*dsqr;     % balance r and l eigenvectors
  34.   else            
  35.     j = 0;                % Osborne's method
  36.     smoffd = sm - diag(diag(sm));           %zero the diagonal
  37.     sqr = smoffd.*smoffd;                
  38.     percd = 1; oldcost = sum(sum(sqr)');
  39.     while j < 30 & percd > 1e-4
  40.         j = j + 1;
  41.         rowsum = sum(sqr') ; colsum = sum(sqr);
  42.         newdsqr = sqrt(colsum./rowsum);      %balance row and
  43.         dsqr = newdsqr.*dsqr;            %    column norms
  44.         newdupdate = ones(nr,1)*newdsqr;
  45.         sqr = newdupdate'.*sqr./newdupdate;
  46.         cost = sum(sum(sqr)');            %cost=frobenius norm
  47.         percd = (oldcost-cost)/oldcost;        %percd= percentage
  48.         oldcost = cost;                %        drop in cost
  49.     end % while j 
  50. end  %if any(opt=='p')
  51.  
  52. d = sqrt(dsqr);
  53. %
  54. % Copyright MUSYN INC 1991,  All Rights Reserved
  55.