MATLAB Interface¶

function [varargout] = primme_eigs(varargin)

primme_eigs() finds a few eigenvalues and their corresponding eigenvectors of a symmetric/Hermitian matrix, A, or of a generalized problem (A,B), by calling PRIMME.

D = primme_eigs(A) returns a vector of A’s 6 largest magnitude eigenvalues.

D = primme_eigs(A,B) returns a vector of the 6 largest magnitude eigenvalues of (A,B)

D = primme_eigs(Afun,Bfun,dim) D = primme_eigs(Afun,dim) accepts a functions Afun and Bfun instead of matrices. Afun and Bfun are function handles. Afun(x) and Bfun(x) return the matrix-vector product A*x and B*x.

D = primme_eigs(...,k) finds the k largest magnitude eigenvalues. k must be less than the dimension of the matrix A.

D = primme_eigs(...,k,target) returns k eigenvalues such that: If target is a real number, it finds the closest eigenvalues to target. If target is

'LA' or 'SA', eigenvalues with the largest or smallest algebraic value.

'LM' or 'SM', eigenvalues with the largest or smallest magnitude if OPTS.targetShifts is empty. If target is a real or complex scalar including 0, primme_eigs() finds the eigenvalues closest to target.

In addition, if some values are provided in OPTS.targetShifts, it finds eigenvalues that are farthest ('LM') or closest ('SM') in absolute value from the given values.

Examples:

k=1, 'LM', OPTS.targetShifts=[] returns the largest magnitude eig(A). k=1, 'SM', OPTS.targetShifts=[] returns the smallest magnitude eig(A). k=3, 'SM', OPTS.targetShifts=[2, 5] returns the closest eigenvalue in absolute sense to 2, and the two closest eigenvalues to 5.

'CLT' or 'CGT', find eigenvalues closest to but less or greater than the given values in OPTS.targetShifts.

D = primme_eigs(...,k,target,OPTS) specifies extra solver parameters. Some default values are indicated in brackets {}:

aNorm: the estimated 2-norm of A {0.0 (estimate the norm internally)}

tol: convergence tolerance: NORM(A*X(:,i)-X(:,i)*D(i,i)) < tol*NORM(A) (see eps) {\(10^4\) times the machine precision}

maxBlockSize: maximum block size (useful for high multiplicities) {1}

reportLevel: reporting level (0-3) (see HIST) {no reporting 0}

display: whether displaying reporting on screen (see HIST) {0 if HIST provided}

isreal: whether A represented by Afun is real or complex {false}

isdouble: whether the class of in/out vectors in Afun are double or single {false}

isgpu: whether the class of in/out vectors in Afun are gpuArray {false}

ishermitian: whether A is Hermitian; otherwise it is considered normal {true}

targetShifts: shifts for interior eigenvalues (see target) {[]}

