function [v,d,flag] = eigs(varargin)
%EIGS Find a few eigenvalues and eigenvectors.
% EIGS solves the eigenvalue problem A*v = lambda*v or the generalized
% eigenvalue problem A*v = lambda*B*v. Only a few selected eigenvalues,
% or eigenvalues and eigenvectors, are computed.
%
% [V,D,FLAG] = EIGS(A)
% [V,D,FLAG] = EIGS('Afun',N)
%
% The first input argument is either a square matrix (which can be
% full or sparse, symmetric or nonsymmetric, real or complex), or a
% string containing the name of an M-file which applies a linear
% operator to the columns of a given matrix. In the latter case,
% the second input argument must be N, the order of the problem.
% For example, EIGS('fft',...) is much faster than EIGS(F,...)
% where F is the explicit FFT matrix.
%
% The remaining input arguments are optional and can be given in
% practically any order:
%
% [V,D,FLAG] = EIGS(A,B,K,SIGMA,OPTIONS)
% [V,D,FLAG] = EIGS('Afun',N,B,K,SIGMA,OPTIONS)
%
% where
%
% B A symmetric positive definite matrix the same size as A.
% K An integer, the number of eigenvalues desired.
% SIGMA A scalar shift or a two letter string.
% OPTIONS A structure containing additional parameters.
%
% With one output argument, D is a vector containing K eigenvalues.
% With two output arguments, D is a K-by-K diagonal matrix and V is
% a matrix with K columns so that A*V = V*D or A*V = B*V*D.
% With three output arguments, FLAG indicates whether or not the
% eigenvalues converged to the desired tolerance. FLAG = 0 indicated
% convergence, FLAG = 1 not. FLAG = 2 indicates that EIGS stagnated
% i.e. two consecutive iterates were the same.
%
% If B is not specified, B = eye(size(A)) is used. Note that if B is
% specified, then its Cholesky factorization is computed.
%
% If K is not specified, K = MIN(N,6) eigenvalues are computed.
%
% If SIGMA is not specified, the K-th eigenvalues largest in magnitude
% are computed. If SIGMA is zero, the K-th eigenvalues smallest in
% magnitude are computed. If SIGMA is a real or complex scalar, the
% "shift", the K-th eigenvalues nearest SIGMA are computed using
% shift-invert. Note that this requires computation of the LU
% factorization of A-SIGMA*B. If SIGMA is one of the following strings,
% it specifies the desired eigenvalues.
%
% SIGMA Specified eigenvalues
%
% 'LM' Largest Magnitude (the default)
% 'SM' Smallest Magnitude (same as sigma = 0)
% 'LR' Largest Real part
% 'SR' Smallest Real part
% 'BE' Both Ends. Computes k/2 eigenvalues
% from each end of the spectrum (one more
% from the high end if k is odd.)
%
% The OPTIONS structure specifies certain parameters in the algorithm.
%
% Field name Parameter Default
%
% OPTIONS.tol Convergence tolerance: 1e-10 (symmetric)
% norm(A*V-V*D,1) <= tol * norm(A,1) 1e-6 (nonsymm.)
% OPTIONS.stagtol Stagnation tolerance: quit when 1e-6
% consecutive iterates are the same.
% OPTIONS.p Dimension of the Arnoldi basis. 2*k
% OPTIONS.maxit Maximum number of iterations. 300
% OPTIONS.disp Number of eigenvalues displayed 20
% at each iteration. Set to 0 for
% no intermediate output.
% OPTIONS.issym Positive if Afun is symmetric. 0
% OPTIONS.cheb Positive if A is a string, 0
% sigma is 'LR','SR', or a shift, and
% polynomial acceleration should be
% applied.
% OPTIONS.v0 Starting vector for the Arnoldi rand(n,1)-.5
% factorization
%
% See also EIG, SVDS.
% Richard J. Radke and Dan Sorensen.
% Copyright (c) 1984-98 by The MathWorks, Inc.
% $Revision: 1.19 $ $Date: 1998/08/14 18:50:27 $
global ChebyStruct % Variables used by accpoly.m
% ====== Check input values and set defaults where needed
A = varargin{1};
if isstr(A)
n = varargin{2};
else
[n,n] = size(A);
if any(size(A) ~= n)
error('Matrix A must be a square matrix or a string.')
end
end
if n == 0
v = [];
d = [];
flag = 0;
return
end
% No need to iterate - also elimiates confusion with
% remaining input arguments (they are all n-by-n).
if n == 1
B = 1;
for j = (2+isstr(A)):nargin
if isequal(size(varargin{j}),[n,n]) & ~isstruct(varargin{j})
B = varargin{j};
if ~isreal(B)
error('B must be real and symmetric positive definite.')
end
break
end
end
if isstr(A)
if nargout < 2
v = feval(A,1) / B;
else
v = 1;
d = feval(A,1) / B;
flag = 0;
end
return
else
if nargout < 2
v = A / B;
else
v = 1;
d = A / B;
flag = 0;
end
return
end
end
B = [];
b2 = 1;
k = [];
sigma = [];
options = [];
for j = (2+isstr(A)):nargin
if isequal(size(varargin{j}),[n,n])
B = varargin{j};
if ~isreal(B)
error('B must be real and symmetric positive definite.')
end
elseif isstruct(varargin{j})
options = varargin{j};
elseif isstr(varargin{j}) & length(varargin{j}) == 2
sigma = varargin{j};
elseif max(size(varargin{j})) == 1
s = varargin{j};
if isempty(k) & isreal(s) & (s == fix(s)) & (s > 0)
k = min(n,s);
else
sigma = s;
end
else
error(['Input argument number ' int2str(j) ' not recognized.'])
end
end
% Set defaults
if isempty(k)
k = min(n,6);
end
if isempty(sigma)
sigma = 'LM';
end
if isempty(options)
fopts = [];
else
fopts = fieldnames(options);
end
% A is the matrix of all zeros (not detectable if A is a string)
if ~isstr(A)
if nnz(A) == 0
if nargout < 2
v = zeros(k,1);
else
v = eye(n,k);
d = zeros(k,k);
flag = 0;
end
return
end
end
% Trick to get faster convergence. The tail end (closest to cut off
% of the sort) will not converge as fast as the leading end of the
% "wanted" Ritz values.
ksave = k;
k = min(n-1,k+3);
% Set issym, specifying a symmetric matrix or operator.
if ~isstr(A)
issym = isequal(A,A');
elseif strmatch('issym',fopts)
issym = options.issym;
else
issym = 0;
end
% Set tol, the convergence tolerance.
if strmatch('tol',fopts);
tol = options.tol;
else
if issym
tol = 1e-10; % Default tol for symmetric problems.
else
tol = 1e-6; % Default tol for nonsymmetric problems.
end
end
% Set stagtol, the stagnation tolerance.
if strmatch('stagtol',fopts);
stagtol = options.stagtol;
else
stagtol = 1e-6;
end
% Set p, the dimension of the Arnoldi basis.
if strmatch('p',fopts);
p = options.p;
else
p = min(max(2*k,20),n);
end
% Set maxit, the maximum number of iterations.
if strmatch('maxit',fopts);
maxit = options.maxit;
else
maxit = max(300,2*n/p);
end
% Set display option.
if strmatch('disp',fopts);
dispn = options.disp;
else
dispn = 20;
end
if strmatch('v0',fopts);
v0 = options.v0;
else
v0 = rand(n,1)-.5;
end
% Check if Chebyshev iteration is requested.
if strmatch('cheb',fopts);
cheb = options.cheb;
else
cheb = 0;
end
if cheb & issym
if isstr(A) & strcmp(sigma,'LR')
sigma = 'LO';
ChebyStruct.sig = 'LR';
ChebyStruct.filtpoly = [];
elseif isstr(A) & strcmp(sigma,'SR')
sigma = 'SO';
ChebyStruct.sig = 'SR';
ChebyStruct.filtpoly = [];
elseif isstr(A) & ~isstr(sigma)
ChebyStruct.sig = sigma;
sigma = 'LO';
ChebyStruct.filtpoly = [];
elseif isstr(A)
warning('Cheb option only available for sym
评论0