Computes the singular value decomposition of a real matrix using a preconditioned Jacobi SVD method.
FORTRAN 77:
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)
C:
lapack_int LAPACKE_<?>gejsv( int matrix_order, char joba, char jobu, char jobv, char jobr, char jobt, char jobp, lapack_int m, lapack_int n, const <datatype>* a, lapack_int lda, <datatype>* sva, <datatype>* u, lapack_int ldu, <datatype>* v, lapack_int ldv, <datatype>* stat, lapack_int* istat );
The routine computes the singular value decomposition (SVD) of a real m-by-n matrix A, where m ≥ n.
The SVD is written as
A = U*Σ*VT,
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 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.
The data types are given for the Fortran interface, except for istat. A <datatype> placeholder, if present, is used for the C interface data types in the C interface section above. See C Interface Conventions for the C interface principal conventions and type definitions.
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 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'.
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.
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.
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.
CHARACTER*1. Must be 'T' or 'N'.
If the matrix is square, the procedure may determine to use a transposed A if A**t seems to be better with respect to convergence. If the matrix is not square, jobt is ignored. This is subject to changes in the future.
The decision is based on two values of entropy over the adjoint orbit of A**t * A. See the descriptions of work(6) and work(7).
If jobt = 'T', the function perfomrs transpositon 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, 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 matrices 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.
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.
INTEGER. The number of rows of the input matrix A; m ≥ 0.
INTEGER. The number of columns in the input matrix A; n ≥ 0.
REAL for sgejsv
DOUBLE PRECISION for dgejsv.
Array a(lda,*) is an array containing the m-by-n matrix A.
The second dimension of a must be at least max(1, n).
work is a workspace array, its dimension max(1, lwork).
sva is a workspace array, its dimension is n.
u is a workspace array, its dimension is (ldu,*); the second dimension of u must be at least max(1, n).
v is a workspace array, its dimension is (ldv,*); the second dimension of u must be at least max(1, n).
INTEGER. The leading dimension of the array a. Must be at least max(1, m).
INTEGER. The leading dimension of the array u; ldu ≥ 1.
jobu = 'U' or 'F' or 'W', ldu ≥ m.
INTEGER. The leading dimension of the array v; ldv ≥ 1.
jobv = 'V' or 'J' or 'W', ldv ≥ n.
INTEGER.
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.
INTEGER. Workspace array, DIMENSION max(3, m+3*n).
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.
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. In that case, v is computed in u as left singular vectors of A**t 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.
On exit:
If jobv = 'V' or 'J', contains the n-by-n matrix of the right singular vectors.
If jobv = 'W' and jobv = 'U', jobt = 'T', and m = n, then v is used as workspace if the procedure replaces A with A**t. In that case, u is computed in v as right singular vectors of A**t 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.
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 provied 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.
INTEGER. On exit,
iwork(1)/istat[0] = the numerical rank determined after the initial QR factorization with pivoting. See the descriptions of joba and jobr.
iwork(2)/istat[1] = the number of the computed nonzero singular value.
iwork(3)/istat[2] = if nonzero, a warning message. If iwork(3)/istat[2]=1, some of the column norms of A were denormalized floats. The requested high accuracy is not warranted by the data.
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.