% Inverse iteration for the eigenvalue of A that is closest to sigma
% A is the n by n input square matrix
% x is the starting eigenvector (n by 1)
% sigma is the shift
% maxit is the max number of iteration
% tol is the error tolerance

%K. Ming Leung, 01/14/03


clear all; format compact; maxit =500; tol=1e-4;
%A = [ 6 2 1; 2 3 1; 1 1 1]; [n,n]=size(A); x=ones(n,1); sigma = 5.1;
%A = [ 1 4; 1 1]; [n,n]=size(A); x=ones(n,1); sigma = 2;
A = [ 3 1; 1 3]; [n,n]=size(A); x=[0; 1]; sigma = 5.5;

[n,n] = size(A); rt = rand(1,n);
AS = A-sigma*eye(n);
error = 1; it = 0; 

disp('       it       eigenvector           eigenvalue ');
while (it < maxit & error > tol)
    y=AS\x; % Should really prefactorize outside loop
    ev = (rt*y)/(rt*x);
    error = norm(y-ev*x);
    x = y/ev;
    out = [ it+1 x' ev]; disp(out);
    it = it+1;
end
evector = x/norm(x),
evalue = sigma+1/ev,
[Mevec,Meval]=eig(A),


