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-processingif (nargout > 3)    error('Too many output arguments')end% Process inputs and do error-checkingif isa(varargin{1},'double')    A = varargin{1};    Amatrix = 1;else	A = fcnchk(varargin{1});    Amatrix = 0;endisrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma)if Amatrix    isrealprob = isreal(A);endissymA = 0;if Amatrix    issymA = isequal(A,A');endif Amatrix    [m,n] = size(A);    if (m ~= n)        error('A must be a square matrix or a function which computes A*x')    endelse    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);    endendBnotthere = 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    endendif Amatrix & ((nargin - ~Bnotthere)>4)    error('Too many inputs')endif (nargin < (4-Amatrix-Bnotthere))    k = min(n,6);else    k = varargin{4-Amatrix-Bnotthere};endkstr = 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)endif issparse(k)    k = full(k);endif (round(k) ~= k)    warning(kstr)    k = round(k);endwhchstr = 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';    endendtol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS)maxit = [];p = [];info = int32(0); % use a random starting vectordisplay = 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 permBif isempty(p)    if isrealprob & ~issymA        p = min(max(2*k+1,20),n);    else        p = min(max(2*k,20),n);    endendif isempty(maxit)    maxit = max(300,ceil(2*n/max(p,1)));endif (info == int32(0))    if isrealprob        resid = zeros(n,1);    else        resid = zeros(2*n,1);