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 ofA
’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 functionsAfun
andBfun
instead of matrices.Afun
andBfun
are function handles.Afun(x)
andBfun(x)
return the matrix-vector productA*x
andB*x
.D = primme_eigs(...,k)
finds thek
largest magnitude eigenvalues.k
must be less than the dimension of the matrixA
.D = primme_eigs(...,k,target)
returnsk
eigenvalues such that: Iftarget
is a real number, it finds the closest eigenvalues totarget
. Iftarget
is'LA'
or'SA'
, eigenvalues with the largest or smallest algebraic value.'LM'
or'SM'
, eigenvalues with the largest or smallest magnitude ifOPTS.targetShifts
is empty. Iftarget
is a real or complex scalar including 0,primme_eigs()
finds the eigenvalues closest totarget
.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 magnitudeeig(A)
.k=1
,'SM'
,OPTS.targetShifts=[]
returns the smallest magnitudeeig(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 inOPTS.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)
(seeeps
) {\(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 byAfun
is real or complex {false}isdouble
: whether the class of in/out vectors inAfun
are double or single {false}isgpu
: whether the class of in/out vectors inAfun
aregpuArray
{false}ishermitian
: whetherA
is Hermitian; otherwise it is considered normal {true}targetShifts
: shifts for interior eigenvalues (seetarget
) {[]}v0
: any number of initial guesses to the eigenvectors (seeinitSize
{[]}orthoConst
: external orthogonalization constraints (seenumOrthoConst
{[]}locking
: 1, hard locking; 0, soft lockingp
: maximum size of the search subspace (seemaxBasisSize
)minRestartSize
: minimum Ritz vectors to keep in restartingmaxMatvecs
: maximum number of matrix vector multiplications {Inf}maxit
: maximum number of outer iterations (seemaxOuterIterations
) {Inf}maxPrevRetain
: number of Ritz vectors from previous iteration that are kept after restart {typically >0}robustShifts
: setting to true may avoid stagnation or misconvergencemaxInnerIterations
: maximum number of inner solver iterationsLeftQ
: use the locked vectors in the left projectorLeftX
: use the approx. eigenvector in the left projectorRightQ
: use the locked vectors in the right projectorRightX
: use the approx. eigenvector in the right projectorSkewQ
: use the preconditioned locked vectors in the right projectorSkewX
: use the preconditioned approx. eigenvector in the right projectorrelTolBase
: a legacy from classical JDQR (not recommended)convTest
: how to stop the inner QMR MethodconvTestFun
: function handler with an alternative convergence criterion. IfFUN(EVAL,EVEC,RNORM)
returns a nonzero value, the pair(EVAL,EVEC)
with residual normRNORM
is considered convergediseed
: random seedreturnUnconverged
: 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 preconditionerP
orP = P1*P2
to accelerate convergence of the method. ApplyingP\x
should approximate(A-sigma*eye(N))\x
, forsigma
near the wanted eigenvalue(s). IfP
is[]
then a preconditioner is not applied.P
may be a function handlePFUN
such thatPFUN(x)
returnsP\x
.[X,D] = primme_eigs(...)
returns a diagonal matrixD
with the eigenvalues and a matrixX
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 astruct
to report statistical information about number of matvecs, elapsed time, and estimates for the largest and smallest algebraic eigenvalues ofA
.[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 matvecsHIST(:,2)
: timeHIST(:,3)
: number of converged/locked pairsHIST(:,4)
: block indexHIST(:,5)
: approximate eigenvalueHIST(:,6)
: residual normHIST(:,7)
: QMR residual norm
OPTS.reportLevel
controls the granularity of the record. IfOPTS.reportLevel == 1
,HIST
has one row per converged eigenpair and only the first three columns together with the fifth and the sixth are reported. IfOPTS.reportLevel == 2
,HIST
has one row per outer iteration and converged value, and only the first six columns are reported. OtherwiseHIST
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 eitherHIST
is not returned orOPTS.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
See also: MATLAB eigs,
primme_svds()