function [X,degX,Y,degY]=diossr(A,degA,B,degB,C,degC,mes)

%  to solve X*A + Y*B = C over polynomial matrices
%  (state-space realization method)                                      
%                                                         
%  A,degA ... nonsingular polynomial matrix                
%  B,degB ... polynomial matrix                
%  C,degC ... polynomial matrix                
%  X,degX ... polynomial matrix   
%  Y,degY ... polynomial matrix   
%  mes (=0,1,2)...message and interrupt level 
%
% call: [X,degX,Y,degY]=diossr(A,degA,B,degB,C,degC,mes)

%  function used: colred, polmul, poladd, rmod, contr,rypol

[D,degD,U,degU,UI,degUI]=colred(A,degA,mes);
[Br,degBr]=polmul(B,degB,UI,degUI);
[Cr,degCr]=polmul(C,degC,UI,degUI);
[Ub,degUb,Vb,degVb]=rmod(D,degD,Br,degBr,mes);
[Uc,degUc,Vc,degVc]=rmod(D,degD,Cr,degCr,mes);

%  Vb*D^{-1}  ... strictly proper part of B*A{-1}
%  Vc*D^{-1}  ... strictly proper part of C*A{-1}  

[F,Hb,Hc]=contr(D,degD,Vb,degVb,Vc,degVc,mes);

%  controlable realizations 
%  F  ... state-feedback matrix corresp. to A^{-1}
%  Hb ... output matrix corresp. to B*A{-1}
%  Hc ... output matrix corresp. to C*A{-1}

[Y,degY]=rypol(F,Hb,Hc,mes);

if isempty(Y), return, end 

%  X from X*A+Y*B=C 

[Ub,degUb]=polmul(Y,degY,Ub,degUb);               
[X,degX]=poladd(-Ub,degUb,Uc,degUc);              
[Vb,degVb]=polmul(Y,degY,Vb,degVb);               
[Ubb,degUbb,Vbb,degVbb]=rmod(D,degD,Vb,degVb,mes);
[X,degX]=poladd(X,degX,-Ubb,degUbb);              

if mes>0
   disp('solution available in: X,degX,Y,degY '); 
   keyboard
end