function [U] = lyapU(R, B);
%
%     Usage:  [U] = lyapU(R,B);
%
%     Solves the Lyapunov equation
%
%         R P + P R' + B B' = 0
%
%     where R is upper triangular square matrix.
%
%     Assumes eigenvalues of R are in left half-plane (strictly)
%     Computes U upper triangular with positive diagonal elements
%     such that
%                P = U U'
%
%     D.C. Sorensen      (sorensen@rice.edu)
%     21 Dec 99
%     modified slightly by:   Y.  Zhou
%     Jan 2000
%
%---------------------------------------------------------------------
%
    [n,n] = size(R);
    U = zeros(n);

    for j = n:-1:2,
       mu1 = sqrt(-2*real(R(j,j)));
       b = B(j,:); mu = norm(b);
       if (mu > 0),
          b = b/mu;
          btmp = B(1:j-1,:)*(b'*mu1) + R(1:j-1,j)*mu/mu1;
          u =  -(R(1:j-1,1:j-1) + (R(j,j)')*eye(j-1)) \ btmp;
          B(1:j-1,:) = B(1:j-1,:) - u*b*mu1;
       else
          u = zeros(j-1,1);
       end
       U(j,j) = mu/mu1;
       U(1:j-1,j) = u;
    end
    U(1,1) = norm(B(1,:))/sqrt(-2*real(R(1,1)));

% end-function-lyapU.m