aeigs60.m

特征结构配置方法（个人）

function  varargout = eigs(varargin)
%EIGS  Find a few eigenvalues and eigenvectors of a matrix using ARPACK.
%   D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues.
%   A must be square and should be large and sparse.
%
%   [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude
%   eigenvalues and a matrix V whose columns are the corresponding eigenvectors.
%
%   [V,D,FLAG] = EIGS(A) also returns a convergence flag.  If FLAG is 0
%   then all the eigenvalues converged; otherwise not all converged.
%
%   EIGS(AFUN,N) accepts the function AFUN instead of the matrix A.
%   Y = AFUN(X) should return Y = A*X.  N is the size of A.  The matrix A
%   represented by AFUN is assumed to be real and nonsymmetric.  In all these
%   calling sequences, EIGS(A,...) may be replaced by EIGS(AFUN,N,...).
%
%   EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D.  B must
%   be symmetric (or Hermitian) positive definite and the same size as A.
%   EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D.
%
%   EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues.
%
%   EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues based on SIGMA:
%      'LM' or 'SM' - Largest or Smallest Magnitude
%   For real symmetric problems, SIGMA may also be:
%      'LA' or 'SA' - Largest or Smallest Algebraic
%      'BE' - Both Ends, one more from high end if K is odd
%   For nonsymmetric and complex problems, SIGMA may also be:
%      'LR' or 'SR' - Largest or Smallest Real part
%      'LI' or 'SI' - Largest or Smallest Imaginary part
%   If SIGMA is a real or complex scalar including 0, EIGS finds the eigenvalues
%   closest to SIGMA.  For scalar SIGMA, and also when SIGMA = 'SM' which uses
%   the same algorithm as SIGMA = 0, B need only be symmetric (or Hermitian)
%   positive semi-definite and the function Y = AFUN(X) must return
%   Y = (A-SIGMA*B)\X, which is Y = A\X when SIGMA = 0 or SIGMA = 'SM'.
%
%   EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options:
%   OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{0} | 1]
%   OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [0 | {1}]
%   OPTS.tol: convergence: abs(d_comp-d_true) < tol*abs(d_comp) [scalar | {eps}]
%   OPTS.maxit: maximum number of iterations [integer | {300}]
%   OPTS.p: number of Lanczos vectors: K+1<p<=N [integer | {2K}]
%   OPTS.v0: starting vector [N-by-1 vector | {randomly generated by ARPACK}]
%   OPTS.disp: diagnostic information display level [0 | {1} | 2]
%   OPTS.cholB: B is actually its Cholesky factor CHOL(B) [{0} | 1]
%   OPTS.permB: sparse B is actually CHOL(B(permB,permB)) [permB | {1:N}]
%
%   EIGS(AFUN,N,K,SIGMA,OPTS,P1,P2,...) and EIGS(AFUN,N,B,K,SIGMA,OPTS,P1,P2,...)
%   provide for additional arguments which are passed to AFUN(X,P1,P2,...).
%
%   Examples:
%      A = delsq(numgrid('C',15));  d1 = eigs(A,5,'SM');
%   Equivalently, if dnRk is the following one-line function:
%      function y = dnRk(x,R,k)
%      y = (delsq(numgrid(R,k))) * x;
%   then pass dnRk's additional arguments, 'C' and 15, to eigs:
%      n = size(A,1);  opts.issym = 1;  d2 = eigs(@dnRk,n,5,'SM',opts,'C',15);
%
%   See also EIG, SVDS, ARPACKC.



cputms = zeros(5,1);
t0 = cputime; % start timing pre-processing

if (nargout > 3)
    error('Too many output arguments')
end

% Process inputs and do error-checking
if isa(varargin{1},'double')
    A = varargin{1};
    Amatrix = 1;
else
	A = fcnchk(varargin{1});
    Amatrix = 0;
end

isrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma)
if Amatrix
    isrealprob = isreal(A);
