function [L] = lyapSL(R,C);
%
%
%     Usage:  [L] = lyapL(R,C);
%
%     Solves the Lyapunov equation
%
%         R' Q + Q R + C' C = 0
%
%     where R is real n by n quasi upper triangular matrix, 
%     and C is an m by n matrix.
%
%     Assumes eigenvalues of R are in open left half-plane.
%     Computes L lower triangular with positive diagonal elements
%     such that
%                Q = L L'
%
%     D.C. Sorensen
%     sorensen@rice.edu
%     8 Jan 00
%
%--------------------------------------------------------------
%
    [n,n] = size(R);
%
%   Be sure R is fully deflated 
%
    for j = 1:n-1,
        if (abs(R(j+1,j)) < 10*eps*max(abs(R(j,j)), abs(R(j+1,j+1)))),
            R(j+1,j) = 0;
        end
    end
%
    L = zeros(n);

%
%  Solve for L in the Schur basis.  
%
%  Uses structure  
%
%       inv(L)*R'*L  is lower triangluar
%
    j = 1;
    while(j < n+1),
       if (j < n), test = R(j+1,j); else test = 0; end
       if (test == 0),
          mu = R(j,j);
          c = C(:,j:n)'*C(:,j);
          L(j:n,j) = -(R(j:n,j:n) + mu*eye(n-j+1))'\c;
          lam = L(j,j); 
          if(lam > 0), 
             slam = sqrt(lam);
             L(j,j) = slam;
             L(j+1:n,j) = L(j+1:n,j)/slam; 
             C(:,j) = C(:,j)/slam;
             C(:,j+1:n) = C(:,j+1:n) - C(:,j)*L(j+1:n,j)';
          end
          j = j+1;
       else           
          [Y,D] = eig(R(j:j+1,j:j+1));
          Y = [real(Y(:,1)) imag(Y(:,1))];
          mu = real(D(1,1)); eta = imag(D(1,1));
          c = C(:,j:n)'*C(:,j:j+1); c = -c*Y;
          M = (R(j:n,j:n) + mu*eye(n-j+1))';
          c = M*c - eta*c*[0 1 ; -1 0];
          L(j:n,j:j+1) = (M*M + eta^2*eye(n-j+1))\c;
          L(j:n,j:j+1) = L(j:n,j:j+1)/Y;  
          L(j+1,j) = (L(j+1,j) + L(j,j+1))/2;
          L(j,j+1) = 0;
          lam = L(j,j); 
          if(lam > 0), 
             slam = sqrt(lam);
             L(j,j) = slam;
             L(j+1:n,j) = L(j+1:n,j)/slam; 
             L(j+1:n,j+1) = L(j+1:n,j+1) - L(j+1:n,j)*L(j+1,j); 
             C(:,j) = C(:,j)/slam;
             C(:,j+1:n) = C(:,j+1:n) - C(:,j)*L(j+1:n,j)';
          end
          j = j+1;
          lam = L(j,j); 
          if(lam > 0), 
             slam = sqrt(lam);
             L(j,j) = slam;
             L(j+1:n,j) = L(j+1:n,j)/slam; 
             C(:,j) = C(:,j)/slam;
             C(:,j+1:n) = C(:,j+1:n) - C(:,j)*L(j+1:n,j)';
          end
          j = j+1;
       end
    end

