?gejsv

Developer Reference for Intel® oneAPI Math Kernel Library for Fortran

Download PDF

ID 766686

Date 12/16/2022

Version

Public

A newer version of this document is available. Customers should click here to go to the newest version.

Visible to Intel only — GUID: GUID-79AE95B6-AC25-42A9-A7BC-811C9BEE825E

View Details

?gejsv

Computes the singular value decomposition using a preconditioned Jacobi SVD method.

Syntax

call sgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, work, lwork, iwork, info)

call dgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, work, lwork, iwork, info)

call cgejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, cwork, lwork, rwork, lrwork, iwork, info )

call zgejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, cwork, lwork, rwork, lrwork, iwork, info )

Include Files

mkl.fi

Description

The routine computes the singular value decomposition (SVD) of a real/complex m-by-n matrix A, where m≥n.

The SVD is written as

A = U*Σ*V^T, for real routines

A = U*Σ*V^H, for complex routines

where Σ is an m-by-n matrix which is zero except for its n diagonal elements, U is an m-by-n (or m-by-m) orthonormal matrix, and V is an n-by-n orthogonal 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 ?gejsv routine can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines.

The routine implements a preconditioned Jacobi SVD algorithm. It uses ?geqp3, ?geqrf, and ?gelqf as preprocessors and preconditioners. Optionally, an additional row pivoting can be used as a preprocessor, which in some cases results in much higher accuracy. An example is matrix A with the structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned diagonal matrices and C is a well-conditioned matrix. In that case, complete pivoting in the first QR factorizations provides accuracy dependent on the condition number of C, and independent of D1, D2. Such higher accuracy is not completely understood theoretically, but it works well in practice.

If A can be written as A = B*D, with well-conditioned B and some diagonal D, then the high accuracy is guaranteed, both theoretically and in software, independent of D. For more details see [Drmac08-1], [Drmac08-2].

The computational range for the singular values can be the full range ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS and LAPACK routines called by ?gejsv are implemented to work in that range. If that is not the case, the restriction for safe computation with the singular values in the range of normalized IEEE numbers is that the spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not overflow. This code (?gejsv) is best used in this restricted range, meaning that singular values of magnitude below ||A||_2 / slamch('O') (for single precision) or ||A||_2 / dlamch('O') (for double precision) are returned as zeros. See jobr for details on this.

This implementation is slower than the one described in [Drmac08-1], [Drmac08-2] due to replacement of some non-LAPACK components, and because the choice of some tuning parameters in the iterative part (?gesvj) is left to the implementer on a particular machine.

The rank revealing QR factorization (in this code: ?geqp3) should be implemented as in [Drmac08-3].

If m is much larger than n, it is obvious that the inital QRF with column pivoting can be preprocessed by the QRF without pivoting. That well known trick is not used in ?gejsv because in some cases heavy row weighting can be treated with complete pivoting. The overhead in cases m much larger than n is then only due to pivoting, but the benefits in accuracy have prevailed. You can incorporate this extra QRF step easily and also improve data movement (matrix transpose, matrix copy, matrix transposed copy) - this implementation of ?gejsv uses only the simplest, naive data movement.

Product and Performance Information
Performance varies by use, configuration and other factors. Learn more at www.Intel.com/PerformanceIndex. Notice revision #20201201

Product and Performance Information

Performance varies by use, configuration and other factors. Learn more at www.Intel.com/PerformanceIndex.

Notice revision #20201201

Input Parameters

joba

CHARACTER*1. Must be 'C', 'E', 'F', 'G', 'A', or 'R'.

Specifies the level of accuracy:

If joba = 'C', high relative accuracy is achieved if A = B*D with well-conditioned B and arbitrary diagonal matrix D. The accuracy cannot be spoiled by column scaling. The accuracy of the computed output depends on the condition of B, and the procedure aims at the best theoretical accuracy. The relative error max_{i=1:N}|d sigma_i| / sigma_i is bounded by f(M,N)*epsilon* cond(B), independent of D. The input matrix is preprocessed with the QRF with column pivoting. This initial preprocessing and preconditioning by a rank revealing QR factorization is common for all values of joba. Additional actions are specified as follows:

If joba = 'E', computation as with 'C' with an additional estimate of the condition number of B. It provides a realistic error bound.

If joba = 'F', accuracy higher than in the 'C' option is achieved, if A = D1*C*D2 with ill-conditioned diagonal scalings D1, D2, and a well-conditioned matrix C. This option is advisable, if the structure of the input matrix is not known and relative accuracy is desirable. The input matrix A is preprocessed with QR factorization with full (row and column) pivoting.

