function [U,J,UI] = lagr(B)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  decomposition of B via Lagrange                 % 
%                B=U'*J*U                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M. Sebek, June 1990                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



[rB,cB]=size(B);
J=B; U = eye(rB); UI=U;
EPS=norm(B,inf)*1e-8;



for i=1:rB
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % reduction of row and column i                  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  V=eye(U);VI=V;

  if abs(J(i,i))>EPS
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % diagonal entry nonzero                       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    V(i,i)=sqrt(abs(J(i,i)));
    VI(i,i)=1/V(i,i);
    for j=i+1:rB
      VI(i,j)=-J(i,j)/J(i,i);
      V(i,j)=J(i,j)*VI(i,i)*sign(J(i,i));
    end
  end

  if abs(J(i,i))<=EPS 

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % diagonal entry zero                          %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    J(i,i)=0;
    if i<rB
      j=i+1;
      while abs(J(i,j))<EPS & j<rB
       j=j+1;
      end
      if abs(J(i,j))>EPS
        c=J(i,j)*2;
        VI(i,i)=(1-J(j,j))/c;
        VI(i,j)=(-1-J(j,j))/c;
        VI(j,i)=1;

        V(i,i)=J(i,j);
        V(i,j)=(1+J(j,j))/2;
        V(j,i)=-J(i,j);
        V(j,j)=(1-J(j,j))/2;
      end
    end
  end

% J=VI'*J*VI;
% U=V*U;
% UI=UI*VI;
  
  J=polmul(VI',0,J,0);
  J=polmul(J,0,VI,0);
  U=polmul(V,0,U,0);
  UI=polmul(UI,0,VI,0);
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % end of reduction of i                          %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% permutation of diagonal J                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

d=diag(J);
K=find(d==0);
z=size(K);
d(K)=ones(z,1)*(-2);
J=diag(d);

[L,K]=sort(-diag(J));
V=zeros(rB,rB);VI=V;
for i=1:rB
  VI(K(i),i)=1;
  V(i,K(i))=1;
end

J=VI'*J*VI;
U=V*U;
UI=UI*VI;

d=diag(J);
d(rB-z+1:rB)=zeros(z,1);
J=diag(round(d));