function U = sylvesterschur(A,B,C)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function U = sylvesterschur(A,B,C)
% Given nxn input matrices A, B, and C, this function solves
% the equation
%       AU + UB = C
% for the nxn matrix U.
%
% The algorithm is to compute a Schur decomposition of A' and of B
% as
%    A = W L W',       B = Y R Y'
% where L is lower triangular, R is upper triangular, and W and Y
% are unitary.
%
% Then we transform the problem to
%    L Us + Us R = Chat
% where Chat = W' C Y.  This can be solved by substitution.
% 
% Finally, we compute U = W Us Y'.
%
% This algorithm is due to Bartels and Stewart.
%
% Its disadvantage is that the Schur decomposition is complex if the
% matrix has any complex eigenvalues, so the returned matrix U may
% have small imaginary part even if all of the input matrices are real.
%
% The solution fails to exist if L(i,i) + R(j,j) = 0 for some value
% of i and j.  If such a value is nearly zero, the function returns
% an error message with U=0.
% sylvesterschurr.m  Dianne O'Leary 08/2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1: compute a Schur decomposition of A' and of B.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   [n,m]=size(C);
   [W,R] = schur(A','complex');
   L = R';
   [Y,R] = schur(B,'complex');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The problem is transformed to L Us + Us R = Chat where Chat = W' C Y.  
% Step 2: Solve the transformed problem by substitution, producing an 
% error message if the divisor becomes too small.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   Us = zeros(n,n);
   Chat = W'*C*Y;
   a = max(abs([diag(R);diag(L)]));
   for i=1:n,
     for j=1:n,
       if (abs(L(i,i)+R(j,j))) > 1000*eps*a
         Us(i,j) = (Chat(i,j)-L(i,1:i-1) *Us(1:i-1,j)...
                             -Us(i,1:j-1)*R(1:j-1,j))/(L(i,i)+R(j,j));
       else
         disp(sprintf('Problem very ill-conditioned or no solution exists.'))
         U = zeros(n,n);
         return
       end
     end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 3: compute U = W Us Y'.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   U = W*Us*Y';