If joba = 'G', computation as with 'F' with an additional estimate of the condition number of B, where A = B*D. If A has heavily weighted rows, using this condition number gives too pessimistic error bound.

If joba = 'A', small singular values are the noise and the matrix is treated as numerically rank defficient. The error in the computed singular values is bounded by f(m,n)*epsilon*||A||. The computed SVD A = U*S*V**t (for real flavors) or A = U*S*V**H (for complex flavors) restores A up to f(m,n)*epsilon*||A||. This enables the procedure to set all singular values below n*epsilon*||A|| to zero.

If joba = 'R', the procedure is similar to the 'A' option. Rank revealing property of the initial QR factorization is used to reveal (using triangular factor) a gap sigma_{r+1} < epsilon * sigma_r, in which case the numerical rank is declared to be r. The SVD is computed with absolute error bounds, but more accurately than with 'A'.

jobu

CHARACTER*1. Must be 'U', 'F', 'W', or 'N'.

Specifies whether to compute the columns of the matrix U:

If jobu = 'U', n columns of U are returned in the array u

If jobu = 'F', a full set of m left singular vectors is returned in the array u.

If jobu = 'W', u may be used as workspace of length m*n. See the description of u.

If jobu = 'N', u is not computed.

jobv

CHARACTER*1. Must be 'V', 'J', 'W', or 'N'.

Specifies whether to compute the matrix V:

If jobv = 'V', n columns of V are returned in the array v; Jacobi rotations are not explicitly accumulated.

If jobv = 'J', n columns of V are returned in the array v but they are computed as the product of Jacobi rotations. This option is allowed only if jobu≠'N'

If jobv = 'W', v may be used as workspace of length n*n. See the description of v.

If jobv = 'N', v is not computed.

jobr

CHARACTER*1. Must be 'N' or 'R'.

Specifies the range for the singular values. If small positive singular values are outside the specified range, they may be set to zero. If A is scaled so that the largest singular value of the scaled matrix is around sqrt(big), big = ?lamch('O'), the function can remove columns of A whose norm in the scaled matrix is less than sqrt(?lamch('S')) (for jobr = 'R'), or less than small = ?lamch('S')/?lamch('E').

If jobr = 'N', the function does not remove small columns of the scaled matrix. This option assumes that BLAS and QR factorizations and triangular solvers are implemented to work in that range. If the condition of A if greater that big, use ?gesvj.

