function f=dio_ts4(am,degA,bn,degB,k)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Test  MATLAB program for solution     %
%%   the Diophantine Equation in polynomial%
%%   matrices by "dioel",  "dioin" and     %
%%   "diossl". The equation is             %
%%             AX + BY = C                 %
%%   where  dimensions of A are (am x am)  %
%%          and degA is degree of A.       %
%%          Dimensions of B are (am x bn), %
%%          degB is degree of B.           %
%%          C is identity matrix (am x am).%
%%                                         %
%%          k is number of experiment      %
%%                                         %
%%   CALL: dio_ts4(am,degA,bn,degB,k)      %
%%   FUNCTION used: dioel,dioin,           %
%%                  diossl                 %
%%   SKIP one method: skip=1               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   S. Pejchova, November 18th  1994      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=int2str(k);
text1=['diodi',s];
text2=['diosa',s];
text2=['save ',text2];
text2=[text2,' A B C degA degB degC ALPHA SJ N1 N2 N3 time1 time2 time3 deg1 deg2 deg3 f1 f2 f3 cl'];
text3=['            ',text1];
text1=['diary ',text1];
eval(text1);
cl=fix(clock);
disp('......................................................');
disp(text3);
disp('......................................................');
disp('            DATE  OF  EXPERIMENT ');
disp(cl(1:3));

A=rand(am,am*(degA+1));
for i=1:am
    for j=1:(am*(degA+1))
        if A(i,j) <= 1e-8
           A(i,j)=0;
        end
     end
end

B=rand(am,bn*(degB+1));
for i=1:am
    for j=1:(bn*(degB+1))
        if B(i,j) <= 1e-8
           B(i,j)=0;
        end
    end
end

C=eye(am);
degC=0;

disp('............Equation is   A X  +  B Y  =  C   ........');
disp('                      ');
disp(sprintf('            MATRIX  A (%g x %g), degA = %g',am,am,degA));
disp(sprintf('            MATRIX  B (%g x %g), degB = %g',am,bn,degB));
disp('......................................................');
No=norm([A B],inf);

[AT,degAT]=transp(A,degA);
[BT,degBT]=transp(B,degB);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   The solution of the Diophantine           %
%   equation by  "diossl"-state-space realiz. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

skip=0;
keyboard;


N2=0;
if skip==0

ms=0;
       flops(0);
       t1=clock;
[X2,degX2,Y2,degY2]=diossl(A,degA,B,degB,C,degC,ms);
       t2=clock;
       f2=flops;
       [A2,degA2]=polmul(A,degA,X2,degX2);
       [B2,degB2]=polmul(B,degB,Y2,degY2);
       [L2,degL2]=poladd(A2,degA2,B2,degB2);
       [LA2,degLA2]=poladd(L2,degL2,-C,degC);
       N2=(norm(LA2,inf))/No;

       time2=etime(t2,t1);
       deg2=max(degX2,degY2);
       disp(sprintf('KRAFFER(SSM-diossl)  -  error=%g,      time=%g,      degree=%g',N2,time2,deg2));
       disp(sprintf('                        flops=%g',f2));
end

if skip~=0
   deg2=max(degA,degB);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   The solution of the Diophantine       %
%   equation  by  "dioin"-interpolation   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

skip=0;
keyboard;


N1=0;
ALPHA=[]; SJ=[];
if skip==0
      kp=am*degA+am*(deg2+1);      %Number of interpolation points.
      ALPHA=rand(am,kp);
      SJ1=(-kp/20:0.1:kp/20);
      SJ=SJ1(1,1:kp);

      flops(0);
      t1=clock;
[XT,degXT,YT,degYT]=dioin(AT,degAT,BT,degBT,C,degC,ALPHA,SJ,deg2);
      t2=clock;
      f1=flops;
      [X1,degX1]=transp(XT,degXT);
      [Y1,degY1]=transp(YT,degYT);

      [A1,degA1]=polmul(A,degA,X1,degX1);
      [B1,degB1]=polmul(B,degB,Y1,degY1);
      [L1,degL1]=poladd(A1,degA1,B1,degB1);
      [LA1,degLA1]=poladd(L1,degL1,-C,degC);
      N1=(norm(LA1,inf))/No;

      time1=etime(t2,t1);
      deg1=max(degX1,degY1);
      disp(sprintf('ANTSAKLIS(PIM-dioin) -  error=%g,      time=%g,      degree=%g',N1,time1,deg1));
      disp(sprintf('                        flops=%g',f1));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    The solution of the Diophantine      %
%    equation by   "dioel"-elementary op. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

skip=0;
keyboard;


N3=0;
if skip==0

      flops(0);
      t1=clock;
[X3,degX3,Y3,degY3,R,degR,S,degS,tm,tn]=dioel(A,degA,B,degB,C,degC);
      t2=clock;
      f3=flops;
      [A3,degA3]=polmul(A,degA,X3,degX3);
      [B3,degB3]=polmul(B,degB,Y3,degY3);
      [L3,degL3]=poladd(A3,degA3,B3,degB3);
      [LA3,degLA3]=poladd(L3,degL3,-C,degC);
      N3=(norm(LA3,inf))/No;

      time3=etime(t2,t1);
      deg3=max(degX3,degY3);
      disp(sprintf('KUCERA(EOM-dioel)    -  error=%g,      time=%g,      degree=%g',N3,time3,deg3));
      disp(sprintf('                        flops=%g',f3));
end
disp('......................................................');

disp('Matrix  A = A0 + A1*s + A2*s^2 + ... + AdegA*s^(degA)');
for i=1:(degA+1)
    disp(sprintf('              A%g',i-1));
    disp(A(:,(i-1)*am+1:(i-1)*am+am));
end

disp('                           ');
disp('Matrix  B = B0 + B1*s + B2*s^2 + ... + BdegB*s^(degB)');
for i=1:(degB+1)
    disp(sprintf('               B%g',i-1));
    disp(B(:,(i-1)*bn+1:(i-1)*bn+bn));
end

disp('                            ');
    disp('                C');
    disp(C);

if (isempty(ALPHA)~=1)|(isempty(SJ)~=1)

disp('                            ');
     disp('Interpolation points');
     disp('                    ');
     disp('                ALPHA');
     disp(ALPHA);
     disp('                    ');
     disp('                SJ');
     disp(SJ);
end
disp('-------------------------------------------------------------------------');
eval(text2);
diary off