function [X] = sylvester(A1, A2, D);
%
%     Usage:  [X] = sylvester(A1, A2, D);
%
%     Solves the Sylvester equation
%
%         A1 X + X A2 + D = 0
%
%     where A1(m by m), A2(n by n), D(m by n) are all real matrices.
%
%     D.C. Sorensen     (sorensen@rice.edu)
%     Y. Zhou           (ykzhou@rice.edu)
%     August, 2001
%
%--------------------------------------------------------------------
%
    [m] = size(A1, 1);
    [n] = size(A2, 1);

    if (m==n & A1 == A2),    %if it is special Sylvester eq. AX + XA + D=0
       [X]=sylv1(A1, D);  return;
    end

    [Q1,R1] = schur(A1);

    if (m==n & A1 == A2')    %if it is Lyapunov eq. AX + XA' + D = 0
                             %get the Schur decomp of A2=A1' from Q1, R1
        dn = [size(A1,1) : -1 : 1];
        Q2 = Q1(:, dn);
        R2 = R1(dn, dn)';
    else                     %have to do Schur decomp of A2 in this case
        [Q2,R2] = schur(A2);
    end

    D = Q1'*D*Q2;

    X = zeros(m,n);  I = eye(m);
    for j = 1:n-1,
      if (abs(R2(j+1,j)) < 10*eps*max(abs(R2(j,j)), abs(R2(j+1,j+1)))),
          R2(j+1,j) = 0;
      end
    end

%
% save R1*R1 to a new matrix Rsq since R1*R1 may be computed
% many times if there are several size-2 diagonal blocks
%
    Rsq = R1 * R1;

%
%  Solve for X in the Real Schur basis.
%
    j = 1;
    while(j < n+1),
       if (j < n), test = R2(j+1,j); else test = 0; end
       if (test ==0),
          mu = R2(j,j);
          b = -( D(:,j) + X(:,1:j-1)*R2(1:j-1,j) );
          X(:,j) = (R1 + mu*I)\b;
          j = j+1;
       else
          r11 = R2(j,j);     r12 =R2(j,j+1);
          r21 = R2(j+1,j);   r22 =R2(j+1,j+1);
          b = -( D(:,j:j+1) + X(:,1:j-1)*R2(1:j-1,j:j+1) );
          btmp =  [R1*b(:,1) + r22*b(:,1) - r21*b(:,2), ...
                   R1*b(:,2) + r11*b(:,2) - r12*b(:,1)];
          X(:,j:j+1) =( Rsq + (r11+r22)*R1 +(r11*r22-r12*r21)*I )\btmp;
          j = j+2;
       end
    end

%
%   convert solution to standard basis
%
    X = Q1*(X*Q2');

% end-function-sylvester.m


%%=========================================================================
%% function sylv1.m is a subroutine of the above Sylvester solver,
%% calling this subroutine in case (A1 == A2) saves the use of R2, Q2
%%=========================================================================

function [X] = sylv1(A, D);
%
%     Usage:  [X] = sylv1(A, D);
%
%     Solves the Sylvester equation
%
%         A X + X A + D  = 0
%
%     where A(n by n), D(n by n) are real matrices.
%
%     D.C. Sorensen     (sorensen@rice.edu)
%     Y. Zhou           (ykzhou@rice.edu)
%     August, 2001
%
%----------------------------------------------------------------------
%
    [Q,R] = schur(A);

    D = Q'*D*Q;
    [n,n] = size(A);
    X = zeros(n);  I = eye(n);
    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

    Rsq = R * R;
%
%  Solve for X in the Real Schur basis.
%
    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);
          b = D(:,j) + X(:,1:j-1)*R(1:j-1,j); b = -b;
          X(:,j) = (R + mu*I)\b;
          j = j+1;
       else
          r11 = R(j,j);     r12 =R(j,j+1);
          r21 = R(j+1,j);   r22 =R(j+1,j+1);
          b = -( D(:,j:j+1) + X(:,1:j-1)*R(1:j-1,j:j+1) );
          btmp = [R*b(:,1) + r22*b(:,1) - r21*b(:,2), ...
                  R*b(:,2) + r11*b(:,2) - r12*b(:,1)];
          X(:,j:j+1) =( Rsq + (r11+r22)*R +(r11*r22-r12*r21)*I )\btmp;
          j = j+2;
       end
    end

%
%   convert solution to standard basis
%
    X = Q*(X*Q');

% end-function-sylv1.m