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-processingflag = 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 arpack@caam.rice.edu'];            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 arpack@caam.rice.edu'];            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    endelse    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    endendif (nargout >= 2) & (mode == 1) & ~isempty(B)    varargout{1} = RB \ varargout{1};endcputms(4) = cputime-t0; % end timing post-processingcputms(5) = sum(cputms(1:4)); % total timeif (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;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%