home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / REIG.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.6 KB  |  61 lines
  1. function [xr,dd] = reig(a,opt)
  2. %
  3. % [XR,DD] = REIG(A,OPT) produces a "real" and "ordered" eigenstructure 
  4. %    decomposition such that    
  5. %                       -1 
  6. %                     XR  * A * XR = DD
  7. %    where 
  8. %               XR is a set of "real" eigenvectors that span the 
  9. %                  same eigenspace as the complex ones.
  10. %               DD is a real block-diagonal matrix with real eigenvalue
  11. %                  (1x1) and complex eigenvalue (2x2) on the main
  12. %                  diagonal.
  13. %   
  14. %    Two options are available: 
  15. %
  16. %                    Opt = 1 - - ordered by real parts (default)
  17. %                    Opt = 2 - - ordered by magnitude
  18.  
  19. % R. Y. Chiang & M. G. Safonov 8/88
  20. % Copyright (c) 1988 by the MathWorks, Inc.
  21. % All Rights Reserved.
  22. % ----------------------------------------------------------------
  23. %
  24.  
  25. if nargin == 1
  26.    opt = 1;
  27. end
  28.  
  29. if norm(imag(a),1)~=0, error('IMAG(A)~=0'),end
  30. [n,n] = size(a);
  31. dd = zeros(n,n); 
  32. xr = zeros(n,n);
  33. [x,d] = eig(a);
  34. dg = diag(d);
  35. if opt == 1
  36.    [dsort,index] = sort(real(dg));
  37. else
  38.    [dsort,index] = sort(abs(dg));
  39. end
  40. dgs = dg(index);
  41. xs = x(:,index);
  42. j=sqrt(-1);
  43. [junk,ixmax] = max(xs);
  44. i=1;
  45. for k=1:n,
  46.    if abs(imag(dgs(k))) <= 1.e-7*abs(real(dgs(k)))+eps,
  47.       dd(i,i) = real(dgs(k));
  48.       xr(:,i) = real(xs(:,k)*exp(-j*angle(xs(ixmax(k),k))));
  49.       i=i+1;
  50.    else 
  51.       if imag(dgs(k)) > 1.e-7*abs(real(dgs(k)))+eps,
  52.          dd(i:i+1,i:i+1) = ...
  53.             [real(dgs(k)) imag(dgs(k));-imag(dgs(k)) real(dgs(k))];
  54.          xr(:,i:i+1) = [real(xs(:,k)) imag(xs(:,k))];
  55.          i=i+2;
  56.       end
  57.    end
  58. end
  59. %
  60. % ------ End of REIG.M % RYC/MGS %
  61.