流形学习_OGLdpf_流形_流形学习_original5uu_目标分类_

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