Computes the singular value decomposition of a real matrix using Jacobi plane rotations.
The routine computes the singular value decomposition (SVD) of a real or complex
The SVD of
Ais written as
for real flavors, or
for complex flavors,
northonormal matrix, and
northogonal/unitary matrix. The diagonal elements of
Σare the singular values of
A; the columns of
Vare the left and right singular vectors of
A, respectively. The matrices
Vare computed and stored in the arrays
v, respectively. The diagonal of
Σis computed and stored in the array
?gesvjroutine can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines.
Vis 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 (
OVERFLOW). In extreme cases, even denormalized singular values can be computed with the corresponding gradual loss of accurate digit.
- Specifies whether matrix storage layout is row major (LAPACK_ROW_MAJOR) or column major (LAPACK_COL_MAJOR).
- 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.
- 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.If,jobu='N'uis not computed. However, see the description ofa.
- 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.
- The number of rows of the input matrixA.1/slamch('E')>form≥0sgesvj.1/dlamch('E')>form≥0dgesvj.
- The number of columns in the input matrixA;m≥n≥0.
- 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
- 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
- Ifjobv='A', the product of Jacobi rotations in?gesvjis applied to the firstmvrows ofv. See the description ofjobv.0.≤mv≤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)
- Array size 6. If,jobu= 'C', wherestat =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ε.
- 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. 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) is small.
- 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).
- Array sizen.If, depending on the valueinfo= 0scale=stat, wherescaleis the scaling factor:
If, the procedureinfo> 0?gesvjdid not converge in the given number of iterations (sweeps) andmay not be accurate.scale*sva(1:n)
- 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.
- 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.