end

issymA = 0;
if Amatrix
    issymA = isequal(A,A');
end

if Amatrix
    [m,n] = size(A);
    if (m ~= n)
        error('A must be a square matrix or a function which computes A*x')
    end
else
    n = varargin{2};
    nstr = sprintf('Size of problem n must be a positive integer');
    if ~isequal(size(n),[1,1]) | ~isreal(n)
        error(nstr)
    end
    if (round(n) ~= n)
        warning(nstr)
        n = round(n);
    end
    if issparse(n)
        n = full(n);
    end
end

Bnotthere = 0;
Bstr = sprintf(['Generalized matrix B must be the same size as A and' ...
        ' either a symmetric positive (semi-)definite matrix or' ...
        ' its Cholesky factor']);
if (nargin < (3-Amatrix-Bnotthere))
    B = [];
    Bnotthere = 1;
else
    Bk = varargin{3-Amatrix-Bnotthere};
    if isempty(Bk) % allow eigs(A,[],k,sigma,opts);
        B = Bk;
    else
        if isequal(size(Bk),[1,1]) & (n ~= 1)
            B = [];
            k = Bk;
            Bnotthere = 1;
        else % eigs(9,8,...) assumes A=9, B=8, ... NOT A=9, k=8, ...
            B = Bk;
            if ~isa(B,'double') | ~isequal(size(B),[n,n])
                error(Bstr)
            end
            isrealprob = isrealprob & isreal(B);
        end
    end
end

if Amatrix & ((nargin - ~Bnotthere)>4)
    error('Too many inputs')
end

if (nargin < (4-Amatrix-Bnotthere))
    k = min(n,6);
else
    k = varargin{4-Amatrix-Bnotthere};
end

kstr = sprintf(['Number of eigenvalues requested k must be a' ...
        ' positive integer <= n']);
if ~isa(k,'double') | ~isequal(size(k),[1,1]) | ~isreal(k) | (k>n)
    error(kstr)
end
if issparse(k)
    k = full(k);
end
if (round(k) ~= k)
    warning(kstr)
    k = round(k);
end

whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string']);
if (nargin < (5-Amatrix-Bnotthere))
    % default: eigs('LM') => ARPACK which='LM', sigma=0
    eigs_sigma = 'LM';
    whch = 'LM';
    sigma = 0;
else
    eigs_sigma = varargin{5-Amatrix-Bnotthere};
    if isstr(eigs_sigma)
        % eigs(string) => ARPACK which=string, sigma=0
        if ~isequal(size(eigs_sigma),[1,2])
            error(whchstr)
        end
        eigs_sigma = upper(eigs_sigma);
        if isequal(eigs_sigma,'SM')
            % eigs('SM') => ARPACK which='LM', sigma=0
            whch = 'LM';
        else
            % eigs(string), where string~='SM' => ARPACK which=string, sigma=0
            whch = eigs_sigma;
        end
        sigma = 0;
    else
        % eigs(scalar) => ARPACK which='LM', sigma=scalar
        if ~isa(eigs_sigma,'double') | ~isequal(size(eigs_sigma),[1,1])
            error('Eigenvalue shift sigma must be a scalar')
        end
        sigma = eigs_sigma;
        if issparse(sigma)
            sigma = full(sigma);
        end
        isrealprob = isrealprob & isreal(sigma);
        whch = 'LM';
    end
end

tol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS)
maxit = [];
p = [];
info = int32(0); % use a random starting vector
display = 1;
cholB = 0;

