Provided by: liblapack-doc_3.12.0-3build1.1_all bug

NAME

       gesvj - gesvj: SVD, Jacobi, low-level

SYNOPSIS

   Functions
       subroutine cgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, cwork, lwork, rwork,
           lrwork, info)
            CGESVJ
       subroutine dgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, work, lwork, info)
           DGESVJ
       subroutine sgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, work, lwork, info)
           SGESVJ
       subroutine zgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, cwork, lwork, rwork,
           lrwork, info)
            ZGESVJ

Detailed Description

Function Documentation

   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)
        CGESVJ

       Purpose:

            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.

       Parameters
           JOBA

                     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.

           JOBU

                     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.

           JOBV

                     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

           M

                     M is INTEGER
                     The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.

           N

                     N is INTEGER
                     The number of columns of the input matrix A.
                     M >= N >= 0.

           A

                     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).

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,M).

           SVA

                     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.

           MV

                     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.

           V

                     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.

           LDV

                     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) .

           CWORK

                     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.

           LWORK

                     LWORK is INTEGER.
                     Length of CWORK, LWORK >= M+N.

           RWORK

                     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.

           LRWORK

                    LRWORK is INTEGER
                    Length of RWORK, LRWORK >= MAX(6,N).

           INFO

                     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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

            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.

       Contributor:

             ============

             Zlatko Drmac (Zagreb, Croatia)

       References:

            [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.

       Bugs, examples and comments:

             ===========================
             Please report all bugs and send interesting test examples and comments to
             drmac@math.hr. Thank you.

   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)
       DGESVJ

       Purpose:

            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.

       Parameters
           JOBA

                     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.

           JOBU

                     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.

           JOBV

                     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

           M

                     M is INTEGER
                     The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.

           N

                     N is INTEGER
                     The number of columns of the input matrix A.
                     M >= N >= 0.

           A

                     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).

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,M).

           SVA

                     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.

           MV

                     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.

           V

                     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.

           LDV

                     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) .

           WORK

                     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.

           LWORK

                     LWORK is INTEGER
                     length of WORK, WORK >= MAX(6,M+N)

           INFO

                     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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             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.

       Contributors:

             ============

             Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

       References:

            [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.

       Bugs, examples and comments:

             ===========================
             Please report all bugs and send interesting test examples and comments to
             drmac@math.hr. Thank you.

   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)
       SGESVJ

       Purpose:

            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.

       Parameters
           JOBA

                     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.

           JOBU

                     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.

           JOBV

                     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

           M

                     M is INTEGER
                     The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.

           N

                     N is INTEGER
                     The number of columns of the input matrix A.
                     M >= N >= 0.

           A

                     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).

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,M).

           SVA

                     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.

           MV

                     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.

           V

                     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.

           LDV

                     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) .

           WORK

                     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.

           LWORK

                     LWORK is INTEGER
                    length of WORK, WORK >= MAX(6,M+N)

           INFO

                     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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:
           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.

       Contributors:
           Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

       References:
           [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.

       Bugs, Examples and Comments:
           Please report all bugs and send interesting test examples and comments to
           drmac@math.hr. Thank you.

   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)
        ZGESVJ

       Purpose:

            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.

       Parameters
           JOBA

                     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.

           JOBU

                     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.

           JOBV

                     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

           M

                     M is INTEGER
                     The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.

           N

                     N is INTEGER
                     The number of columns of the input matrix A.
                     M >= N >= 0.

           A

                     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).

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,M).

           SVA

                     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.

           MV

                     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.

           V

                     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.

           LDV

                     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) .

           CWORK

                     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.

           LWORK

                     LWORK is INTEGER.
                     Length of CWORK, LWORK >= M+N.

           RWORK

                     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.

           LRWORK

                    LRWORK is INTEGER
                    Length of RWORK, LRWORK >= MAX(6,N).

           INFO

                     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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

            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.

       Contributor:

             ============

             Zlatko Drmac (Zagreb, Croatia)

       References:

            [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.

       Bugs, examples and comments:

             ===========================
             Please report all bugs and send interesting test examples and comments to
             drmac@math.hr. Thank you.

Author

       Generated automatically by Doxygen for LAPACK from the source code.