function [L0,G0,U0] = pb(L0,G0,U0,X,e)
%
%  Perturbation algorithm PB (new version)
%
   tol = .005;
   [L,M] = size(U0);
   X = sqrt(e)*X;

   ETA = U0' * X;
   YSQ = real(ETA).^2 + imag(ETA).^2;
   Y = sqrt(YSQ);
   U0 = U0 * diag(ETA./Y);

   XS = U0 * Y;   
   XN = X - XS;
   Y(M+1) = norm(XN);
   U0(:,M+1) = XN / Y(M+1); 

   L0(M+1) = G0;
   for I = 1:M+1
      B(I,I) = 1;
      for J = I+1:M+1
%        B(J,I) = Y(I)*Y(J) / (L0(I)-L0(J));
         B(J,I) = Y(I)*Y(J) / max(tol*L0(1),L0(I)-L0(J));
         B(I,J) = -B(J,I);
      end
   end

   L0 = (1-e)*L0(1:M) + YSQ;
   G0 = (1-e)*G0 + Y(M+1)^2 / (L-M);
   U0 = U0 * B(:,1:M); 
   for I = 1:M
      U0(:,I) = gorth(U0(:,I));
   end

%   ordering:
%   [L0,I] = sort(L0);
%   L0 = flipud(L0);
%   U0 = U0(:,I);
%   U0 = fliplr(U0);
