my_spec_cut.m

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

        zworkd(1:2:end-1) = real(workd);
        zworkd(2:2:end) = imag(workd);
        arpackc( 'znaupd', ido, ...
            bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
            ldv, iparam, ipntr, zworkd, workl, lworkl, ...
            rwork, info );
        workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]);
    end

    if (info < 0)
        es = sprintf('Error with ARPACK routine %saupd: info = %d',...
           prefix,full(info));
        error(es)
    end
     
    cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd
    t0 = cputime; % start timing MATLAB OP(X)
    
    % Compute which columns of workd ipntr references
    
    
    
    
    
    %[row,col1] = ind2sub([n,3],double(ipntr(1)));
    %tim
    row = mod(double(ipntr(1))-1,n) + 1 ;
    col1 = floor((double(ipntr(1))-1)/n) + 1;
    
    
    if (row ~= 1)
        str = sprintf(['ipntr(1)=%d does not refer to the start of a' ...
                ' column of the %d-by-3 array workd.'],full(ipntr(1)),n);
        error(str)
    end
    
    
    
    %[row,col2] = ind2sub([n,3],double(ipntr(2)));
    %tim
    row = mod(double(ipntr(2))-1,n) + 1 ;
    col2 = floor((double(ipntr(2))-1)/n) + 1;

    
    
    if (row ~= 1)
        str = sprintf(['ipntr(2)=%d does not refer to the start of a' ...
                ' column of the %d-by-3 array workd.'],full(ipntr(2)),n);
        error(str)
    end
    if ~isempty(B) & (mode == 3) & (ido == 1)
        [row,col3] = ind2sub([n,3],double(ipntr(3)));
        if (row ~= 1)
            str = sprintf(['ipntr(3)=%d does not refer to the start of a' ...
                    ' column of the %d-by-3 array workd.'],full(ipntr(3)),n);
            error(str)
        end
    end

    switch (ido)
        case {-1,1}
        if Amatrix
            if (mode == 1)
                if isempty(B)
                    % OP = A*x
                    workd(:,col2) = A * workd(:,col1);
                else
                    % OP = R'\(A*(R\x))
                    if issparse(B)
                        workd(permB,col2) = RB \ workd(:,col1);
                        workd(:,col2) = A * workd(:,col2);
                        workd(:,col2) = RBT \ workd(permB,col2);
                    else
                        workd(:,col2) = RBT \ (A * (RB \ workd(:,col1)));
                    end
                end
            elseif (mode == 3)
                if isempty(B)
                    if issparse(A)
                        workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
                    else
                        workd(:,col2) = U \ (L \ (P * workd(:,col1)));
                    end
                else % B is not empty
                    if (ido == -1)
                        if cholB
                            workd(:,col2) = RBT * (RB * workd(:,col1));
                        else
                            workd(:,col2) = B * workd(:,col1);
                        end
                        if issparse(A) & issparse(B)
                            workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
                        else
                            workd(:,col2) = U \ (L \ (P * workd(:,col1)));
                        end
                    elseif (ido == 1)
                        if issparse(A) & issparse(B)
                            workd(perm,col2) = U \ (L \ (P * workd(:,col3)));
                        else
                            workd(:,col2) = U \ (L \ (P * workd(:,col3)));
                        end
                    end
                end
            else % mode is not 1 or 3
                error(['Unknown mode returned from ' prefix 'aupd.'])
            end
        else % A is not a matrix
            if (mode == 1)
                if isempty(B)
                    % OP = A*x
                    %workd(:,col2) = feval(A,workd(:,col1),varargin{7-Amatrix-Bnotthere:end});                    
                    
                    
                    
                    
                    
                    nombre_A_times_X = nombre_A_times_X + 1;
                    
                    
                    
                    %pause(0); %voir
                    
                    if nbArgumentsAfun == 1
                        workd(:,col2) = feval(A,workd(:,col1),arguments_Afun);
                        %workd(:,col2) = max(workd(:,col2),0);
                    elseif nbArgumentsAfun == 2                        
                        workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2);
                    elseif nbArgumentsAfun == 3                        
                        workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2,arguments_Afun3);
                    else
                        workd(:,col2) = feval(A,workd(:,col1),varargin{indexArgumentsAfun});                    
                    end                      
                    %workd(:,col2) = tim_w_times_x_c(workd(:,col1),arguments_Afun); %slower
                   
                else
                    % OP = R'\(A*(R\x))
                    if issparse(B)
                        workd(permB,col2) = RB \ workd(:,col1);
                        workd(:,col2) = feval(A,workd(:,col2),arguments_Afun);
                        workd(:,col2) = RBT \ workd(permB,col2);

                    else
                        workd(:,col2) = RBT \ feval(A,(RB\workd(:,col1)),arguments_Afun);
                    end
                end
            elseif (mode == 3)
                if isempty(B)
                    workd(:,col2) = feval(A,workd(:,col1), arguments_Afun);
                else
                    if (ido == -1)
                        if cholB
                            workd(:,col2) = RBT * (RB * workd(:,col1));
                        else
                            workd(:,col2) = B * workd(:,col1);
                        end
                        workd(:,col2) = feval(A,workd(:,col2), arguments_Afun);
                    elseif (ido == 1)
                        workd(:,col2) = feval(A,workd(:,col3), arguments_Afun);
                    end
                end
            else % mode is not 1 or 3
                error(['Unknown mode returned from ' prefix 'aupd.'])
            end
        end % if Amatrix
    case 2
        if (mode == 3)
            if cholB
                workd(:,col2) = RBT * (RB * workd(:,col1));
            else
                workd(:,col2) = B * workd(:,col1);
            end
        else
            error(['Unknown mode returned from ' prefix 'aupd.'])
        end
    case 3
        % setting iparam(1) = ishift = 1 ensures this never happens
        warning('MATLAB:eigs:WorklShiftsUnsupported', ...
            ['EIGS does not yet support computing the shifts in workl.' ...
            ' Setting reverse communication parameter to 99 and returning.'])
        ido = int32(99);
    case 99
    otherwise
        error(['Unknown value of reverse communication parameter' ...
                ' returned from ' prefix 'aupd.'])      
    end
    
    cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X)
    
    %tim
    if nombreIterations ~= double(ipntr(15))
        nombreIterations = double(ipntr(15));
    end

    if display >= 1 && display <=2 
        iter = double(ipntr(15));
        if (iter > ariter) & (ido ~= 99)
            ariter = iter;
            ds = sprintf(['Iteration %d: a few Ritz values of the' ...
                    ' %d-by-%d matrix:'],iter,p,p);
            disp(ds)
            if isrealprob
                if issymA
                    dispvec = [workl(double(ipntr(6))+(0:p-1))];
                    if isequal(whch,'BE')
                        % roughly k Large eigenvalues and k Small eigenvalues
                        disp(dispvec(max(end-2*k+1,1):end))
                    else
                        % k eigenvalues
                        disp(dispvec(max(end-k+1,1):end))
                    end
                else
                    dispvec = [complex(workl(double(ipntr(6))+(0:p-1)), ...
                            workl(double(ipntr(7))+(0:p-1)))];
                    % k+1 eigenvalues (keep complex conjugate pairs together)
                    disp(dispvec(max(end-k,1):end))
                end
            else
                dispvec = [complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ...
                        workl(2*double(ipntr(6))+(0:2:2*(p-1))))];
                disp(dispvec(max(end-k+1,1):end))
            end
        end
    end
    
