function [lam,x] = inv_power(A, x0, mu, nsteps,displ)
%
%  call: [lam,x] = inv_power(A, x0, mu, nsteps,displ)
%
%  Inverse power iteration with A, starting with vector x0, with shift mu.
%  Returns approximation to eigenvalue and eigenvector.
%
x = x0;
I = eye(length(x0));
[L,U] = lu(A - mu*I);
for i = 1:nsteps
   y = U\(L\x);         % y = (A - mu*I)\x would cost O(n^3) flops
   x = y/norm(y);
   lam = x'*A*x;
   if displ
      lam
      pause
   end
end
