function [evector, evalue] = InverseIterationF(A,x,maxit,tol)
% Inverse iteration for the least dominant eigenvalue of A
% A is the n by n input square matrix
% x is the starting eigenvector (n by 1)
% maxit is the max number of iteration
% tol is the error tolerance

%K. Ming Leung, 01/11/03
[n,n] = size(A); rt = rand(1,n);
error = 1; it = 0; evalueOld = eps;
disp('       it            eigenvector          eigenvalue');
[L,U] = LU_factor(A);
while (it < maxit & error > tol)
    y = LU_solve(L,U,x);
    evalue = max(y); evalue = (rt*y)/(rt*x);
    x = y/evalue;
    out = [ it+1 x' evalue]; disp(out);
    error = abs(evalue-evalueOld);
    evalueOld = evalue;
    it = it+1;
end
evector = x;
evalue = 1/evalue;

function [L,U] = LU_factor(A)
% LU factorization of an n by n matrix A
% using Gaussian elimination without pivoting
% A is an n by n matrix
% L is a lower triangylar matrix
% U is an upper triangular matrix
[n,m] = size(A); L = eye(n); U = A;
for j = 1:n
    for i = j+1:n
        L(i,j) = U(i,j)/U(j,j);
        U(i,:) = U(i,:)-L(i,j)*U(j,:);
    end
end

function x = LU_solve(L,U,b)
% Solves the equation LU x = b for x
% using forward and backward substitution
% L is a lower triangular matrix
% U is an upper triangular matrix
% b is a n by 1 vector
[n,m] = size(L); z = zeros(n,1); x = z;
% Solves L z = b using forward substitution
z(1) = b(1);
for i = 2:n
    z(i) = b(i)-L(i,1:i-1)*z(1:i-1);
end
% Solves U x = z using backward substitution
x(n) = z(n)/U(n,n);
for i = n-1:-1:1
    x(i) = (z(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);
end

 