end % while (ido ~= 99)

t0 = cputime; % start timing post-processing

flag = 0;
if (info < 0)
    es = sprintf('Error with ARPACK routine %saupd: info = %d',prefix,full(info));
    error(es)
else
    if (nargout >= 2)
        rvec = int32(1); % compute eigenvectors
    else
        rvec = int32(0); % do not compute eigenvectors
    end
    
    if isrealprob
        if issymA
            arpackc( 'dseupd', rvec, 'A', select, ...
                d, v, ldv, sigma, ...
                bmat, int32(n), whch, nev, tol, resid, ncv, ...
                v, ldv, iparam, ipntr, workd, workl, lworkl, info );
            if isequal(whch,'LM') | isequal(whch,'LA')
                d = flipud(d);
                if (rvec == 1)
                    v(:,1:k) = v(:,k:-1:1);
                end
            end
            if ((isequal(whch,'SM') | isequal(whch,'SA')) & (rvec == 0))
                d = flipud(d);
            end
        else
            arpackc( 'dneupd', rvec, 'A', select, ...
                d, di, v, ldv, sigmar, sigmai, workev, ...
                bmat, int32(n), whch, nev, tol, resid, ncv, ...
                v, ldv, iparam, ipntr, workd, workl, lworkl, info );
            d = complex(d,di);
            if rvec
                d(k+1) = [];
            else
                zind = find(d == 0);
                if isempty(zind)
                    d = d(k+1:-1:2);
                else
                    d(max(zind)) = [];
                    d = flipud(d);
                end
            end
        end
    else
        zsigma = [real(sigma); imag(sigma)];
        arpackc( 'zneupd', rvec, 'A', select, ...
            zd, zv, ldv, zsigma, workev, ...
            bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
            ldv, iparam, ipntr, zworkd, workl, lworkl, ...
            rwork, info );
        if issymA
            d = zd(1:2:end-1);
        else
            d = complex(zd(1:2:end-1),zd(2:2:end));
        end
        v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]);
    end
    
    if (info ~= 0)
        es = ['Error with ARPACK routine ' prefix 'eupd: '];
        switch double(info)
            case 2
            ss = sum(select);
            if (ss < k)
            es = [es ...
                    '  The logical variable select was only set with ' int2str(ss) ...
                    ' 1''s instead of nconv=' int2str(double(iparam(5))) ...
                    ' (k=' int2str(k) ').' ...
                    ' Please contact the ARPACK authors at [e-mail].'];
            else
            es = [es ...
                    'The LAPACK reordering routine ' prefix(1) ...
                    'trsen did not return all ' int2str(k) ' eigenvalues.']
            end
            case 1
            es = [es ...
                    'The Schur form could not be reordered by the LAPACK routine ' ...
                    prefix(1) 'trsen.' ...
                    ' Please contact the ARPACK authors at [e-mail].'];
            case -14
            es = [es prefix ...
                    'aupd did not find any eigenvalues to sufficient accuracy.'];
            otherwise
            es = [es sprintf('info = %d',full(info))];
        end
        error(es)
    else
        nconv = double(iparam(5));
        if (nconv == 0)
            if (nargout < 3)
                warning('MATLAB:eigs:NoEigsConverged', ...
                    'None of the %d requested eigenvalues converged.',k)
            else
                flag = 1;
            end
        elseif (nconv < k)
            if (nargout < 3)
                warning('MATLAB:eigs:NotAllEigsConverged', ...
                    'Only %d of the %d requested eigenvalues converged.',nconv,k)
            else
                flag = 1;
            end
        end
    end % if (info ~= 0)