if (nargin >= (6-Amatrix-Bnotthere))
    opts = varargin{6-Amatrix-Bnotthere};
    if ~isa(opts,'struct')
        error('Options argument must be a structure')
    end
    
    if isfield(opts,'issym') & ~Amatrix
        issymA = opts.issym;
        if (issymA ~= 0) & (issymA ~= 1)
            error('opts.issym must be 0 or 1')
        end
    end
    
    if isfield(opts,'isreal') & ~Amatrix
        if (opts.isreal ~= 0) & (opts.isreal ~= 1)
            error('opts.isreal must be 0 or 1')
        end
        isrealprob = isrealprob & opts.isreal;
    end
    
    if ~isempty(B) & (isfield(opts,'cholB') | isfield(opts,'permB'))
        if isfield(opts,'cholB')
            cholB = opts.cholB;
            if (cholB ~= 0) & (cholB ~= 1)
                error('opts.cholB must be 0 or 1')
            end
            if isfield(opts,'permB')
                if issparse(B) & cholB
                    permB = opts.permB;
                    if ~isequal(sort(permB),(1:n)) & ...
                            ~isequal(sort(permB),(1:n)')
                        error('opts.permB must be a permutation of 1:n')
                    end
                else
                    warning(['Ignoring opts.permB since B is not its sparse' ...
                            ' Cholesky factor'])
                end
            else
                permB = 1:n;
            end
        end
    end
    
    if isfield(opts,'tol')
        if ~isequal(size(opts.tol),[1,1]) | ~isreal(opts.tol) | (opts.tol<=0)
            error(['Convergence tolerance opts.tol must be a strictly' ...
                    ' positive real scalar'])
        else
            tol = full(opts.tol);
        end
    end
    
    if isfield(opts,'p')
        p = opts.p;
        pstr = sprintf(['Number of basis vectors opts.p must be a positive' ...
                ' integer <= n']);
        if ~isequal(size(p),[1,1]) | ~isreal(p) | (p<=0) | (p>n)
            error(pstr)
        end
        if issparse(p)
            p = full(p);
        end
        if (round(p) ~= p)
            warning(pstr)
            p = round(p);
        end
    end
    
    if isfield(opts,'maxit')
        maxit = opts.maxit;
        str = sprintf(['Maximum number of iterations opts.maxit must be' ...
                ' a positive integer']);
        if ~isequal(size(maxit),[1,1]) | ~isreal(maxit) | (maxit<=0)
            error(str)
        end
        if issparse(maxit)
            maxit = full(maxit);
        end
        if (round(maxit) ~= maxit)
            warning(str)
            maxit = round(maxit);
        end
    end
    
    if isfield(opts,'v0')
        if ~isequal(size(opts.v0),[n,1])
            error('Start vector opts.v0 must be n-by-1')
        end
        if isrealprob
            if ~isreal(opts.v0)
                error(['Start vector opts.v0 must be real for real problems'])
            end
            resid = full(opts.v0);
        else
            resid(1:2:(2*n-1),1) = full(real(opts.v0));
            resid(2:2:2*n,1) = full(imag(opts.v0));
        end
        info = int32(1); % use resid as starting vector
    end
    
    if isfield(opts,'disp')
        display = opts.disp;
        dispsyt = sprintf('Diagnostic level opts.disp must be an integer');
        if (~isequal(size(display),[1,1])) | (~isreal(display)) | (display<0)
            error(dispstr)
        end
        if (round(display) ~= display)
            warning(dispstr)
            display = round(display);
        end
    end
    
    if isfield(opts,'cheb')
        warning(['Ignoring polynomial acceleration opts.cheb' ...
                ' (no longer an option)']);
    end
    if isfield(opts,'stagtol')
        warning(['Ignoring stagnation tolerance opts.stagtol' ...
                ' (no longer an option)']);
    end
    
end

% Now we know issymA, isrealprob, cholB, and permB

if isempty(p)
    if isrealprob & ~issymA
        p = min(max(2*k+1,20),n);
    else
        p = min(max(2*k,20),n);
    end
end
if isempty(maxit)
    maxit = max(300,ceil(2*n/max(p,1)));
end
if (info == int32(0))
    if isrealprob
        resid = zeros(n,1);
    else
        resid = zeros(2*n,1);