function [mu,x,M,X]=invpow93(A,x0,Nmax,TOL);
%%function [mu,x,M,X]=invpow93(A,x0,Nmax,TOL);
%% The inverse power method for finding an eigenvalue and corresponding 
%% eigenvector of the nxn matrix A nearest q= x0'Ax0/x0'x0.
%% INPUTS
%%   A - the given matrix
%%   x0 - a starting vector as a column vector
%%   Nmax - the maximum number of iterations
%%   TOL - tolerence in the stopping criterion
%% OUTPUTS
%%   mu - the approximate eigenvalue
%%   x - the approximate coresponding eigenvector
%%   M  - matrix of successive approximations to the eigenvalue
%%   X  - matrix of successive approximations to the eigenvectors as rows

x=x0;q=x0'*A*x0/(x0'*x0);
n=length(x0);
k=1;flag=0;
I=find(abs(x)==norm(x,inf));
x=x/x(I(1));
[L,U,P]=lu_gescp(A-q*eye(n));
while (k<=Nmax & flag==0)
         y=lup_solv(L,U,P,x);
    mu =y(I(1));
    I=find(abs(y)==norm(y,inf));
    err= norm(x-y/y(I(1)),inf);
    x=y/y(I(1));
    if err<TOL
       flag=1;
    end;
    mu=(1/mu)+q;
    k=k+1;M=[M;mu]; X=[X;x'];
end;
