function [eig1,x,k] = my_power_shift(A,tau,tol, x0)

% Input:  A : an n by n matrix;
%        tau: Shift number
%        tol: tolerance 
%        x0 : non-zero initial guess
%
% Output: eig1: the eigenvalue of A that is closest to tau
%          x : the approximation of a unit eigenvector
%          k: the number of iterations

x = x0/norm(x0);
eig0 = 0; eig1 =1; k=0;
[m,n] = size(A); B = A - tau*eye(m); 
[l u p ] = lu(B);
while ((abs(eig0-eig1)) >= tol)
   k = k +1;
   eig0 = eig1; 
   b = p*x; y1 = l\b;  y=u\y1;
   x = y/norm(y);
   eig1 = x'*A*x;
end