v0: any number of initial guesses to the eigenvectors (see initSize {[]}

orthoConst: external orthogonalization constraints (see numOrthoConst {[]}

locking: 1, hard locking; 0, soft locking

p: maximum size of the search subspace (see maxBasisSize)

minRestartSize: minimum Ritz vectors to keep in restarting

maxMatvecs: maximum number of matrix vector multiplications {Inf}

maxit: maximum number of outer iterations (see maxOuterIterations) {Inf}

maxPrevRetain: number of Ritz vectors from previous iteration that are kept after restart {typically >0}

robustShifts: setting to true may avoid stagnation or misconvergence

maxInnerIterations: maximum number of inner solver iterations

LeftQ: use the locked vectors in the left projector

LeftX: use the approx. eigenvector in the left projector

RightQ: use the locked vectors in the right projector

RightX: use the approx. eigenvector in the right projector

SkewQ: use the preconditioned locked vectors in the right projector

SkewX: use the preconditioned approx. eigenvector in the right projector

relTolBase: a legacy from classical JDQR (not recommended)

convTest: how to stop the inner QMR Method

convTestFun: function handler with an alternative convergence criterion. If FUN(EVAL,EVEC,RNORM) returns a nonzero value, the pair (EVAL,EVEC) with residual norm RNORM is considered converged

iseed: random seed

returnUnconverged: whether to return unconverged pairs if maximum iterations or matvecs is reached.

D = primme_eigs(A,k,target,OPTS,METHOD) specifies the eigensolver method. METHOD can be one of the next strings:

‘PRIMME_DYNAMIC’, (default) switches dynamically to the best method

‘PRIMME_DEFAULT_MIN_TIME’, best method for low-cost matrix-vector product

‘PRIMME_DEFAULT_MIN_MATVECS’, best method for heavy matvec/preconditioner

‘PRIMME_Arnoldi’, Arnoldi not implemented efficiently

‘PRIMME_GD’, classical block Generalized Davidson

‘PRIMME_GD_plusK’, GD+k block GD with recurrence restarting

‘PRIMME_GD_Olsen_plusK’, GD+k with approximate Olsen precond.

‘PRIMME_JD_Olsen_plusK’, GD+k, exact Olsen (two precond per step)

‘PRIMME_RQI’, Rayleigh Quotient Iteration. Also INVIT, but for INVIT provide OPTS.targetShifts

‘PRIMME_JDQR’, Original block, Jacobi Davidson

‘PRIMME_JDQMR’, Our block JDQMR method (similar to JDCG)

‘PRIMME_JDQMR_ETol’, Slight, but efficient JDQMR modification

‘PRIMME_STEEPEST_DESCENT’, equivalent to GD(block,2*block)

‘PRIMME_LOBPCG_OrthoBasis’, equivalent to GD(nev,3*nev)+nev

‘PRIMME_LOBPCG_OrthoBasis_Window’ equivalent to GD(block,3*block)+block nev>block

D = primme_eigs(A,k,target,OPTS,METHOD,P)

D = primme_eigs(A,k,target,OPTS,METHOD,P1,P2) uses preconditioner P or P = P1*P2 to accelerate convergence of the method. Applying P\x should approximate (A-sigma*eye(N))\x, for sigma near the wanted eigenvalue(s). If P is [] then a preconditioner is not applied. P may be a function handle PFUN such that PFUN(x) returns P\x.

[X,D] = primme_eigs(...) returns a diagonal matrix D with the eigenvalues and a matrix X whose columns are the corresponding eigenvectors.

[X,D,R] = primme_eigs(...) also returns an array of the residual norms of the computed eigenpairs.

[X,D,R,STATS] = primme_eigs(...) returns a struct to report statistical information about number of matvecs, elapsed time, and estimates for the largest and smallest algebraic eigenvalues of A.

[X,D,R,STATS,HIST] = primme_eigs(...) it returns the convergence history, instead of printing it. Every row is a record, and the columns report:

HIST(:,1): number of matvecs

HIST(:,2): time

HIST(:,3): number of converged/locked pairs

HIST(:,4): block index

HIST(:,5): approximate eigenvalue

HIST(:,6): residual norm

HIST(:,7): QMR residual norm

OPTS.reportLevel controls the granularity of the record. If OPTS.reportLevel == 1, HIST has one row per converged eigenpair and only the first three columns together with the fifth and the sixth are reported. If OPTS.reportLevel == 2, HIST has one row per outer iteration and converged value, and only the first six columns are reported. Otherwise HIST has one row per QMR iteration, outer iteration and converged value, and all columns are reported.

The convergence history is displayed if OPTS.reportLevel > 0 and either HIST is not returned or OPTS.display == 1.

Examples:

A = diag(1:100);

d = primme_eigs(A,10) % the 10 largest magnitude eigenvalues

d = primme_eigs(A,10,'SM') % the 10 smallest magnitude eigenvalues

d = primme_eigs(A,10,25.0) % the 10 closest eigenvalues to 25.0

opts.targetShifts = [2 20];
d = primme_eigs(A,10,'SM',opts) % 1 eigenvalue closest to 2 and
                                % 9 eigenvalues closest to 20

B = diag(100:-1:1);
d = primme_eigs(A,B,10,'SM') % the 10 smallest magnitude eigenvalues

opts = struct();
opts.tol = 1e-4; % set tolerance
opts.maxBlockSize = 2; % set block size
[x,d] = primme_eigs(A,10,'SA',opts,'DEFAULT_MIN_TIME')

opts.orthoConst = x;
[d,rnorms] = primme_eigs(A,10,'SA',opts) % find another 10 with the default method

% Compute the 6 eigenvalues closest to 30.5 using ILU(0) as a preconditioner
% by passing the matrices L and U.
A = sparse(diag(1:50) + diag(ones(49,1), 1) + diag(ones(49,1), -1));
[L,U] = ilu(A, struct('type', 'nofill'));
d = primme_eigs(A, k, 30.5, [], [], L, U);

% Compute the 6 eigenvalues closest to 30.5 using Jacobi preconditioner
% by passing a function.
Pfun = @(x)(diag(A) - 30.5)\x;
d = primme_eigs(A,6,30.5,[],[],Pfun) % find the closest 5 to 30.5