function  [C,degC,K] = polsf2(A,degA,mes)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB FUNCTION polsf2 FOR SPECTRAL FACTORIZATION OF A   %
% POLYNOMIAL MATRIX                                        %
%                                                          %
%    A=C~*K*C                                              %
%                                                          %
% mes (=0,1,2) ... message  and interrupt level            %
%                                                          %
% version with singular part extraction                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% functions used: adj, coljoin, dispoldi, poladd, polcheck,
%		  poldiagr, poldiv, polen, polgld,polmatcl, 
%                 polmul,
%		  polnr, polprenr, polseten, polsize, 
%		  poltrf, poluf, rowjoin;

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



  disp('..................................................');
  disp('Start of spectral factorization');
  disp('..................................................');


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

[rA,n]=polsize(A,degA);
if rA~=n
  error('polsf: matrix is not square');
end
rank=n;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% zero matrix                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if degA <= 0 & A(:,1:rA)==zeros(rA,rA)
  K=A;
  C=eye(rA);
  degC=0;
  return
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% greatest common left divisor                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  disp('..................................................');
	           disp('Start of singular part extraction');
  disp('..................................................');

  [AP,degAP,UI,degUI,U,degU]=polgld(A,degA,mes-1);
  [AB,degAB]=polmul(adj(UI,degUI),degUI,AP,degAP);
  [AB,degAB,DEGAB]=polmatcl(AB,degAB);
  [rAB,rank]=polsize(AB,degAB);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% singular part extraction                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if rank<n

  while max(DEGAB(rAB,:))<0
     DEGAB(rAB,:)=[];
     AB(rAB,:)=[];
     rAB=rAB-1;
  end
           
                Atest=[AB;zeros(n-rank,rank*(degAB+1))];
                Atest=coljoin(Atest,degAB,zeros(n,n-rank),0);
	e=polcheck(A,degA,adj(U,degU),degU,Atest,degAB,U,degU);
disp(sprintf('Singular part extracted with rel. error %g',e));

                               if mes >0
                                    AB,degAB,keyboard
				 if mes>1, U,keyboard,end
                               end
  disp('..................................................');


else

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A is non-singular                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   AB=A; degAB=degA;
  [rAB,rank]=polsize(A,degA);
  U=eye(n);degU=0;
  UI=U;degUI=0;
  disp('Matrix non-singular - extraction cancelled');
  disp('..................................................');

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% diagonal form, AB=V*D*W, WI=W^(-1)                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  

[D,degD,V,degV,W,degW,VI,degVI,WI,degWI]=poldiagr(AB,degAB,mes-1);

                         Dif=poladd(D,degD,-adj(D,degD),degD);
                               e=norm(Dif,inf)/norm(D,inf);

   disp(sprintf('Relative error of Hermitian symmetry %g',e));
  disp('..................................................');

                                

  				      if mes>0
                                         keyboard
					 dispoldi(D,degD);
					 keyboard
                                         if mes>1
		                	    V,W,keyboard
				            WI,keyboard
                                         end
                                      end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% scalar spectral factorizations                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  disp('..................................................');
                      disp('Start of scalar factorizations');

F=zeros(rank,rank); degF=0;
c=zeros(1,rank);
for i=1:rank
                                       disp(sprintf('%g)',i));
  [d,degd]=polen(D,degD,i,i);
                               if mes>1,i,d,degd,keyboard,end
  if degd>=0

    [d,degd,t]=polprenr(d,degd);
			       if mes>1,d,t,keyboard,end
    [f,degf,sig]=polnr(d,degd,mes-1);
                             if mes>1,f,degf,sig,keyboard,end
    [f,degf]=poltrf(f,degf,1/t);
    sc=f(degf+1);
    c(i)=sig*sc^2;
    f=f/sc;
    [F,degF]=polseten(F,degF,f,degf,i,i);

  else

    F(i,i)=1;

  end

		  if mes>1, polen(F,degF,i,i),keyboard,end
end

V=polmul(V,degV,diag(c),0);

                                           if mes>0
		                              dispoldi(F,degF); 
                                           end

[FS,degFS]=adj(F,degF);

                                           if mes>1
		                             FS,keyboard
		                             W,WI,keyboard
		                           end

                   e=polcheck(D,degD,FS,degFS,diag(c),0,F,degF);
 disp(sprintf('Relative error of diagonal factorization %g',e));
     disp('..................................................');

			                 if mes>0,keyboard,end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X=(W~)^(-1)*V                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

WIS=adj(WI,degWI);
                                   if mes>1,WIS,keyboard,end

[X,degX]=polmul(WIS,degWI,V,degV);
[X,degX]=polmatcl(X,degX);

                             if mes>0
	                        X,degX,keyboard
	                     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Y=(F~)^(-1)*X*F~                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Y=zeros(rank,rank);
degY=0;

for i=1:rank
  for j=i:rank
    [x,degx]=polen(X,degX,i,j);
                                 if mes>1,i,j,x,keyboard,end
    if i==j
      y=x;
      degy=degx;
      [Y,degY]=polseten(Y,degY,y,degy,i,j);
    else
      [fsj,degfsj]=polen(FS,degFS,j,j);
      [x,degx]=polmul(x,degx,fsj,degfsj);
      [fsi,degfsi]=polen(FS,degFS,i,i);
      [y,degy,u,degu]=poldiv(x,degx,fsi,degfsi);
                                    if mes>1,y,u,keyboard,end
      if degu>=0
	error('polsf: not divisible');
      end
      [Y,degY]=polseten(Y,degY,y,degy,i,j);
      [Y,degY]=polseten(Y,degY,adj(y,degy),degy,j,i);
    end
  end
end

                      	                  Ytest=adj(Y,degY);
                             Dif=poladd(Y,degY,-Ytest,degY);
                                e=norm(Dif,inf)/norm(Y,inf);

       disp(sprintf('relative error of symmetry of Y %g',e));

                            if mes>0
		               Y,degY,keyboard
		            end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% factorization of unimodular matrix Y=Z~*J*Z              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  

[J,Z,degZ]=poluf(Y,degY,mes-1);

				if mes>0
				   J,Z,keyboard
                                end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spectral factor of nonsingular part CB=Z*F*W             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[CB,degCB]=polmul(Z,degZ,F,degF);
[CB,degCB]=polmul(CB,degCB,W,degW);
[CB,degCB]=polmatcl(CB,degCB);
  
  disp('..................................................');
         e=polcheck(AB,degAB,adj(CB,degCB),degCB,J,0,CB,degCB);
disp(sprintf('Non-singular spectral factor found with relative error %g',e));

                               if mes >0
                                            CB,degCB,keyboard
                                                   J,keyboard
                               end
  disp('..................................................');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spectral factor  C=(CB (+) I) * UI                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=n-rank;
degC=degCB;

if m==0
  C=CB;
  K=J;
else
  C=coljoin(CB,degCB,zeros(rank,m),0);
  C=rowjoin(C,degC,coljoin(zeros(m,rank),0,eye(m),0),0);
  K=[J zeros(rank,m);zeros(m,n)];
end

[C,degC]=polmul(C,degC,U,degU);
[C,degC]=polmatcl(C,degC);

  
  disp('..................................................');
              e=polcheck(A,degA,adj(C,degC),degC,K,0,C,degC);
disp(sprintf('Spectral factor found with relative error %g',e));
   disp('..................................................');