function [evalest evecest iters]=powereig(A,shift,maxiter)

% Shift and Invert Method for estimating the eigenvalue of A closest to 
% the shift, and the corresponding
% eigenvector.  

n=size(A,1);
if n~=size(A,2)
    error('Matrix must be square')
end

if ~exist('shift','var')
    shift=0;
end

if ~exist('maxiter','var')
    maxiter=10*n;
end

iters.x=zeros(n,maxiter);
iters.x(:,1)=randn(n,1);
iters.x(:,1)=randn(n,1)/norm(iters.x(:,1));

iters.theta=zeros(maxiter,1);
iters.theta(1)=iters.x(:,1)'*A*iters.x(:,1);
iters.normx=zeros(maxiter+1,1);
iters.normx(1)=1;

[L U]=lu(A-shift*eye(n));   % Do this once for the linear solves in each loop below, cost of 2/3n^3

for i=1:maxiter
    y=U\(L\iters.x(:,i));   % cost of 2n^2 each iteration
    iters.theta(i)=iters.x(:,i)'*y;   % Rayleigh Quotient (x'(A-shift*I)^{-1}x)/(x'x)

    iters.x(:,i+1)=y/norm(y);
    
    % This gives estimates too, but now you must perform another matrix
    % vector multiply every iter, and these thetas are direct estimates of
    % the eigs of A
    % iters.theta(i+1)=iters.x(:,i+1)'*A*iters.x(:,i+1);   
    iters.normx(i+1)=norm(iters.x(:,i+1)-iters.x(:,i));
end
    

evalest=1/iters.theta(end) + shift;
evecest=iters.x(:,end);


figure(1)
% This plot needs to be changed to just plot theta if you use theta = x'Ax above
plot(abs(1./iters.theta + shift))
title('1/\theta_i+\sigma')

figure(2)
semilogy(abs(diff(1./iters.theta)))
title('Check convergence: |\theta_{i+1}-\theta_i|')

figure(3)
semilogy(iters.normx)
title('Check convergence: ||x_{i+1}-x_i||')

% For comparisons
eigA=eig(A)