home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / HQR10.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.8 KB  |  67 lines
  1. function [qn,tn,m,swap] = hqr10(a)
  2. %
  3. % [QN,TN,M,SWAP] = HQR10(A) produces an ordered Complex Schur decom-
  4. %          position such that the stable part is on the top of the 
  5. %          the diagonal of Tn, and the unstable part at the bottom,
  6. %     where
  7. %             Tn = Qn' * A * Qn 
  8. %              m = no. of stable poles.
  9. %           swap = total no. of swaps.
  10. %
  11. %         The complex Givens rotation is used to swap the two adjacent
  12. %    diagonal terms.
  13. %
  14.  
  15. % R. Y. Chiang & M. G. Safonov 9/11/87
  16. % Copyright (c) 1988 by the MathWorks, Inc.
  17. % All Rights Reserved.
  18. %----------------------------------------------------------------------
  19. %
  20. [ma,na] = size(a);
  21. % 
  22. % -------- Regular Schur Decomposition :
  23. %
  24. [q,t] = schur(a); 
  25. %
  26. [mt,nt] = size(t); 
  27. dt = diag(t);
  28. %
  29. % -------- Find the no. of stable roots :
  30. %
  31. idm = find(real(dt) < zeros(mt,1));
  32. [m,dum] = size(idm);
  33. if (m==0) | (m==mt), qn = eye(mt); tn = t; swap = 0; return, end
  34. %
  35. % -------- Assign I.C. for the loop counter :
  36. %
  37. pcnt = 0.; kcnt = 0.; gcnt = 0.;
  38. %
  39. % -------- Begin major loop :
  40. %
  41. for p = 1 : (nt-m)*m
  42.     pcnt = pcnt + 1;
  43.     for k = 1 : (nt-1)
  44.         kcnt = kcnt + 1;
  45.         tkk = t(k,k); tkp1 = t(k+1,k+1);
  46.         if (tkk>0.) & (tkp1<0.)
  47.            gcnt = gcnt + 1;
  48.            jb = eye(nt);
  49.            jb(k:k+1,k:k+1) = givens(t(k,k+1),t(k,k)-t(k+1,k+1));
  50.            t = jb' * t * jb;
  51.            q = q*jb;
  52.         end
  53.     end
  54.     dt = diag(t); ndt = dt(1:m,1);
  55.     flag = find(real(ndt)<zeros(m,1)); [fm,dum] = size(flag);
  56.     if fm == m & sum(size(q)) ~= 0
  57.        qn = q; tn = t; swap = gcnt; break
  58.     end
  59.     if sum(size(q)) == 0 | sum(size(t)) == 0
  60.        error('WARNING: the problem is ill-conditioned, no convergence!')
  61.     end   
  62. end
  63. %
  64. % -------- End of HQR10.M ---- RYC/MGS %
  65.  
  66.  
  67.