my_spec_cut.m

利用加速随机游走计算初始点确定的分类算法

    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('MATLAB:eigs:IgnoredOptionPermB', ...
                            ['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 = ['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('MATLAB:eigs:NonIntegerVecQty',['%s\n         ' ...
                    'Rounding number of basis vectors.'],pstr)
            p = round(p);
        end
    end
    
    if isfield(opts,'maxit')
        maxit = opts.maxit;
        str = ['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('MATLAB:eigs:NonIntegerIterationQty',['%s\n         ' ...
                    'Rounding number of iterations.'],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;
        dispstr = '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('MATLAB:eigs:NonIntegerDiagnosticLevel', ...
                    '%s\n         Rounding diagnostic level.',dispstr)
            display = round(display);
        end
    end
    
    if isfield(opts,'cheb')
        warning('MATLAB:eigs:ObsoleteOptionCheb', ...
                ['Ignoring polynomial acceleration opts.cheb' ...
                ' (no longer an option).']);
    end
    if isfield(opts,'stagtol')
        warning('MATLAB:eigs:ObsoleteOptionStagtol', ...
                ['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);
    end
end

if ~isempty(B) % B must be symmetric (Hermitian) positive (semi-)definite
    if cholB
        if ~isequal(triu(B),B)
            error(Bstr)
        end
    else
        if ~isequal(B,B')
            error(Bstr)
        end
    end
end

useeig = 0;
if isrealprob & issymA
    knstr = sprintf(['For real symmetric problems, must have number' ...
            ' of eigenvalues k < n.\n']);
else
    knstr = sprintf(['For nonsymmetric and complex problems, must have' ...
            ' number of eigenvalues k < n-1.\n']);
end
if isempty(B)
    knstr = [knstr sprintf(['Using eig(full(A)) instead.'])];
else
    knstr = [knstr sprintf(['Using eig(full(A),full(B)) instead.'])];
end
if (k == 0)
    useeig = 1;
end
if isrealprob & issymA
    if (k > n-1)
        if (n >= 6)
            warning('MATLAB:eigs:TooManyRequestedEigsForRealSym', ...
                '%s',knstr)
        end
        useeig = 1;
    end
else
    if (k > n-2)
        if (n >= 7)
            warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ...
                '%s',knstr)
        end
        useeig = 1;
    end
end

if isrealprob & issymA
    if ~isreal(sigma)
        error(['For real symmetric problems, eigenvalue shift sigma must' ...
                ' be real.'])
    end
else
    if ~isrealprob & issymA & ~isreal(sigma)
        warning('MATLAB:eigs:ComplexShiftForRealProblem', ...
            ['Complex eigenvalue shift sigma on a Hermitian problem' ...
            ' (all real eigenvalues).'])
    end
end

if isrealprob & issymA
    if strcmp(whch,'LR')
        whch = 'LA';
        warning('MATLAB:eigs:SigmaChangedToLA', ... 
            ['For real symmetric problems, sigma value ''LR''' ...
            ' (Largest Real) is now ''LA'' (Largest Algebraic).'])
    end
    if strcmp(whch,'SR')
        whch = 'SA';
        warning('MATLAB:eigs:SigmaChangedToSA', ...
            ['For real symmetric problems, sigma value ''SR''' ...
            ' (Smallest Real) is now ''SA'' (Smallest Algebraic).'])
    end
    %if (~strcmp(whch,'LM')) && (~strcmp(whch,'SM')) && (~strcmp(whch,'LA')) && (~strcmp(whch,'SA')) && (~strcmp(whch,'BE'))
    if (~strcmp(whch,'LM')) & (~strcmp(whch,'SM')) & (~strcmp(whch,'LA')) & (~strcmp(whch,'SA')) & (~strcmp(whch,'BE'))
    %if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'})
        whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
                'LM','SM','LA','SA','BE')];
        error(whchstr)
    end
else
    if strcmp(whch,'BE')
        warning('MATLAB:eigs:SigmaChangedToLM', ...
            ['Sigma value ''BE'' is now only available for real' ...
            ' symmetric problems.  Computing ''LM'' eigenvalues instead.'])
        whch = 'LM';
    end
    if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'})
         whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ...
                 'LM','SM','LR','SR','LI','SI')];
        error(whchstr)
    end
end

% Now have enough information to do early return on cases eigs does not handle
if useeig
    if (nargout <= 1)
        varargout{1} = eigs2(A,n,B,k,whch,sigma,cholB, ...
            varargin{7-Amatrix-Bnotthere:end});
    else
        [varargout{1},varargout{2}] = eigs2(A,n,B,k,whch,sigma,cholB, ...
            varargin{7-Amatrix-Bnotthere:end});
    end
    if (nargout == 3)
        varargout{3} = 0; % flag indicates "convergence"
    end
    return
end

if isrealprob & ~issymA
    sigmar = real(sigma);
    sigmai = imag(sigma);
end

if isrealprob & issymA
    if (p <= k)
        error(['For real symmetric problems, must have number of' ...
                ' basis vectors opts.p > k.'])
    end
else
    if (p <= k+1)
        error(['For nonsymmetric and complex problems, must have number of' ...
                ' basis vectors opts.p > k+1.'])
    end
end