If jobr = 'R', restricted range for singular values of the scaled matrix A is [sqrt(?lamch('S'), sqrt(big)], roughly as described above. This option is recommended.

For computing the singular values in the full range [?lamch('S'),big], use ?gesvj.

jobt

CHARACTER*1. Must be 'T' or 'N'.

If the matrix is square, the procedure may determine to use a transposed A if A^T (for real flavors) or A^H (for complex flavors) seems to be better with respect to convergence. If the matrix is not square, jobt is ignored.

The decision is based on two values of entropy over the adjoint orbit of A^T * A (for real flavors) or A^H * A (for complex flavors). See the descriptions of work(6) and work(7).

If jobt = 'T', the function performs transposition if the entropy test indicates possibly faster convergence of the Jacobi process, if A is taken as input. If A is replaced with A^T or A^H, the row pivoting is included automatically.

If jobt = 'N', the functions attempts no speculations. This option can be used to compute only the singular values, or the full SVD (u, sigma, and v). For only one set of singular vectors (u or v), the caller should provide both u and v, as one of the arrays is used as workspace if the matrix A is transposed. The implementer can easily remove this constraint and make the code more complicated. See the descriptions of u and v.

CAUTION:

The jobt = 'T' option is experimental and its effect might not be the same in subsequent releases. Consider using the jobt = 'N' instead.

jobp

CHARACTER*1. Must be 'P' or 'N'.

Enables structured perturbations of denormalized numbers. This option should be active if the denormals are poorly implemented, causing slow computation, especially in cases of fast convergence. For details, see [Drmac08-1], [Drmac08-2] . For simplicity, such perturbations are included only when the full SVD or only the singular values are requested. You can add the perturbation for the cases of computing one set of singular vectors.

If jobp = 'P', the function introduces perturbation.

If jobp = 'N', the function introduces no perturbation.

m

INTEGER. The number of rows of the input matrix A; m≥ 0.

n

INTEGER. The number of columns in the input matrix A; m≥n≥ 0.

a, u, v

REAL for sgejsv

DOUBLE PRECISION for dgejsv.

COMPLEX for cgejsv

DOUBLE COMPLEX for zgejsv

Array a(lda,n) is an array containing the m-by-n matrix A.

u is a workspace array, its size is (ldu,*); the second dimension of u must be m if jobu = 'F' and n otherwise, . When jobt = 'T' and m = n, u must be provided even though jobu = 'N'.

v is a workspace array, its size is (ldv,n). When jobt = 'T' and m = n, v must be provided even though jobv = 'N'.

lda

INTEGER. The leading dimension of the array a. Must be at least max(1, m) .

sva

REAL for sgejsv

DOUBLE PRECISION for dgejsv.

REAL for cgejsv

DOUBLE PRECISION for zgejsv

sva is a workspace array, its size is n.

ldu

INTEGER. The leading dimension of the array u; ldu≥ 1.

jobu = 'U' or 'F' or 'W', ldu≥m for column major layout.

ldv

INTEGER. The leading dimension of the array v; ldv≥ 1.

jobv = 'V' or 'J' or 'W', ldv≥n.

work

REAL for sgejsv

DOUBLE PRECISION for dgejsv.

work is a workspace array, its dimension max(7, lwork).

lwork

INTEGER.

For real flavors:

Length of work to confirm proper allocation of work space. lwork depends on the task performed:

If only sigma is needed (jobu = 'N', jobv = 'N') and

no scaled condition estimate is required, then lwork≥ max(2*m+n,4*n+1,7). This is the minimal requirement. For optimal performance (blocked code) the optimal value is lwork≥ max(2*m+n,3*n+(n+1)*nb,7). Here nb is the optimal block size for ?geqp3/?geqrf.

In general, the optimal length lwork is computed as

lwork≥ max(2*m+n,n+lwork(sgeqp3),n+lwork(sgeqrf),7) for sgejsv

lwork≥ max(2*m+n,n+lwork(dgeqp3),n+lwork(dgeqrf),7) for dgejsv
... an estimate of the scaled condition number of A is required (joba = 'E', 'G'). In this case, lwork is the maximum of the above and n*n+4*n, that is, lwork≥ max(2*m+n,n*n+4*n,7). For optimal performance (blocked code) the optimal value is lwork≥ max(2*m+n,3*n+(n+1)*nb, n*n+4*n, 7).

In general, the optimal length lwork is computed as

lwork≥ max(2*m+n,n+lwork(sgeqp3),n+lwork(sgeqrf),n+n*n+lwork(spocon, 7) for sgejsv

lwork≥ max(2*m+n,n+lwork(dgeqp3),n+lwork(dgeqrf),n+n*n+lwork(dpocon, 7) for dgejsv

If sigma and the right singular vectors are needed (jobv = 'V'),

the minimal requirement is lwork≥ max(2*m+n,4*n+1,7).
for optimal performance, lwork≥ max(2*m+n,3*n+(n+1)*nb,7), where nb is the optimal block size for ?geqp3, ?geqrf, ?gelqf, ?ormlq. In general, the optimal length lwork is computed as

lwork≥ max(2*m+n, n+lwork(sgeqp3), n+lwork(spocon), n+lwork(sgelqf), 2*n+lwork(sgeqrf), n+lwork(sormlq) for sgejsv

lwork≥ max(2*m+n, n+lwork(dgeqp3), n+lwork(dpocon), n+lwork(dgelqf), 2*n+lwork(dgeqrf), n+lwork(dormlq) for dgejsv

If sigma and the left singular vectors are needed

the minimal requirement is lwork≥ max(2*n+m,4*n+1,7).
for optimal performance,

if jobu = 'U' :: lwork≥ max(2*m+n,3*n+(n+1)*nb, 7),

if jobu = 'F' :: lwork≥ max(2*m+n,3*n+(n+1)*nb, n+m*nb, 7),

where nb is the optimal block size for ?geqp3, ?geqrf, ?ormlq . In general, the optimal length lwork is computed as

lwork≥ max(2*m+n, n+lwork(sgeqp3), n+lwork(spocon), 2*n+lwork(sgeqrf), n+lwork(sormlq) for sgejsv

lwork≥ max(2*m+n, n+lwork(dgeqp3), n+lwork(dpocon), 2*n+lwork(dgeqrf), n+lwork(dormlq) for dgejsv

Here lwork(?ormlq) equals n*nb (for jobu = 'U') or m*nb (for jobu = 'F')

If full SVD is needed (jobu = 'U' or 'F') and

if jobv = 'V',

the minimal requirement is lwork≥ max(2*m+n, 6*n+2*n*n)
if jobv = 'J',

the minimal requirement is lwork≥ max(2*m+n, 4*n+n*n, 2*n+n*n+6)
For optimal performance, lwork should be additionally larger than n+m*nb, where nb is the optimal block size for ?ormlq.

For complex flavors:

Length of cwork to confirm proper allocation of workspace. The value of lwork depends on the job:

If only sigma is needed ( jobu.EQ.'N', jobv.EQ.'N' ) and
- no scaled condition estimate is required: lwork≥ 2*n+1. This is the minimal requirement. For optimal performance (blocked code) the optimal value is lwork≥n + (n+1)*nb. Here nb is the optimal block size for ?geqp3 and ?geqrf. In general, optimal lwork is computed as lwork≥ max(n+lwork(?geqp3), n+lwork(?geqrf), lwork(?gesvj)).
- an estimate of the scaled condition number of a is required (joba='E' or 'G'). In this case, the minimal requirement is lwork≥n*n + 2*n. For optimal performance (blocked code) the optimal value is lwork≥ max(n+(n+1)*nb, n*n+3*n) = n² + 2*n. In general, the optimal length lwork is computed as lwork≥ max(n+lwork(?geqp3),n+lwork(?geqrf), lwork(?gesvj), n+n*n+lwork(?pocon)).
If sigma and the right singular vectors are needed (jobv.EQ.'V' or jobu.EQ.'N') and
- no scaled condition estimate is requested (jobe .EQ.'N'), then the minimal requirement is lwork≥ 3*n. For optimal performance, lwork≥ max(n+(n+1)*nb, 2*n+n*nb) = 2*n+n*nb, where nb is the optimal block size for ?geqp3, ?geqrf, ?gelq, ?unmlq.
  
  In general, the optimal length lwork is computed as lwork≥ max(n+lwork(?geqp3), n+lwork(?gesvj), n+lwork(?gelqf), 2*n+lwork(?geqrf), n+lwork(?unmlq)).
- an estimate of the scaled condition number of a is required (joba='E' or 'G'), then the minimal requirement is lwork≥ 3*n. For optimal performance, lwork≥ max(n+(n+1)*nb, 2*n, 2*n+n*nb) = 2*n+n*nb, where nb is the optimal block size for ?geqp3, ?geqrf, ?gelq, ?unmlq
  
  In general, the optimal length lwork is computed as lwork≥ max(n+lwork(?geqp3), n+lwork(?pocon), n+lwork(?gesvj), n+lwork(?gelqf), 2*n+lwork(?geqrf), n+lwork(?unmlq))
If sigma and the left singular vectors are needed and
- no scaled condition estimate is requested (jobe .EQ.'N'), then the minimal requirement is lwork≥ 3*n.
  
  For optimal performance: if jobu.EQ.'U' :: lwork≥ max(3*n, n+(n+1)*nb, 2*n+n*nb) = 2*n+n*nb, where nb is the optimal block size for ?geqp3, ?geqrf, ?unmqr. In general, the optimal length lwork is computed as lwork≥ max(n+lwork(?geqp3),2*n+lwork(?geqrf), n+lwork(?unmqr)).
- an estimate of the scaled condition number of a is required (joba='E' or 'G'), then the minimal requirement is lwork≥ 3*n.
  
  For optimal performance: if jobu.EQ.'U' :: lwork≥ max(3*n, n+(n+1)*nb, 2*n+n*nb) = 2*n+n*nb, where nb is the optimal block size for ?geqp3, ?geqrf, ?unmqr. In general, the optimal length lwork is computed as lwork≥ max(n+lwork(?geqp3),n+lwork(?pocon), 2*n+lwork(?geqrf), n+lwork(?unmqr)).
then the minimal requirement is lwork≥ 3*n. For optimal performance: if jobu.EQ.'U' :: lwork≥ max(3*n, n+(n+1)*nb, 2*n+n*nb), where nb is the optimal block size for ?geqp3, ?geqrf, ?unmqr. In general, the optimal length lwork is computed as lwork≥ max(n+lwork(?geqp3),n+lwork(?pocon), 2*n+lwork(?geqrf), n+lwork(?unmqr)).
If the full SVD is needed: (jobu.EQ.'U' or jobu.EQ.'F') and
- if jobv.EQ.'V' the minimal requirement is lwork≥ 5*n+2*n*n.
- if jobv.EQ.'J' the minimal requirement is lwork≥ 4*n+n*n. In both cases, the allocated cwork can accommodate blocked runs of ?geqp3, ?geqrf, ?gelqf, ?unmqr, ?unmlq.

If the call to ?gejsv is a workspace query (indicated by lwork = -1 or lrwork = -1), then on exit cwork(1) contains the required length of cwork for the job parameters used in the call.

cwork

COMPLEX for cgejsv

DOUBLE COMPLEX for zgejsv

cwork is a workspace array of size max(2, lwork).

If the call to ?gejsv is a workspace query (indicated by lwork = -1 or lrwork = -1), then on exit cwork(1) contains the required length of cwork for the job parameters used in the call.

rwork

REAL for cgejsv

DOUBLE PRECISION for zgejsv

rwork is an array of size at least max(7, lrwork) for real flavors and at least max(7, lwork) for complex flavors.

lrwork

INTEGER. Length of rwork to confirm proper allocation of workspace. lrwork depends on the job:

1. If only singular values are requested i.e. if lsame(jobu,'N') .AND. lsame(jobv,'N') then:

If lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then lrwork = max( 7, 2 * m ).
Otherwise, lrwork = max( 7, n ).

2. If singular values with the right singular vectors are requested i.e. if (lsame(jobv,'V').OR.lsame(jobv,'J')) .AND. .NOT.(lsame(jobu,'U').OR.lsame(jobu,'F')) then:

If lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then lrwork = max( 7, 2 * m ).
Otherwise, lrwork = max( 7, n ).

3. If singular values with the left singular vectors are requested, i.e. if (lsame(jobu,'U').OR.lsame(jobu,'F')) .AND. .NOT.(lsame(jobv,'V').OR.lsame(jobv,'J')) then:

If lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then lrwork = max( 7, 2 * m ).
Otherwise, lrwork = max( 7, n ).

4. If singular values with both the left and the right singular vectors are requested, i.e. if (lsame(jobu,'U').OR.lsame(jobu,'F')) .AND. (lsame(jobv,'V').OR.lsame(jobv,'J')) then:

If lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then lrwork = max( 7, 2 * m ).
Otherwise, lrwork = max( 7, n ).

For complex flavors, if the call to ?gejsv is a workspace query (indicated by lwork = -1 or lrwork = -1), then on exit rwork(1) contains the required length of rwork for the job parameters used in the call.

iwork

INTEGER. Workspace array, of size

For real flavors:

).

max(3, m+3*nmax( 3, 2 * n + m ).

For complex flavors, the size depends on the job but is at least 4:

If only the singular values are requested and lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then the length of iwork is n + m; otherwise the length of iwork is n.
If the singular values and the right singular vectors are requested and lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then the length of iwork is n + m; otherwise the length of iwork is n.
If the singular values and the left singular vectors are requested and lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then the length of iwork is n + m; otherwise the length of iwork is n.
If the singular values and both the left and the right singular vectors are requested and
- if lsame(jobv,'J'), if lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then the length of iwork is n + m; otherwise the length of iwork is n.
- if lsame(jobv,'V') , if lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then the length of iwork is 2*n + m; otherwise the length of iwork is 2*n.
and lsame(jobt,'T') .OR. lsame(joba,'F') .OR. lsame(joba,'G'), then the length of iwork is n + m; otherwise the length of iwork is n.

Output Parameters

sva

On exit:

For work(1)/work(2) = one: the singular values of A. During the computation sva contains Euclidean column norms of the iterated matrices in the array a.

For work(1)≠work(2): the singular values of A are (work(1)/work(2)) * sva(1:n). This factored form is used if sigma_max(A) overflows or if small singular values have been saved from underflow by scaling the input matrix A.

jobr = 'R', some of the singular values may be returned as exact zeros obtained by 'setting to zero' because they are below the numerical rank threshold or are denormalized numbers.

u

On exit:

If jobu = 'U', contains the m-by-n matrix of the left singular vectors.

If jobu = 'F', contains the m-by-m matrix of the left singular vectors, including an orthonormal basis of the orthogonal complement of the range of A.

If jobu = 'W' and jobv = 'V', jobt = 'T', and m = n, then u is used as workspace if the procedure replaces A with A^T (for real flavors) or A^H (for complex flavors). In that case, v is computed in u as left singular vectors of A^T or A^H and copied back to the v array. This 'W' option is just a reminder to the caller that in this case u is reserved as workspace of length n*n.

If jobu = 'N', u is not referenced.

v

On exit:

If jobv = 'V' or 'J', contains the n-by-n matrix of the right singular vectors.

If jobv = 'W' and jobu = 'U', jobt = 'T', and m = n, then v is used as workspace if the procedure replaces A with A^T (for real flavors) or A^H (for complex flavors). In that case, u is computed in v as right singular vectors of A^T or A^H and copied back to the u array. This 'W' option is just a reminder to the caller that in this case v is reserved as workspace of length n*n.

If jobv = 'N', v is not referenced.

work

On exit,

work(1) = scale = work(2)/work(1) is the scaling factor such that scale*sva(1:n) are the computed singular values of A. See the description of sva().

work(2) = see the description of work(1).

work(3) = sconda is an estimate for the condition number of column equilibrated A. If joba = 'E' or 'G', sconda is an estimate of sqrt(||(R**t * R)**(-1)||_1). It is computed using ?pocon. It holds n**(-1/4) * sconda≤ ||R**(-1)||_2 ≤n**(1/4) * sconda, where R is the triangular factor from the QRF of A. However, if R is truncated and the numerical rank is determined to be strictly smaller than n, sconda is returned as -1, indicating that the smallest singular values might be lost.

If full SVD is needed, the following two condition numbers are useful for the analysis of the algorithm. They are provided for a user who is familiar with the details of the method.

work(4) = an estimate of the scaled condition number of the triangular factor in the first QR factorization.

work(5) = an estimate of the scaled condition number of the triangular factor in the second QR factorization.

The following two parameters are computed if jobt = 'T'. They are provided for a user who is familiar with the details of the method.

work(6) = the entropy of A**t*A :: this is the Shannon entropy of diag(A**t*A) / Trace(A**t*A) taken as point in the probability simplex.

work(7) = the entropy of A*A**t.

rwork

On exit,

rwork(1) determines the scaling factor scale = rwork(2) / rwork(1) such that scale*sva(1:n) are the computed singular values of a. (See the description of sva().)

rwork(2) = see the description of rwork(1).

rwork(3) = sconda is an estimate for the condition number of column equilibrated A. If joba = 'E' or 'G', sconda is an estimate of SQRT(||(R^* * R)^(-1)||_1). It is computed using ?pocon. It holds n^(-1/4) * sconda≤ ||R^(-1)||_2 ≤n^(1/4) * sconda where R is the triangular factor from the QRF of A. However, if R is truncated and the numerical rank is determined to be strictly smaller than n, sconda is returned as -1, thus indicating that the smallest singular values might be lost.

If full SVD is needed, the following two condition numbers are useful for the analysis of the algorithm. They are provided for a user who is familiar with the details of the method.

rwork(4) = an estimate of the scaled condition number of the triangular factor in the first QR factorization.

rwork(5) = an estimate of the scaled condition number of the triangular factor in the second QR factorization.

The following two parameters are computed if jobt = 'T'. They are provided for a user who is familiar with the details of the method.

rwork(6) = the entropy of A^* * A :: this is the Shannon entropy of diag(A^* * A) / Trace(A^* * A) taken as point in the probability simplex.

rwork(7) = the entropy of A * A^*. (See the description of rwork(6).)

iwork

INTEGER. On exit,

iwork(1) = the numerical rank determined after the initial QR factorization with pivoting. See the descriptions of joba and jobr.

iwork(2) = the number of the computed nonzero singular value.

iwork(3) = if nonzero, a warning message. If iwork(3)=1, some of the column norms of A were denormalized floats. The requested high accuracy is not warranted by the data.

For complex flavors, iwork(4) = 1 or -1. If iwork(4) = 1, then the procedure used A^H to do the job as specified by the job parameters.

For complex flavors, if the call to ?gejsv is a workspace query (indicated by lwork = -1 or lrwork = -1), then on exit iwork(1) contains the required length of iwork for the job parameters used in the call.

info

INTEGER.

If info = 0, the execution is successful.

If info = -i, the i-th parameter had an illegal value.

If info > 0, the function did not converge in the maximal number of sweeps. The computed values may be inaccurate.

Parent topic: Singular Value Decomposition: LAPACK Driver Routines

Select Your Language

Using Intel.com Search

Quick Links

Recent Searches

Advanced Search

Only search in

Developer Reference for Intel® oneAPI Math Kernel Library for Fortran

?gejsv

See Also