aeigs60.m

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

    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(knstr)
        end
        useeig = 1;
    end
else
    if (k > n-2)
        if (n >= 7)
            warning(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(['Complex eigenvalue shift sigma on a Hermitian problem' ...
                ' (all real eigenvalues)'])
    end
end

if 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)
    end
else
    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)
    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
        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)
    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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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)));