function [X,E]=vermatreqn(A,B,C,D,F)
%    VERMATREQN       Verified solution of the matrix equation A*X*B+C*X*D=F
%                     (in particular, of the Sylvester or Lyapunov equation).
%
%    This is an INTLAB file. It requires to have INTLAB installed under
%    MATLAB to function properly.
%
%    For interval (or real) m-by-m matrices A, C, p-by-p matrices B, D, and
%    an m-by-p matrix F,
%        [X,E]=vermatreqn(A,B,C,D,F)
%    returns a verified m-by-p solution X of the matrix equation
%        A*X*B+C*X*D=F,
%    or yields no verified result.
%
%    Possible outcomes:
%
%    ~isnan(X.inf(1,1)) :  it is verified that for each Ao in A, Bo in B,
%                          Co in C, Do in D, and Fo in F, the equation
%                             Ao*Xo*Bo+Co*Xo*Do=Fo
%                          possesses a unique matrix solution Xo which is
%                          verified to be contained in X,
%     isnan(X.inf(1,1)) :  no verified result (the interval matrix X
%                          consists of NaN's).
%
%    The structure E explains reasons for NaN output.
%
%    For a verified solution of the SYLVESTER equation
%        A*X+X*D=F
%    (A, D square of possibly different sizes), use
%        [X,E]=vermatreqn(A,eye(size(D)),eye(size(A)),D,F).
%
%    For a verified solution of the LYAPUNOV equation
%        A*X+X*A'=F
%    (A, F square of the same size), use
%        [X,E]=vermatreqn(A,eye(size(A')),eye(size(A)),A',F).
%
%    See also VERIFYLSS.

%    Copyright 2008 Jiri Rohn
%
%    Based on solving the square system of linear equations
%        (kron(B',A)+kron(D',C))*vec(X)=vec(F),
%    where vec(X)=X(:) and vec(F)=F(:); see R. A. Horn and C. R. Johnson,
%    Topics in Matrix Analysis, Cambridge University Press, Cambridge 1991,
%    Lemma 4.3.1.
%
%    This work was done during author's employment at the Anglo-American
%    University in Prague, Czech Republic.
%
%    WARRANTY
%
%    Because the program is licensed free of charge, there is
%    no warranty for the program, to the extent permitted by applicable
%    law. Except when otherwise stated in writing the copyright holder
%    and/or other parties provide the program "as is" without warranty
%    of any kind, either expressed or implied, including, but not
%    limited to, the implied warranties of merchantability and fitness
%    for a particular purpose. The entire risk as to the quality and
%    performance of the program is with you. Should the program prove
%    defective, you assume the cost of all necessary servicing, repair
%    or correction.
%
%    History
%
%    2008-05-16   first version
%    2008-05-27   version for posting
%
gr=getround;
setround(0);
[m,n]=size(A); [p,q]=size(B);
% defaults
X=repmat(infsup(NaN,NaN),n,p);
E.error='vermatreqn: none';
E.where='NaN';
E.value='NaN';
% data check
if ~(nargin==5&&nargout<=2&&m==n&&p==q&&all(size(A)==size(C))&&all(size(B)==size(D))&&all(size(F)==[m,q]))
    E.error='vermatreqn: wrong data';
    setround(gr); return
end
% equation:  A*X*B+C*X*D=F % A, C: mxm, B, D: pxp
% solved as: (kron(B',A)+kron(D',C))*vec(X)=vec(F) % see Horn and Johnson, Topics ..., Lemma 4.3.1
% the system matrix is mpxmp
AK=kronmod(B',A)+kronmod(D',C);
if issparse(AK)
    AK=full(AK); % sparse matrices not implemented in verifylss
end
bK=F(:);
xK=verifylss(AK,bK); % matrix in vector form
X=reshape(xK,n,p);   % reshaped back
if isnan(X.inf(1,1)) % no result
    E.error='vermatreqn: verifylss failed';
end
setround(gr);
%
%%%%%%%%%%%%%%%%%%%%%
% Subfunction kronmod
%%%%%%%%%%%%%%%%%%%%%
%
function K=kronmod(A,B)
% Kronecker product of matrices of type intval
if ~isintval(A)
    A=infsup(A,A);
end
if ~isintval(B)
    B=infsup(B,B);
end
% the following five lines (construction of the Kronecker product)
% copied from KRON.M, Copyright 1984-2004 The MathWorks, Inc.
% (but the last product is performed in interval arithmetic)
[ma,na]=size(A);
[mb,nb]=size(B);
[ia,ib]=meshgrid(1:ma,1:mb);
[ja,jb]=meshgrid(1:na,1:nb);
K=A(ia,ja).*B(ib,jb);