?gesvj
?gesvj
Computes the singular value decomposition of a real matrix using Jacobi plane rotations.
Syntax
lapack_int
LAPACKE_sgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
float
*
a
,
lapack_int
lda
,
float
*
sva
,
lapack_int
mv
,
float
*
v
,
lapack_int
ldv
,
float
*
stat
);
lapack_int
LAPACKE_dgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
double
*
a
,
lapack_int
lda
,
double
*
sva
,
lapack_int
mv
,
double
*
v
,
lapack_int
ldv
,
double
*
stat
);
lapack_int
LAPACKE_cgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
lapack_complex_float
*
a
,
lapack_int
lda
,
float
*
sva
,
lapack_int
mv
,
lapack_complex_float
*
v
,
lapack_int
ldv
,
float
*
stat
);
lapack_int
LAPACKE_zgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
lapack_complex_double
*
a
,
lapack_int
lda
,
double
*
sva
,
lapack_int
mv
,
lapack_complex_double
*
v
,
lapack_int
ldv
,
double
*
stat
);
Include Files
- mkl.h
Description
The routine computes the singular value decomposition (SVD) of a real or complex
m
-by-n
matrix A
, where m
≥
n
. The SVD of
A
is written asA
= U
*Σ
*V
T
A
= U
*Σ
*V
H
where
Σ
is an m
-by-n
diagonal matrix, U
is an m
-by-n
orthonormal matrix, and V
is an n
-by-n
orthogonal/unitary matrix. The diagonal elements of Σ
are the singular values of A
; the columns of U
and V
are the left and right singular vectors of A
, respectively. The matrices U
and V
are computed and stored in the arrays u
and v
, respectively. The diagonal of Σ
is computed and
stored in the array sva
.The
?gesvj
routine can sometimes compute tiny singular values and their
singular vectors much more accurately than other SVD
routines.The
n
-by-n
orthogonal matrix V
is obtained as a product of Jacobi plane rotations. The rotations
are implemented as fast scaled rotations of Anda and Park [AndaPark94]. In the case of underflow of the Jacobi angle, a modified Jacobi transformation of Drmac ([Drmac08-4]) is used. Pivot strategy uses column interchanges of de Rijk ([deRijk98]). The relative accuracy of the computed singular values and the accuracy of the computed singular vectors (in angle metric) is as guaranteed by the theory of Demmel and Veselic [Demmel92]. The condition number that determines the accuracy in the full rank case is essentially
The computational range for the nonzero singular values is the machine number interval (
UNDERFLOW
,OVERFLOW
). In extreme cases, even denormalized singular values can be
computed with the corresponding gradual loss of accurate digit.Input Parameters
- matrix_layout
- Specifies whether matrix storage layout is row major (LAPACK_ROW_MAJOR) or column major (LAPACK_COL_MAJOR).
- joba
- Must be'L','U'or'G'.Specifies the structure ofA:If, the input matrixjoba='L'Ais lower triangular.If, the input matrixjoba='U'Ais upper triangular.If, the input matrixjoba='G'Ais a generalm-by-n,m≥n.
- jobu
- Must be'U','C'or'N'.Specifies whether to compute the left singular vectors (columns ofU):If, the left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns ofjobu='U'A. See more details in the description ofa. The default numerical orthogonality threshold is set to approximatelyTOL=CTOL*EPS, CTOL=sqrt(m), EPS =?lamch('E')If, analogous tojobu='C', except that you can control the level of numerical orthogonality of the computed left singular vectors.jobu='U'TOLcan be set toTOL=CTOL*EPS, whereCTOLis given on input in the arraystat. NoCTOLsmaller than ONE is allowed.CTOLgreater than1 / EPSis meaningless. The option'C'can be used ifis satisfactory orthogonality of the computed left singular vectors, som*EPSCTOL=could save a few sweeps of Jacobi rotations. See the descriptions ofmaandstat[0].If,jobu='N'uis not computed. However, see the description ofa.
- jobv
- Must be'V','A'or'N'.Specifies whether to compute the right singular vectors, that is, the matrixV:If, the matrixjobv='V'Vis computed and returned in the arrayv.If, the Jacobi rotations are applied to thejobv='A'mv-bynarrayv. In other words, the right singular vector matrixVis not computed explicitly, instead it is applied to anmv-bynmatrix initially stored in the firstmvrows ofV.If, the matrixjobv='N'Vis not computed and the arrayvis not referenced.
- m
- The number of rows of the input matrixA.1/slamch('E')>form≥0sgesvj.1/dlamch('E')>form≥0dgesvj.
- n
- The number of columns in the input matrixA;m≥n≥0.
- a,v
- Arraya(size at leastis an array containing thelda*nfor column major layout andlda*mfor row major layout)m-by-nmatrixA.Arrayv(size at least max(1,.ldv*n)) contains, ifthejobv= 'A'mv-by-nmatrix to be post-multiplied by Jacobi rotations
- lda
- The leading dimension of the arraya. Must be at leastmax(1,m)for column major layout and at least max(1,.n) for row major layout
- mv
- Ifjobv='A', the product of Jacobi rotations in?gesvjis applied to the firstmvrows ofv. See the description ofjobv.0.≤mv≤ldv
- ldv
- The leading dimension of the arrayv;.ldv≥1jobv='V',.ldv≥max(1,n)jobv='A',ldv≥max(1,mv)for column major layout and.for row major layoutldv≥max(1,n)
- stat
- Array size 6. If,jobu= 'C', wherestat[0] =CTOLCTOLdefines the threshold for convergence. The process stops if all columns ofAare mutually orthogonal up toCTOL*EPS, whereEPS=. It is required that?lamch('E')CTOL≥1 - that is, it is not allowed to force the routine to obtain orthogonality belowε.
Output Parameters
- a
- On exit:Iforjobu='U':jobu='C'
- if, the leading columns ofinfo= 0Acontain left singular vectors corresponding to the computed singular values ofathat are above the underflow threshold?lamch('S'), that is, non-zero singular values. The number of the computed non-zero singular values is returned instat[1]. Also see the descriptions ofsvaandstat. The computed columns ofuare mutually numerically orthogonal up to approximatelyTOL=sqrt((default); orm)*EPSTOL=CTOL*EPS, see the description ofjobu='C'jobu.
- if, the procedureinfo> 0?gesvjdid not converge in the given number of iterations (sweeps). In that case, the computed columns ofumay not be orthogonal up toTOL. The outputu(stored ina),sigma(given by the computed singular values in) andsva(1:n)vis still a decomposition of the input matrixAin the sense that the residual||for real flavors orA-scale*U*sigma*VT||2/ ||A||2||for complex flavors (whereA-scale*U*sigma*VH||2/ ||A||2scale =stat[0]) is small.
If:jobu='N'- if, note that the left singular vectors are 'for free' in the one-sided Jacobi SVD algorithm. However, if only the singular values are needed, the level of numerical orthogonality ofinfo= 0uis not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately. Thus, on exit,m*EPSacontains the columns ofuscaled with the corresponding singular values.
- if, the procedureinfo> 0?gesvjdid not converge in the given number of iterations (sweeps).
- sva
- Array sizen.If, depending on the valueinfo= 0scale=stat[0], wherescaleis the scaling factor:
- if,scale= 1sva[0:contains the computed singular values ofn- 1]a.
- if, the singular values ofscale≠1aare, and this factored representation is due to the fact that some of the singular values ofscale*sva(1:n)amight underflow or overflow.
If, the procedureinfo> 0?gesvjdid not converge in the given number of iterations (sweeps) andmay not be accurate.scale*sva(1:n) - v
- On exit:If, contains thejobv='V'n-by-nmatrix of the right singular vectors.If, thenjobv='A'vcontains the product of the computed right singular vector matrix and the initial matrix in the arrayv.If,jobv='N'vis not referenced.
- stat
- On exit,stat[0]=scaleis the scaling factor such thatare the computed singular values ofscale*sva(1:n)A. See the description ofsva.stat[1]is the number of the computed nonzero singular values.stat[2]is the number of the computed singular values that are larger than the underflow threshold.stat[3]is the number of sweeps of Jacobi rotations needed for numerical convergence.stat[4]=max_{iin the last sweep. This is useful information in cases when≠j} |COS(A(:,i),A(:,j))|?gesvjdid not converge, as it can be used to estimate whether the output is still useful and for post festum analysis.stat[5]is the largest absolute value over all sines of the Jacobi rotation angles in the last sweep. It can be useful in a post festum analysis.
Return Values
This function returns a value
info
.If , the execution is successful.
info
=0If , the
info
= -i
i
-th parameter had an illegal value.If , the function did not converge in the maximal number (30) of sweeps. The output may still be useful. See the description of
info
> 0stat
.