gsylvester(A,B,C,D,E,F,version)

function [X, Y, info] = gsylvester(A,B,C,D,E,F,version)
%  PURPOSE: solves the generalized Sylvester equations:
%
%              A * R   -   L * B            = scale * C            (1)         
%              D * R   -   L * E            = scale * F
%
%              A * X * E  +  D * X * B      = scale * C            (2)
%
%       where R and L are unknown m-by-n matrices, (A, D), (B, E) and
%       (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
%       respectively, with real entries. 0 <= scale <= 1 is an output
%       scaling factor chosen to avoid overflow.
%
%  USAGE: [X,Y,info] = gsylvester_schur(A,B,C,D,E,F,version)
%  where: 
%         A,B,C,D,E,F               matrices
%         version                   equation to solve
%                                   = 1 - solves equation (1)
%                                   = 2 - solves equation (2), matrix F is
%                                           not reffered
%
%         X, Y                      matrices forming solution, if version = 2
%                                   then Y = 0
%         info                      = 0: successful exit
%                                   <0: If info = -i, the i-th argument had an illegal value.
%                                   1: (A, D) and (B, E) have common or close eigenvalues.
%                                   2:
%                                   3: system is ill-conditioned, solution
%                                   may be inaccurate
%
%   COMMENTS:
%
% Copyright  (c) Pawel Kowal (2006)
% All rights reserved
% LREM_SOLVE toolbox is available free for noncommercial academic use only.
% [e-mail]

% check size
m       = size(A,1);
n       = size(B,1);
if size(A,2)~=m || size(D,1)~=m ||size(D,2)~=m || size(C,1)~=m
    error('Inconsistent matrix size');
end
if size(B,2)~=n || size(E,1)~=n ||size(E,2)~=n || size(C,1)~=m
    error('Inconsistent matrix size');
end
if version==1
    if size(F,1)~=m || size(F,2)~=n
        error('Inconsistent matrix size');
    end
end

switch version
    case 1
        [P,Q,TA,TD]             = gschur(A,D,0);
        [U,V,TB,TE]             = gschur(B,E,0);

        [X, Y, scale, info]     = gsylvester_schur(TA,TB,P'*C*V,TD,TE,P'*F*V);
        if scale==0
            Y                   = P*Y*U';
            X                   = Q*X*V';
            info                = 2;
        else
            Y                   = P*Y*U'/scale;
            X                   = Q*X*V'/scale;
        end
    case 2
        [P,Q,TA,TD]             = gschur(A,-D,0);
        [U,V,TB,TE,AR,AI,BETA]  = gschur(B,E,0);
        F                       = zeros(size(D,1),size(E,2));

        [X, Y, scale, info]     = gsylvester_schur(TA,TB,P'*C*V,TD,TE,F);
        Y                       = P*Y*U';
        X                       = Q*X*V';
        
        eig1                    = min(abs(AR+i*AI));
        eig2                    = min(abs(BETA));
        
        if eig1<eig2
            tol                 = 100 *max(size(E)')*norm(E,1)*eps;
            if eig2<tol
                info            = 3;
            end                
            X                   = (E'\(X)')';
            Y                   = 0;
        else
            tol                 = 100 *max(size(D)')*norm(D,1)*eps;
            if eig1<tol
                info            = 3;
            end                            
            X                   = -(D\Y);
            Y                   = 0;            
        end
        if scale==0
            info                = 2;
        else
            X                   = X/scale;
        end        
end