MATLAB Interface¶
-
function [varargout] = primme_svds(varargin)
primme_svds()
finds a few singular values and vectors of a matrixA
by calling PRIMME.A
is typically large and sparse.S = primme_svds(A)
returns a vector with the 6 largest singular values ofA
.S = primme_svds(AFUN,M,N)
accepts the function handleAFUN
to perform the matrix vector products with an M-by-N matrixA
.AFUN(X,'notransp')
returnsA*X
whileAFUN(X,'transp')
returnsA’*X
. In all the following,A
can be replaced byAFUN,M,N
.S = primme_svds(A,k)
computes thek
largest singular values ofA
.S = primme_svds(A,k,sigma)
computes thek
singular values closest to the scalar shiftsigma
.If
sigma
is a vector, find the singular valueS(i)
closest to eachsigma(i)
, fori<=k
.If
sigma
is'L'
, it computes the largest singular values.if
sigma
is'S'
, it computes the smallest singular values.
S = primme_svds(A,k,sigma,OPTIONS)
specifies extra solver parameters. Some default values are indicated in brackets {}:aNorm
: estimation of the 2-norm ofA
{0.0 (estimate the norm internally)}tol
: convergence toleranceNORM([A*V-U*S;A'*U-V*S]) <= tol * NORM(A)
(seeeps
) {1e-10
for double precision and1e-3
for single precision}maxit
: maximum number of matvecs withA
andA'
(seemaxMatvecs
) {inf}p
: maximum basis size (seemaxBasisSize
)reportLevel
: reporting level (0-3) (see HIST) {no reporting 0}display
: whether displaying reporting on screen (see HIST) {0 if HIST provided}isreal
: if 0, the matrix is complex; else it’s real {0: complex}isdouble
: if 0, the matrix is single; else it’s double {1: double}method
: which equivalent eigenproblem to solve‘
primme_svds_normalequations
’:A'*A
orA*A'
‘
primme_svds_augmented
’:[0 A';A 0]
‘
primme_svds_hybrid
’: first normal equations and then augmented (default)
u0
: initial guesses to the left singular vectors (seeinitSize
) {[]}v0
: initial guesses to the right singular vectors {[]}orthoConst
: external orthogonalization constraints (seenumOrthoConst
) {[]}locking
: 1, hard locking; 0, soft lockingmaxBlockSize
: maximum block sizeiseed
: random seedprimme
: options for first stage solverprimmeStage2
: options for second stage solverconvTestFun
: function handler with an alternative convergence criterion. IfFUN(SVAL,LSVEC,RSVEC,RNORM)
returns a nonzero value, the triplet(SVAL,LSVEC,RSVEC)
with residual normRNORM
is considered converged.
The available options for
OPTIONS.primme
andprimmeStage2
are the same asprimme_eigs()
, plus the option'method'
.S = primme_svds(A,k,sigma,OPTIONS,P)
applies a preconditionerP
as follows:If
P
is a matrix it appliesP\X
andP'\X
to approximateA\X
andA'\X
.If
P
is a function handle,PFUN
,PFUN(X,'notransp')
returnsP\X
andPFUN(X,'transp')
returnsP’\X
, approximatingA\X
andA'\X
respectively.- If
P
is astruct
, it can have one or more of the following fields: P.AHA\X
orP.AHA(X)
returns an approximation of(A'*A)\X
,P.AAH\X
orP.AAH(X)
returns an approximation of(A*A')\X
,P.aug\X
orP.aug(X)
returns an approximation of[zeros(N,N) A';A zeros(M,M)]\X
.
- If
If
P
is[]
then no preconditioner is applied.
S = primme_svds(A,k,sigma,OPTIONS,P1,P2
) applies a factorized preconditioner:If both
P1
andP2
are nonempty, apply(P1*P2)\X
to approximateA\X
.If
P1
is[]
andP2
is nonempty, then(P2'*P2)\X
approximatesA'*A
.P2
can be the R factor of an (incomplete) QR factorization ofA
or the L factor of an (incomplete) LL’ factorization ofA'*A
(RIF).If both
P1
andP2
are[]
then no preconditioner is applied.
[U,S,V] = primme_svds(...)
returns also the corresponding singular vectors. IfA
is M-by-N andk
singular triplets are computed, thenU
is M-by-k with orthonormal columns,S
is k-by-k diagonal, andV
is N-by-k with orthonormal columns.[S,R] = primme_svds(...)
[U,S,V,R] = primme_svds(...)
returns the residual norm of eachk
triplet,NORM([A*V(:,i)-S(i,i)*U(:,i); A'*U(:,i)-S(i,i)*V(:,i)])
.[U,S,V,R,STATS] = primme_svds(...)
returns how many timesA
andP
were used and elapsed time. The application ofA
is counted independently from the application ofA'
.[U,S,V,R,STATS,HIST] = primme_svds(...)
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 tripletsHIST(:,4)
: stageHIST(:,5)
: block indexHIST(:,6)
: approximate singular valueHIST(:,7)
: residual normHIST(:,8)
: 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:50); A(200,1) = 0; % rectangular matrix of size 200x50 s = primme_svds(A,10) % the 10 largest singular values s = primme_svds(A,10,'S') % the 10 smallest singular values s = primme_svds(A,10,25) % the 10 closest singular values to 25 opts = struct(); opts.tol = 1e-4; % set tolerance opts.method = 'primme_svds_normalequations' % set svd solver method opts.primme.method = 'DEFAULT_MIN_TIME' % set first stage eigensolver method opts.primme.maxBlockSize = 2; % set block size for first stage [u,s,v] = primme_svds(A,10,'S',opts); % find 10 smallest svd triplets opts.orthoConst = {u,v}; [s,rnorms] = primme_svds(A,10,'S',opts) % find another 10 % Compute the 5 smallest singular values of a rectangular matrix using % Jacobi preconditioner on (A'*A) A = sparse(diag(1:50) + diag(ones(49,1), 1)); A(200,50) = 1; % size(A)=[200 50] P = diag(sum(abs(A).^2)); precond.AHA = @(x)P\x; s = primme_svds(A,5,'S',[],precond) % find the 5 smallest values % Estimation of the smallest singular value A = diag([1 repmat(2,1,1000) 3:100]); [~,sval,~,rnorm] = primme_svds(A,1,'S',struct('convTestFun',@(s,u,v,r)r<s*.1)); sval - rnorm % approximate smallest singular value
See also: MATLAB svds,
primme_eigs()