function [U] = lyapUR(R,B);
%
%     Usage:  [U] = lyapUR(R,B);
%
%     Solves the Lyapunov equation
%
%         R P + P R' + B B' = 0
%
%     where R is real n by n quasi upper triangular matrix,
%     and B is an n by p real matrix.
%
%     Assumes eigenvalues of R are in open left half-plane.
%     Computes U upper triangular with positive diagonal elements
%     such that
%                P = U U'
%
%     D.C. Sorensen     (sorensen@rice.edu)
%     Y.   Zhou         (ykzhou@rice.edu)
%     August  2001
%
%--------------------------------------------------------------
%
    [n,n] = size(R);
    U = zeros(n);

    j = n;
    while(j > 0),
       if (j > 1), test = R(j,j-1); else test = 0; end
       if (test == 0),
          mu = R(j,j);
          b = B(1:j,:)*B(j,:)';
          U(1:j,j) = -(R(1:j,1:j) + mu*eye(j))\b;
          lam = U(j,j);
          if(lam > 0),
             slam = sqrt(lam);
             U(j,j) = slam;
             U(1:j-1,j) = U(1:j-1,j)/slam;
             B(j,:) = B(j,:)/slam;
             B(1:j-1,:) = B(1:j-1,:) - U(1:j-1,j)*B(j,:);
          end
          j = j-1;
       else
          r11 = R(j-1,j-1);    r21 = R(j-1,j);
          r12 = R(j,j-1);      r22 = R(j,j);
          b = B(1:j,:)*B(j-1:j,:)'; b = -b;
          M = R(1:j,1:j);  I =eye(j);
          M2=M*M;
          btmp = [M*b(:,1) + r22*b(:,1) - r21*b(:,2), ...
                  M*b(:,2) + r11*b(:,2) - r12*b(:,1)];
          U(1:j,j-1:j) = (M2+ (r11+r22)*M +(r11*r22-r12*r21)*I ) \btmp;
          U(j-1,j) = (U(j-1,j) + U(j,j-1))/2;
          U(j,j-1) = 0;
          lam = U(j,j);
          if(lam > 0),
             slam = sqrt(lam);
             U(j,j) = slam;
             U(1:j-1,j) = U(1:j-1,j)/slam;
             U(1:j-1,j-1) = U(1:j-1,j-1) - U(1:j-1,j)*U(j-1,j);
             B(j,:) = B(j,:)/slam;
             B(1:j-1,:) = B(1:j-1,:) - U(1:j-1,j)*B(j,:);
          end
          j = j-1;
          lam = U(j,j);
          if(lam > 0),
             slam = sqrt(lam);
             U(j,j) = slam;
             U(1:j-1,j) = U(1:j-1,j)/slam;
             B(j,:) = B(j,:)/slam;
             B(1:j-1,:) = B(1:j-1,:) - U(1:j-1,j)*B(j,:);
          end
          j = j-1;
       end
    end

% end-function-lyapUR.m