endendif ~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    endenduseeig = 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']);endif isempty(B)    knstr = [knstr sprintf(['Using eig(full(A)) instead.'])];else    knstr = [knstr sprintf(['Using eig(full(A),full(B)) instead.'])];endif (k == 0)    useeig = 1;endif isrealprob & issymA    if (k > n-1)        if (n >= 6)            warning(knstr)        end        useeig = 1;    endelse    if (k > n-2)        if (n >= 7)            warning(knstr)        end        useeig = 1;    endendif isrealprob & issymA    if ~isreal(sigma)        error(['For real symmetric problems, eigenvalue shift sigma must' ...                ' be real'])    endelse    if ~isrealprob & issymA & ~isreal(sigma)        warning(['Complex eigenvalue shift sigma on a Hermitian problem' ...                ' (all real eigenvalues)'])    endendif isrealprob & issymA    if strcmp(whch,'LR')        whch = 'LA';        warning(['For real symmetric problems, sigma value ''LR''' ...                ' (Largest Real) is now ''LA'' (Largest Algebraic)'])    end    if strcmp(whch,'SR')        whch = 'SA';        warning(['For real symmetric problems, sigma value ''SR''' ...                ' (Smallest Real) is now ''SA'' (Smallest Algebraic)'])    end    if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'})        error(whchstr)    endelse    if strcmp(whch,'BE')        warning(['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'})        error(whchstr)    endend% Now have enough information to do early return on cases eigs does not handleif 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    returnendif isrealprob & ~issymA    sigmar = real(sigma);    sigmai = imag(sigma);endif isrealprob & issymA    if (p <= k)        error(['For real symmetric problems, must have number of' ...                ' basis vectors opts.p > k'])    endelse    if (p <= k+1)        error(['For nonsymmetric and complex problems, must have number of' ...                ' basis vectors opts.p > k+1'])    endendif 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;endif cholB    pB = 0;    RB = B;    RBT = B';endif (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        perm = colamd(AsB);        [L,U,P] = lu(AsB(:,perm));    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 = double(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'])];        disp(ds)        pause(2)    endendif (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)    endend% Allocate outputs and ARPACK work variablesif 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);    endelse % 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);endido = int32(0); % reverse communication parameterif isempty(B) | (~isempty(B) & (mode == 1))    bmat = 'I'; % standard eigenvalue problemelse    bmat = 'G'; % generalized eigenvalue problemendnev = int32(k); % number of eigenvalues requestedncv = int32(p); % number of Lanczos vectorsiparam = int32(zeros(11,1));iparam([1 3 7]) = int32([1 maxit mode]);select = int32(zeros(p,1));cputms(1) = cputime - t0; % end timing pre-processingiter = 0;ariter = 0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added by Tom Wright, October 2002 for EigTool interaction%guiiter = 0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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        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,double(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)));    if (row ~= 1)        str = sprintf(['ipntr(1)=%d does not refer to the start of a' ...                ' column of the %d-by-3 array workd'],double(ipntr(1)),n);        error(str)    end    [row,col2] = ind2sub([n,3],double(ipntr(2)));    if (row ~= 1)        str = sprintf(['ipntr(2)=%d does not refer to the start of a' ...                ' column of the %d-by-3 array workd'],double(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'],double(ipntr(3)),n);            error(str)        end    end        if ((ido == -1) | (ido == 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) = RBT \ (A * (RB \ workd(permB,col1)));                    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)));