home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / ORSF.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.4 KB  |  93 lines
  1. % function [u,t,k]=orsf(u,a,fl,bord)
  2. %
  3. %  *****  UNTESTED SUBROTUINE  *****
  4. %
  5. % Given a matrix A in real Schur form, ORSF finds 
  6. %  T = U'*A*U with U unitary and the eigenvalues of 
  7. %  T are ordered to have increasing real parts, when 
  8. %  FL = 'o', or if FL = 's' sorted into two groups 
  9. %  with the eigenvalues of the first group having 
  10. %  real part < BORD. The default value of BORD is zero.
  11. %  K gives the number of poles with real part < BORD.
  12. %
  13. %  See Also: OCSF
  14.  
  15. function [u,t,k]=orsf(u,a,fl,bord)
  16.  
  17. jay=sqrt(-1);
  18. [n,m]=size(a);
  19. if n~=m
  20.   disp('A is not square')
  21.   return
  22. end
  23. if (nargin>2)&(fl~='o')&(fl~='s'), 
  24.   disp(' usage:  [u,t]=orsf(u,a,fl,bord) with fl=''o'' pr ''s'' ')
  25.   return
  26. end
  27. if nargin<=3, bord=0;end
  28. [u,t]=rsf2csf(u,a);
  29. i=1; while i<n, if imag(t(i,i))==0,i=i+1;
  30.                else t(i+1,i+1)=t(i,i)';i=i+2;end;end %ensures roots exact conj.
  31. if fl=='o',
  32.   ia=1;id=1;ib=n-1;il=1;
  33.   while il~=0,
  34.     il=0;
  35.     for i=ia:id:ib,
  36.       tii=t(i,i);tioio=t(i+1,i+1);
  37.       if real(tii)>real(tioio),
  38.         [s,r]=qr([t(i,i+1);tioio-tii]);
  39.         t(1:i+1,i:i+1)=t(1:i+1,i:i+1)*s;
  40.         t(i:i+1,i:n)=s'*t(i:i+1,i:n);
  41.         u(1:n,i:i+1)=u(1:n,i:i+1)*s;
  42.         t(i:i+1,i:i+1)=[tioio,t(i,i+1);0 tii];
  43.         il=i;
  44.       end;
  45.     end;
  46.     ib=ia;id=-id;ia=il+id;
  47.   end;
  48. end; %if fl=='o'
  49.  
  50. if fl=='s',
  51.   ia=1;id=1;ib=n-1;il=1;
  52.   while il~=0,
  53.     il=0;
  54.     for i=ia:id:ib,
  55.       tii=t(i,i);tioio=t(i+1,i+1);
  56.       if (real(tii)>=bord) & (real(tioio)<bord),
  57.         [s,r]=qr([t(i,i+1);tioio-tii]);
  58.         t(1:i+1,i:i+1)=t(1:i+1,i:i+1)*s;
  59.         t(i:i+1,i:n)=s'*t(i:i+1,i:n);
  60.         u(1:n,i:i+1)=u(1:n,i:i+1)*s;
  61.         t(i:i+1,i:i+1)=[tioio,t(i,i+1);0 tii];
  62.         il=i;
  63.       end;
  64.     end;
  65.     ib=ia;id=-id;ia=il+id;
  66.   end;
  67. end; %if fl=='s'
  68.  
  69. % find k
  70. k= max(find(real(diag(t))<bord));
  71. % now regain real form.
  72. i=1;
  73. while i<=n,
  74.   if imag(t(i,i))==0,
  75.      [q,r,e]=qr([real(u(:,i)) imag(u(:,i))]');
  76.      s=jay*q(1,2)+q(2,2);
  77.      u(:,i)=u(:,i)*s;
  78.      if i>1, t(1:i-1,i)=t(1:i-1,i)*s;end
  79.      if i<n, t(i,i+1:n)=s'*t(i,i+1:n);end
  80.      i=i+1;
  81.     else
  82.      [q,r,e]=qr([real(u(:,i:i+1)) imag(u(:,i:i+1))]');
  83.      s=jay*q(1:2,3:4)+q(3:4,3:4);
  84.      t(1:i+1,i:i+1)=t(1:i+1,i:i+1)*s;
  85.      t(i:i+1,i:n)=s'*t(i:i+1,i:n);
  86.      u(1:n,i:i+1)=u(1:n,i:i+1)*s;
  87.      i=i+2;
  88.    end %if imag(t(i,i))==0,
  89. end %while i<=n
  90. u=real(u);t=real(t);
  91. %
  92. % Copyright MUSYN INC 1991,  All Rights Reserved
  93.