.TH "SRC/cgejsv.f" 3 "Version 3.12.0" "LAPACK" \" -*- nroff -*- .ad l .nh .SH NAME SRC/cgejsv.f .SH SYNOPSIS .br .PP .SS "Functions/Subroutines" .in +1c .ti -1c .RI "subroutine \fBcgejsv\fP (joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, cwork, lwork, rwork, lrwork, iwork, info)" .br .RI "\fBCGEJSV\fP " .in -1c .SH "Function/Subroutine Documentation" .PP .SS "subroutine cgejsv (character*1 joba, character*1 jobu, character*1 jobv, character*1 jobr, character*1 jobt, character*1 jobp, integer m, integer n, complex, dimension( lda, * ) a, integer lda, real, dimension( n ) sva, complex, dimension( ldu, * ) u, integer ldu, complex, dimension( ldv, * ) v, integer ldv, complex, dimension( lwork ) cwork, integer lwork, real, dimension( lrwork ) rwork, integer lrwork, integer, dimension( * ) iwork, integer info)" .PP \fBCGEJSV\fP .PP \fBPurpose:\fP .RS 4 .PP .nf CGEJSV computes the singular value decomposition (SVD) of a complex M-by-N matrix [A], where M >= N\&. The SVD of [A] is written as [A] = [U] * [SIGMA] * [V]^*, where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and [V] is an N-by-N unitary matrix\&. The diagonal elements of [SIGMA] are the singular values of [A]\&. The columns of [U] and [V] are the left and the 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 [SIGMA] is computed and stored in the array SVA\&. .fi .PP .RE .PP \fBParameters\fP .RS 4 \fIJOBA\fP .PP .nf JOBA is CHARACTER*1 Specifies the level of accuracy: = 'C': This option works well (high relative accuracy) 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: = 'E': Computation as with 'C' with an additional estimate of the condition number of B\&. It provides a realistic error bound\&. = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings D1, D2, and well-conditioned matrix C, this option gives higher accuracy than the 'C' option\&. If the structure of the input matrix is not known, and relative accuracy is desirable, then this option is advisable\&. The input matrix A is preprocessed with QR factorization with FULL (row and column) pivoting\&. = '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, then using this condition number gives too pessimistic error bound\&. = 'A': Small singular values are not well determined by the data and are considered as noisy; the matrix is treated as numerically rank deficient\&. The error in the computed singular values is bounded by f(m,n)*epsilon*||A||\&. The computed SVD A = U * S * V^* restores A up to f(m,n)*epsilon*||A||\&. This gives the procedure the licence to discard (set to zero) all singular values below N*epsilon*||A||\&. = 'R': Similar as in 'A'\&. Rank revealing property of the initial QR factorization is used do 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'\&. .fi .PP .br \fIJOBU\fP .PP .nf JOBU is CHARACTER*1 Specifies whether to compute the columns of U: = 'U': N columns of U are returned in the array U\&. = 'F': full set of M left sing\&. vectors is returned in the array U\&. = 'W': U may be used as workspace of length M*N\&. See the description of U\&. = 'N': U is not computed\&. .fi .PP .br \fIJOBV\fP .PP .nf JOBV is CHARACTER*1 Specifies whether to compute the matrix V: = 'V': N columns of V are returned in the array V; Jacobi rotations are not explicitly accumulated\&. = 'J': N columns of V are returned in the array V, but they are computed as the product of Jacobi rotations, if JOBT = 'N'\&. = 'W': V may be used as workspace of length N*N\&. See the description of V\&. = 'N': V is not computed\&. .fi .PP .br \fIJOBR\fP .PP .nf JOBR is CHARACTER*1 Specifies the RANGE for the singular values\&. Issues the licence to set to zero small positive singular values if they are outside specified range\&. If A \&.NE\&. 0 is scaled so that the largest singular value of c*A is around SQRT(BIG), BIG=SLAMCH('O'), then JOBR issues the licence to kill columns of A whose norm in c*A is less than SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN, where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E')\&. = 'N': Do not kill small columns of c*A\&. This option assumes that BLAS and QR factorizations and triangular solvers are implemented to work in that range\&. If the condition of A is greater than BIG, use CGESVJ\&. = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)] (roughly, as described above)\&. This option is recommended\&. =========================== For computing the singular values in the FULL range [SFMIN,BIG] use CGESVJ\&. .fi .PP .br \fIJOBT\fP .PP .nf JOBT is CHARACTER*1 If the matrix is square then the procedure may determine to use transposed A if A^* 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^* * A\&. See the descriptions of RWORK(6) and RWORK(7)\&. = 'T': transpose if entropy test indicates possibly faster convergence of Jacobi process if A^* is taken as input\&. If A is replaced with A^*, then the row pivoting is included automatically\&. = 'N': do not speculate\&. The option 'T' 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\&. In general, this option is considered experimental, and 'N'; should be preferred\&. This is subject to changes in the future\&. .fi .PP .br \fIJOBP\fP .PP .nf JOBP is CHARACTER*1 Issues the licence to introduce structured perturbations to drown denormalized numbers\&. This licence should be active if the denormals are poorly implemented, causing slow computation, especially in cases of fast convergence (!)\&. For details see [1,2]\&. For the sake of simplicity, this perturbations are included only when the full SVD or only the singular values are requested\&. The implementer/user can easily add the perturbation for the cases of computing one set of singular vectors\&. = 'P': introduce perturbation = 'N': do not perturb .fi .PP .br \fIM\fP .PP .nf M is INTEGER The number of rows of the input matrix A\&. M >= 0\&. .fi .PP .br \fIN\fP .PP .nf N is INTEGER The number of columns of the input matrix A\&. M >= N >= 0\&. .fi .PP .br \fIA\fP .PP .nf A is COMPLEX array, dimension (LDA,N) On entry, the M-by-N matrix A\&. .fi .PP .br \fILDA\fP .PP .nf LDA is INTEGER The leading dimension of the array A\&. LDA >= max(1,M)\&. .fi .PP .br \fISVA\fP .PP .nf SVA is REAL array, dimension (N) On exit, - For RWORK(1)/RWORK(2) = ONE: The singular values of A\&. During the computation SVA contains Euclidean column norms of the iterated matrices in the array A\&. - For RWORK(1) \&.NE\&. RWORK(2): The singular values of A are (RWORK(1)/RWORK(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\&. - If JOBR='R' then some of the singular values may be returned as exact zeros obtained by 'set to zero' because they are below the numerical rank threshold or are denormalized numbers\&. .fi .PP .br \fIU\fP .PP .nf U is COMPLEX array, dimension ( LDU, N ) or ( LDU, M ) If JOBU = 'U', then U contains on exit the M-by-N matrix of the left singular vectors\&. If JOBU = 'F', then U contains on exit the M-by-M matrix of the left singular vectors, including an ONB of the orthogonal complement of the Range(A)\&. If JOBU = 'W' \&.AND\&. (JOBV = 'V' \&.AND\&. JOBT = 'T' \&.AND\&. M = N), then U is used as workspace if the procedure replaces A with A^*\&. In that case, [V] is computed in U as left singular vectors of A^* and then 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, unless JOBT='T'\&. .fi .PP .br \fILDU\fP .PP .nf LDU is INTEGER The leading dimension of the array U, LDU >= 1\&. IF JOBU = 'U' or 'F' or 'W', then LDU >= M\&. .fi .PP .br \fIV\fP .PP .nf V is COMPLEX array, dimension ( LDV, N ) If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of the right singular vectors; If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N), then V is used as workspace if the procedure replaces A with A^*\&. In that case, [U] is computed in V as right singular vectors of A^* and then 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, unless JOBT='T'\&. .fi .PP .br \fILDV\fP .PP .nf LDV is INTEGER The leading dimension of the array V, LDV >= 1\&. If JOBV = 'V' or 'J' or 'W', then LDV >= N\&. .fi .PP .br \fICWORK\fP .PP .nf CWORK is COMPLEX array, dimension (MAX(2,LWORK)) If the call to CGEJSV 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\&. .fi .PP .br \fILWORK\fP .PP .nf LWORK is INTEGER Length of CWORK to confirm proper allocation of workspace\&. LWORK depends on the job: 1\&. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and 1\&.1 \&.\&. no scaled condition estimate required (JOBA\&.NE\&.'E'\&.AND\&.JOBA\&.NE\&.'G'): 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 CGEQP3 and CGEQRF\&. In general, optimal LWORK is computed as LWORK >= max(N+LWORK(CGEQP3),N+LWORK(CGEQRF), LWORK(CGESVJ))\&. 1\&.2\&. \&.\&. an estimate of the scaled condition number of A is required (JOBA='E', or 'G')\&. In this case, LWORK 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+2*N)=N**2+2*N\&. In general, the optimal length LWORK is computed as LWORK >= max(N+LWORK(CGEQP3),N+LWORK(CGEQRF), LWORK(CGESVJ), N*N+LWORK(CPOCON))\&. 2\&. If SIGMA and the right singular vectors are needed (JOBV = 'V'), (JOBU = 'N') 2\&.1 \&.\&. no scaled condition estimate requested (JOBE = 'N'): -> 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 CGEQP3, CGEQRF, CGELQF, CUNMLQ\&. In general, the optimal length LWORK is computed as LWORK >= max(N+LWORK(CGEQP3), N+LWORK(CGESVJ), N+LWORK(CGELQF), 2*N+LWORK(CGEQRF), N+LWORK(CUNMLQ))\&. 2\&.2 \&.\&. an estimate of the scaled condition number of A is required (JOBA='E', or 'G')\&. -> 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 CGEQP3, CGEQRF, CGELQF, CUNMLQ\&. In general, the optimal length LWORK is computed as LWORK >= max(N+LWORK(CGEQP3), LWORK(CPOCON), N+LWORK(CGESVJ), N+LWORK(CGELQF), 2*N+LWORK(CGEQRF), N+LWORK(CUNMLQ))\&. 3\&. If SIGMA and the left singular vectors are needed 3\&.1 \&.\&. no scaled condition estimate requested (JOBE = 'N'): -> the minimal requirement is LWORK >= 3*N\&. -> For optimal performance: if JOBU = '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 CGEQP3, CGEQRF, CUNMQR\&. In general, the optimal length LWORK is computed as LWORK >= max(N+LWORK(CGEQP3), 2*N+LWORK(CGEQRF), N+LWORK(CUNMQR))\&. 3\&.2 \&.\&. an estimate of the scaled condition number of A is required (JOBA='E', or 'G')\&. -> the minimal requirement is LWORK >= 3*N\&. -> For optimal performance: if JOBU = '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 CGEQP3, CGEQRF, CUNMQR\&. In general, the optimal length LWORK is computed as LWORK >= max(N+LWORK(CGEQP3),N+LWORK(CPOCON), 2*N+LWORK(CGEQRF), N+LWORK(CUNMQR))\&. 4\&. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and 4\&.1\&. if JOBV = 'V' the minimal requirement is LWORK >= 5*N+2*N*N\&. 4\&.2\&. if JOBV = 'J' the minimal requirement is LWORK >= 4*N+N*N\&. In both cases, the allocated CWORK can accommodate blocked runs of CGEQP3, CGEQRF, CGELQF, CUNMQR, CUNMLQ\&. If the call to CGEJSV is a workspace query (indicated by LWORK=-1 or LRWORK=-1), then on exit CWORK(1) contains the optimal and CWORK(2) contains the minimal length of CWORK for the job parameters used in the call\&. .fi .PP .br \fIRWORK\fP .PP .nf RWORK is REAL array, dimension (MAX(7,LRWORK)) 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 CPOCON\&. 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 developer/implementer 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 developer/implementer 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)\&.) If the call to CGEJSV 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\&. .fi .PP .br \fILRWORK\fP .PP .nf LRWORK is INTEGER Length of RWORK to confirm proper allocation of workspace\&. LRWORK depends on the job: 1\&. If only the singular values are requested i\&.e\&. if LSAME(JOBU,'N') \&.AND\&. LSAME(JOBV,'N') then: 1\&.1\&. If LSAME(JOBT,'T') \&.OR\&. LSAME(JOBA,'F') \&.OR\&. LSAME(JOBA,'G'), then: LRWORK = max( 7, 2 * M )\&. 1\&.2\&. 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: 2\&.1\&. If LSAME(JOBT,'T') \&.OR\&. LSAME(JOBA,'F') \&.OR\&. LSAME(JOBA,'G'), then LRWORK = max( 7, 2 * M )\&. 2\&.2\&. 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: 3\&.1\&. If LSAME(JOBT,'T') \&.OR\&. LSAME(JOBA,'F') \&.OR\&. LSAME(JOBA,'G'), then LRWORK = max( 7, 2 * M )\&. 3\&.2\&. 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: 4\&.1\&. If LSAME(JOBT,'T') \&.OR\&. LSAME(JOBA,'F') \&.OR\&. LSAME(JOBA,'G'), then LRWORK = max( 7, 2 * M )\&. 4\&.2\&. Otherwise, LRWORK = max( 7, N )\&. If, on entry, LRWORK = -1 or LWORK=-1, a workspace query is assumed and the length of RWORK is returned in RWORK(1)\&. .fi .PP .br \fIIWORK\fP .PP .nf IWORK is INTEGER array, of dimension at least 4, that further depends on the job: 1\&. If only the singular values are requested then: 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\&. 2\&. If the singular values and the right singular vectors are requested then: 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\&. 3\&. If the singular values and the left singular vectors are requested then: 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\&. 4\&. If the singular values with both the left and the right singular vectors are requested, then: 4\&.1\&. If LSAME(JOBV,'J') the length of IWORK is determined as follows: 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\&. 4\&.2\&. If LSAME(JOBV,'V') the length of IWORK is determined as follows: 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\&. 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 values IWORK(3) = if nonzero, a warning message: If IWORK(3) = 1 then some of the column norms of A were denormalized floats\&. The requested high accuracy is not warranted by the data\&. IWORK(4) = 1 or -1\&. If IWORK(4) = 1, then the procedure used A^* to do the job as specified by the JOB parameters\&. If the call to CGEJSV is a workspace query (indicated by LWORK = -1 and LRWORK = -1), then on exit IWORK(1) contains the required length of IWORK for the job parameters used in the call\&. .fi .PP .br \fIINFO\fP .PP .nf INFO is INTEGER < 0: if INFO = -i, then the i-th argument had an illegal value\&. = 0: successful exit; > 0: CGEJSV did not converge in the maximal allowed number of sweeps\&. The computed values may be inaccurate\&. .fi .PP .RE .PP \fBAuthor\fP .RS 4 Univ\&. of Tennessee .PP Univ\&. of California Berkeley .PP Univ\&. of Colorado Denver .PP NAG Ltd\&. .RE .PP \fBFurther Details:\fP .RS 4 .PP .nf CGEJSV implements a preconditioned Jacobi SVD algorithm\&. It uses CGEQP3, CGEQRF, and CGELQF 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 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\&. Further, 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 [1], [2]\&. The computational range for the singular values can be the full range ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS & LAPACK routines called by CGEJSV are implemented to work in that range\&. If that is not the case, then 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 (CGEJSV) is best used in this restricted range, meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are returned as zeros\&. See JOBR for details on this\&. Further, this implementation is somewhat slower than the one described in [1,2] due to replacement of some non-LAPACK components, and because the choice of some tuning parameters in the iterative part (CGESVJ) is left to the implementer on a particular machine\&. The rank revealing QR factorization (in this code: CGEQP3) should be implemented as in [3]\&. We have a new version of CGEQP3 under development that is more robust than the current one in LAPACK, with a cleaner cut in rank deficient cases\&. It will be available in the SIGMA library [4]\&. If M is much larger than N, it is obvious that the initial QRF with column pivoting can be preprocessed by the QRF without pivoting\&. That well known trick is not used in CGEJSV 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 terms of accuracy have prevailed\&. The implementer/user can incorporate this extra QRF step easily\&. The implementer can also improve data movement (matrix transpose, matrix copy, matrix transposed copy) - this implementation of CGEJSV uses only the simplest, naive data movement\&. .fi .PP .RE .PP \fBContributor:\fP .RS 4 Zlatko Drmac (Zagreb, Croatia) .RE .PP \fBReferences:\fP .RS 4 .PP .nf [1] Z\&. Drmac and K\&. Veselic: New fast and accurate Jacobi SVD algorithm I\&. SIAM J\&. Matrix Anal\&. Appl\&. Vol\&. 35, No\&. 2 (2008), pp\&. 1322-1342\&. LAPACK Working note 169\&. [2] Z\&. Drmac and K\&. Veselic: New fast and accurate Jacobi SVD algorithm II\&. SIAM J\&. Matrix Anal\&. Appl\&. Vol\&. 35, No\&. 2 (2008), pp\&. 1343-1362\&. LAPACK Working note 170\&. [3] Z\&. Drmac and Z\&. Bujanovic: On the failure of rank-revealing QR factorization software - a case study\&. ACM Trans\&. Math\&. Softw\&. Vol\&. 35, No 2 (2008), pp\&. 1-28\&. LAPACK Working note 176\&. [4] Z\&. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations\&. Department of Mathematics, University of Zagreb, 2008, 2016\&. .fi .PP .RE .PP \fBBugs, examples and comments:\fP .RS 4 Please report all bugs and send interesting examples and/or comments to drmac@math.hr\&. Thank you\&. .RE .PP .PP Definition at line \fB565\fP of file \fBcgejsv\&.f\fP\&. .SH "Author" .PP Generated automatically by Doxygen for LAPACK from the source code\&.