function [D,degD,U,degU,V,degV,UI,degUI,VI,degVI]=poldiagr(A,degA,mes)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB function to find a diagonal form D of a polynomial    %
% matrix A along with corresponding unimodular matrices        %
% U and V:                                                     %  
%                   A=U*D*V                                    %
%                                                              %
% UI=U^-1, VI=V^-1,                                            %
% mes (=1,2,3) ... message and interrupt level                 %
% row-column pivotting                                         % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M. Sebek, September 1990                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% at the moment works just for a square matrix 

% function used: polen, polcheck, polmatcl, polmul, polseten,
%		 polsize.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% preliminaries                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('..................................................');
disp('Start of diagonalization');
disp('..................................................');

[rD,cD]=polsize(A,degA);
m=min(rD,cD);

[D,degD,DEGD,LCD]=polmatcl(A,degA);

			    if mes>0,D,DEGD,LCD,keyboard,end

U=eye(rD); UI=U; degU=0;degUI=0;
V=eye(cD); VI=V; degV=0; degVI=0;

ner=0;                       % number of elementary reductions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% diagonalization                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:m-1
			    if mes>0,i,end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% zeroing the i-th row and column                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  while max(DEGD(i,i+1:cD))>=0 | max(DEGD(i+1:rD,i))>=0

    ner=ner+1;
    TU=eye(rD); TUI=TU; degTU=0;degTUI=0;
    TV=eye(cD); TVI=TV; degTV=0; degTVI=0;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % finding the "pivot element" in row or column i        %                         
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    
    d=abs(DEGD(i,i));
    if d>max(DEGD(i,i+1:cD)) & d>max(DEGD(i+1:rD,i))
      d=d-1;
    end

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % row i                                                 %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    L=find(abs(DEGD(i,i:cD))<=d);

    if ~isempty(L)
      L=L+(i-1)*ones(1,length(L));
      [ML,LL]=sort(abs(LCD(i,L)));
      ML=ML(length(ML));
      l=L(LL(length(LL)));
    
      if DEGD(i,i)==DEGD(i,l) & abs(LCD(i,i))==ML
         l=i;
      end
    else
      ML=0;
    end 
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % column i                                              %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    K=find(abs(DEGD(i:rD,i))<=d);
    
    if ~isempty(K)
      K=K+(i-1)*ones(length(K),1);
      [MK,KK]=sort(abs(LCD(K,i)));
      MK=MK(length(MK));
      k=K(KK(length(KK)));
    
      if DEGD(i,i)==DEGD(k,i) & abs(LCD(i,i))==MK
        k=i;
      end
    else
      MK=0;
    end
    
    
    
    if MK>ML 

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % interchanging rows i and k                           %
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      if k~=i 
       PU=eye(rD);
       PU(i,i)=0; PU(k,k)=0; PU(i,k)=1; PU(k,i)=1;

       D=PU*D;  LCD=PU*LCD;
       b=DEGD(i,:);DEGD(i,:)=DEGD(k,:);DEGD(k,:)=b;
       U=polmul(U,degU,PU,0); UI=PU*UI;
       ner=ner+1;
      end

    else

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % interchanging columns i and l                       %
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      if l~=i
       PV=eye(cD);
       PV(i,i)=0; PV(l,l)=0; PV(i,l)=1; PV(l,i)=1;

       D=polmul(D,degD,PV,0);  LCD=LCD*PV;
       b=DEGD(:,i);DEGD(:,i)=DEGD(:,l);DEGD(:,l)=b;
       V=PV*V; VI=polmul(VI,degVI,PV,0);
       ner=ner+1;
      end

    end				

				  if mes>0,DEGD,keyboard,end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % reduction of row i                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for j=i+1:cD
      if DEGD(i,j)>=0
	c=LCD(i,j)/LCD(i,i);
	p=DEGD(i,j)-DEGD(i,i);
			 if mes>1,i,j,c,p,keyboard,end
	if p==0
	  s=c;
        else
	  s=[zeros(1,p) c];
        end
	[TV,degTV]=polseten(TV,degTV,s,p,i,j);
	[TVI,degTVI]=polseten(TVI,degTVI,-s,p,i,j);
                         if mes>1,   TV,TVI,keyboard,end
        ner=ner+1;
      end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % reduction of column i                              %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for j=i+1:rD
      if DEGD(j,i)>=0
	c=LCD(j,i)/LCD(i,i);
	p=DEGD(j,i)-DEGD(i,i);
			  if mes>1,i,j,c,p,keyboard,end
	if p==0
	  s=c;
        else
	  s=[zeros(1,p) c];
        end
	[TU,degTU]=polseten(TU,degTU,s,p,j,i);
	[TUI,degTUI]=polseten(TUI,degTUI,-s,p,j,i);
                              if mes>1,TU,keyboard,end
        ner=ner+1;
      end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % application of T...                                %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [D,degD]=polmul(TUI,degTUI,D,degD);
    [D,degD]=polmul(D,degD,TVI,degTVI);
    [D,degD,DEGD,LCD]=polmatcl(D,degD);
    
    [U,degU]=polmul(U,degU,TU,degTU);
    [UI,degUI]=polmul(TUI,degTUI,UI,degUI);

    [V,degV]=polmul(TV,degTV,V,degV);
    [VI,degVI]=polmul(VI,degVI,TVI,degTVI);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % end of one step of reductions                      %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      [U,degU]=polmatcl(U,degU);
      [UI,degUI]=polmatcl(UI,degUI);
      [V,degV]=polmatcl(V,degV);
      [VI,degVI]=polmatcl(VI,degVI);
			   
   if mes>0
     DEGD
     e=polcheck(A,degA,U,degU,D,degD,V,degV);
     disp(sprintf('relative error after %g reductions is %g',ner,e));
     keyboard
   end
  
  
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % end of reduction of row and column i                 %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% last row                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if cD>m

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % "pivoting"                                          %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  d=abs(DEGD(rD,rD));
  if d>max(DEGD(rD,rD+1:cD))
    d=d-1;
  end
  L=find(DEGD(rD+1:cD)<=d);
  L=L+rD*ones(1,length(L));
  [a,LL]=sort(abs(LCD(rD,L)));
  l=L(LL(length(LL)));

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % reduction of the last row                           %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  % to be completed


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% last column                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if rD>m
  % to be completed
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('..................................................');
disp(sprintf('Diagonalization finished after %g reductions',ner));
e=polcheck(A,degA,U,degU,D,degD,V,degV);
disp(sprintf('                relative error %g ',e));
disp('..................................................');