home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / CSCHUR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.8 KB  |  102 lines
  1. function [v,t,m,swap] = cschur(a,Type)
  2. %
  3. % [V,T,M,SWAP] = CSCHUR(A,TYPE) produces an ordered Complex Schur form
  4. %         
  5. %              V' * A * V  = T = |T1  T12|
  6. %                                | 0   T2|,
  7. %
  8. %              m    = no. of stable poles (s or z plane)
  9. %              swap = total no. of swaps.
  10. %
  11. %       The unordered Schur form produced by "schur" is ordered using
  12. %       complex Givens rotation to swap adjacent diagonal terms.
  13. %
  14. %       Six types of ordering can be selected:
  15. %       
  16. %         Type = 1 --- Re(eig(T1)) < 0, Re(eig(T2)) > 0.
  17. %         Type = 2 --- Re(eig(T1)) > 0, Re(eig(T2)) < 0.
  18. %         Type = 3 --- eigenvalue real parts in descending order.
  19. %         Type = 4 --- eigenvalue real parts in ascending order.
  20. %         Type = 5 --- modulus of eigenvalues in descending order.
  21. %         Type = 6 --- modulus of eigenvalues in ascending order.
  22. %
  23.  
  24. % R. Y. Chiang & M. G. Safonov 4/4/86
  25. % Copyright (c) 1988 by the MathWorks, Inc.
  26. % All Rights Reserved.
  27. %---------------------------------------------------------------------
  28. [ma,na] = size(a);
  29. %
  30. % -------- Regular Complex Schur Decomposition :
  31. %
  32. [v,t] = schur(a); 
  33. [v,t] = rsf2csf(v,t);
  34. [mt,nt] = size(t); 
  35. dt = diag(t);
  36. if Type < 5
  37.    idm = find(real(dt)<zeros(mt,1));
  38. else
  39.    idm = find(abs(dt)<ones(mt,1));
  40. end
  41. [m,dum] = size(idm);  
  42. %
  43. % -------- Assign I.C. for the loop counter :
  44. %
  45. pcnt = 0.; kcnt = 0.; gcnt = 0.;
  46.  
  47. % The Chiang and Safonov routine ends here.  The rest has
  48. %  been replaced by optimized code.
  49.  
  50. % Charles R. Denham, MathWorks, 1989.
  51. % Copyright (c) 1989 by the MathWorks, Inc.
  52. % All Rights Reserved.
  53.  
  54. % Predetermine the order of sorting, to avoid possible
  55. %  infinite loop due to round-off errors.
  56.  
  57. if Type > 0 & Type < 7
  58.  
  59.    if Type == 1
  60.       d = sign(real(diag(t)));
  61.      elseif Type == 2
  62.       d = -sign(real(diag(t)));
  63.      elseif Type == 3
  64.       d = -real(diag(t));
  65.      elseif Type == 4
  66.       d = real(diag(t));
  67.      elseif Type == 5
  68.       d = -abs(diag(t));
  69.      elseif Type == 6
  70.       d = abs(diag(t));
  71.    end
  72.  
  73. % Simple bubble sorting via Given's rotations.  Matrix t
  74. %  is re-ordered in the same manner as an ascending sort
  75. %  of matrix d.  Then, matrix v is made compatible.  The
  76. %  bubble sort code is reminiscent of the SCHORD function
  77. %  in the Control_Toolbox, but the order of interchanges
  78. %  is different.
  79.  
  80.    q = eye(nt,nt);
  81.    swap = 0;
  82.    okay = 1;
  83.    while okay
  84.       okay = 0;
  85.       for k = 1:nt-1
  86.          if d(k) > d(k+1)
  87.             okay = 1;
  88.             swap = swap + 1;
  89.             g = givens(t(k,k+1), t(k,k)-t(k+1,k+1));
  90.             t(:,k:k+1) = t(:,k:k+1) * g;
  91.             t(k:k+1,:) = g' * t(k:k+1,:);
  92.             q(:,k:k+1) = q(:,k:k+1) * g;
  93.             temp = d(k); d(k) = d(k+1); d(k+1) = temp;
  94.          end
  95.       end
  96.    end
  97.  
  98.    v = v * q;
  99.  
  100. end
  101.  
  102. % End of C.R. Denham code.