aeigs60.m

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

                        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});
                else
                    % OP = R'\(A*(R\x))
                    if issparse(B)
                        workd(permB,col2) = RBT \ ...
                            feval(A,(RB\workd(permB,col1)), ...
                            varargin{7-Amatrix-Bnotthere:end});
                    else
                        workd(:,col2) = RBT \ ...
                            feval(A,(RB\workd(:,col1)), ...
                            varargin{7-Amatrix-Bnotthere:end});
                    end
                end
            elseif (mode == 3)
                if isempty(B)
                    workd(:,col2) = feval(A,workd(:,col1), ...
                        varargin{7-Amatrix-Bnotthere:end});
                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), ...
                            varargin{7-Amatrix-Bnotthere:end});
                    elseif (ido == 1)
                        workd(:,col2) = feval(A,workd(:,col3), ...
                            varargin{7-Amatrix-Bnotthere:end});
                    end
                end
            else % mode is not 1 or 3
                error(['Unknown mode returned from ' prefix 'aupd'])
            end
        end % if Amatrix
    elseif (ido == 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
    elseif (ido == 3)
        % setting iparam(1) = ishift = 1 ensures this never happens
        warning(['eigs does not yet support computing the shifts in workl.' ...
                ' Setting reverse communication parameter to 99 and returning'])
        ido = int32(99);
    elseif (ido ~= 99)
        error(['Unknown value of reverse communication parameter' ...
                ' returned from ' prefix 'aupd'])      
    end
    
    cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X)

    if display
        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(end-2*k+1:end))
                    else
                        % k eigenvalues
                        disp(dispvec(end-k+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(end-k: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(end-k+1:end))
            end
        end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Added by Tom Wright, October 2002 for EigTool interaction
%
%% Call the commands to update EigTool, and return if STOP has been pressed
    iter = double(ipntr(15));
    if (iter > guiiter) & (ido ~= 99)
      guiiter = iter;
      arpackgui_update;
      if return_now==1, varargout{1} = 'stopped'; return; end;
    end;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end % while (ido ~= 99)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Added by Tom Wright, October 2002 for EigTool interaction
%
% Duplicate the variables for later as v may get overwritten by the Ritz
% vectors; 'modify' an entry to ensure they are duplicated here, and not
% overwritten inside the arpackc mex file below (otherwise both v and gui_v
% will be destroyed).
if exist('v','var'), gui_v = v; gui_v(1) = gui_v(1)*1; end;
if exist('zv','var'), gui_zv = zv; gui_zv(1) = gui_zv(1)*1; end;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t0 = cputime; % start timing post-processing

flag = 0;
if (info < 0)
    es = sprintf('Error with ARPACK routine %saupd: info = %d',prefix,double(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',double(info))];
        end
        error(es)
    else
        nconv = double(iparam(5));
        if (nconv == 0)
            if (nargout < 3)
                ws = sprintf(['None of the %d requested eigenvalues' ...
                        ' converged'],k);
                warning(ws)
            else
                flag = 1;
            end
        elseif (nconv < k)
            if (nargout < 3)
                ws = sprintf(['Only %d of the %d requested eigenvalues' ...
                        ' converged'],nconv,k);
                warning(ws)
            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)
    varargout{1} = RB \ varargout{1};
end

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

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

if (display == 2)
    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)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Added by Tom Wright, October 2002 for EigTool interaction
%
setupagui;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%