end % if (info < 0)

if (issymA) | (~isrealprob)
    if (nargout <= 1)
        if isrealprob
            varargout{1} = d;
        else
            varargout{1} = d(k:-1:1,1);
        end
    else
        varargout{1} = v(:,1:k);
        varargout{2} = diag(d(1:k,1));
        if (nargout >= 3)
            varargout{3} = flag;
        end
    end
else
    if (nargout <= 1)
        varargout{1} = d;
    else
        cplxd = find(di ~= 0);
        % complex conjugate pairs of eigenvalues occur together
        cplxd = cplxd(1:2:end);
        v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ...
                complex(v(:,cplxd),-v(:,cplxd+1))];
        varargout{1} = v(:,1:k);
        varargout{2} = diag(d);
        if (nargout >= 3)
            varargout{3} = flag;
        end
    end
end

if (nargout >= 2) & (mode == 1) & ~isempty(B)
    if issparse(B)
        varargout{1}(permB,:) = RB \ varargout{1};
    else
        varargout{1} = RB \ varargout{1};
    end
end

cputms(4) = cputime-t0; % end timing post-processing

cputms(5) = sum(cputms(1:4)); % total time

if (display >= 2) %tim
    if (mode == 1)
        innerstr = sprintf(['Compute A*X:' ...
                '                               %f\n'],cputms(3));
    elseif (mode == 2)
        innerstr = sprintf(['Compute A*X and solve B*X=Y for X:' ...
                '         %f\n'],cputms(3));
    elseif (mode == 3)
        if isempty(B)
            innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ...
                    '               %f\n'],cputms(3));
        else
            innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ...
                    '             %f\n'],cputms(3));
        end
    end
    if ((mode == 3) & (Amatrix))
        if isempty(B)
            prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ...
                    '   %f\n'],cputms(1));
        else
            prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ...
                    '   %f\n'],cputms(1));
        end
    elseif ((mode == 2) & (~cholB))
        prepstr = sprintf(['Pre-processing, including chol(B):' ...
                '         %f\n'],cputms(1));
    else
        prepstr = sprintf(['Pre-processing:' ...
                '                            %f\n'],cputms(1));
    end
    sstr = sprintf(['***********CPU Timing Results in seconds***********']);
    ds = sprintf(['\n' sstr '\n' ...
            prepstr ...
            'ARPACK''s %saupd:                           %f\n' ...
            innerstr ...
            'Post-processing with ARPACK''s %seupd:      %f\n' ...
            '***************************************************\n' ...
            'Total:                                     %f\n' ...
            sstr '\n'], ...
        prefix,cputms(2),prefix,cputms(4),cputms(5));
    disp(ds)
    if nombre_A_times_X > 0      %tim
        disp(sprintf('Number of iterations : %f\n',nombreIterations));
        disp(sprintf('Number of times A*X was computed : %f\n',nombre_A_times_X));
        disp(sprintf('Average time for A*X : %f\n',cputms(3)/nombre_A_times_X));          
    end
end