home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / SYSCL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.0 KB  |  92 lines
  1. % function sysout = syscl(sys)
  2. %
  3. %  SYSCL scales a state space realization in order to 
  4. %  reduce errors in later calculations, using diagonal 
  5. %  similarity transformations, and deleting any rows 
  6. %  and columns which are zero (apart from the diagonal of A).
  7.  
  8. function sysout=syscl(sys)
  9.  
  10. [a,b,c,d]=unpck(sys);
  11. [n,n1]=size(a);
  12. if n1~=n
  13.   disp('A is not square')
  14.   return
  15. end
  16. [n1,m]=size(b)
  17. if n1~=n
  18.   disp('dimensions do not match')
  19.   return
  20. end
  21. [p,n1]=size(c)
  22. if n1~=n
  23.   disp('dimensions do not match')
  24.   return
  25. end
  26. ib=16;bb=256;il=1;ih=n;
  27.  
  28. % remove and save the diagonal of a.
  29.  
  30. adiag=diag(a); a=a-diag(adiag,0);
  31.  
  32. % find zero rows in a & b and zero columns in a & c
  33.  
  34. perm=ones(1,n);for i=1:n,perm(i)=i;end;
  35. j=ih;
  36. while j>0,
  37.   nozra=any((a(perm,perm))');
  38.   nozb=any([b(perm,:),zeros(n,1)]')|nozra;
  39.   nozca=any(a(perm,perm));
  40.   nozc=any([c(:,perm);zeros(1,n)])|nozca;
  41.   if ~nozb(j)|~nozc(j),
  42.     for i=j+1:n,perm(i-1)=perm(i);end;
  43.     n=n-1;ih=ih-1;j=ih+1;
  44.     perm=perm(1:n);
  45.   elseif j==1,break;
  46.   end;
  47.   j=j-1;
  48. end;
  49. j=ih;
  50. while j>0,
  51.   if ~nozra(j),t=perm(j);perm(j)=perm(ih);perm(ih)=t;
  52.     t=nozra(j);nozra(j)=nozra(ih);nozra(ih)=t;
  53.     ih=ih-1;j=ih+1;
  54.   elseif j==1,break;
  55.   end;
  56.   j=j-1;
  57. end;
  58. j=il;
  59. while j<=ih,
  60.   if ~nozca(j),t=perm(j);perm(j)=perm(il);perm(il)=t;
  61.     t=nozca(j);nozca(j)=nozca(il);nozca(il)=t;
  62.     il=il+1;j=il-1;
  63.   elseif j==ih,break;
  64.   end;
  65.   j=j+1;
  66. end;
  67. perm=perm(1:n);a=a(perm,perm);b=b(perm,:);c=c(:,perm);
  68. adiag=adiag(perm);
  69. %n,a,b,c,il,ih
  70. % scaling
  71.  
  72. fail=1;
  73. while fail,fail=0;
  74.   for i=1:n,
  75.     cosum=sum([abs(a(:,i));abs(c(:,i))]);
  76.     rosum=sum([abs(a(i,:)), abs(b(i,:))]);
  77.     f=1;g=rosum/ib;s=cosum+rosum;
  78.     while cosum<g,f=f*ib;cosum=cosum*bb;end;
  79.     g=rosum*ib;
  80.     while cosum>=g,f=f/ib;cosum=cosum/bb;end;
  81.     if(cosum+rosum)/f<0.95*s,
  82.       g=1/f;fail=1;
  83.       a(i,:)=a(i,:)*g;a(:,i)=a(:,i)*f;
  84.       b(i,:)=b(i,:)*g;c(:,i)=c(:,i)*f;
  85.     end;
  86.   end;
  87. end;
  88. a=a+diag(adiag,0);
  89. sysout=pck(a,b,c,d);
  90. %
  91. % Copyright MUSYN INC 1991,  All Rights Reserved
  92.