invshift.m

sareli<matlab科学计算>配套程序

function [lambda,x,iter]=invshift(A,mu,tol,nmax,x0)
%INVSHIFT     Numerically evaluate one eigenvalue of a matrix.
%   LAMBDA = INVSHIFT(A) compute with the inverse power method the 
%   eigenvalue of A of minimum modulus.
%   LAMBDA = INVSHIFT(A,MU) compute the eigenvalue
%   of A closest to the given number (real or complex) MU.
%   LAMBDA = INVSHIFT(A,MU,TOL,NMAX,X0) uses an absolute error tolerance of TOL 
%   instead of the default, which is 1.e-6 and a maximum number of iterations
%   of NMAX (the default is 100), starting from the initial vector X0.
%   [LAMBDA,V,ITER] = INVSHIFT(A,MU,TOL,NMAX,X0) also returns the eigenvector
%   V such that A*V=LAMBDA*V and the iteration number at which V was computed.
[n,m]=size(A);
if n ~= m, error('Only for square matrices'); end
if nargin == 1
  x0 = rand(n,1); nmax = 100; tol = 1.e-06; mu = 0;
elseif nargin == 2
  x0 = rand(n,1); nmax = 100; tol = 1.e-06;
end
[L,U]=lu(A-mu*eye(n)); 
if norm(x0) == 0 
   x0 = rand(n,1); 
end 
x0=x0/norm(x0); 
z0=L\x0; 
pro=U\z0; 
lambda=x0'*pro; 
err=tol+1; 
iter=0;
while err > tol*abs(lambda) & abs(lambda) ~= 0 & iter <= nmax
   x = pro; x = x/norm(x); 
   z=L\x; 
   pro=U\z; 
   lambdanew = x'*pro;   
   err = abs(lambdanew - lambda); 
   lambda = lambdanew;   
   iter = iter + 1;
end
lambda = 1/lambda + mu;
return