home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / SJH6.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.4 KB  |  44 lines
  1. % function u = sjh6(a,c)
  2. %
  3. %   Solves the Lyapunov equation  a'*x + x*a + c'*c = 0
  4. %   the solution returned is the Cholesky factor u where 
  5. %   u is upper triangular and x = u'*u
  6. %
  7. %   modified by K. Glover  June 1987.
  8. %   it is assumed that a is in upper Schur form on entry.
  9.  
  10. function u = sjh6(a,c)
  11. %  check if input has any complex inputs
  12.   if any(any(imag(a)))| any(any(imag(c))) realfl=0; 
  13.       else realfl=1; end
  14.   [na,ma]=size(a);
  15.   [mc,nc]=size(c);
  16.   n = na; m = mc;
  17. % perform complex Schur form on a
  18.   if realfl [ua,a] = rsf2csf(eye(n),a); end
  19. % qr decomposition of c
  20.   if realfl c=c*ua;end
  21.   c = triu(qr(c));
  22. % now calculate the solution from equations (5.13) - (5.16).
  23.   u = zeros(n);  id=eye(n); y=zeros(1,n-1);
  24.   for i=1:(n-1)
  25.     f = sqrt(-2*real(a(i,i)));
  26.     u(i,i) = abs(c(1,1))/f;
  27.     e = sign(c(1,1))*f; if c(1,1)==0 e=f;end
  28.     u(i,(i+1):n) = (-e'*c(1,2:(n-i+1)) -u(i,i)'*a(i,(i+1):n))/...
  29.         (a(i,i)'*id(1:(n-i),1:(n-i)) + a((i+1):n,(i+1):n));
  30.     y(1,1:(n-i)) = c(1,2:(n-i+1)) - e*u(i,(i+1):n);
  31.     if m>1  c(1:m,1:(n-i)) = triu(qr([c(2:m,2:(n-i+1));y(1,1:(n-i))]));
  32.               else c(1:m,1:(n-i)) = y(1:m,1:(n-i)); end
  33.     end
  34.   f = sqrt(-2*real(a(n,n)));
  35.   u(n,n) = abs(c(1,1))/f;
  36. % now obtain u in upper triangular form
  37.   if ~realfl u = triu(qr(u));, end;
  38.   if realfl u=u*ua'; u(1:n,1:n) = qrsc(u);end;
  39.  end
  40.  
  41.  
  42. %
  43. % Copyright MUSYN INC 1991,  All Rights Reserved
  44.