function [L0,G0,U0] = pc(L0,G0,U0,X,e)
%
%  Perturbation algorithm PC (new version)
%

   [L,M] = size(U0);
   X = sqrt(e)*X;

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

   V1(:,1) = X - V0(:,1);
   for m = 2:M
      V1(:,m) = V1(:,m-1) - V0(:,m);
   end 

   XN = V1(:,M);
   Y(M+1) = norm(XN);

   W1 = zeros(L,M);
   for m = 2:M
      W1(:,m) = W1(:,m-1) + W0(:,m-1);
   end 

   L0 = (1-e)*L0 + YSQ;
   G0 = (1-e)*G0 + Y(M+1)^2/(L-M);
   U0 = U0 + V1 * diag(Z) - W1 * diag(Y(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);