if isequal(whch,'LM') & ~isequal(eigs_sigma,'LM')
    % A*x = lambda*M*x, M symmetric (positive) semi-definite
    % => OP = inv(A - sigma*M)*M and B = M
    % => shift-and-invert mode
    mode = 3;
elseif isempty(B)
    % A*x = lambda*x
    % => OP = A and B = I
    mode = 1;
else % B is not empty
    % Do not use mode=2.
    % Use mode = 1 with OP = R'\(A*(R\x)) and B = I
    % where R is B's upper triangular Cholesky factor: B = R'*R.
    % Finally, V = R\V returns the actual generalized eigenvectors of A and B.
    mode = 1;
end

if cholB
    pB = 0;
    RB = B;
    RBT = B';
end

if (mode == 3) & Amatrix % need lu(A-sigma*B)
    if issparse(A) & (isempty(B) | issparse(B))
        if isempty(B)
            AsB = A - sigma * speye(n);
        else
            if cholB
                AsB = A - sigma * RBT * RB;
            else
                AsB = A - sigma * B;
            end
        end
        [L,U,P,Q] = lu(AsB);
        [perm,dummy] = find(Q);
    else
        if isempty(B)
            AsB = A - sigma * eye(n);
        else
            if cholB
                AsB = A - sigma * RBT * RB;
            else
                AsB = A - sigma * B;
            end
        end
        [L,U,P] = lu(AsB);
    end
    dU = diag(U);
    rcondestU = full(min(abs(dU)) / max(abs(dU)));
    if (rcondestU < eps)
        if isempty(B)
            ds = sprintf(['(A-sigma*I) has small reciprocal condition' ...
                    ' estimate: %f\n'],rcondestU);
        else
            ds = sprintf(['(A-sigma*B) has small reciprocal condition' ...
                    ' estimate: %f\n'],rcondestU);
        end
        ds = [ds sprintf(['indicating that sigma is near an exact' ...
                    ' eigenvalue. The\nalgorithm may not converge unless' ...
                    ' you try a new value for sigma.\n'])];
        warning('MATLAB:eigs:SigmaNearExactEig',ds)
    end
end

if (mode == 1) & ~isempty(B) & ~cholB % need chol(B)
    if issparse(B)
        permB = symamd(B);
        [RB,pB] = chol(B(permB,permB));
    else
        [RB,pB] = chol(B);
    end
    if (pB == 0)
        RBT = RB';
    else
        error(Bstr)
    end
end

% Allocate outputs and ARPACK work variables
if isrealprob
    if issymA % real and symmetric
        prefix = 'ds';
        v = zeros(n,p);
        ldv = int32(size(v,1));
        ipntr = int32(zeros(15,1));
        workd = zeros(n,3);
        lworkl = p*(p+8);
        workl = zeros(lworkl,1);
        lworkl = int32(lworkl);
        d = zeros(k,1);
    else % real but not symmetric
        prefix = 'dn';
        v = zeros(n,p);
        ldv = int32(size(v,1));
        ipntr = int32(zeros(15,1));
        workd = zeros(n,3);
        lworkl = 3*p*(p+2);
        workl = zeros(lworkl,1);
        lworkl = int32(lworkl);
        workev = zeros(3*p,1);
        d = zeros(k+1,1);
        di = zeros(k+1,1);
    end
else % complex
    prefix = 'zn';
    zv = zeros(2*n*p,1);
    ldv = int32(n);
    ipntr = int32(zeros(15,1));
    workd = complex(zeros(n,3));
    zworkd = zeros(2*prod(size(workd)),1);
    lworkl = 3*p^2+5*p;
    workl = zeros(2*lworkl,1);
    lworkl = int32(lworkl);
    workev = zeros(2*2*p,1);
    zd = zeros(2*(k+1),1);
    rwork = zeros(p,1);
end

ido = int32(0); % reverse communication parameter
if isempty(B) | (~isempty(B) & (mode == 1))
    bmat = 'I'; % standard eigenvalue problem
else
    bmat = 'G'; % generalized eigenvalue problem
end
nev = int32(k); % number of eigenvalues requested
ncv = int32(p); % number of Lanczos vectors
iparam = int32(zeros(11,1));
iparam([1 3 7]) = int32([1 maxit mode]);
select = int32(zeros(p,1));

cputms(1) = cputime - t0; % end timing pre-processing

iter = 0;
ariter = 0;


%tim


indexArgumentsAfun = 7-Amatrix-Bnotthere:length(varargin);
nbArgumentsAfun = length(indexArgumentsAfun);
if nbArgumentsAfun >=1
    arguments_Afun = varargin{7-Amatrix-Bnotthere};
end
if nbArgumentsAfun >=2
    arguments_Afun2 = varargin{7-Amatrix-Bnotthere+1};
end
if nbArgumentsAfun >=3
    arguments_Afun3 = varargin{7-Amatrix-Bnotthere+2};
end
%fin tim



while (ido ~= 99)
        
    t0 = cputime; % start timing ARPACK calls **aupd
        
    if isrealprob
        arpackc( [prefix 'aupd'], ido, ...
            bmat, int32(n), whch, nev, tol, resid, ncv, ...
            v, ldv, iparam, ipntr, workd, workl, lworkl, info);
    else