function [L0,U0] = pa(L0,U0,X,e)
%
%  Perturbation algorithm PA
%
   tol = .005;
   L = length(X);
   X = sqrt(e)*X;

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

   for I = 1:L
      B(I,I) = 1;
      for J = I+1:L
%        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 + YSQ;
   U0 = U0 * B; 
   for I = 1:L
      U0(:,I) = gorth(U0(:,I));
   end

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