.TH "gesvj" 3 "Wed Feb 7 2024 11:30:40" "Version 3.12.0" "LAPACK" \" -*- nroff -*- .ad l .nh .SH NAME gesvj \- gesvj: SVD, Jacobi, low-level .SH SYNOPSIS .br .PP .SS "Functions" .in +1c .ti -1c .RI "subroutine \fBcgesvj\fP (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, cwork, lwork, rwork, lrwork, info)" .br .RI "\fB CGESVJ \fP " .ti -1c .RI "subroutine \fBdgesvj\fP (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, work, lwork, info)" .br .RI "\fBDGESVJ\fP " .ti -1c .RI "subroutine \fBsgesvj\fP (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, work, lwork, info)" .br .RI "\fBSGESVJ\fP " .ti -1c .RI "subroutine \fBzgesvj\fP (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, cwork, lwork, rwork, lrwork, info)" .br .RI "\fB ZGESVJ \fP " .in -1c .SH "Detailed Description" .PP .SH "Function Documentation" .PP .SS "subroutine cgesvj (character*1 joba, character*1 jobu, character*1 jobv, integer m, integer n, complex, dimension( lda, * ) a, integer lda, real, dimension( n ) sva, integer mv, complex, dimension( ldv, * ) v, integer ldv, complex, dimension( lwork ) cwork, integer lwork, real, dimension( lrwork ) rwork, integer lrwork, integer info)" .PP \fB CGESVJ \fP .PP \fBPurpose:\fP .RS 4 .PP .nf CGESVJ computes the singular value decomposition (SVD) of a complex M-by-N matrix A, where M >= N\&. The SVD of A is written as [++] [xx] [x0] [xx] A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx] [++] [xx] where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal 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\&. .fi .PP .RE .PP \fBParameters\fP .RS 4 \fIJOBA\fP .PP .nf JOBA is CHARACTER*1 Specifies the structure of A\&. = 'L': The input matrix A is lower triangular; = 'U': The input matrix A is upper triangular; = 'G': The input matrix A is general M-by-N matrix, M >= N\&. .fi .PP .br \fIJOBU\fP .PP .nf JOBU is CHARACTER*1 Specifies whether to compute the left singular vectors (columns of U): = 'U' or 'F': The left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of A\&. See more details in the description of A\&. The default numerical orthogonality threshold is set to approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E')\&. = 'C': Analogous to JOBU='U', except that user can control the level of numerical orthogonality of the computed left singular vectors\&. TOL can be set to TOL = CTOL*EPS, where CTOL is given on input in the array WORK\&. No CTOL smaller than ONE is allowed\&. CTOL greater than 1 / EPS is meaningless\&. The option 'C' can be used if M*EPS is satisfactory orthogonality of the computed left singular vectors, so CTOL=M could save few sweeps of Jacobi rotations\&. See the descriptions of A and WORK(1)\&. = 'N': The matrix U is not computed\&. However, see the description of A\&. .fi .PP .br \fIJOBV\fP .PP .nf JOBV is CHARACTER*1 Specifies whether to compute the right singular vectors, that is, the matrix V: = 'V' or 'J': the matrix V is computed and returned in the array V = 'A': the Jacobi rotations are applied to the MV-by-N array V\&. In other words, the right singular vector matrix V is not computed explicitly; instead it is applied to an MV-by-N matrix initially stored in the first MV rows of V\&. = 'N': the matrix V is not computed and the array V is not referenced .fi .PP .br \fIM\fP .PP .nf M is INTEGER The number of rows of the input matrix A\&. 1/SLAMCH('E') > 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\&. On exit, If JOBU = 'U' \&.OR\&. JOBU = 'C': If INFO = 0 : RANKA orthonormal columns of U are returned in the leading RANKA columns of the array A\&. Here RANKA <= N is the number of computed singular values of A that are above the underflow threshold SLAMCH('S')\&. The singular vectors corresponding to underflowed or zero singular values are not computed\&. The value of RANKA is returned in the array RWORK as RANKA=NINT(RWORK(2))\&. Also see the descriptions of SVA and RWORK\&. The computed columns of U are mutually numerically orthogonal up to approximately TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), see the description of JOBU\&. If INFO > 0, the procedure CGESVJ did not converge in the given number of iterations (sweeps)\&. In that case, the computed columns of U may not be orthogonal up to TOL\&. The output U (stored in A), SIGMA (given by the computed singular values in SVA(1:N)) and V is still a decomposition of the input matrix A in the sense that the residual || A - SCALE * U * SIGMA * V^* ||_2 / ||A||_2 is small\&. If JOBU = 'N': If INFO = 0 : 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 of U is not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately M*EPS\&. Thus, on exit, A contains the columns of U scaled with the corresponding singular values\&. If INFO > 0 : the procedure CGESVJ did not converge in the given number of iterations (sweeps)\&. .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, If INFO = 0 : depending on the value SCALE = RWORK(1), we have: If SCALE = ONE: SVA(1:N) contains the computed singular values of A\&. During the computation SVA contains the Euclidean column norms of the iterated matrices in the array A\&. If SCALE \&.NE\&. ONE: The singular values of A are SCALE*SVA(1:N), and this factored representation is due to the fact that some of the singular values of A might underflow or overflow\&. If INFO > 0 : the procedure CGESVJ did not converge in the given number of iterations (sweeps) and SCALE*SVA(1:N) may not be accurate\&. .fi .PP .br \fIMV\fP .PP .nf MV is INTEGER If JOBV = 'A', then the product of Jacobi rotations in CGESVJ is applied to the first MV rows of V\&. See the description of JOBV\&. .fi .PP .br \fIV\fP .PP .nf V is COMPLEX array, dimension (LDV,N) If JOBV = 'V', then V contains on exit the N-by-N matrix of the right singular vectors; If JOBV = 'A', then V contains the product of the computed right singular vector matrix and the initial matrix in the array V\&. If JOBV = 'N', then V is not referenced\&. .fi .PP .br \fILDV\fP .PP .nf LDV is INTEGER The leading dimension of the array V, LDV >= 1\&. If JOBV = 'V', then LDV >= max(1,N)\&. If JOBV = 'A', then LDV >= max(1,MV) \&. .fi .PP .br \fICWORK\fP .PP .nf CWORK is COMPLEX array, dimension (max(1,LWORK)) Used as workspace\&. If on entry LWORK = -1, then a workspace query is assumed and no computation is done; CWORK(1) is set to the minial (and optimal) length of CWORK\&. .fi .PP .br \fILWORK\fP .PP .nf LWORK is INTEGER\&. Length of CWORK, LWORK >= M+N\&. .fi .PP .br \fIRWORK\fP .PP .nf RWORK is REAL array, dimension (max(6,LRWORK)) On entry, If JOBU = 'C' : RWORK(1) = CTOL, where CTOL defines the threshold for convergence\&. The process stops if all columns of A are mutually orthogonal up to CTOL*EPS, EPS=SLAMCH('E')\&. It is required that CTOL >= ONE, i\&.e\&. it is not allowed to force the routine to obtain orthogonality below EPSILON\&. On exit, RWORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A\&. (See description of SVA()\&.) RWORK(2) = NINT(RWORK(2)) is the number of the computed nonzero singular values\&. RWORK(3) = NINT(RWORK(3)) is the number of the computed singular values that are larger than the underflow threshold\&. RWORK(4) = NINT(RWORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence\&. RWORK(5) = max_{i\&.NE\&.j} |COS(A(:,i),A(:,j))| in the last sweep\&. This is useful information in cases when CGESVJ did not converge, as it can be used to estimate whether the output is still useful and for post festum analysis\&. RWORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep\&. It can be useful for a post festum analysis\&. If on entry LRWORK = -1, then a workspace query is assumed and no computation is done; RWORK(1) is set to the minial (and optimal) length of RWORK\&. .fi .PP .br \fILRWORK\fP .PP .nf LRWORK is INTEGER Length of RWORK, LRWORK >= MAX(6,N)\&. .fi .PP .br \fIINFO\fP .PP .nf INFO is INTEGER = 0: successful exit\&. < 0: if INFO = -i, then the i-th argument had an illegal value > 0: CGESVJ did not converge in the maximal allowed number (NSWEEP=30) of sweeps\&. The output may still be useful\&. See the description of RWORK\&. .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 The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations\&. In the case of underflow of the tangent of the Jacobi angle, a modified Jacobi transformation of Drmac [3] is used\&. Pivot strategy uses column interchanges of de Rijk [1]\&. 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 [2]\&. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(\&.) is the spectral condition number\&. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [4,5], and it is the kernel routine in the SIGMA library [6]\&. Some tuning parameters (marked with [TP]) are available for the implementer\&. 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 digits\&. .fi .PP .RE .PP \fBContributor:\fP .RS 4 .PP .nf ============ Zlatko Drmac (Zagreb, Croatia) .fi .PP .RE .PP \fBReferences:\fP .RS 4 .PP .nf [1] P\&. P\&. M\&. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer\&. SIAM J\&. Sci\&. Stat\&. Comp\&., Vol\&. 10 (1998), pp\&. 359-371\&. [2] J\&. Demmel and K\&. Veselic: Jacobi method is more accurate than QR\&. [3] Z\&. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic\&. SIAM J\&. Sci\&. Comp\&., Vol\&. 18 (1997), pp\&. 1200-1222\&. [4] 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\&. [5] 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\&. [6] Z\&. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations\&. Department of Mathematics, University of Zagreb, 2008, 2015\&. .fi .PP .RE .PP \fBBugs, examples and comments:\fP .RS 4 .PP .nf =========================== Please report all bugs and send interesting test examples and comments to drmac@math\&.hr\&. Thank you\&. .fi .PP .RE .PP .SS "subroutine dgesvj (character*1 joba, character*1 jobu, character*1 jobv, integer m, integer n, double precision, dimension( lda, * ) a, integer lda, double precision, dimension( n ) sva, integer mv, double precision, dimension( ldv, * ) v, integer ldv, double precision, dimension( lwork ) work, integer lwork, integer info)" .PP \fBDGESVJ\fP .PP \fBPurpose:\fP .RS 4 .PP .nf DGESVJ computes the singular value decomposition (SVD) of a real M-by-N matrix A, where M >= N\&. The SVD of A is written as [++] [xx] [x0] [xx] A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx] [++] [xx] where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal matrix, and V is an N-by-N orthogonal 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\&. DGESVJ can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines, see below under Further Details\&. .fi .PP .RE .PP \fBParameters\fP .RS 4 \fIJOBA\fP .PP .nf JOBA is CHARACTER*1 Specifies the structure of A\&. = 'L': The input matrix A is lower triangular; = 'U': The input matrix A is upper triangular; = 'G': The input matrix A is general M-by-N matrix, M >= N\&. .fi .PP .br \fIJOBU\fP .PP .nf JOBU is CHARACTER*1 Specifies whether to compute the left singular vectors (columns of U): = 'U': The left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of A\&. See more details in the description of A\&. The default numerical orthogonality threshold is set to approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E')\&. = 'C': Analogous to JOBU='U', except that user can control the level of numerical orthogonality of the computed left singular vectors\&. TOL can be set to TOL = CTOL*EPS, where CTOL is given on input in the array WORK\&. No CTOL smaller than ONE is allowed\&. CTOL greater than 1 / EPS is meaningless\&. The option 'C' can be used if M*EPS is satisfactory orthogonality of the computed left singular vectors, so CTOL=M could save few sweeps of Jacobi rotations\&. See the descriptions of A and WORK(1)\&. = 'N': The matrix U is not computed\&. However, see the description of A\&. .fi .PP .br \fIJOBV\fP .PP .nf JOBV is CHARACTER*1 Specifies whether to compute the right singular vectors, that is, the matrix V: = 'V': the matrix V is computed and returned in the array V = 'A': the Jacobi rotations are applied to the MV-by-N array V\&. In other words, the right singular vector matrix V is not computed explicitly, instead it is applied to an MV-by-N matrix initially stored in the first MV rows of V\&. = 'N': the matrix V is not computed and the array V is not referenced .fi .PP .br \fIM\fP .PP .nf M is INTEGER The number of rows of the input matrix A\&. 1/DLAMCH('E') > 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 DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A\&. On exit : If JOBU = 'U' \&.OR\&. JOBU = 'C' : If INFO = 0 : RANKA orthonormal columns of U are returned in the leading RANKA columns of the array A\&. Here RANKA <= N is the number of computed singular values of A that are above the underflow threshold DLAMCH('S')\&. The singular vectors corresponding to underflowed or zero singular values are not computed\&. The value of RANKA is returned in the array WORK as RANKA=NINT(WORK(2))\&. Also see the descriptions of SVA and WORK\&. The computed columns of U are mutually numerically orthogonal up to approximately TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), see the description of JOBU\&. If INFO > 0 : the procedure DGESVJ did not converge in the given number of iterations (sweeps)\&. In that case, the computed columns of U may not be orthogonal up to TOL\&. The output U (stored in A), SIGMA (given by the computed singular values in SVA(1:N)) and V is still a decomposition of the input matrix A in the sense that the residual ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small\&. If JOBU = 'N' : If INFO = 0 : 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 of U is not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately M*EPS\&. Thus, on exit, A contains the columns of U scaled with the corresponding singular values\&. If INFO > 0 : the procedure DGESVJ did not converge in the given number of iterations (sweeps)\&. .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 DOUBLE PRECISION array, dimension (N) On exit : If INFO = 0 : depending on the value SCALE = WORK(1), we have: If SCALE = ONE : SVA(1:N) contains the computed singular values of A\&. During the computation SVA contains the Euclidean column norms of the iterated matrices in the array A\&. If SCALE \&.NE\&. ONE : The singular values of A are SCALE*SVA(1:N), and this factored representation is due to the fact that some of the singular values of A might underflow or overflow\&. If INFO > 0 : the procedure DGESVJ did not converge in the given number of iterations (sweeps) and SCALE*SVA(1:N) may not be accurate\&. .fi .PP .br \fIMV\fP .PP .nf MV is INTEGER If JOBV = 'A', then the product of Jacobi rotations in DGESVJ is applied to the first MV rows of V\&. See the description of JOBV\&. .fi .PP .br \fIV\fP .PP .nf V is DOUBLE PRECISION array, dimension (LDV,N) If JOBV = 'V', then V contains on exit the N-by-N matrix of the right singular vectors; If JOBV = 'A', then V contains the product of the computed right singular vector matrix and the initial matrix in the array V\&. If JOBV = 'N', then V is not referenced\&. .fi .PP .br \fILDV\fP .PP .nf LDV is INTEGER The leading dimension of the array V, LDV >= 1\&. If JOBV = 'V', then LDV >= max(1,N)\&. If JOBV = 'A', then LDV >= max(1,MV) \&. .fi .PP .br \fIWORK\fP .PP .nf WORK is DOUBLE PRECISION array, dimension (LWORK) On entry : If JOBU = 'C' : WORK(1) = CTOL, where CTOL defines the threshold for convergence\&. The process stops if all columns of A are mutually orthogonal up to CTOL*EPS, EPS=DLAMCH('E')\&. It is required that CTOL >= ONE, i\&.e\&. it is not allowed to force the routine to obtain orthogonality below EPS\&. On exit : WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A\&. (See description of SVA()\&.) WORK(2) = NINT(WORK(2)) is the number of the computed nonzero singular values\&. WORK(3) = NINT(WORK(3)) is the number of the computed singular values that are larger than the underflow threshold\&. WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence\&. WORK(5) = max_{i\&.NE\&.j} |COS(A(:,i),A(:,j))| in the last sweep\&. This is useful information in cases when DGESVJ did not converge, as it can be used to estimate whether the output is still useful and for post festum analysis\&. WORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep\&. It can be useful for a post festum analysis\&. .fi .PP .br \fILWORK\fP .PP .nf LWORK is INTEGER length of WORK, WORK >= MAX(6,M+N) .fi .PP .br \fIINFO\fP .PP .nf INFO is INTEGER = 0: successful exit\&. < 0: if INFO = -i, then the i-th argument had an illegal value > 0: DGESVJ did not converge in the maximal allowed number (30) of sweeps\&. The output may still be useful\&. See the description of WORK\&. .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 The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations\&. The rotations are implemented as fast scaled rotations of Anda and Park [1]\&. In the case of underflow of the Jacobi angle, a modified Jacobi transformation of Drmac [4] is used\&. Pivot strategy uses column interchanges of de Rijk [2]\&. 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 [3]\&. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(\&.) is the spectral condition number\&. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [5,6], and it is the kernel routine in the SIGMA library [7]\&. Some tuning parameters (marked with [TP]) are available for the implementer\&. 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 digits\&. .fi .PP .RE .PP \fBContributors:\fP .RS 4 .PP .nf ============ Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) .fi .PP .RE .PP \fBReferences:\fP .RS 4 .PP .nf [1] A\&. A\&. Anda and H\&. Park: Fast plane rotations with dynamic scaling\&. SIAM J\&. matrix Anal\&. Appl\&., Vol\&. 15 (1994), pp\&. 162-174\&. [2] P\&. P\&. M\&. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer\&. SIAM J\&. Sci\&. Stat\&. Comp\&., Vol\&. 10 (1998), pp\&. 359-371\&. [3] J\&. Demmel and K\&. Veselic: Jacobi method is more accurate than QR\&. [4] Z\&. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic\&. SIAM J\&. Sci\&. Comp\&., Vol\&. 18 (1997), pp\&. 1200-1222\&. [5] 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\&. [6] 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\&. [7] Z\&. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations\&. Department of Mathematics, University of Zagreb, 2008\&. .fi .PP .RE .PP \fBBugs, examples and comments:\fP .RS 4 .PP .nf =========================== Please report all bugs and send interesting test examples and comments to drmac@math\&.hr\&. Thank you\&. .fi .PP .RE .PP .SS "subroutine sgesvj (character*1 joba, character*1 jobu, character*1 jobv, integer m, integer n, real, dimension( lda, * ) a, integer lda, real, dimension( n ) sva, integer mv, real, dimension( ldv, * ) v, integer ldv, real, dimension( lwork ) work, integer lwork, integer info)" .PP \fBSGESVJ\fP .PP \fBPurpose:\fP .RS 4 .PP .nf SGESVJ computes the singular value decomposition (SVD) of a real M-by-N matrix A, where M >= N\&. The SVD of A is written as [++] [xx] [x0] [xx] A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx] [++] [xx] where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal matrix, and V is an N-by-N orthogonal 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\&. SGESVJ can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines, see below under Further Details\&. .fi .PP .RE .PP \fBParameters\fP .RS 4 \fIJOBA\fP .PP .nf JOBA is CHARACTER*1 Specifies the structure of A\&. = 'L': The input matrix A is lower triangular; = 'U': The input matrix A is upper triangular; = 'G': The input matrix A is general M-by-N matrix, M >= N\&. .fi .PP .br \fIJOBU\fP .PP .nf JOBU is CHARACTER*1 Specifies whether to compute the left singular vectors (columns of U): = 'U': The left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of A\&. See more details in the description of A\&. The default numerical orthogonality threshold is set to approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E')\&. = 'C': Analogous to JOBU='U', except that user can control the level of numerical orthogonality of the computed left singular vectors\&. TOL can be set to TOL = CTOL*EPS, where CTOL is given on input in the array WORK\&. No CTOL smaller than ONE is allowed\&. CTOL greater than 1 / EPS is meaningless\&. The option 'C' can be used if M*EPS is satisfactory orthogonality of the computed left singular vectors, so CTOL=M could save few sweeps of Jacobi rotations\&. See the descriptions of A and WORK(1)\&. = 'N': The matrix U is not computed\&. However, see the description of A\&. .fi .PP .br \fIJOBV\fP .PP .nf JOBV is CHARACTER*1 Specifies whether to compute the right singular vectors, that is, the matrix V: = 'V': the matrix V is computed and returned in the array V = 'A': the Jacobi rotations are applied to the MV-by-N array V\&. In other words, the right singular vector matrix V is not computed explicitly; instead it is applied to an MV-by-N matrix initially stored in the first MV rows of V\&. = 'N': the matrix V is not computed and the array V is not referenced .fi .PP .br \fIM\fP .PP .nf M is INTEGER The number of rows of the input matrix A\&. 1/SLAMCH('E') > 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 REAL array, dimension (LDA,N) On entry, the M-by-N matrix A\&. On exit, If JOBU = 'U' \&.OR\&. JOBU = 'C': If INFO = 0: RANKA orthonormal columns of U are returned in the leading RANKA columns of the array A\&. Here RANKA <= N is the number of computed singular values of A that are above the underflow threshold SLAMCH('S')\&. The singular vectors corresponding to underflowed or zero singular values are not computed\&. The value of RANKA is returned in the array WORK as RANKA=NINT(WORK(2))\&. Also see the descriptions of SVA and WORK\&. The computed columns of U are mutually numerically orthogonal up to approximately TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), see the description of JOBU\&. If INFO > 0, the procedure SGESVJ did not converge in the given number of iterations (sweeps)\&. In that case, the computed columns of U may not be orthogonal up to TOL\&. The output U (stored in A), SIGMA (given by the computed singular values in SVA(1:N)) and V is still a decomposition of the input matrix A in the sense that the residual ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small\&. If JOBU = 'N': If INFO = 0: 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 of U is not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately M*EPS\&. Thus, on exit, A contains the columns of U scaled with the corresponding singular values\&. If INFO > 0: the procedure SGESVJ did not converge in the given number of iterations (sweeps)\&. .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, If INFO = 0 : depending on the value SCALE = WORK(1), we have: If SCALE = ONE: SVA(1:N) contains the computed singular values of A\&. During the computation SVA contains the Euclidean column norms of the iterated matrices in the array A\&. If SCALE \&.NE\&. ONE: The singular values of A are SCALE*SVA(1:N), and this factored representation is due to the fact that some of the singular values of A might underflow or overflow\&. If INFO > 0 : the procedure SGESVJ did not converge in the given number of iterations (sweeps) and SCALE*SVA(1:N) may not be accurate\&. .fi .PP .br \fIMV\fP .PP .nf MV is INTEGER If JOBV = 'A', then the product of Jacobi rotations in SGESVJ is applied to the first MV rows of V\&. See the description of JOBV\&. .fi .PP .br \fIV\fP .PP .nf V is REAL array, dimension (LDV,N) If JOBV = 'V', then V contains on exit the N-by-N matrix of the right singular vectors; If JOBV = 'A', then V contains the product of the computed right singular vector matrix and the initial matrix in the array V\&. If JOBV = 'N', then V is not referenced\&. .fi .PP .br \fILDV\fP .PP .nf LDV is INTEGER The leading dimension of the array V, LDV >= 1\&. If JOBV = 'V', then LDV >= max(1,N)\&. If JOBV = 'A', then LDV >= max(1,MV) \&. .fi .PP .br \fIWORK\fP .PP .nf WORK is REAL array, dimension (LWORK) On entry, If JOBU = 'C' : WORK(1) = CTOL, where CTOL defines the threshold for convergence\&. The process stops if all columns of A are mutually orthogonal up to CTOL*EPS, EPS=SLAMCH('E')\&. It is required that CTOL >= ONE, i\&.e\&. it is not allowed to force the routine to obtain orthogonality below EPSILON\&. On exit, WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular vcalues of A\&. (See description of SVA()\&.) WORK(2) = NINT(WORK(2)) is the number of the computed nonzero singular values\&. WORK(3) = NINT(WORK(3)) is the number of the computed singular values that are larger than the underflow threshold\&. WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence\&. WORK(5) = max_{i\&.NE\&.j} |COS(A(:,i),A(:,j))| in the last sweep\&. This is useful information in cases when SGESVJ did not converge, as it can be used to estimate whether the output is still useful and for post festum analysis\&. WORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep\&. It can be useful for a post festum analysis\&. .fi .PP .br \fILWORK\fP .PP .nf LWORK is INTEGER length of WORK, WORK >= MAX(6,M+N) .fi .PP .br \fIINFO\fP .PP .nf INFO is INTEGER = 0: successful exit\&. < 0: if INFO = -i, then the i-th argument had an illegal value > 0: SGESVJ did not converge in the maximal allowed number (30) of sweeps\&. The output may still be useful\&. See the description of WORK\&. .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 The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations\&. The rotations are implemented as fast scaled rotations of Anda and Park [1]\&. In the case of underflow of the Jacobi angle, a modified Jacobi transformation of Drmac [4] is used\&. Pivot strategy uses column interchanges of de Rijk [2]\&. 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 [3]\&. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(\&.) is the spectral condition number\&. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [5,6], and it is the kernel routine in the SIGMA library [7]\&. Some tuning parameters (marked with [TP]) are available for the implementer\&. .br 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 digits\&. .RE .PP \fBContributors:\fP .RS 4 Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) .RE .PP \fBReferences:\fP .RS 4 [1] A\&. A\&. Anda and H\&. Park: Fast plane rotations with dynamic scaling\&. .br SIAM J\&. matrix Anal\&. Appl\&., Vol\&. 15 (1994), pp\&. 162-174\&. .br .br [2] P\&. P\&. M\&. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer\&. .br SIAM J\&. Sci\&. Stat\&. Comp\&., Vol\&. 10 (1998), pp\&. 359-371\&. .br .br [3] J\&. Demmel and K\&. Veselic: Jacobi method is more accurate than QR\&. .br [4] Z\&. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic\&. .br SIAM J\&. Sci\&. Comp\&., Vol\&. 18 (1997), pp\&. 1200-1222\&. .br .br [5] Z\&. Drmac and K\&. Veselic: New fast and accurate Jacobi SVD algorithm I\&. .br SIAM J\&. Matrix Anal\&. Appl\&. Vol\&. 35, No\&. 2 (2008), pp\&. 1322-1342\&. .br LAPACK Working note 169\&. .br .br [6] Z\&. Drmac and K\&. Veselic: New fast and accurate Jacobi SVD algorithm II\&. .br SIAM J\&. Matrix Anal\&. Appl\&. Vol\&. 35, No\&. 2 (2008), pp\&. 1343-1362\&. .br LAPACK Working note 170\&. .br .br [7] Z\&. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations\&. .br Department of Mathematics, University of Zagreb, 2008\&. .RE .PP \fBBugs, Examples and Comments:\fP .RS 4 Please report all bugs and send interesting test examples and comments to drmac@math.hr\&. Thank you\&. .RE .PP .SS "subroutine zgesvj (character*1 joba, character*1 jobu, character*1 jobv, integer m, integer n, complex*16, dimension( lda, * ) a, integer lda, double precision, dimension( n ) sva, integer mv, complex*16, dimension( ldv, * ) v, integer ldv, complex*16, dimension( lwork ) cwork, integer lwork, double precision, dimension( lrwork ) rwork, integer lrwork, integer info)" .PP \fB ZGESVJ \fP .PP \fBPurpose:\fP .RS 4 .PP .nf ZGESVJ computes the singular value decomposition (SVD) of a complex M-by-N matrix A, where M >= N\&. The SVD of A is written as [++] [xx] [x0] [xx] A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx] [++] [xx] where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal 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\&. .fi .PP .RE .PP \fBParameters\fP .RS 4 \fIJOBA\fP .PP .nf JOBA is CHARACTER*1 Specifies the structure of A\&. = 'L': The input matrix A is lower triangular; = 'U': The input matrix A is upper triangular; = 'G': The input matrix A is general M-by-N matrix, M >= N\&. .fi .PP .br \fIJOBU\fP .PP .nf JOBU is CHARACTER*1 Specifies whether to compute the left singular vectors (columns of U): = 'U' or 'F': The left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of A\&. See more details in the description of A\&. The default numerical orthogonality threshold is set to approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=DLAMCH('E')\&. = 'C': Analogous to JOBU='U', except that user can control the level of numerical orthogonality of the computed left singular vectors\&. TOL can be set to TOL = CTOL*EPS, where CTOL is given on input in the array WORK\&. No CTOL smaller than ONE is allowed\&. CTOL greater than 1 / EPS is meaningless\&. The option 'C' can be used if M*EPS is satisfactory orthogonality of the computed left singular vectors, so CTOL=M could save few sweeps of Jacobi rotations\&. See the descriptions of A and WORK(1)\&. = 'N': The matrix U is not computed\&. However, see the description of A\&. .fi .PP .br \fIJOBV\fP .PP .nf JOBV is CHARACTER*1 Specifies whether to compute the right singular vectors, that is, the matrix V: = 'V' or 'J': the matrix V is computed and returned in the array V = 'A': the Jacobi rotations are applied to the MV-by-N array V\&. In other words, the right singular vector matrix V is not computed explicitly; instead it is applied to an MV-by-N matrix initially stored in the first MV rows of V\&. = 'N': the matrix V is not computed and the array V is not referenced .fi .PP .br \fIM\fP .PP .nf M is INTEGER The number of rows of the input matrix A\&. 1/DLAMCH('E') > 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*16 array, dimension (LDA,N) On entry, the M-by-N matrix A\&. On exit, If JOBU = 'U' \&.OR\&. JOBU = 'C': If INFO = 0 : RANKA orthonormal columns of U are returned in the leading RANKA columns of the array A\&. Here RANKA <= N is the number of computed singular values of A that are above the underflow threshold DLAMCH('S')\&. The singular vectors corresponding to underflowed or zero singular values are not computed\&. The value of RANKA is returned in the array RWORK as RANKA=NINT(RWORK(2))\&. Also see the descriptions of SVA and RWORK\&. The computed columns of U are mutually numerically orthogonal up to approximately TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), see the description of JOBU\&. If INFO > 0, the procedure ZGESVJ did not converge in the given number of iterations (sweeps)\&. In that case, the computed columns of U may not be orthogonal up to TOL\&. The output U (stored in A), SIGMA (given by the computed singular values in SVA(1:N)) and V is still a decomposition of the input matrix A in the sense that the residual || A - SCALE * U * SIGMA * V^* ||_2 / ||A||_2 is small\&. If JOBU = 'N': If INFO = 0 : 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 of U is not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately M*EPS\&. Thus, on exit, A contains the columns of U scaled with the corresponding singular values\&. If INFO > 0: the procedure ZGESVJ did not converge in the given number of iterations (sweeps)\&. .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 DOUBLE PRECISION array, dimension (N) On exit, If INFO = 0 : depending on the value SCALE = RWORK(1), we have: If SCALE = ONE: SVA(1:N) contains the computed singular values of A\&. During the computation SVA contains the Euclidean column norms of the iterated matrices in the array A\&. If SCALE \&.NE\&. ONE: The singular values of A are SCALE*SVA(1:N), and this factored representation is due to the fact that some of the singular values of A might underflow or overflow\&. If INFO > 0: the procedure ZGESVJ did not converge in the given number of iterations (sweeps) and SCALE*SVA(1:N) may not be accurate\&. .fi .PP .br \fIMV\fP .PP .nf MV is INTEGER If JOBV = 'A', then the product of Jacobi rotations in ZGESVJ is applied to the first MV rows of V\&. See the description of JOBV\&. .fi .PP .br \fIV\fP .PP .nf V is COMPLEX*16 array, dimension (LDV,N) If JOBV = 'V', then V contains on exit the N-by-N matrix of the right singular vectors; If JOBV = 'A', then V contains the product of the computed right singular vector matrix and the initial matrix in the array V\&. If JOBV = 'N', then V is not referenced\&. .fi .PP .br \fILDV\fP .PP .nf LDV is INTEGER The leading dimension of the array V, LDV >= 1\&. If JOBV = 'V', then LDV >= max(1,N)\&. If JOBV = 'A', then LDV >= max(1,MV) \&. .fi .PP .br \fICWORK\fP .PP .nf CWORK is COMPLEX*16 array, dimension (max(1,LWORK)) Used as workspace\&. If on entry LWORK = -1, then a workspace query is assumed and no computation is done; CWORK(1) is set to the minial (and optimal) length of CWORK\&. .fi .PP .br \fILWORK\fP .PP .nf LWORK is INTEGER\&. Length of CWORK, LWORK >= M+N\&. .fi .PP .br \fIRWORK\fP .PP .nf RWORK is DOUBLE PRECISION array, dimension (max(6,LRWORK)) On entry, If JOBU = 'C' : RWORK(1) = CTOL, where CTOL defines the threshold for convergence\&. The process stops if all columns of A are mutually orthogonal up to CTOL*EPS, EPS=DLAMCH('E')\&. It is required that CTOL >= ONE, i\&.e\&. it is not allowed to force the routine to obtain orthogonality below EPSILON\&. On exit, RWORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A\&. (See description of SVA()\&.) RWORK(2) = NINT(RWORK(2)) is the number of the computed nonzero singular values\&. RWORK(3) = NINT(RWORK(3)) is the number of the computed singular values that are larger than the underflow threshold\&. RWORK(4) = NINT(RWORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence\&. RWORK(5) = max_{i\&.NE\&.j} |COS(A(:,i),A(:,j))| in the last sweep\&. This is useful information in cases when ZGESVJ did not converge, as it can be used to estimate whether the output is still useful and for post festum analysis\&. RWORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep\&. It can be useful for a post festum analysis\&. If on entry LRWORK = -1, then a workspace query is assumed and no computation is done; RWORK(1) is set to the minial (and optimal) length of RWORK\&. .fi .PP .br \fILRWORK\fP .PP .nf LRWORK is INTEGER Length of RWORK, LRWORK >= MAX(6,N)\&. .fi .PP .br \fIINFO\fP .PP .nf INFO is INTEGER = 0: successful exit\&. < 0: if INFO = -i, then the i-th argument had an illegal value > 0: ZGESVJ did not converge in the maximal allowed number (NSWEEP=30) of sweeps\&. The output may still be useful\&. See the description of RWORK\&. .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 The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations\&. In the case of underflow of the tangent of the Jacobi angle, a modified Jacobi transformation of Drmac [3] is used\&. Pivot strategy uses column interchanges of de Rijk [1]\&. 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 [2]\&. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(\&.) is the spectral condition number\&. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [4,5], and it is the kernel routine in the SIGMA library [6]\&. Some tuning parameters (marked with [TP]) are available for the implementer\&. 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 digits\&. .fi .PP .RE .PP \fBContributor:\fP .RS 4 .PP .nf ============ Zlatko Drmac (Zagreb, Croatia) .fi .PP .RE .PP \fBReferences:\fP .RS 4 .PP .nf [1] P\&. P\&. M\&. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer\&. SIAM J\&. Sci\&. Stat\&. Comp\&., Vol\&. 10 (1998), pp\&. 359-371\&. [2] J\&. Demmel and K\&. Veselic: Jacobi method is more accurate than QR\&. [3] Z\&. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic\&. SIAM J\&. Sci\&. Comp\&., Vol\&. 18 (1997), pp\&. 1200-1222\&. [4] 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\&. [5] 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\&. [6] Z\&. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations\&. Department of Mathematics, University of Zagreb, 2008, 2015\&. .fi .PP .RE .PP \fBBugs, examples and comments:\fP .RS 4 .PP .nf =========================== Please report all bugs and send interesting test examples and comments to drmac@math\&.hr\&. Thank you\&. .fi .PP .RE .PP .SH "Author" .PP Generated automatically by Doxygen for LAPACK from the source code\&.