rpl/lapack/lapack/zgejsv.f - diff

Return to zgejsv.f CVS log

Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/zgejsv.f between versions 1.4 and 1.9

-version 1.4, 2017/06/17 10:54:09
+version 1.9, 2023/08/07 08:39:17
  Line 1
  *> \brief \b ZGEJSV
  *
  *  =========== DOCUMENTATION ===========
  *
  * Online html documentation available at
  *            http://www.netlib.org/lapack/explore-html/
  *
  *> \htmlonly
  *> Download ZGEJSV + dependencies
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgejsv.f">
  *> [TGZ]</a>
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgejsv.f">
  *> [ZIP]</a>
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgejsv.f">
  *> [TXT]</a>
  *> \endhtmlonly
  *
  *  Definition:
  *  ===========
  *
  *     SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
  *                         M, N, A, LDA, SVA, U, LDU, V, LDV,
  *                         CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
  *
  *     .. Scalar Arguments ..
  *     IMPLICIT    NONE
  *     INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
  *     ..
  *     .. Array Arguments ..
  *     COMPLEX*16     A( LDA, * ),  U( LDU, * ), V( LDV, * ), CWORK( LWORK )
  *     DOUBLE PRECISION   SVA( N ), RWORK( LRWORK )
  *     INTEGER     IWORK( * )
  *     CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
  *       ..
  *
  *
  *> \par Purpose:
  *  =============
  *>
  *> \verbatim
  *>
  *> ZGEJSV computes the singular value decomposition (SVD) of a complex M-by-N
  *> matrix [A], where M >= N. The SVD of [A] is written as
  *>
  *>              [A] = [U] * [SIGMA] * [V]^*,
  *>
  *> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
  *> diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and
  *> [V] is an N-by-N unitary matrix. The diagonal elements of [SIGMA] are
  *> the singular values of [A]. The columns of [U] and [V] are the left and
  *> the right singular vectors of [A], respectively. The matrices [U] and [V]
  *> are computed and stored in the arrays U and V, respectively. The diagonal
  *> of [SIGMA] is computed and stored in the array SVA.
  *> \endverbatim
  *>
  *>  Arguments:
  *>  ==========
  *>
  *> \param[in] JOBA
  *> \verbatim
  *>          JOBA is CHARACTER*1
  *>         Specifies the level of accuracy:
  *>       = 'C': This option works well (high relative accuracy) if A = B * D,
  *>              with well-conditioned B and arbitrary diagonal matrix D.
  *>              The accuracy cannot be spoiled by COLUMN scaling. The
  *>              accuracy of the computed output depends on the condition of
  *>              B, and the procedure aims at the best theoretical accuracy.
  *>              The relative error max_{i=1:N}|d sigma_i| / sigma_i is
  *>              bounded by f(M,N)*epsilon* cond(B), independent of D.
  *>              The input matrix is preprocessed with the QRF with column
  *>              pivoting. This initial preprocessing and preconditioning by
  *>              a rank revealing QR factorization is common for all values of
  *>              JOBA. Additional actions are specified as follows:
  *>       = 'E': Computation as with 'C' with an additional estimate of the
  *>              condition number of B. It provides a realistic error bound.
  *>       = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
  *>              D1, D2, and well-conditioned matrix C, this option gives
  *>              higher accuracy than the 'C' option. If the structure of the
  *>              input matrix is not known, and relative accuracy is
  *>              desirable, then this option is advisable. The input matrix A
  *>              is preprocessed with QR factorization with FULL (row and
  *>              column) pivoting.
- *>       = 'G'  Computation as with 'F' with an additional estimate of the
+ *>       = 'G': Computation as with 'F' with an additional estimate of the
  *>              condition number of B, where A=B*D. If A has heavily weighted
  *>              rows, then using this condition number gives too pessimistic
  *>              error bound.
  *>       = 'A': Small singular values are not well determined by the data
  *>              and are considered as noisy; the matrix is treated as
- *>              numerically rank defficient. The error in the computed
+ *>              numerically rank deficient. The error in the computed
  *>              singular values is bounded by f(m,n)*epsilon*||A||.
  *>              The computed SVD A = U * S * V^* restores A up to
  *>              f(m,n)*epsilon*||A||.
  *>              This gives the procedure the licence to discard (set to zero)
  *>              all singular values below N*epsilon*||A||.
  *>       = 'R': Similar as in 'A'. Rank revealing property of the initial
  *>              QR factorization is used do reveal (using triangular factor)
  *>              a gap sigma_{r+1} < epsilon * sigma_r in which case the
  *>              numerical RANK is declared to be r. The SVD is computed with
  *>              absolute error bounds, but more accurately than with 'A'.
  *> \endverbatim
  *>
  *> \param[in] JOBU
  *> \verbatim
  *>          JOBU is CHARACTER*1
  *>         Specifies whether to compute the columns of U:
  *>       = 'U': N columns of U are returned in the array U.
  *>       = 'F': full set of M left sing. vectors is returned in the array U.
  *>       = 'W': U may be used as workspace of length M*N. See the description
  *>              of U.
  *>       = 'N': U is not computed.
  *> \endverbatim
  *>
  *> \param[in] JOBV
  *> \verbatim
  *>          JOBV is CHARACTER*1
  *>         Specifies whether to compute the matrix V:
  *>       = 'V': N columns of V are returned in the array V; Jacobi rotations
  *>              are not explicitly accumulated.
  *>       = 'J': N columns of V are returned in the array V, but they are
- *>              computed as the product of Jacobi rotations, if JOBT .EQ. 'N'.
+ *>              computed as the product of Jacobi rotations, if JOBT = 'N'.
  *>       = 'W': V may be used as workspace of length N*N. See the description
  *>              of V.
  *>       = 'N': V is not computed.
  *> \endverbatim
  *>
  *> \param[in] JOBR
  *> \verbatim
  *>          JOBR is CHARACTER*1
  *>         Specifies the RANGE for the singular values. Issues the licence to
  *>         set to zero small positive singular values if they are outside
  *>         specified range. If A .NE. 0 is scaled so that the largest singular
  *>         value of c*A is around SQRT(BIG), BIG=DLAMCH('O'), then JOBR issues
  *>         the licence to kill columns of A whose norm in c*A is less than
- *>         SQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
+ *>         SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
  *>         where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('E').
  *>       = 'N': Do not kill small columns of c*A. This option assumes that
  *>              BLAS and QR factorizations and triangular solvers are
  *>              implemented to work in that range. If the condition of A
  *>              is greater than BIG, use ZGESVJ.
  *>       = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)]
  *>              (roughly, as described above). This option is recommended.
  *>                                             ===========================
  *>         For computing the singular values in the FULL range [SFMIN,BIG]
  *>         use ZGESVJ.
  *> \endverbatim
  *>
  *> \param[in] JOBT
  *> \verbatim
  *>          JOBT is CHARACTER*1
  *>         If the matrix is square then the procedure may determine to use
  *>         transposed A if A^* seems to be better with respect to convergence.
  *>         If the matrix is not square, JOBT is ignored.
  *>         The decision is based on two values of entropy over the adjoint
  *>         orbit of A^* * A. See the descriptions of WORK(6) and WORK(7).
  *>       = 'T': transpose if entropy test indicates possibly faster
  *>         convergence of Jacobi process if A^* is taken as input. If A is
  *>         replaced with A^*, then the row pivoting is included automatically.
  *>       = 'N': do not speculate.
  *>         The option 'T' can be used to compute only the singular values, or
  *>         the full SVD (U, SIGMA and V). For only one set of singular vectors
  *>         (U or V), the caller should provide both U and V, as one of the
  *>         matrices is used as workspace if the matrix A is transposed.
  *>         The implementer can easily remove this constraint and make the
  *>         code more complicated. See the descriptions of U and V.
  *>         In general, this option is considered experimental, and 'N'; should
  *>         be preferred. This is subject to changes in the future.
  *> \endverbatim
  *>
  *> \param[in] JOBP
  *> \verbatim
  *>          JOBP is CHARACTER*1
  *>         Issues the licence to introduce structured perturbations to drown
  *>         denormalized numbers. This licence should be active if the
  *>         denormals are poorly implemented, causing slow computation,
  *>         especially in cases of fast convergence (!). For details see [1,2].
  *>         For the sake of simplicity, this perturbations are included only
  *>         when the full SVD or only the singular values are requested. The
  *>         implementer/user can easily add the perturbation for the cases of
  *>         computing one set of singular vectors.
  *>       = 'P': introduce perturbation
  *>       = 'N': do not perturb
  *> \endverbatim
  *>
  *> \param[in] M
  *> \verbatim
  *>          M is INTEGER
  *>         The number of rows of the input matrix A.  M >= 0.
  *> \endverbatim
  *>
  *> \param[in] N
  *> \verbatim
  *>          N is INTEGER
  *>         The number of columns of the input matrix A. M >= N >= 0.
  *> \endverbatim
  *>
  *> \param[in,out] A
  *> \verbatim
  *>          A is COMPLEX*16 array, dimension (LDA,N)
  *>          On entry, the M-by-N matrix A.
  *> \endverbatim
  *>
  *> \param[in] LDA
  *> \verbatim
  *>          LDA is INTEGER
  *>          The leading dimension of the array A.  LDA >= max(1,M).
  *> \endverbatim
  *>
  *> \param[out] SVA
  *> \verbatim
  *>          SVA is DOUBLE PRECISION array, dimension (N)
  *>          On exit,
  *>          - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
  *>            computation SVA contains Euclidean column norms of the
  *>            iterated matrices in the array A.
  *>          - For WORK(1) .NE. WORK(2): The singular values of A are
  *>            (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
  *>            sigma_max(A) overflows or if small singular values have been
  *>            saved from underflow by scaling the input matrix A.
  *>          - If JOBR='R' then some of the singular values may be returned
  *>            as exact zeros obtained by "set to zero" because they are
  *>            below the numerical rank threshold or are denormalized numbers.
  *> \endverbatim
  *>
  *> \param[out] U
  *> \verbatim
  *>          U is COMPLEX*16 array, dimension ( LDU, N )
  *>          If JOBU = 'U', then U contains on exit the M-by-N matrix of
  *>                         the left singular vectors.
  *>          If JOBU = 'F', then U contains on exit the M-by-M matrix of
  *>                         the left singular vectors, including an ONB
  *>                         of the orthogonal complement of the Range(A).
- *>          If JOBU = 'W'  .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
+ *>          If JOBU = 'W'  .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
  *>                         then U is used as workspace if the procedure
  *>                         replaces A with A^*. In that case, [V] is computed
  *>                         in U as left singular vectors of A^* and then
  *>                         copied back to the V array. This 'W' option is just
  *>                         a reminder to the caller that in this case U is
  *>                         reserved as workspace of length N*N.
  *>          If JOBU = 'N'  U is not referenced, unless JOBT='T'.
  *> \endverbatim
  *>
  *> \param[in] LDU
  *> \verbatim
  *>          LDU is INTEGER
  *>          The leading dimension of the array U,  LDU >= 1.
  *>          IF  JOBU = 'U' or 'F' or 'W',  then LDU >= M.
  *> \endverbatim
  *>
  *> \param[out] V
  *> \verbatim
  *>          V is COMPLEX*16 array, dimension ( LDV, N )
  *>          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
  *>                         the right singular vectors;
- *>          If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
+ *>          If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
  *>                         then V is used as workspace if the pprocedure
  *>                         replaces A with A^*. In that case, [U] is computed
  *>                         in V as right singular vectors of A^* and then
  *>                         copied back to the U array. This 'W' option is just
  *>                         a reminder to the caller that in this case V is
  *>                         reserved as workspace of length N*N.
  *>          If JOBV = 'N'  V is not referenced, unless JOBT='T'.
  *> \endverbatim
  *>
  *> \param[in] LDV
  *> \verbatim
  *>          LDV is INTEGER
  *>          The leading dimension of the array V,  LDV >= 1.
  *>          If JOBV = 'V' or 'J' or 'W', then LDV >= N.
  *> \endverbatim
  *>
  *> \param[out] CWORK
  *> \verbatim
  *>          CWORK is COMPLEX*16 array, dimension (MAX(2,LWORK))
  *>          If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
  *>          LRWORK=-1), then on exit CWORK(1) contains the required length of
  *>          CWORK for the job parameters used in the call.
  *> \endverbatim
  *>
  *> \param[in] LWORK
  *> \verbatim
  *>          LWORK is INTEGER
  *>          Length of CWORK to confirm proper allocation of workspace.
  *>          LWORK depends on the job:
  *>
- *>          1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
+ *>          1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and
  *>            1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'):
  *>               LWORK >= 2*N+1. This is the minimal requirement.
  *>               ->> For optimal performance (blocked code) the optimal value
  *>               is LWORK >= N + (N+1)*NB. Here NB is the optimal
  *>               block size for ZGEQP3 and ZGEQRF.
  *>               In general, optimal LWORK is computed as
  *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ)).
  *>            1.2. .. an estimate of the scaled condition number of A is
  *>               required (JOBA='E', or 'G'). In this case, LWORK the minimal
  *>               requirement is LWORK >= N*N + 2*N.
  *>               ->> For optimal performance (blocked code) the optimal value
  *>               is LWORK >= max(N+(N+1)*NB, N*N+2*N)=N**2+2*N.
  *>               In general, the optimal length LWORK is computed as
  *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ),
  *>                            N*N+LWORK(ZPOCON)).
- *>          2. If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
+ *>          2. If SIGMA and the right singular vectors are needed (JOBV = 'V'),
- *>             (JOBU.EQ.'N')
+ *>             (JOBU = 'N')
- *>            2.1   .. no scaled condition estimate requested (JOBE.EQ.'N'):
+ *>            2.1   .. no scaled condition estimate requested (JOBE = 'N'):
  *>            -> the minimal requirement is LWORK >= 3*N.
  *>            -> For optimal performance,
  *>               LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
  *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,
  *>               ZUNMLQ. In general, the optimal length LWORK is computed as
  *>               LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZGESVJ),
  *>                       N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).
  *>            2.2 .. an estimate of the scaled condition number of A is
  *>               required (JOBA='E', or 'G').
  *>            -> the minimal requirement is LWORK >= 3*N.
  *>            -> For optimal performance,
  *>               LWORK >= max(N+(N+1)*NB, 2*N,2*N+N*NB)=2*N+N*NB,
  *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,
  *>               ZUNMLQ. In general, the optimal length LWORK is computed as
  *>               LWORK >= max(N+LWORK(ZGEQP3), LWORK(ZPOCON), N+LWORK(ZGESVJ),
  *>                       N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).
  *>          3. If SIGMA and the left singular vectors are needed
- *>            3.1  .. no scaled condition estimate requested (JOBE.EQ.'N'):
+ *>            3.1  .. no scaled condition estimate requested (JOBE = 'N'):
  *>            -> the minimal requirement is LWORK >= 3*N.
  *>            -> For optimal performance:
- *>               if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
+ *>               if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
  *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
  *>               In general, the optimal length LWORK is computed as
  *>               LWORK >= max(N+LWORK(ZGEQP3), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).
  *>            3.2  .. an estimate of the scaled condition number of A is
  *>               required (JOBA='E', or 'G').
  *>            -> the minimal requirement is LWORK >= 3*N.
  *>            -> For optimal performance:
- *>               if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
+ *>               if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
  *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
  *>               In general, the optimal length LWORK is computed as
  *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON),
  *>                        2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).
- *>          4. If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and
+ *>          4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
- *>            4.1. if JOBV.EQ.'V'
+ *>            4.1. if JOBV = 'V'
  *>               the minimal requirement is LWORK >= 5*N+2*N*N.
- *>            4.2. if JOBV.EQ.'J' the minimal requirement is
+ *>            4.2. if JOBV = 'J' the minimal requirement is
  *>               LWORK >= 4*N+N*N.
- *>            In both cases, the allocated CWORK can accomodate blocked runs
+ *>            In both cases, the allocated CWORK can accommodate blocked runs
  *>            of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, ZUNMLQ.
  *>
  *>          If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
  *>          LRWORK=-1), then on exit CWORK(1) contains the optimal and CWORK(2) contains the
  *>          minimal length of CWORK for the job parameters used in the call.
  *> \endverbatim
  *>
  *> \param[out] RWORK
  *> \verbatim
  *>          RWORK is DOUBLE PRECISION array, dimension (MAX(7,LWORK))
  *>          On exit,
  *>          RWORK(1) = Determines the scaling factor SCALE = RWORK(2) / RWORK(1)
  *>                    such that SCALE*SVA(1:N) are the computed singular values
  *>                    of A. (See the description of SVA().)
  *>          RWORK(2) = See the description of RWORK(1).
  *>          RWORK(3) = SCONDA is an estimate for the condition number of
- *>                    column equilibrated A. (If JOBA .EQ. 'E' or 'G')
+ *>                    column equilibrated A. (If JOBA = 'E' or 'G')
  *>                    SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
- *>                    It is computed using SPOCON. It holds
+ *>                    It is computed using ZPOCON. It holds
  *>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
  *>                    where R is the triangular factor from the QRF of A.
  *>                    However, if R is truncated and the numerical rank is
  *>                    determined to be strictly smaller than N, SCONDA is
  *>                    returned as -1, thus indicating that the smallest
  *>                    singular values might be lost.
  *>
  *>          If full SVD is needed, the following two condition numbers are
- *>          useful for the analysis of the algorithm. They are provied for
+ *>          useful for the analysis of the algorithm. They are provided for
  *>          a developer/implementer who is familiar with the details of
  *>          the method.
  *>
  *>          RWORK(4) = an estimate of the scaled condition number of the
  *>                    triangular factor in the first QR factorization.
  *>          RWORK(5) = an estimate of the scaled condition number of the
  *>                    triangular factor in the second QR factorization.
- *>          The following two parameters are computed if JOBT .EQ. 'T'.
+ *>          The following two parameters are computed if JOBT = 'T'.
  *>          They are provided for a developer/implementer who is familiar
  *>          with the details of the method.
  *>          RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy
  *>                    of diag(A^* * A) / Trace(A^* * A) taken as point in the
  *>                    probability simplex.
  *>          RWORK(7) = the entropy of A * A^*. (See the description of RWORK(6).)
  *>          If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
  *>          LRWORK=-1), then on exit RWORK(1) contains the required length of
  *>          RWORK for the job parameters used in the call.
  *> \endverbatim
  *>
  *> \param[in] LRWORK
  *> \verbatim
  *>          LRWORK is INTEGER
  *>          Length of RWORK to confirm proper allocation of workspace.
  *>          LRWORK depends on the job:
  *>
  *>       1. If only the singular values are requested i.e. if
  *>          LSAME(JOBU,'N') .AND. LSAME(JOBV,'N')
  *>          then:
  *>          1.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
  *>               then: LRWORK = max( 7, 2 * M ).
  *>          1.2. Otherwise, LRWORK  = max( 7,  N ).
  *>       2. If singular values with the right singular vectors are requested
  *>          i.e. if
  *>          (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) .AND.
  *>          .NOT.(LSAME(JOBU,'U').OR.LSAME(JOBU,'F'))
  *>          then:
  *>          2.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
  *>          then LRWORK = max( 7, 2 * M ).
  *>          2.2. Otherwise, LRWORK  = max( 7,  N ).
  *>       3. If singular values with the left singular vectors are requested, i.e. if
  *>          (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND.
  *>          .NOT.(LSAME(JOBV,'V').OR.LSAME(JOBV,'J'))
  *>          then:
  *>          3.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
  *>          then LRWORK = max( 7, 2 * M ).
  *>          3.2. Otherwise, LRWORK  = max( 7,  N ).
  *>       4. If singular values with both the left and the right singular vectors
  *>          are requested, i.e. if
  *>          (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND.
  *>          (LSAME(JOBV,'V').OR.LSAME(JOBV,'J'))
  *>          then:
  *>          4.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
  *>          then LRWORK = max( 7, 2 * M ).
  *>          4.2. Otherwise, LRWORK  = max( 7, N ).
  *>
- *>          If, on entry, LRWORK = -1 ot LWORK=-1, a workspace query is assumed and
+ *>          If, on entry, LRWORK = -1 or LWORK=-1, a workspace query is assumed and
  *>          the length of RWORK is returned in RWORK(1).
  *> \endverbatim
  *>
  *> \param[out] IWORK
  *> \verbatim
  *>          IWORK is INTEGER array, of dimension at least 4, that further depends
  *>          on the job:
  *>
  *>          1. If only the singular values are requested then:
  *>             If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
  *>             then the length of IWORK is N+M; otherwise the length of IWORK is N.
  *>          2. If the singular values and the right singular vectors are requested then:
  *>             If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
  *>             then the length of IWORK is N+M; otherwise the length of IWORK is N.
  *>          3. If the singular values and the left singular vectors are requested then:
  *>             If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
  *>             then the length of IWORK is N+M; otherwise the length of IWORK is N.
  *>          4. If the singular values with both the left and the right singular vectors
  *>             are requested, then:
  *>             4.1. If LSAME(JOBV,'J') the length of IWORK is determined as follows:
  *>                  If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
  *>                  then the length of IWORK is N+M; otherwise the length of IWORK is N.
  *>             4.2. If LSAME(JOBV,'V') the length of IWORK is determined as follows:
  *>                  If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
  *>                  then the length of IWORK is 2*N+M; otherwise the length of IWORK is 2*N.
  *>
  *>          On exit,
  *>          IWORK(1) = the numerical rank determined after the initial
  *>                     QR factorization with pivoting. See the descriptions
  *>                     of JOBA and JOBR.
  *>          IWORK(2) = the number of the computed nonzero singular values
  *>          IWORK(3) = if nonzero, a warning message:
- *>                     If IWORK(3).EQ.1 then some of the column norms of A
+ *>                     If IWORK(3) = 1 then some of the column norms of A
  *>                     were denormalized floats. The requested high accuracy
  *>                     is not warranted by the data.
- *>          IWORK(4) = 1 or -1. If IWORK(4) .EQ. 1, then the procedure used A^* to
+ *>          IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to
  *>                     do the job as specified by the JOB parameters.
- *>          If the call to ZGEJSV is a workspace query (indicated by LWORK .EQ. -1 or
+ *>          If the call to ZGEJSV is a workspace query (indicated by LWORK = -1 or
- *>          LRWORK .EQ. -1), then on exit IWORK(1) contains the required length of
+ *>          LRWORK = -1), then on exit IWORK(1) contains the required length of
  *>          IWORK for the job parameters used in the call.
  *> \endverbatim
  *>
  *> \param[out] INFO
  *> \verbatim
  *>          INFO is INTEGER
- *>           < 0  : if INFO = -i, then the i-th argument had an illegal value.
+ *>           < 0:  if INFO = -i, then the i-th argument had an illegal value.
- *>           = 0 :  successful exit;
+ *>           = 0:  successful exit;
- *>           > 0 :  ZGEJSV  did not converge in the maximal allowed number
+ *>           > 0:  ZGEJSV  did not converge in the maximal allowed number
- *>                  of sweeps. The computed values may be inaccurate.
+ *>                 of sweeps. The computed values may be inaccurate.
  *> \endverbatim
  *
  *  Authors:
  *  ========
  *
  *> \author Univ. of Tennessee
  *> \author Univ. of California Berkeley
  *> \author Univ. of Colorado Denver
  *> \author NAG Ltd.
  *
- *> \date June 2016
+ *> \ingroup complex16GEsing
  *
- *> \ingroup complex16GEsing
+ *> \par Further Details:
- *
+ *  =====================
- *> \par Further Details:
+ *>
- *  =====================
+ *> \verbatim
  *>
- *> \verbatim
+ *>  ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3,
- *>
+ *>  ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an
- *>  ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3,
+ *>  additional row pivoting can be used as a preprocessor, which in some
- *>  ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an
+ *>  cases results in much higher accuracy. An example is matrix A with the
- *>  additional row pivoting can be used as a preprocessor, which in some
+ *>  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
- *>  cases results in much higher accuracy. An example is matrix A with the
+ *>  diagonal matrices and C is well-conditioned matrix. In that case, complete
- *>  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
+ *>  pivoting in the first QR factorizations provides accuracy dependent on the
- *>  diagonal matrices and C is well-conditioned matrix. In that case, complete
+ *>  condition number of C, and independent of D1, D2. Such higher accuracy is
- *>  pivoting in the first QR factorizations provides accuracy dependent on the
+ *>  not completely understood theoretically, but it works well in practice.
- *>  condition number of C, and independent of D1, D2. Such higher accuracy is
+ *>  Further, if A can be written as A = B*D, with well-conditioned B and some
- *>  not completely understood theoretically, but it works well in practice.
+ *>  diagonal D, then the high accuracy is guaranteed, both theoretically and
- *>  Further, if A can be written as A = B*D, with well-conditioned B and some
+ *>  in software, independent of D. For more details see [1], [2].
- *>  diagonal D, then the high accuracy is guaranteed, both theoretically and
+ *>     The computational range for the singular values can be the full range
- *>  in software, independent of D. For more details see [1], [2].
+ *>  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
- *>     The computational range for the singular values can be the full range
+ *>  & LAPACK routines called by ZGEJSV are implemented to work in that range.
- *>  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
+ *>  If that is not the case, then the restriction for safe computation with
- *>  & LAPACK routines called by ZGEJSV are implemented to work in that range.
+ *>  the singular values in the range of normalized IEEE numbers is that the
- *>  If that is not the case, then the restriction for safe computation with
+ *>  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
- *>  the singular values in the range of normalized IEEE numbers is that the
+ *>  overflow. This code (ZGEJSV) is best used in this restricted range,
- *>  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
+ *>  meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
- *>  overflow. This code (ZGEJSV) is best used in this restricted range,
+ *>  returned as zeros. See JOBR for details on this.
- *>  meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
+ *>     Further, this implementation is somewhat slower than the one described
- *>  returned as zeros. See JOBR for details on this.
+ *>  in [1,2] due to replacement of some non-LAPACK components, and because
- *>     Further, this implementation is somewhat slower than the one described
+ *>  the choice of some tuning parameters in the iterative part (ZGESVJ) is
- *>  in [1,2] due to replacement of some non-LAPACK components, and because
+ *>  left to the implementer on a particular machine.
- *>  the choice of some tuning parameters in the iterative part (ZGESVJ) is
+ *>     The rank revealing QR factorization (in this code: ZGEQP3) should be
- *>  left to the implementer on a particular machine.
+ *>  implemented as in [3]. We have a new version of ZGEQP3 under development
- *>     The rank revealing QR factorization (in this code: ZGEQP3) should be
+ *>  that is more robust than the current one in LAPACK, with a cleaner cut in
- *>  implemented as in [3]. We have a new version of ZGEQP3 under development
+ *>  rank deficient cases. It will be available in the SIGMA library [4].
- *>  that is more robust than the current one in LAPACK, with a cleaner cut in
+ *>  If M is much larger than N, it is obvious that the initial QRF with
- *>  rank deficient cases. It will be available in the SIGMA library [4].
+ *>  column pivoting can be preprocessed by the QRF without pivoting. That
- *>  If M is much larger than N, it is obvious that the initial QRF with
+ *>  well known trick is not used in ZGEJSV because in some cases heavy row
- *>  column pivoting can be preprocessed by the QRF without pivoting. That
+ *>  weighting can be treated with complete pivoting. The overhead in cases
- *>  well known trick is not used in ZGEJSV because in some cases heavy row
+ *>  M much larger than N is then only due to pivoting, but the benefits in
- *>  weighting can be treated with complete pivoting. The overhead in cases
+ *>  terms of accuracy have prevailed. The implementer/user can incorporate
- *>  M much larger than N is then only due to pivoting, but the benefits in
+ *>  this extra QRF step easily. The implementer can also improve data movement
- *>  terms of accuracy have prevailed. The implementer/user can incorporate
+ *>  (matrix transpose, matrix copy, matrix transposed copy) - this
- *>  this extra QRF step easily. The implementer can also improve data movement
+ *>  implementation of ZGEJSV uses only the simplest, naive data movement.
- *>  (matrix transpose, matrix copy, matrix transposed copy) - this
+ *> \endverbatim
- *>  implementation of ZGEJSV uses only the simplest, naive data movement.
+ *
- *> \endverbatim
+ *> \par Contributor:
- *
+ *  ==================
- *> \par Contributor:
+ *>
- *  ==================
+ *>  Zlatko Drmac, Department of Mathematics, Faculty of Science,
- *>
+ *>  University of Zagreb (Zagreb, Croatia); drmac@math.hr
- *>  Zlatko Drmac, Department of Mathematics, Faculty of Science,
+ *
- *>  University of Zagreb (Zagreb, Croatia); drmac@math.hr
+ *> \par References:
- *
+ *  ================
- *> \par References:
+ *>
- *  ================
+ *> \verbatim
  *>
- *> \verbatim
+ *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
- *>
+ *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
- *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
+ *>     LAPACK Working note 169.
- *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
+ *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
- *>     LAPACK Working note 169.
+ *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
- *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
+ *>     LAPACK Working note 170.
- *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
+ *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
- *>     LAPACK Working note 170.
+ *>     factorization software - a case study.
- *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
+ *>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
- *>     factorization software - a case study.
+ *>     LAPACK Working note 176.
- *>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
+ *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
- *>     LAPACK Working note 176.
+ *>     QSVD, (H,K)-SVD computations.
- *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
+ *>     Department of Mathematics, University of Zagreb, 2008, 2016.
- *>     QSVD, (H,K)-SVD computations.
+ *> \endverbatim
- *>     Department of Mathematics, University of Zagreb, 2008, 2016.
+ *
- *> \endverbatim
+ *>  \par Bugs, examples and comments:
- *
+ *   =================================
- *>  \par Bugs, examples and comments:
+ *>
- *   =================================
+ *>  Please report all bugs and send interesting examples and/or comments to
- *>
+ *>  drmac@math.hr. Thank you.
- *>  Please report all bugs and send interesting examples and/or comments to
+ *>
- *>  drmac@math.hr. Thank you.
+ *  =====================================================================
- *>
+       SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
- *  =====================================================================
+      $                   M, N, A, LDA, SVA, U, LDU, V, LDV,
-       SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
+      $                   CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
-      $                   M, N, A, LDA, SVA, U, LDU, V, LDV,
+ *
-      $                   CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
+ *  -- LAPACK computational routine --
- *
+ *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
- *  -- LAPACK computational routine (version 3.7.0) --
+ *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
- *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+ *
- *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+ *     .. Scalar Arguments ..
- *     December 2016
+       IMPLICIT    NONE
- *
+       INTEGER     INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N
- *     .. Scalar Arguments ..
+ *     ..
-       IMPLICIT    NONE
+ *     .. Array Arguments ..
-       INTEGER     INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N
+       COMPLEX*16       A( LDA, * ), U( LDU, * ), V( LDV, * ),
- *     ..
+      $                 CWORK( LWORK )
- *     .. Array Arguments ..
+       DOUBLE PRECISION SVA( N ), RWORK( LRWORK )
-       COMPLEX*16       A( LDA, * ), U( LDU, * ), V( LDV, * ),
+       INTEGER          IWORK( * )
-      $                 CWORK( LWORK )
+       CHARACTER*1      JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
-       DOUBLE PRECISION SVA( N ), RWORK( LRWORK )
+ *     ..
-       INTEGER          IWORK( * )
+ *
-       CHARACTER*1      JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
+ *  ===========================================================================
- *     ..
+ *
- *
+ *     .. Local Parameters ..
- *  ===========================================================================
+       DOUBLE PRECISION ZERO, ONE
- *
+       PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
- *     .. Local Parameters ..
+       COMPLEX*16 CZERO, CONE
-       DOUBLE PRECISION ZERO, ONE
+       PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) )
-       PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+ *     ..
-       COMPLEX*16 CZERO, CONE
+ *     .. Local Scalars ..
-       PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) )
+       COMPLEX*16       CTEMP
- *     ..
+       DOUBLE PRECISION AAPP,    AAQQ,   AATMAX, AATMIN, BIG,    BIG1,
- *     .. Local Scalars ..
+      $                 COND_OK, CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,
-       COMPLEX*16       CTEMP
+      $                 MAXPRJ,  SCALEM, SCONDA, SFMIN,  SMALL,  TEMP1,
-       DOUBLE PRECISION AAPP,    AAQQ,   AATMAX, AATMIN, BIG,    BIG1,
+      $                 USCAL1,  USCAL2, XSC
-      $                 COND_OK, CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,
+       INTEGER IERR,   N1,     NR,     NUMRANK,        p, q,   WARNING
-      $                 MAXPRJ,  SCALEM, SCONDA, SFMIN,  SMALL,  TEMP1,
+       LOGICAL ALMORT, DEFR,   ERREST, GOSCAL,  JRACC,  KILL,   LQUERY,
-      $                 USCAL1,  USCAL2, XSC
+      $        LSVEC,  L2ABER, L2KILL, L2PERT,  L2RANK, L2TRAN, NOSCAL,
-       INTEGER IERR,   N1,     NR,     NUMRANK,        p, q,   WARNING
+      $        ROWPIV, RSVEC,  TRANSP
-       LOGICAL ALMORT, DEFR,   ERREST, GOSCAL,  JRACC,  KILL,   LQUERY,
+ *
-      $        LSVEC,  L2ABER, L2KILL, L2PERT,  L2RANK, L2TRAN, NOSCAL,
+       INTEGER OPTWRK, MINWRK, MINRWRK, MINIWRK
-      $        ROWPIV, RSVEC,  TRANSP
+       INTEGER LWCON,  LWLQF, LWQP3, LWQRF, LWUNMLQ, LWUNMQR, LWUNMQRM,
- *
+      $        LWSVDJ, LWSVDJV, LRWQP3, LRWCON, LRWSVDJ, IWOFF
-       INTEGER OPTWRK, MINWRK, MINRWRK, MINIWRK
+       INTEGER LWRK_ZGELQF, LWRK_ZGEQP3,  LWRK_ZGEQP3N, LWRK_ZGEQRF,
-       INTEGER LWCON,  LWLQF, LWQP3, LWQRF, LWUNMLQ, LWUNMQR, LWUNMQRM,
+      $        LWRK_ZGESVJ, LWRK_ZGESVJV, LWRK_ZGESVJU, LWRK_ZUNMLQ,
-      $        LWSVDJ, LWSVDJV, LRWQP3, LRWCON, LRWSVDJ, IWOFF
+      $        LWRK_ZUNMQR, LWRK_ZUNMQRM
-       INTEGER LWRK_ZGELQF, LWRK_ZGEQP3,  LWRK_ZGEQP3N, LWRK_ZGEQRF,
+ *     ..
-      $        LWRK_ZGESVJ, LWRK_ZGESVJV, LWRK_ZGESVJU, LWRK_ZUNMLQ,
+ *     .. Local Arrays
-      $        LWRK_ZUNMQR, LWRK_ZUNMQRM
+       COMPLEX*16         CDUMMY(1)
- *     ..
+       DOUBLE PRECISION   RDUMMY(1)
- *     .. Local Arrays
+ *
-       COMPLEX*16         CDUMMY(1)
+ *     .. Intrinsic Functions ..
-       DOUBLE PRECISION   RDUMMY(1)
+       INTRINSIC ABS, DCMPLX, CONJG, DLOG, MAX, MIN, DBLE, NINT, SQRT
- *
+ *     ..
- *     .. Intrinsic Functions ..
+ *     .. External Functions ..
-       INTRINSIC ABS, DCMPLX, CONJG, DLOG, MAX, MIN, DBLE, NINT, SQRT
+       DOUBLE PRECISION      DLAMCH, DZNRM2
- *     ..
+       INTEGER   IDAMAX, IZAMAX
- *     .. External Functions ..
+       LOGICAL   LSAME
-       DOUBLE PRECISION      DLAMCH, DZNRM2
+       EXTERNAL  IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2
-       INTEGER   IDAMAX, IZAMAX
+ *     ..
-       LOGICAL   LSAME
+ *     .. External Subroutines ..
-       EXTERNAL  IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2
+       EXTERNAL  DLASSQ, ZCOPY,  ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLAPMR,
- *     ..
+      $          ZLASCL, DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,
- *     .. External Subroutines ..
+      $          ZUNMQR, ZPOCON, DSCAL,  ZDSCAL, ZSWAP,  ZTRSM,  ZLACGV,
-       EXTERNAL  DLASSQ, ZCOPY,  ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLAPMR,
+      $          XERBLA
-      $          ZLASCL, DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,
+ *
-      $          ZUNMQR, ZPOCON, DSCAL,  ZDSCAL, ZSWAP,  ZTRSM,  ZLACGV,
+       EXTERNAL  ZGESVJ
-      $          XERBLA
+ *     ..
  *
-       EXTERNAL  ZGESVJ
+ *     Test the input arguments
- *     ..
+ *
- *
+       LSVEC  = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
- *     Test the input arguments
+       JRACC  = LSAME( JOBV, 'J' )
- *
+       RSVEC  = LSAME( JOBV, 'V' ) .OR. JRACC
-       LSVEC  = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
+       ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
-       JRACC  = LSAME( JOBV, 'J' )
+       L2RANK = LSAME( JOBA, 'R' )
-       RSVEC  = LSAME( JOBV, 'V' ) .OR. JRACC
+       L2ABER = LSAME( JOBA, 'A' )
-       ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
+       ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
-       L2RANK = LSAME( JOBA, 'R' )
+       L2TRAN = LSAME( JOBT, 'T' ) .AND. ( M .EQ. N )
-       L2ABER = LSAME( JOBA, 'A' )
+       L2KILL = LSAME( JOBR, 'R' )
-       ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
+       DEFR   = LSAME( JOBR, 'N' )
-       L2TRAN = LSAME( JOBT, 'T' ) .AND. ( M .EQ. N )
+       L2PERT = LSAME( JOBP, 'P' )
-       L2KILL = LSAME( JOBR, 'R' )
+ *
-       DEFR   = LSAME( JOBR, 'N' )
+       LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
-       L2PERT = LSAME( JOBP, 'P' )
+ *
- *
+       IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
-       LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
+      $     ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
- *
+          INFO = - 1
-       IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
+       ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR.
-      $     ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
+      $   ( LSAME( JOBU, 'W' ) .AND. RSVEC .AND. L2TRAN ) ) ) THEN
-          INFO = - 1
+          INFO = - 2
-       ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR.
+       ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
-      $   ( LSAME( JOBU, 'W' ) .AND. RSVEC .AND. L2TRAN ) ) ) THEN
+      $   ( LSAME( JOBV, 'W' ) .AND. LSVEC .AND. L2TRAN ) ) ) THEN
-          INFO = - 2
+          INFO = - 3
-       ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
+       ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) )    THEN
-      $   ( LSAME( JOBV, 'W' ) .AND. LSVEC .AND. L2TRAN ) ) ) THEN
+          INFO = - 4
-          INFO = - 3
+       ELSE IF ( .NOT. ( LSAME(JOBT,'T') .OR. LSAME(JOBT,'N') ) ) THEN
-       ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) )    THEN
+          INFO = - 5
-          INFO = - 4
+       ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
-       ELSE IF ( .NOT. ( LSAME(JOBT,'T') .OR. LSAME(JOBT,'N') ) ) THEN
+          INFO = - 6
-          INFO = - 5
+       ELSE IF ( M .LT. 0 ) THEN
-       ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
+          INFO = - 7
-          INFO = - 6
+       ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
-       ELSE IF ( M .LT. 0 ) THEN
+          INFO = - 8
-          INFO = - 7
+       ELSE IF ( LDA .LT. M ) THEN
-       ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
+          INFO = - 10
-          INFO = - 8
+       ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
-       ELSE IF ( LDA .LT. M ) THEN
+          INFO = - 13
-          INFO = - 10
+       ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
-       ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
+          INFO = - 15
-          INFO = - 13
+       ELSE
-       ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
+ *        #:)
-          INFO = - 15
+          INFO = 0
-       ELSE
+       END IF
- *        #:)
+ *
-          INFO = 0
+       IF ( INFO .EQ. 0 ) THEN
-       END IF
+ *         .. compute the minimal and the optimal workspace lengths
- *
+ *         [[The expressions for computing the minimal and the optimal
-       IF ( INFO .EQ. 0 ) THEN
+ *         values of LCWORK, LRWORK are written with a lot of redundancy and
- *         .. compute the minimal and the optimal workspace lengths
+ *         can be simplified. However, this verbose form is useful for
- *         [[The expressions for computing the minimal and the optimal
+ *         maintenance and modifications of the code.]]
- *         values of LCWORK, LRWORK are written with a lot of redundancy and
+ *
- *         can be simplified. However, this verbose form is useful for
+ *        .. minimal workspace length for ZGEQP3 of an M x N matrix,
- *         maintenance and modifications of the code.]]
+ *         ZGEQRF of an N x N matrix, ZGELQF of an N x N matrix,
- *
+ *         ZUNMLQ for computing N x N matrix, ZUNMQR for computing N x N
- *        .. minimal workspace length for ZGEQP3 of an M x N matrix,
+ *         matrix, ZUNMQR for computing M x N matrix, respectively.
- *         ZGEQRF of an N x N matrix, ZGELQF of an N x N matrix,
+           LWQP3 = N+1
- *         ZUNMLQ for computing N x N matrix, ZUNMQR for computing N x N
+           LWQRF = MAX( 1, N )
- *         matrix, ZUNMQR for computing M x N matrix, respectively.
+           LWLQF = MAX( 1, N )
-           LWQP3 = N+1
+           LWUNMLQ  = MAX( 1, N )
-           LWQRF = MAX( 1, N )
+           LWUNMQR  = MAX( 1, N )
-           LWLQF = MAX( 1, N )
+           LWUNMQRM = MAX( 1, M )
-           LWUNMLQ  = MAX( 1, N )
+ *        .. minimal workspace length for ZPOCON of an N x N matrix
-           LWUNMQR  = MAX( 1, N )
+           LWCON = 2 * N
-           LWUNMQRM = MAX( 1, M )
+ *        .. minimal workspace length for ZGESVJ of an N x N matrix,
- *        .. minimal workspace length for ZPOCON of an N x N matrix
+ *         without and with explicit accumulation of Jacobi rotations
-           LWCON = 2 * N
+           LWSVDJ  = MAX( 2 * N, 1 )
- *        .. minimal workspace length for ZGESVJ of an N x N matrix,
+           LWSVDJV = MAX( 2 * N, 1 )
- *         without and with explicit accumulation of Jacobi rotations
+ *         .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ
-           LWSVDJ  = MAX( 2 * N, 1 )
+           LRWQP3  = 2 * N
-           LWSVDJV = MAX( 2 * N, 1 )
+           LRWCON  = N
- *         .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ
+           LRWSVDJ = N
-           LRWQP3  = N
+           IF ( LQUERY ) THEN
-           LRWCON  = N
+               CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1,
-           LRWSVDJ = N
+      $             RDUMMY, IERR )
-           IF ( LQUERY ) THEN
+               LWRK_ZGEQP3 = INT( CDUMMY(1) )
-               CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1,
+               CALL ZGEQRF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
-      $             RDUMMY, IERR )
+               LWRK_ZGEQRF = INT( CDUMMY(1) )
-               LWRK_ZGEQP3 = CDUMMY(1)
+               CALL ZGELQF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
-               CALL ZGEQRF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
+               LWRK_ZGELQF = INT( CDUMMY(1) )
-               LWRK_ZGEQRF = CDUMMY(1)
+           END IF
-               CALL ZGELQF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
+           MINWRK  = 2
-               LWRK_ZGELQF = CDUMMY(1)
+           OPTWRK  = 2
-           END IF
+           MINIWRK = N
-           MINWRK  = 2
+           IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN
-           OPTWRK  = 2
+ *             .. minimal and optimal sizes of the complex workspace if
-           MINIWRK = N
+ *             only the singular values are requested
-           IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN
+               IF ( ERREST ) THEN
- *             .. minimal and optimal sizes of the complex workspace if
+                   MINWRK = MAX( N+LWQP3, N**2+LWCON, N+LWQRF, LWSVDJ )
- *             only the singular values are requested
+               ELSE
-               IF ( ERREST ) THEN
+                   MINWRK = MAX( N+LWQP3, N+LWQRF, LWSVDJ )
-                   MINWRK = MAX( N+LWQP3, N**2+LWCON, N+LWQRF, LWSVDJ )
+               END IF
-               ELSE
+               IF ( LQUERY ) THEN
-                   MINWRK = MAX( N+LWQP3, N+LWQRF, LWSVDJ )
+                   CALL ZGESVJ( 'L', 'N', 'N', N, N, A, LDA, SVA, N, V,
-               END IF
+      $                 LDV, CDUMMY, -1, RDUMMY, -1, IERR )
-               IF ( LQUERY ) THEN
+                   LWRK_ZGESVJ = INT( CDUMMY(1) )
-                   CALL ZGESVJ( 'L', 'N', 'N', N, N, A, LDA, SVA, N, V,
+                   IF ( ERREST ) THEN
-      $                 LDV, CDUMMY, -1, RDUMMY, -1, IERR )
+                       OPTWRK = MAX( N+LWRK_ZGEQP3, N**2+LWCON,
-                   LWRK_ZGESVJ = CDUMMY(1)
+      $                              N+LWRK_ZGEQRF, LWRK_ZGESVJ )
-                   IF ( ERREST ) THEN
+                   ELSE
-                       OPTWRK = MAX( N+LWRK_ZGEQP3, N**2+LWCON,
+                       OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWRK_ZGEQRF,
-      $                              N+LWRK_ZGEQRF, LWRK_ZGESVJ )
+      $                              LWRK_ZGESVJ )
-                   ELSE
+                   END IF
-                       OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWRK_ZGEQRF,
+               END IF
-      $                              LWRK_ZGESVJ )
+               IF ( L2TRAN .OR. ROWPIV ) THEN
-                   END IF
+                   IF ( ERREST ) THEN
-               END IF
+                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWCON, LRWSVDJ )
-               IF ( L2TRAN .OR. ROWPIV ) THEN
+                   ELSE
-                   IF ( ERREST ) THEN
+                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ )
-                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWCON, LRWSVDJ )
+                   END IF
-                   ELSE
+               ELSE
-                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ )
+                   IF ( ERREST ) THEN
-                   END IF
+                      MINRWRK = MAX( 7, LRWQP3, LRWCON, LRWSVDJ )
-               ELSE
+                   ELSE
-                   IF ( ERREST ) THEN
+                      MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
-                      MINRWRK = MAX( 7, LRWQP3, LRWCON, LRWSVDJ )
+                   END IF
-                   ELSE
+               END IF
-                      MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
+               IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
-                   END IF
+           ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
-               END IF
+ *            .. minimal and optimal sizes of the complex workspace if the
-               IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
+ *            singular values and the right singular vectors are requested
-           ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
+              IF ( ERREST ) THEN
- *            .. minimal and optimal sizes of the complex workspace if the
+                  MINWRK = MAX( N+LWQP3, LWCON, LWSVDJ, N+LWLQF,
- *            singular values and the right singular vectors are requested
+      $                         2*N+LWQRF, N+LWSVDJ, N+LWUNMLQ )
-              IF ( ERREST ) THEN
+              ELSE
-                  MINWRK = MAX( N+LWQP3, LWCON, LWSVDJ, N+LWLQF,
+                  MINWRK = MAX( N+LWQP3, LWSVDJ, N+LWLQF, 2*N+LWQRF,
-      $                         2*N+LWQRF, N+LWSVDJ, N+LWUNMLQ )
+      $                         N+LWSVDJ, N+LWUNMLQ )
-              ELSE
+              END IF
-                  MINWRK = MAX( N+LWQP3, LWSVDJ, N+LWLQF, 2*N+LWQRF,
+              IF ( LQUERY ) THEN
-      $                         N+LWSVDJ, N+LWUNMLQ )
+                  CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
-              END IF
+      $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )
-              IF ( LQUERY ) THEN
+                  LWRK_ZGESVJ = INT( CDUMMY(1) )
-                  CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
+                  CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
-      $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )
+      $                V, LDV, CDUMMY, -1, IERR )
-                  LWRK_ZGESVJ = CDUMMY(1)
+                  LWRK_ZUNMLQ = INT( CDUMMY(1) )
-                  CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
+                  IF ( ERREST ) THEN
-      $                V, LDV, CDUMMY, -1, IERR )
+                  OPTWRK = MAX( N+LWRK_ZGEQP3, LWCON, LWRK_ZGESVJ,
-                  LWRK_ZUNMLQ = CDUMMY(1)
+      $                         N+LWRK_ZGELQF, 2*N+LWRK_ZGEQRF,
-                  IF ( ERREST ) THEN
+      $                         N+LWRK_ZGESVJ,  N+LWRK_ZUNMLQ )
-                  OPTWRK = MAX( N+LWRK_ZGEQP3, LWCON, LWRK_ZGESVJ,
+                  ELSE
-      $                         N+LWRK_ZGELQF, 2*N+LWRK_ZGEQRF,
+                  OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVJ,N+LWRK_ZGELQF,
-      $                         N+LWRK_ZGESVJ,  N+LWRK_ZUNMLQ )
+      $                         2*N+LWRK_ZGEQRF, N+LWRK_ZGESVJ,
-                  ELSE
+      $                         N+LWRK_ZUNMLQ )
-                  OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVJ,N+LWRK_ZGELQF,
+                  END IF
-      $                         2*N+LWRK_ZGEQRF, N+LWRK_ZGESVJ,
+              END IF
-      $                         N+LWRK_ZUNMLQ )
+              IF ( L2TRAN .OR. ROWPIV ) THEN
-                  END IF
+                   IF ( ERREST ) THEN
-              END IF
+                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ, LRWCON )
-              IF ( L2TRAN .OR. ROWPIV ) THEN
+                   ELSE
-                   IF ( ERREST ) THEN
+                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ )
-                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ, LRWCON )
+                   END IF
-                   ELSE
+              ELSE
-                      MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ )
+                   IF ( ERREST ) THEN
-                   END IF
+                      MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
-              ELSE
+                   ELSE
-                   IF ( ERREST ) THEN
+                      MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
-                      MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
+                   END IF
-                   ELSE
+              END IF
-                      MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
+              IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
-                   END IF
+           ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
-              END IF
+ *            .. minimal and optimal sizes of the complex workspace if the
-              IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
+ *            singular values and the left singular vectors are requested
-           ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
+              IF ( ERREST ) THEN
- *            .. minimal and optimal sizes of the complex workspace if the
+                  MINWRK = N + MAX( LWQP3,LWCON,N+LWQRF,LWSVDJ,LWUNMQRM )
- *            singular values and the left singular vectors are requested
+              ELSE
-              IF ( ERREST ) THEN
+                  MINWRK = N + MAX( LWQP3, N+LWQRF, LWSVDJ, LWUNMQRM )
-                  MINWRK = N + MAX( LWQP3,LWCON,N+LWQRF,LWSVDJ,LWUNMQRM )
+              END IF
-              ELSE
+              IF ( LQUERY ) THEN
-                  MINWRK = N + MAX( LWQP3, N+LWQRF, LWSVDJ, LWUNMQRM )
+                  CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
-              END IF
+      $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )
-              IF ( LQUERY ) THEN
+                  LWRK_ZGESVJ = INT( CDUMMY(1) )
-                  CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
+                  CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
-      $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )
+      $               LDU, CDUMMY, -1, IERR )
-                  LWRK_ZGESVJ = CDUMMY(1)
+                  LWRK_ZUNMQRM = INT( CDUMMY(1) )
-                  CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
+                  IF ( ERREST ) THEN
-      $               LDU, CDUMMY, -1, IERR )
+                  OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, N+LWRK_ZGEQRF,
-                  LWRK_ZUNMQRM = CDUMMY(1)
+      $                             LWRK_ZGESVJ, LWRK_ZUNMQRM )
-                  IF ( ERREST ) THEN
+                  ELSE
-                  OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, N+LWRK_ZGEQRF,
+                  OPTWRK = N + MAX( LWRK_ZGEQP3, N+LWRK_ZGEQRF,
       $                             LWRK_ZGESVJ, LWRK_ZUNMQRM )
-                  ELSE
+                  END IF
-                  OPTWRK = N + MAX( LWRK_ZGEQP3, N+LWRK_ZGEQRF,
+              END IF
-      $                             LWRK_ZGESVJ, LWRK_ZUNMQRM )
+              IF ( L2TRAN .OR. ROWPIV ) THEN
-                  END IF
+                  IF ( ERREST ) THEN
-              END IF
+                     MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ, LRWCON )
-              IF ( L2TRAN .OR. ROWPIV ) THEN
+                  ELSE
-                  IF ( ERREST ) THEN
+                     MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ )
-                     MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ, LRWCON )
+                  END IF
-                  ELSE
+              ELSE
-                     MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ )
+                  IF ( ERREST ) THEN
-                  END IF
+                     MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
-              ELSE
+                  ELSE
-                  IF ( ERREST ) THEN
+                     MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
-                     MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
+                  END IF
-                  ELSE
+              END IF
-                     MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
+              IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
-                  END IF
+           ELSE
-              END IF
+ *            .. minimal and optimal sizes of the complex workspace if the
-              IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
+ *            full SVD is requested
-           ELSE
+              IF ( .NOT. JRACC ) THEN
- *            .. minimal and optimal sizes of the complex workspace if the
+                  IF ( ERREST ) THEN
- *            full SVD is requested
+                     MINWRK = MAX( N+LWQP3, N+LWCON,  2*N+N**2+LWCON,
-              IF ( .NOT. JRACC ) THEN
+      $                         2*N+LWQRF,         2*N+LWQP3,
-                  IF ( ERREST ) THEN
+      $                         2*N+N**2+N+LWLQF,  2*N+N**2+N+N**2+LWCON,
-                     MINWRK = MAX( N+LWQP3, N+LWCON,  2*N+N**2+LWCON,
+      $                         2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
-      $                         2*N+LWQRF,         2*N+LWQP3,
+      $                         2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
-      $                         2*N+N**2+N+LWLQF,  2*N+N**2+N+N**2+LWCON,
+      $                         N+N**2+LWSVDJ,   N+LWUNMQRM )
-      $                         2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
+                  ELSE
-      $                         2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
+                     MINWRK = MAX( N+LWQP3,        2*N+N**2+LWCON,
-      $                         N+N**2+LWSVDJ,   N+LWUNMQRM )
+      $                         2*N+LWQRF,         2*N+LWQP3,
-                  ELSE
+      $                         2*N+N**2+N+LWLQF,  2*N+N**2+N+N**2+LWCON,
-                     MINWRK = MAX( N+LWQP3,        2*N+N**2+LWCON,
+      $                         2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
-      $                         2*N+LWQRF,         2*N+LWQP3,
+      $                         2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
-      $                         2*N+N**2+N+LWLQF,  2*N+N**2+N+N**2+LWCON,
+      $                         N+N**2+LWSVDJ,      N+LWUNMQRM )
-      $                         2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
+                  END IF
-      $                         2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
+                  MINIWRK = MINIWRK + N
-      $                         N+N**2+LWSVDJ,      N+LWUNMQRM )
+                  IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
-                  END IF
+              ELSE
-                  MINIWRK = MINIWRK + N
+                  IF ( ERREST ) THEN
-                  IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
+                     MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+LWQRF,
-              ELSE
+      $                         2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
-                  IF ( ERREST ) THEN
+      $                         N+LWUNMQRM )
-                     MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+LWQRF,
+                  ELSE
-      $                         2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
+                     MINWRK = MAX( N+LWQP3, 2*N+LWQRF,
-      $                         N+LWUNMQRM )
+      $                         2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
-                  ELSE
+      $                         N+LWUNMQRM )
-                     MINWRK = MAX( N+LWQP3, 2*N+LWQRF,
+                  END IF
-      $                         2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
+                  IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
-      $                         N+LWUNMQRM )
+              END IF
-                  END IF
+              IF ( LQUERY ) THEN
-                  IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
+                  CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
-              END IF
+      $                LDU, CDUMMY, -1, IERR )
-              IF ( LQUERY ) THEN
+                  LWRK_ZUNMQRM = INT( CDUMMY(1) )
-                  CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
+                  CALL ZUNMQR( 'L', 'N', N, N, N, A, LDA, CDUMMY, U,
       $                LDU, CDUMMY, -1, IERR )
-                  LWRK_ZUNMQRM = CDUMMY(1)
+                  LWRK_ZUNMQR = INT( CDUMMY(1) )
-                  CALL ZUNMQR( 'L', 'N', N, N, N, A, LDA, CDUMMY, U,
+                  IF ( .NOT. JRACC ) THEN
-      $                LDU, CDUMMY, -1, IERR )
+                      CALL ZGEQP3( N,N, A, LDA, IWORK, CDUMMY,CDUMMY, -1,
-                  LWRK_ZUNMQR = CDUMMY(1)
+      $                    RDUMMY, IERR )
-                  IF ( .NOT. JRACC ) THEN
+                      LWRK_ZGEQP3N = INT( CDUMMY(1) )
-                      CALL ZGEQP3( N,N, A, LDA, IWORK, CDUMMY,CDUMMY, -1,
+                      CALL ZGESVJ( 'L', 'U', 'N', N, N, U, LDU, SVA,
-      $                    RDUMMY, IERR )
+      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
-                      LWRK_ZGEQP3N = CDUMMY(1)
+                      LWRK_ZGESVJ = INT( CDUMMY(1) )
-                      CALL ZGESVJ( 'L', 'U', 'N', N, N, U, LDU, SVA,
+                      CALL ZGESVJ( 'U', 'U', 'N', N, N, U, LDU, SVA,
       $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
-                      LWRK_ZGESVJ = CDUMMY(1)
+                      LWRK_ZGESVJU = INT( CDUMMY(1) )
-                      CALL ZGESVJ( 'U', 'U', 'N', N, N, U, LDU, SVA,
+                      CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
       $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
-                      LWRK_ZGESVJU = CDUMMY(1)
+                      LWRK_ZGESVJV = INT( CDUMMY(1) )
-                      CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
+                      CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
-      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
+      $                    V, LDV, CDUMMY, -1, IERR )
-                      LWRK_ZGESVJV = CDUMMY(1)
+                      LWRK_ZUNMLQ = INT( CDUMMY(1) )
-                      CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
+                      IF ( ERREST ) THEN
-      $                    V, LDV, CDUMMY, -1, IERR )
+                        OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
-                      LWRK_ZUNMLQ = CDUMMY(1)
+      $                          2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
-                      IF ( ERREST ) THEN
+      $                          2*N+LWRK_ZGEQP3N,
-                        OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
+      $                          2*N+N**2+N+LWRK_ZGELQF,
-      $                          2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
+      $                          2*N+N**2+N+N**2+LWCON,
-      $                          2*N+LWRK_ZGEQP3N,
+      $                          2*N+N**2+N+LWRK_ZGESVJ,
-      $                          2*N+N**2+N+LWRK_ZGELQF,
+      $                          2*N+N**2+N+LWRK_ZGESVJV,
-      $                          2*N+N**2+N+N**2+LWCON,
+      $                          2*N+N**2+N+LWRK_ZUNMQR,
-      $                          2*N+N**2+N+LWRK_ZGESVJ,
+      $                          2*N+N**2+N+LWRK_ZUNMLQ,
-      $                          2*N+N**2+N+LWRK_ZGESVJV,
+      $                          N+N**2+LWRK_ZGESVJU,
-      $                          2*N+N**2+N+LWRK_ZUNMQR,
+      $                          N+LWRK_ZUNMQRM )
-      $                          2*N+N**2+N+LWRK_ZUNMLQ,
+                      ELSE
-      $                          N+N**2+LWRK_ZGESVJU,
+                        OPTWRK = MAX( N+LWRK_ZGEQP3,
-      $                          N+LWRK_ZUNMQRM )
+      $                          2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
-                      ELSE
+      $                          2*N+LWRK_ZGEQP3N,
-                        OPTWRK = MAX( N+LWRK_ZGEQP3,
+      $                          2*N+N**2+N+LWRK_ZGELQF,
-      $                          2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
+      $                          2*N+N**2+N+N**2+LWCON,
-      $                          2*N+LWRK_ZGEQP3N,
+      $                          2*N+N**2+N+LWRK_ZGESVJ,
-      $                          2*N+N**2+N+LWRK_ZGELQF,
+      $                          2*N+N**2+N+LWRK_ZGESVJV,
-      $                          2*N+N**2+N+N**2+LWCON,
+      $                          2*N+N**2+N+LWRK_ZUNMQR,
-      $                          2*N+N**2+N+LWRK_ZGESVJ,
+      $                          2*N+N**2+N+LWRK_ZUNMLQ,
-      $                          2*N+N**2+N+LWRK_ZGESVJV,
+      $                          N+N**2+LWRK_ZGESVJU,
-      $                          2*N+N**2+N+LWRK_ZUNMQR,
+      $                          N+LWRK_ZUNMQRM )
-      $                          2*N+N**2+N+LWRK_ZUNMLQ,
+                      END IF
-      $                          N+N**2+LWRK_ZGESVJU,
+                  ELSE
-      $                          N+LWRK_ZUNMQRM )
+                      CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
-                      END IF
+      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
-                  ELSE
+                      LWRK_ZGESVJV = INT( CDUMMY(1) )
-                      CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
+                      CALL ZUNMQR( 'L', 'N', N, N, N, CDUMMY, N, CDUMMY,
-      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
+      $                    V, LDV, CDUMMY, -1, IERR )
-                      LWRK_ZGESVJV = CDUMMY(1)
+                      LWRK_ZUNMQR = INT( CDUMMY(1) )
-                      CALL ZUNMQR( 'L', 'N', N, N, N, CDUMMY, N, CDUMMY,
+                      CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
-      $                    V, LDV, CDUMMY, -1, IERR )
+      $                    LDU, CDUMMY, -1, IERR )
-                      LWRK_ZUNMQR = CDUMMY(1)
+                      LWRK_ZUNMQRM = INT( CDUMMY(1) )
-                      CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
+                      IF ( ERREST ) THEN
-      $                    LDU, CDUMMY, -1, IERR )
+                         OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
-                      LWRK_ZUNMQRM = CDUMMY(1)
+      $                           2*N+LWRK_ZGEQRF, 2*N+N**2,
-                      IF ( ERREST ) THEN
+      $                           2*N+N**2+LWRK_ZGESVJV,
-                         OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
+      $                           2*N+N**2+N+LWRK_ZUNMQR,N+LWRK_ZUNMQRM )
-      $                           2*N+LWRK_ZGEQRF, 2*N+N**2,
+                      ELSE
-      $                           2*N+N**2+LWRK_ZGESVJV,
+                         OPTWRK = MAX( N+LWRK_ZGEQP3, 2*N+LWRK_ZGEQRF,
-      $                           2*N+N**2+N+LWRK_ZUNMQR,N+LWRK_ZUNMQRM )
+      $                           2*N+N**2, 2*N+N**2+LWRK_ZGESVJV,
-                      ELSE
+      $                           2*N+N**2+N+LWRK_ZUNMQR,
-                         OPTWRK = MAX( N+LWRK_ZGEQP3, 2*N+LWRK_ZGEQRF,
+      $                           N+LWRK_ZUNMQRM )
-      $                           2*N+N**2, 2*N+N**2+LWRK_ZGESVJV,
+                      END IF
-      $                           2*N+N**2+N+LWRK_ZUNMQR,
+                  END IF
-      $                           N+LWRK_ZUNMQRM )
+              END IF
-                      END IF
+              IF ( L2TRAN .OR. ROWPIV ) THEN
-                  END IF
+                  MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ, LRWCON )
-              END IF
+              ELSE
-              IF ( L2TRAN .OR. ROWPIV ) THEN
+                  MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
-                  MINRWRK = MAX( 7, 2*M,  LRWQP3, LRWSVDJ, LRWCON )
+              END IF
-              ELSE
+           END IF
-                  MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
+           MINWRK = MAX( 2, MINWRK )
-              END IF
+           OPTWRK = MAX( MINWRK, OPTWRK )
-           END IF
+           IF ( LWORK  .LT. MINWRK  .AND. (.NOT.LQUERY) ) INFO = - 17
-           MINWRK = MAX( 2, MINWRK )
+           IF ( LRWORK .LT. MINRWRK .AND. (.NOT.LQUERY) ) INFO = - 19
-           OPTWRK = MAX( 2, OPTWRK )
+       END IF
-           IF ( LWORK  .LT. MINWRK  .AND. (.NOT.LQUERY) ) INFO = - 17
+ *
-           IF ( LRWORK .LT. MINRWRK .AND. (.NOT.LQUERY) ) INFO = - 19
+       IF ( INFO .NE. 0 ) THEN
-       END IF
+ *       #:(
- *
+          CALL XERBLA( 'ZGEJSV', - INFO )
-       IF ( INFO .NE. 0 ) THEN
+          RETURN
- *       #:(
+       ELSE IF ( LQUERY ) THEN
-          CALL XERBLA( 'ZGEJSV', - INFO )
+           CWORK(1) = OPTWRK
-          RETURN
+           CWORK(2) = MINWRK
-       ELSE IF ( LQUERY ) THEN
+           RWORK(1) = MINRWRK
-           CWORK(1) = OPTWRK
+           IWORK(1) = MAX( 4, MINIWRK )
-           CWORK(2) = MINWRK
+           RETURN
-           RWORK(1) = MINRWRK
+       END IF
-           IWORK(1) = MAX( 4, MINIWRK )
+ *
-           RETURN
+ *     Quick return for void matrix (Y3K safe)
-       END IF
+ * #:)
- *
+       IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
- *     Quick return for void matrix (Y3K safe)
+          IWORK(1:4) = 0
- * #:)
+          RWORK(1:7) = 0
-       IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
+          RETURN
-          IWORK(1:4) = 0
+       ENDIF
-          RWORK(1:7) = 0
+ *
-          RETURN
+ *     Determine whether the matrix U should be M x N or M x M
-       ENDIF
+ *
- *
+       IF ( LSVEC ) THEN
- *     Determine whether the matrix U should be M x N or M x M
+          N1 = N
- *
+          IF ( LSAME( JOBU, 'F' ) ) N1 = M
-       IF ( LSVEC ) THEN
+       END IF
-          N1 = N
+ *
-          IF ( LSAME( JOBU, 'F' ) ) N1 = M
+ *     Set numerical parameters
-       END IF
+ *
- *
+ *!    NOTE: Make sure DLAMCH() does not fail on the target architecture.
- *     Set numerical parameters
+ *
- *
+       EPSLN = DLAMCH('Epsilon')
- *!    NOTE: Make sure DLAMCH() does not fail on the target architecture.
+       SFMIN = DLAMCH('SafeMinimum')
- *
+       SMALL = SFMIN / EPSLN
-       EPSLN = DLAMCH('Epsilon')
+       BIG   = DLAMCH('O')
-       SFMIN = DLAMCH('SafeMinimum')
+ *     BIG   = ONE / SFMIN
-       SMALL = SFMIN / EPSLN
+ *
-       BIG   = DLAMCH('O')
+ *     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
- *     BIG   = ONE / SFMIN
+ *
- *
+ *(!)  If necessary, scale SVA() to protect the largest norm from
- *     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
+ *     overflow. It is possible that this scaling pushes the smallest
- *
+ *     column norm left from the underflow threshold (extreme case).
- *(!)  If necessary, scale SVA() to protect the largest norm from
+ *
- *     overflow. It is possible that this scaling pushes the smallest
+       SCALEM  = ONE / SQRT(DBLE(M)*DBLE(N))
- *     column norm left from the underflow threshold (extreme case).
+       NOSCAL  = .TRUE.
- *
+       GOSCAL  = .TRUE.
-       SCALEM  = ONE / SQRT(DBLE(M)*DBLE(N))
+       DO 1874 p = 1, N
-       NOSCAL  = .TRUE.
+          AAPP = ZERO
-       GOSCAL  = .TRUE.
+          AAQQ = ONE
-       DO 1874 p = 1, N
+          CALL ZLASSQ( M, A(1,p), 1, AAPP, AAQQ )
-          AAPP = ZERO
+          IF ( AAPP .GT. BIG ) THEN
-          AAQQ = ONE
+             INFO = - 9
-          CALL ZLASSQ( M, A(1,p), 1, AAPP, AAQQ )
+             CALL XERBLA( 'ZGEJSV', -INFO )
-          IF ( AAPP .GT. BIG ) THEN
+             RETURN
-             INFO = - 9
+          END IF
-             CALL XERBLA( 'ZGEJSV', -INFO )
+          AAQQ = SQRT(AAQQ)
-             RETURN
+          IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL  ) THEN
-          END IF
+             SVA(p)  = AAPP * AAQQ
-          AAQQ = SQRT(AAQQ)
+          ELSE
-          IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL  ) THEN
+             NOSCAL  = .FALSE.
-             SVA(p)  = AAPP * AAQQ
+             SVA(p)  = AAPP * ( AAQQ * SCALEM )
-          ELSE
+             IF ( GOSCAL ) THEN
-             NOSCAL  = .FALSE.
+                GOSCAL = .FALSE.
-             SVA(p)  = AAPP * ( AAQQ * SCALEM )
+                CALL DSCAL( p-1, SCALEM, SVA, 1 )
-             IF ( GOSCAL ) THEN
+             END IF
-                GOSCAL = .FALSE.
+          END IF
-                CALL DSCAL( p-1, SCALEM, SVA, 1 )
+CONTINUE
-             END IF
+ *
-          END IF
+       IF ( NOSCAL ) SCALEM = ONE
-CONTINUE
+ *
- *
+       AAPP = ZERO
-       IF ( NOSCAL ) SCALEM = ONE
+       AAQQ = BIG
- *
+       DO 4781 p = 1, N
-       AAPP = ZERO
+          AAPP = MAX( AAPP, SVA(p) )
-       AAQQ = BIG
+          IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) )
-       DO 4781 p = 1, N
+CONTINUE
-          AAPP = MAX( AAPP, SVA(p) )
+ *
-          IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) )
+ *     Quick return for zero M x N matrix
-CONTINUE
+ * #:)
- *
+       IF ( AAPP .EQ. ZERO ) THEN
- *     Quick return for zero M x N matrix
+          IF ( LSVEC ) CALL ZLASET( 'G', M, N1, CZERO, CONE, U, LDU )
- * #:)
+          IF ( RSVEC ) CALL ZLASET( 'G', N, N,  CZERO, CONE, V, LDV )
-       IF ( AAPP .EQ. ZERO ) THEN
+          RWORK(1) = ONE
-          IF ( LSVEC ) CALL ZLASET( 'G', M, N1, CZERO, CONE, U, LDU )
+          RWORK(2) = ONE
-          IF ( RSVEC ) CALL ZLASET( 'G', N, N,  CZERO, CONE, V, LDV )
+          IF ( ERREST ) RWORK(3) = ONE
-          RWORK(1) = ONE
+          IF ( LSVEC .AND. RSVEC ) THEN
-          RWORK(2) = ONE
+             RWORK(4) = ONE
-          IF ( ERREST ) RWORK(3) = ONE
+             RWORK(5) = ONE
-          IF ( LSVEC .AND. RSVEC ) THEN
+          END IF
-             RWORK(4) = ONE
+          IF ( L2TRAN ) THEN
-             RWORK(5) = ONE
+             RWORK(6) = ZERO
-          END IF
+             RWORK(7) = ZERO
-          IF ( L2TRAN ) THEN
+          END IF
-             RWORK(6) = ZERO
+          IWORK(1) = 0
-             RWORK(7) = ZERO
+          IWORK(2) = 0
-          END IF
+          IWORK(3) = 0
-          IWORK(1) = 0
+          IWORK(4) = -1
-          IWORK(2) = 0
+          RETURN
-          IWORK(3) = 0
+       END IF
-          IWORK(4) = -1
+ *
-          RETURN
+ *     Issue warning if denormalized column norms detected. Override the
-       END IF
+ *     high relative accuracy request. Issue licence to kill nonzero columns
- *
+ *     (set them to zero) whose norm is less than sigma_max / BIG (roughly).
- *     Issue warning if denormalized column norms detected. Override the
+ * #:(
- *     high relative accuracy request. Issue licence to kill nonzero columns
+       WARNING = 0
- *     (set them to zero) whose norm is less than sigma_max / BIG (roughly).
+       IF ( AAQQ .LE. SFMIN ) THEN
- * #:(
+          L2RANK = .TRUE.
-       WARNING = 0
+          L2KILL = .TRUE.
-       IF ( AAQQ .LE. SFMIN ) THEN
+          WARNING = 1
-          L2RANK = .TRUE.
+       END IF
-          L2KILL = .TRUE.
+ *
-          WARNING = 1
+ *     Quick return for one-column matrix
-       END IF
+ * #:)
- *
+       IF ( N .EQ. 1 ) THEN
- *     Quick return for one-column matrix
+ *
- * #:)
+          IF ( LSVEC ) THEN
-       IF ( N .EQ. 1 ) THEN
+             CALL ZLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
- *
+             CALL ZLACPY( 'A', M, 1, A, LDA, U, LDU )
-          IF ( LSVEC ) THEN
+ *           computing all M left singular vectors of the M x 1 matrix
-             CALL ZLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
+             IF ( N1 .NE. N  ) THEN
-             CALL ZLACPY( 'A', M, 1, A, LDA, U, LDU )
+               CALL ZGEQRF( M, N, U,LDU, CWORK, CWORK(N+1),LWORK-N,IERR )
- *           computing all M left singular vectors of the M x 1 matrix
+               CALL ZUNGQR( M,N1,1, U,LDU,CWORK,CWORK(N+1),LWORK-N,IERR )
-             IF ( N1 .NE. N  ) THEN
+               CALL ZCOPY( M, A(1,1), 1, U(1,1), 1 )
-               CALL ZGEQRF( M, N, U,LDU, CWORK, CWORK(N+1),LWORK-N,IERR )
+             END IF
-               CALL ZUNGQR( M,N1,1, U,LDU,CWORK,CWORK(N+1),LWORK-N,IERR )
+          END IF
-               CALL ZCOPY( M, A(1,1), 1, U(1,1), 1 )
+          IF ( RSVEC ) THEN
-             END IF
+              V(1,1) = CONE
           END IF
-          IF ( RSVEC ) THEN
+          IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
-              V(1,1) = CONE
+             SVA(1)  = SVA(1) / SCALEM
-          END IF
+             SCALEM  = ONE
-          IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
+          END IF
-             SVA(1)  = SVA(1) / SCALEM
+          RWORK(1) = ONE / SCALEM
-             SCALEM  = ONE
+          RWORK(2) = ONE
-          END IF
+          IF ( SVA(1) .NE. ZERO ) THEN
-          RWORK(1) = ONE / SCALEM
+             IWORK(1) = 1
-          RWORK(2) = ONE
+             IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
-          IF ( SVA(1) .NE. ZERO ) THEN
+                IWORK(2) = 1
-             IWORK(1) = 1
+             ELSE
-             IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
+                IWORK(2) = 0
-                IWORK(2) = 1
+             END IF
-             ELSE
+          ELSE
-                IWORK(2) = 0
+             IWORK(1) = 0
-             END IF
+             IWORK(2) = 0
-          ELSE
+          END IF
-             IWORK(1) = 0
+          IWORK(3) = 0
-             IWORK(2) = 0
+          IWORK(4) = -1
-          END IF
+          IF ( ERREST ) RWORK(3) = ONE
-          IWORK(3) = 0
+          IF ( LSVEC .AND. RSVEC ) THEN
-          IWORK(4) = -1
+             RWORK(4) = ONE
-          IF ( ERREST ) RWORK(3) = ONE
+             RWORK(5) = ONE
-          IF ( LSVEC .AND. RSVEC ) THEN
+          END IF
-             RWORK(4) = ONE
+          IF ( L2TRAN ) THEN
-             RWORK(5) = ONE
+             RWORK(6) = ZERO
-          END IF
+             RWORK(7) = ZERO
-          IF ( L2TRAN ) THEN
+          END IF
-             RWORK(6) = ZERO
+          RETURN
-             RWORK(7) = ZERO
+ *
-          END IF
+       END IF
-          RETURN
+ *
- *
+       TRANSP = .FALSE.
-       END IF
+ *
- *
+       AATMAX = -ONE
-       TRANSP = .FALSE.
+       AATMIN =  BIG
- *
+       IF ( ROWPIV .OR. L2TRAN ) THEN
-       AATMAX = -ONE
+ *
-       AATMIN =  BIG
+ *     Compute the row norms, needed to determine row pivoting sequence
-       IF ( ROWPIV .OR. L2TRAN ) THEN
+ *     (in the case of heavily row weighted A, row pivoting is strongly
- *
+ *     advised) and to collect information needed to compare the
- *     Compute the row norms, needed to determine row pivoting sequence
+ *     structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.).
- *     (in the case of heavily row weighted A, row pivoting is strongly
+ *
- *     advised) and to collect information needed to compare the
+          IF ( L2TRAN ) THEN
- *     structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.).
+             DO 1950 p = 1, M
- *
+                XSC   = ZERO
-          IF ( L2TRAN ) THEN
+                TEMP1 = ONE
-             DO 1950 p = 1, M
+                CALL ZLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
-                XSC   = ZERO
+ *              ZLASSQ gets both the ell_2 and the ell_infinity norm
-                TEMP1 = ONE
+ *              in one pass through the vector
-                CALL ZLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
+                RWORK(M+p)  = XSC * SCALEM
- *              ZLASSQ gets both the ell_2 and the ell_infinity norm
+                RWORK(p)    = XSC * (SCALEM*SQRT(TEMP1))
- *              in one pass through the vector
+                AATMAX = MAX( AATMAX, RWORK(p) )
-                RWORK(M+p)  = XSC * SCALEM
+                IF (RWORK(p) .NE. ZERO)
-                RWORK(p)    = XSC * (SCALEM*SQRT(TEMP1))
+      $            AATMIN = MIN(AATMIN,RWORK(p))
-                AATMAX = MAX( AATMAX, RWORK(p) )
+      CONTINUE
-                IF (RWORK(p) .NE. ZERO)
+          ELSE
-      $            AATMIN = MIN(AATMIN,RWORK(p))
+             DO 1904 p = 1, M
-      CONTINUE
+                RWORK(M+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
-          ELSE
+                AATMAX = MAX( AATMAX, RWORK(M+p) )
-             DO 1904 p = 1, M
+                AATMIN = MIN( AATMIN, RWORK(M+p) )
-                RWORK(M+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
+      CONTINUE
-                AATMAX = MAX( AATMAX, RWORK(M+p) )
+          END IF
-                AATMIN = MIN( AATMIN, RWORK(M+p) )
+ *
-      CONTINUE
+       END IF
-          END IF
+ *
- *
+ *     For square matrix A try to determine whether A^*  would be better
-       END IF
+ *     input for the preconditioned Jacobi SVD, with faster convergence.
- *
+ *     The decision is based on an O(N) function of the vector of column
- *     For square matrix A try to determine whether A^*  would be better
+ *     and row norms of A, based on the Shannon entropy. This should give
- *     input for the preconditioned Jacobi SVD, with faster convergence.
+ *     the right choice in most cases when the difference actually matters.
- *     The decision is based on an O(N) function of the vector of column
+ *     It may fail and pick the slower converging side.
- *     and row norms of A, based on the Shannon entropy. This should give
+ *
- *     the right choice in most cases when the difference actually matters.
+       ENTRA  = ZERO
- *     It may fail and pick the slower converging side.
+       ENTRAT = ZERO
- *
+       IF ( L2TRAN ) THEN
-       ENTRA  = ZERO
+ *
-       ENTRAT = ZERO
+          XSC   = ZERO
-       IF ( L2TRAN ) THEN
+          TEMP1 = ONE
- *
+          CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
-          XSC   = ZERO
+          TEMP1 = ONE / TEMP1
-          TEMP1 = ONE
+ *
-          CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
+          ENTRA = ZERO
-          TEMP1 = ONE / TEMP1
+          DO 1113 p = 1, N
- *
+             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
-          ENTRA = ZERO
+             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
-          DO 1113 p = 1, N
+   CONTINUE
-             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
+          ENTRA = - ENTRA / DLOG(DBLE(N))
-             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
+ *
-   CONTINUE
+ *        Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex.
-          ENTRA = - ENTRA / DLOG(DBLE(N))
+ *        It is derived from the diagonal of  A^* * A.  Do the same with the
- *
+ *        diagonal of A * A^*, compute the entropy of the corresponding
- *        Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex.
+ *        probability distribution. Note that A * A^* and A^* * A have the
- *        It is derived from the diagonal of  A^* * A.  Do the same with the
+ *        same trace.
- *        diagonal of A * A^*, compute the entropy of the corresponding
+ *
- *        probability distribution. Note that A * A^* and A^* * A have the
+          ENTRAT = ZERO
- *        same trace.
+          DO 1114 p = 1, M
- *
+             BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1
-          ENTRAT = ZERO
+             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
-          DO 1114 p = 1, M
+   CONTINUE
-             BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1
+          ENTRAT = - ENTRAT / DLOG(DBLE(M))
-             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
+ *
-   CONTINUE
+ *        Analyze the entropies and decide A or A^*. Smaller entropy
-          ENTRAT = - ENTRAT / DLOG(DBLE(M))
+ *        usually means better input for the algorithm.
  *
- *        Analyze the entropies and decide A or A^*. Smaller entropy
+          TRANSP = ( ENTRAT .LT. ENTRA )
- *        usually means better input for the algorithm.
+ *
- *
+ *        If A^* is better than A, take the adjoint of A. This is allowed
-          TRANSP = ( ENTRAT .LT. ENTRA )
+ *        only for square matrices, M=N.
- *
+          IF ( TRANSP ) THEN
- *        If A^* is better than A, take the adjoint of A. This is allowed
+ *           In an optimal implementation, this trivial transpose
- *        only for square matrices, M=N.
+ *           should be replaced with faster transpose.
-          IF ( TRANSP ) THEN
+             DO 1115 p = 1, N - 1
- *           In an optimal implementation, this trivial transpose
+                A(p,p) = CONJG(A(p,p))
- *           should be replaced with faster transpose.
+                DO 1116 q = p + 1, N
-             DO 1115 p = 1, N - 1
+                    CTEMP = CONJG(A(q,p))
-                A(p,p) = CONJG(A(p,p))
+                   A(q,p) = CONJG(A(p,q))
-                DO 1116 q = p + 1, N
+                   A(p,q) = CTEMP
-                    CTEMP = CONJG(A(q,p))
+         CONTINUE
-                   A(q,p) = CONJG(A(p,q))
+      CONTINUE
-                   A(p,q) = CTEMP
+             A(N,N) = CONJG(A(N,N))
-         CONTINUE
+             DO 1117 p = 1, N
-      CONTINUE
+                RWORK(M+p) = SVA(p)
-             A(N,N) = CONJG(A(N,N))
+                SVA(p)     = RWORK(p)
-             DO 1117 p = 1, N
+ *              previously computed row 2-norms are now column 2-norms
-                RWORK(M+p) = SVA(p)
+ *              of the transposed matrix
-                SVA(p)     = RWORK(p)
+      CONTINUE
- *              previously computed row 2-norms are now column 2-norms
+             TEMP1  = AAPP
- *              of the transposed matrix
+             AAPP   = AATMAX
-      CONTINUE
+             AATMAX = TEMP1
-             TEMP1  = AAPP
+             TEMP1  = AAQQ
-             AAPP   = AATMAX
+             AAQQ   = AATMIN
-             AATMAX = TEMP1
+             AATMIN = TEMP1
-             TEMP1  = AAQQ
+             KILL   = LSVEC
-             AAQQ   = AATMIN
+             LSVEC  = RSVEC
-             AATMIN = TEMP1
+             RSVEC  = KILL
-             KILL   = LSVEC
+             IF ( LSVEC ) N1 = N
-             LSVEC  = RSVEC
+ *
-             RSVEC  = KILL
+             ROWPIV = .TRUE.
-             IF ( LSVEC ) N1 = N
+          END IF
  *
-             ROWPIV = .TRUE.
+       END IF
-          END IF
+ *     END IF L2TRAN
  *
-       END IF
+ *     Scale the matrix so that its maximal singular value remains less
- *     END IF L2TRAN
+ *     than SQRT(BIG) -- the matrix is scaled so that its maximal column
- *
+ *     has Euclidean norm equal to SQRT(BIG/N). The only reason to keep
- *     Scale the matrix so that its maximal singular value remains less
+ *     SQRT(BIG) instead of BIG is the fact that ZGEJSV uses LAPACK and
- *     than SQRT(BIG) -- the matrix is scaled so that its maximal column
+ *     BLAS routines that, in some implementations, are not capable of
- *     has Euclidean norm equal to SQRT(BIG/N). The only reason to keep
+ *     working in the full interval [SFMIN,BIG] and that they may provoke
- *     SQRT(BIG) instead of BIG is the fact that ZGEJSV uses LAPACK and
+ *     overflows in the intermediate results. If the singular values spread
- *     BLAS routines that, in some implementations, are not capable of
+ *     from SFMIN to BIG, then ZGESVJ will compute them. So, in that case,
- *     working in the full interval [SFMIN,BIG] and that they may provoke
+ *     one should use ZGESVJ instead of ZGEJSV.
- *     overflows in the intermediate results. If the singular values spread
+ *     >> change in the April 2016 update: allow bigger range, i.e. the
- *     from SFMIN to BIG, then ZGESVJ will compute them. So, in that case,
+ *     largest column is allowed up to BIG/N and ZGESVJ will do the rest.
- *     one should use ZGESVJ instead of ZGEJSV.
+       BIG1   = SQRT( BIG )
- *     >> change in the April 2016 update: allow bigger range, i.e. the
+       TEMP1  = SQRT( BIG / DBLE(N) )
- *     largest column is allowed up to BIG/N and ZGESVJ will do the rest.
+ *      TEMP1  = BIG/DBLE(N)
-       BIG1   = SQRT( BIG )
+ *
-       TEMP1  = SQRT( BIG / DBLE(N) )
+       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
- *      TEMP1  = BIG/DBLE(N)
+       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
- *
+           AAQQ = ( AAQQ / AAPP ) * TEMP1
-       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
+       ELSE
-       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
+           AAQQ = ( AAQQ * TEMP1 ) / AAPP
-           AAQQ = ( AAQQ / AAPP ) * TEMP1
+       END IF
-       ELSE
+       TEMP1 = TEMP1 * SCALEM
-           AAQQ = ( AAQQ * TEMP1 ) / AAPP
+       CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
-       END IF
+ *
-       TEMP1 = TEMP1 * SCALEM
+ *     To undo scaling at the end of this procedure, multiply the
-       CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
+ *     computed singular values with USCAL2 / USCAL1.
  *
- *     To undo scaling at the end of this procedure, multiply the
+       USCAL1 = TEMP1
- *     computed singular values with USCAL2 / USCAL1.
+       USCAL2 = AAPP
  *
-       USCAL1 = TEMP1
+       IF ( L2KILL ) THEN
-       USCAL2 = AAPP
+ *        L2KILL enforces computation of nonzero singular values in
- *
+ *        the restricted range of condition number of the initial A,
-       IF ( L2KILL ) THEN
+ *        sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN).
- *        L2KILL enforces computation of nonzero singular values in
+          XSC = SQRT( SFMIN )
- *        the restricted range of condition number of the initial A,
+       ELSE
- *        sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN).
+          XSC = SMALL
-          XSC = SQRT( SFMIN )
+ *
-       ELSE
+ *        Now, if the condition number of A is too big,
-          XSC = SMALL
+ *        sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN,
- *
+ *        as a precaution measure, the full SVD is computed using ZGESVJ
- *        Now, if the condition number of A is too big,
+ *        with accumulated Jacobi rotations. This provides numerically
- *        sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN,
+ *        more robust computation, at the cost of slightly increased run
- *        as a precaution measure, the full SVD is computed using ZGESVJ
+ *        time. Depending on the concrete implementation of BLAS and LAPACK
- *        with accumulated Jacobi rotations. This provides numerically
+ *        (i.e. how they behave in presence of extreme ill-conditioning) the
- *        more robust computation, at the cost of slightly increased run
+ *        implementor may decide to remove this switch.
- *        time. Depending on the concrete implementation of BLAS and LAPACK
+          IF ( ( AAQQ.LT.SQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
- *        (i.e. how they behave in presence of extreme ill-conditioning) the
+             JRACC = .TRUE.
- *        implementor may decide to remove this switch.
+          END IF
-          IF ( ( AAQQ.LT.SQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
+ *
-             JRACC = .TRUE.
+       END IF
-          END IF
+       IF ( AAQQ .LT. XSC ) THEN
- *
+          DO 700 p = 1, N
-       END IF
+             IF ( SVA(p) .LT. XSC ) THEN
-       IF ( AAQQ .LT. XSC ) THEN
+                CALL ZLASET( 'A', M, 1, CZERO, CZERO, A(1,p), LDA )
-          DO 700 p = 1, N
+                SVA(p) = ZERO
-             IF ( SVA(p) .LT. XSC ) THEN
+             END IF
-                CALL ZLASET( 'A', M, 1, CZERO, CZERO, A(1,p), LDA )
+    CONTINUE
-                SVA(p) = ZERO
+       END IF
-             END IF
+ *
-    CONTINUE
+ *     Preconditioning using QR factorization with pivoting
-       END IF
+ *
- *
+       IF ( ROWPIV ) THEN
- *     Preconditioning using QR factorization with pivoting
+ *        Optional row permutation (Bjoerck row pivoting):
- *
+ *        A result by Cox and Higham shows that the Bjoerck's
-       IF ( ROWPIV ) THEN
+ *        row pivoting combined with standard column pivoting
- *        Optional row permutation (Bjoerck row pivoting):
+ *        has similar effect as Powell-Reid complete pivoting.
- *        A result by Cox and Higham shows that the Bjoerck's
+ *        The ell-infinity norms of A are made nonincreasing.
- *        row pivoting combined with standard column pivoting
+          IF ( ( LSVEC .AND. RSVEC ) .AND. .NOT.( JRACC ) ) THEN
- *        has similar effect as Powell-Reid complete pivoting.
+               IWOFF = 2*N
- *        The ell-infinity norms of A are made nonincreasing.
+          ELSE
-          IF ( ( LSVEC .AND. RSVEC ) .AND. .NOT.( JRACC ) ) THEN
+               IWOFF = N
-               IWOFF = 2*N
+          END IF
-          ELSE
+          DO 1952 p = 1, M - 1
-               IWOFF = N
+             q = IDAMAX( M-p+1, RWORK(M+p), 1 ) + p - 1
-          END IF
+             IWORK(IWOFF+p) = q
-          DO 1952 p = 1, M - 1
+             IF ( p .NE. q ) THEN
-             q = IDAMAX( M-p+1, RWORK(M+p), 1 ) + p - 1
+                TEMP1      = RWORK(M+p)
-             IWORK(IWOFF+p) = q
+                RWORK(M+p) = RWORK(M+q)
-             IF ( p .NE. q ) THEN
+                RWORK(M+q) = TEMP1
-                TEMP1      = RWORK(M+p)
+             END IF
-                RWORK(M+p) = RWORK(M+q)
+   CONTINUE
-                RWORK(M+q) = TEMP1
+          CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(IWOFF+1), 1 )
-             END IF
+       END IF
-   CONTINUE
+ *
-          CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(IWOFF+1), 1 )
+ *     End of the preparation phase (scaling, optional sorting and
-       END IF
+ *     transposing, optional flushing of small columns).
  *
- *     End of the preparation phase (scaling, optional sorting and
+ *     Preconditioning
- *     transposing, optional flushing of small columns).
+ *
- *
+ *     If the full SVD is needed, the right singular vectors are computed
- *     Preconditioning
+ *     from a matrix equation, and for that we need theoretical analysis
- *
+ *     of the Businger-Golub pivoting. So we use ZGEQP3 as the first RR QRF.
- *     If the full SVD is needed, the right singular vectors are computed
+ *     In all other cases the first RR QRF can be chosen by other criteria
- *     from a matrix equation, and for that we need theoretical analysis
+ *     (eg speed by replacing global with restricted window pivoting, such
- *     of the Businger-Golub pivoting. So we use ZGEQP3 as the first RR QRF.
+ *     as in xGEQPX from TOMS # 782). Good results will be obtained using
- *     In all other cases the first RR QRF can be chosen by other criteria
+ *     xGEQPX with properly (!) chosen numerical parameters.
- *     (eg speed by replacing global with restricted window pivoting, such
+ *     Any improvement of ZGEQP3 improves overall performance of ZGEJSV.
- *     as in xGEQPX from TOMS # 782). Good results will be obtained using
+ *
- *     xGEQPX with properly (!) chosen numerical parameters.
+ *     A * P1 = Q1 * [ R1^* 0]^*:
- *     Any improvement of ZGEQP3 improves overal performance of ZGEJSV.
+       DO 1963 p = 1, N
- *
+ *        .. all columns are free columns
- *     A * P1 = Q1 * [ R1^* 0]^*:
+          IWORK(p) = 0
-       DO 1963 p = 1, N
+CONTINUE
- *        .. all columns are free columns
+       CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LWORK-N,
-          IWORK(p) = 0
+      $             RWORK, IERR )
-CONTINUE
+ *
-       CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LWORK-N,
+ *     The upper triangular matrix R1 from the first QRF is inspected for
-      $             RWORK, IERR )
+ *     rank deficiency and possibilities for deflation, or possible
- *
+ *     ill-conditioning. Depending on the user specified flag L2RANK,
- *     The upper triangular matrix R1 from the first QRF is inspected for
+ *     the procedure explores possibilities to reduce the numerical
- *     rank deficiency and possibilities for deflation, or possible
+ *     rank by inspecting the computed upper triangular factor. If
- *     ill-conditioning. Depending on the user specified flag L2RANK,
+ *     L2RANK or L2ABER are up, then ZGEJSV will compute the SVD of
- *     the procedure explores possibilities to reduce the numerical
+ *     A + dA, where ||dA|| <= f(M,N)*EPSLN.
- *     rank by inspecting the computed upper triangular factor. If
+ *
- *     L2RANK or L2ABER are up, then ZGEJSV will compute the SVD of
+       NR = 1
- *     A + dA, where ||dA|| <= f(M,N)*EPSLN.
+       IF ( L2ABER ) THEN
- *
+ *        Standard absolute error bound suffices. All sigma_i with
-       NR = 1
+ *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
-       IF ( L2ABER ) THEN
+ *        aggressive enforcement of lower numerical rank by introducing a
- *        Standard absolute error bound suffices. All sigma_i with
+ *        backward error of the order of N*EPSLN*||A||.
- *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
+          TEMP1 = SQRT(DBLE(N))*EPSLN
- *        agressive enforcement of lower numerical rank by introducing a
+          DO 3001 p = 2, N
- *        backward error of the order of N*EPSLN*||A||.
+             IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN
-          TEMP1 = SQRT(DBLE(N))*EPSLN
+                NR = NR + 1
-          DO 3001 p = 2, N
+             ELSE
-             IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN
+                GO TO 3002
-                NR = NR + 1
+             END IF
-             ELSE
+   CONTINUE
-                GO TO 3002
+   CONTINUE
-             END IF
+       ELSE IF ( L2RANK ) THEN
-   CONTINUE
+ *        .. similarly as above, only slightly more gentle (less aggressive).
-   CONTINUE
+ *        Sudden drop on the diagonal of R1 is used as the criterion for
-       ELSE IF ( L2RANK ) THEN
+ *        close-to-rank-deficient.
- *        .. similarly as above, only slightly more gentle (less agressive).
+          TEMP1 = SQRT(SFMIN)
- *        Sudden drop on the diagonal of R1 is used as the criterion for
+          DO 3401 p = 2, N
- *        close-to-rank-deficient.
+             IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
-          TEMP1 = SQRT(SFMIN)
+      $           ( ABS(A(p,p)) .LT. SMALL ) .OR.
-          DO 3401 p = 2, N
+      $           ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
-             IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
+             NR = NR + 1
-      $           ( ABS(A(p,p)) .LT. SMALL ) .OR.
+   CONTINUE
-      $           ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
+   CONTINUE
-             NR = NR + 1
+ *
-   CONTINUE
+       ELSE
-   CONTINUE
+ *        The goal is high relative accuracy. However, if the matrix
- *
+ *        has high scaled condition number the relative accuracy is in
-       ELSE
+ *        general not feasible. Later on, a condition number estimator
- *        The goal is high relative accuracy. However, if the matrix
+ *        will be deployed to estimate the scaled condition number.
- *        has high scaled condition number the relative accuracy is in
+ *        Here we just remove the underflowed part of the triangular
- *        general not feasible. Later on, a condition number estimator
+ *        factor. This prevents the situation in which the code is
- *        will be deployed to estimate the scaled condition number.
+ *        working hard to get the accuracy not warranted by the data.
- *        Here we just remove the underflowed part of the triangular
+          TEMP1  = SQRT(SFMIN)
- *        factor. This prevents the situation in which the code is
+          DO 3301 p = 2, N
- *        working hard to get the accuracy not warranted by the data.
+             IF ( ( ABS(A(p,p)) .LT. SMALL ) .OR.
-          TEMP1  = SQRT(SFMIN)
+      $           ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
-          DO 3301 p = 2, N
+             NR = NR + 1
-             IF ( ( ABS(A(p,p)) .LT. SMALL ) .OR.
+   CONTINUE
-      $           ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
+   CONTINUE
-             NR = NR + 1
+ *
-   CONTINUE
+       END IF
-   CONTINUE
+ *
- *
+       ALMORT = .FALSE.
-       END IF
+       IF ( NR .EQ. N ) THEN
- *
+          MAXPRJ = ONE
-       ALMORT = .FALSE.
+          DO 3051 p = 2, N
-       IF ( NR .EQ. N ) THEN
+             TEMP1  = ABS(A(p,p)) / SVA(IWORK(p))
-          MAXPRJ = ONE
+             MAXPRJ = MIN( MAXPRJ, TEMP1 )
-          DO 3051 p = 2, N
+   CONTINUE
-             TEMP1  = ABS(A(p,p)) / SVA(IWORK(p))
+          IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
-             MAXPRJ = MIN( MAXPRJ, TEMP1 )
+       END IF
-   CONTINUE
+ *
-          IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
+ *
-       END IF
+       SCONDA = - ONE
- *
+       CONDR1 = - ONE
- *
+       CONDR2 = - ONE
-       SCONDA = - ONE
+ *
-       CONDR1 = - ONE
+       IF ( ERREST ) THEN
-       CONDR2 = - ONE
+          IF ( N .EQ. NR ) THEN
- *
+             IF ( RSVEC ) THEN
-       IF ( ERREST ) THEN
+ *              .. V is available as workspace
-          IF ( N .EQ. NR ) THEN
+                CALL ZLACPY( 'U', N, N, A, LDA, V, LDV )
-             IF ( RSVEC ) THEN
+                DO 3053 p = 1, N
- *              .. V is available as workspace
+                   TEMP1 = SVA(IWORK(p))
-                CALL ZLACPY( 'U', N, N, A, LDA, V, LDV )
+                   CALL ZDSCAL( p, ONE/TEMP1, V(1,p), 1 )
-                DO 3053 p = 1, N
+         CONTINUE
-                   TEMP1 = SVA(IWORK(p))
+                IF ( LSVEC )THEN
-                   CALL ZDSCAL( p, ONE/TEMP1, V(1,p), 1 )
+                    CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1,
-         CONTINUE
+      $                  CWORK(N+1), RWORK, IERR )
-                IF ( LSVEC )THEN
+                ELSE
                     CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1,
-      $                  CWORK(N+1), RWORK, IERR )
+      $                  CWORK, RWORK, IERR )
-                ELSE
+                END IF
-                    CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1,
+ *
-      $                  CWORK, RWORK, IERR )
+             ELSE IF ( LSVEC ) THEN
-                END IF
+ *              .. U is available as workspace
- *
+                CALL ZLACPY( 'U', N, N, A, LDA, U, LDU )
-             ELSE IF ( LSVEC ) THEN
+                DO 3054 p = 1, N
- *              .. U is available as workspace
+                   TEMP1 = SVA(IWORK(p))
-                CALL ZLACPY( 'U', N, N, A, LDA, U, LDU )
+                   CALL ZDSCAL( p, ONE/TEMP1, U(1,p), 1 )
-                DO 3054 p = 1, N
+         CONTINUE
-                   TEMP1 = SVA(IWORK(p))
+                CALL ZPOCON( 'U', N, U, LDU, ONE, TEMP1,
-                   CALL ZDSCAL( p, ONE/TEMP1, U(1,p), 1 )
+      $              CWORK(N+1), RWORK, IERR )
-         CONTINUE
+             ELSE
-                CALL ZPOCON( 'U', N, U, LDU, ONE, TEMP1,
+                CALL ZLACPY( 'U', N, N, A, LDA, CWORK, N )
-      $              CWORK(N+1), RWORK, IERR )
+ *[]            CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
-             ELSE
+ *              Change: here index shifted by N to the left, CWORK(1:N)
-                CALL ZLACPY( 'U', N, N, A, LDA, CWORK, N )
+ *              not needed for SIGMA only computation
- *[]            CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
+                DO 3052 p = 1, N
- *              Change: here index shifted by N to the left, CWORK(1:N)
+                   TEMP1 = SVA(IWORK(p))
- *              not needed for SIGMA only computation
+ *[]               CALL ZDSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 )
-                DO 3052 p = 1, N
+                   CALL ZDSCAL( p, ONE/TEMP1, CWORK((p-1)*N+1), 1 )
-                   TEMP1 = SVA(IWORK(p))
+         CONTINUE
- *[]               CALL ZDSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 )
+ *           .. the columns of R are scaled to have unit Euclidean lengths.
-                   CALL ZDSCAL( p, ONE/TEMP1, CWORK((p-1)*N+1), 1 )
+ *[]               CALL ZPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1,
-         CONTINUE
+ *[]     $              CWORK(N+N*N+1), RWORK, IERR )
- *           .. the columns of R are scaled to have unit Euclidean lengths.
+                CALL ZPOCON( 'U', N, CWORK, N, ONE, TEMP1,
- *[]               CALL ZPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1,
+      $              CWORK(N*N+1), RWORK, IERR )
- *[]     $              CWORK(N+N*N+1), RWORK, IERR )
+ *
-                CALL ZPOCON( 'U', N, CWORK, N, ONE, TEMP1,
+             END IF
-      $              CWORK(N*N+1), RWORK, IERR )
+             IF ( TEMP1 .NE. ZERO ) THEN
- *
+                SCONDA = ONE / SQRT(TEMP1)
-             END IF
+             ELSE
-             IF ( TEMP1 .NE. ZERO ) THEN
+                SCONDA = - ONE
-                SCONDA = ONE / SQRT(TEMP1)
+             END IF
-             ELSE
+ *           SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
-                SCONDA = - ONE
+ *           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
-             END IF
+          ELSE
- *           SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
+             SCONDA = - ONE
- *           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
+          END IF
-          ELSE
+       END IF
-             SCONDA = - ONE
+ *
-          END IF
+       L2PERT = L2PERT .AND. ( ABS( A(1,1)/A(NR,NR) ) .GT. SQRT(BIG1) )
-       END IF
+ *     If there is no violent scaling, artificial perturbation is not needed.
  *
-       L2PERT = L2PERT .AND. ( ABS( A(1,1)/A(NR,NR) ) .GT. SQRT(BIG1) )
+ *     Phase 3:
- *     If there is no violent scaling, artificial perturbation is not needed.
+ *
- *
+       IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
- *     Phase 3:
+ *
- *
+ *         Singular Values only
-       IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
+ *
- *
+ *         .. transpose A(1:NR,1:N)
- *         Singular Values only
+          DO 1946 p = 1, MIN( N-1, NR )
- *
+             CALL ZCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
- *         .. transpose A(1:NR,1:N)
+             CALL ZLACGV( N-p+1, A(p,p), 1 )
-          DO 1946 p = 1, MIN( N-1, NR )
+   CONTINUE
-             CALL ZCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
+          IF ( NR .EQ. N ) A(N,N) = CONJG(A(N,N))
-             CALL ZLACGV( N-p+1, A(p,p), 1 )
+ *
-   CONTINUE
+ *        The following two DO-loops introduce small relative perturbation
-          IF ( NR .EQ. N ) A(N,N) = CONJG(A(N,N))
+ *        into the strict upper triangle of the lower triangular matrix.
- *
+ *        Small entries below the main diagonal are also changed.
- *        The following two DO-loops introduce small relative perturbation
+ *        This modification is useful if the computing environment does not
- *        into the strict upper triangle of the lower triangular matrix.
+ *        provide/allow FLUSH TO ZERO underflow, for it prevents many
- *        Small entries below the main diagonal are also changed.
+ *        annoying denormalized numbers in case of strongly scaled matrices.
- *        This modification is useful if the computing environment does not
+ *        The perturbation is structured so that it does not introduce any
- *        provide/allow FLUSH TO ZERO underflow, for it prevents many
+ *        new perturbation of the singular values, and it does not destroy
- *        annoying denormalized numbers in case of strongly scaled matrices.
+ *        the job done by the preconditioner.
- *        The perturbation is structured so that it does not introduce any
+ *        The licence for this perturbation is in the variable L2PERT, which
- *        new perturbation of the singular values, and it does not destroy
+ *        should be .FALSE. if FLUSH TO ZERO underflow is active.
- *        the job done by the preconditioner.
+ *
- *        The licence for this perturbation is in the variable L2PERT, which
+          IF ( .NOT. ALMORT ) THEN
- *        should be .FALSE. if FLUSH TO ZERO underflow is active.
+ *
- *
+             IF ( L2PERT ) THEN
-          IF ( .NOT. ALMORT ) THEN
+ *              XSC = SQRT(SMALL)
- *
+                XSC = EPSLN / DBLE(N)
-             IF ( L2PERT ) THEN
+                DO 4947 q = 1, NR
- *              XSC = SQRT(SMALL)
+                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
-                XSC = EPSLN / DBLE(N)
+                   DO 4949 p = 1, N
-                DO 4947 q = 1, NR
+                      IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
-                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
+      $                    .OR. ( p .LT. q ) )
-                   DO 4949 p = 1, N
+ *     $                     A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
-                      IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
+      $                     A(p,q) = CTEMP
-      $                    .OR. ( p .LT. q ) )
+            CONTINUE
- *     $                     A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
+         CONTINUE
-      $                     A(p,q) = CTEMP
+             ELSE
-            CONTINUE
+                CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, A(1,2),LDA )
-         CONTINUE
+             END IF
-             ELSE
+ *
-                CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, A(1,2),LDA )
+ *            .. second preconditioning using the QR factorization
-             END IF
+ *
- *
+             CALL ZGEQRF( N,NR, A,LDA, CWORK, CWORK(N+1),LWORK-N, IERR )
- *            .. second preconditioning using the QR factorization
+ *
- *
+ *           .. and transpose upper to lower triangular
-             CALL ZGEQRF( N,NR, A,LDA, CWORK, CWORK(N+1),LWORK-N, IERR )
+             DO 1948 p = 1, NR - 1
- *
+                CALL ZCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
- *           .. and transpose upper to lower triangular
+                CALL ZLACGV( NR-p+1, A(p,p), 1 )
-             DO 1948 p = 1, NR - 1
+      CONTINUE
-                CALL ZCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
+ *
-                CALL ZLACGV( NR-p+1, A(p,p), 1 )
+       END IF
-      CONTINUE
+ *
- *
+ *           Row-cyclic Jacobi SVD algorithm with column pivoting
-       END IF
+ *
- *
+ *           .. again some perturbation (a "background noise") is added
- *           Row-cyclic Jacobi SVD algorithm with column pivoting
+ *           to drown denormals
- *
+             IF ( L2PERT ) THEN
- *           .. again some perturbation (a "background noise") is added
+ *              XSC = SQRT(SMALL)
- *           to drown denormals
+                XSC = EPSLN / DBLE(N)
-             IF ( L2PERT ) THEN
+                DO 1947 q = 1, NR
- *              XSC = SQRT(SMALL)
+                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
-                XSC = EPSLN / DBLE(N)
+                   DO 1949 p = 1, NR
-                DO 1947 q = 1, NR
+                      IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
-                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
+      $                       .OR. ( p .LT. q ) )
-                   DO 1949 p = 1, NR
+ *     $                   A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
-                      IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
+      $                   A(p,q) = CTEMP
-      $                       .OR. ( p .LT. q ) )
+            CONTINUE
- *     $                   A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
+         CONTINUE
-      $                   A(p,q) = CTEMP
+             ELSE
-            CONTINUE
+                CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, A(1,2), LDA )
-         CONTINUE
+             END IF
-             ELSE
+ *
-                CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, A(1,2), LDA )
+ *           .. and one-sided Jacobi rotations are started on a lower
-             END IF
+ *           triangular matrix (plus perturbation which is ignored in
- *
+ *           the part which destroys triangular form (confusing?!))
- *           .. and one-sided Jacobi rotations are started on a lower
+ *
- *           triangular matrix (plus perturbation which is ignored in
+             CALL ZGESVJ( 'L', 'N', 'N', NR, NR, A, LDA, SVA,
- *           the part which destroys triangular form (confusing?!))
+      $                N, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
  *
-             CALL ZGESVJ( 'L', 'N', 'N', NR, NR, A, LDA, SVA,
+             SCALEM  = RWORK(1)
-      $                N, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
+             NUMRANK = NINT(RWORK(2))
  *
-             SCALEM  = RWORK(1)
+ *
-             NUMRANK = NINT(RWORK(2))
+       ELSE IF ( ( RSVEC .AND. ( .NOT. LSVEC ) .AND. ( .NOT. JRACC ) )
- *
+      $       .OR.
- *
+      $   ( JRACC .AND. ( .NOT. LSVEC ) .AND. ( NR .NE. N ) ) ) THEN
-       ELSE IF ( ( RSVEC .AND. ( .NOT. LSVEC ) .AND. ( .NOT. JRACC ) )
+ *
-      $       .OR.
+ *        -> Singular Values and Right Singular Vectors <-
-      $   ( JRACC .AND. ( .NOT. LSVEC ) .AND. ( NR .NE. N ) ) ) THEN
+ *
- *
+          IF ( ALMORT ) THEN
- *        -> Singular Values and Right Singular Vectors <-
+ *
- *
+ *           .. in this case NR equals N
-          IF ( ALMORT ) THEN
+             DO 1998 p = 1, NR
- *
+                CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
- *           .. in this case NR equals N
+                CALL ZLACGV( N-p+1, V(p,p), 1 )
-             DO 1998 p = 1, NR
+      CONTINUE
-                CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
+             CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
-                CALL ZLACGV( N-p+1, V(p,p), 1 )
+ *
-      CONTINUE
+             CALL ZGESVJ( 'L','U','N', N, NR, V, LDV, SVA, NR, A, LDA,
-             CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
+      $                  CWORK, LWORK, RWORK, LRWORK, INFO )
- *
+             SCALEM  = RWORK(1)
-             CALL ZGESVJ( 'L','U','N', N, NR, V, LDV, SVA, NR, A, LDA,
+             NUMRANK = NINT(RWORK(2))
-      $                  CWORK, LWORK, RWORK, LRWORK, INFO )
-             SCALEM  = RWORK(1)
+          ELSE
-             NUMRANK = NINT(RWORK(2))
+ *
+ *        .. two more QR factorizations ( one QRF is not enough, two require
-          ELSE
+ *        accumulated product of Jacobi rotations, three are perfect )
  *
- *        .. two more QR factorizations ( one QRF is not enough, two require
+             CALL ZLASET( 'L', NR-1,NR-1, CZERO, CZERO, A(2,1), LDA )
- *        accumulated product of Jacobi rotations, three are perfect )
+             CALL ZGELQF( NR,N, A, LDA, CWORK, CWORK(N+1), LWORK-N, IERR)
- *
+             CALL ZLACPY( 'L', NR, NR, A, LDA, V, LDV )
-             CALL ZLASET( 'L', NR-1,NR-1, CZERO, CZERO, A(2,1), LDA )
+             CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
-             CALL ZGELQF( NR,N, A, LDA, CWORK, CWORK(N+1), LWORK-N, IERR)
+             CALL ZGEQRF( NR, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
-             CALL ZLACPY( 'L', NR, NR, A, LDA, V, LDV )
+      $                   LWORK-2*N, IERR )
-             CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
+             DO 8998 p = 1, NR
-             CALL ZGEQRF( NR, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
+                CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
-      $                   LWORK-2*N, IERR )
+                CALL ZLACGV( NR-p+1, V(p,p), 1 )
-             DO 8998 p = 1, NR
+      CONTINUE
-                CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
+             CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV)
-                CALL ZLACGV( NR-p+1, V(p,p), 1 )
+ *
-      CONTINUE
+             CALL ZGESVJ( 'L', 'U','N', NR, NR, V,LDV, SVA, NR, U,
-             CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV)
+      $                  LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
- *
+             SCALEM  = RWORK(1)
-             CALL ZGESVJ( 'L', 'U','N', NR, NR, V,LDV, SVA, NR, U,
+             NUMRANK = NINT(RWORK(2))
-      $                  LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
+             IF ( NR .LT. N ) THEN
-             SCALEM  = RWORK(1)
+                CALL ZLASET( 'A',N-NR, NR, CZERO,CZERO, V(NR+1,1),  LDV )
-             NUMRANK = NINT(RWORK(2))
+                CALL ZLASET( 'A',NR, N-NR, CZERO,CZERO, V(1,NR+1),  LDV )
-             IF ( NR .LT. N ) THEN
+                CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE, V(NR+1,NR+1),LDV )
-                CALL ZLASET( 'A',N-NR, NR, CZERO,CZERO, V(NR+1,1),  LDV )
+             END IF
-                CALL ZLASET( 'A',NR, N-NR, CZERO,CZERO, V(1,NR+1),  LDV )
+ *
-                CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE, V(NR+1,NR+1),LDV )
+          CALL ZUNMLQ( 'L', 'C', N, N, NR, A, LDA, CWORK,
-             END IF
+      $               V, LDV, CWORK(N+1), LWORK-N, IERR )
  *
-          CALL ZUNMLQ( 'L', 'C', N, N, NR, A, LDA, CWORK,
+          END IF
-      $               V, LDV, CWORK(N+1), LWORK-N, IERR )
+ *         .. permute the rows of V
- *
+ *         DO 8991 p = 1, N
-          END IF
+ *            CALL ZCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
- *         .. permute the rows of V
+ * 8991    CONTINUE
- *         DO 8991 p = 1, N
+ *         CALL ZLACPY( 'All', N, N, A, LDA, V, LDV )
- *            CALL ZCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
+          CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
- * 8991    CONTINUE
+ *
- *         CALL ZLACPY( 'All', N, N, A, LDA, V, LDV )
+           IF ( TRANSP ) THEN
-          CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
+             CALL ZLACPY( 'A', N, N, V, LDV, U, LDU )
- *
+           END IF
-           IF ( TRANSP ) THEN
+ *
-             CALL ZLACPY( 'A', N, N, V, LDV, U, LDU )
+       ELSE IF ( JRACC .AND. (.NOT. LSVEC) .AND. ( NR.EQ. N ) ) THEN
-           END IF
+ *
- *
+          CALL ZLASET( 'L', N-1,N-1, CZERO, CZERO, A(2,1), LDA )
-       ELSE IF ( JRACC .AND. (.NOT. LSVEC) .AND. ( NR.EQ. N ) ) THEN
+ *
- *
+          CALL ZGESVJ( 'U','N','V', N, N, A, LDA, SVA, N, V, LDV,
-          CALL ZLASET( 'L', N-1,N-1, CZERO, CZERO, A(2,1), LDA )
+      $               CWORK, LWORK, RWORK, LRWORK, INFO )
- *
+           SCALEM  = RWORK(1)
-          CALL ZGESVJ( 'U','N','V', N, N, A, LDA, SVA, N, V, LDV,
+           NUMRANK = NINT(RWORK(2))
-      $               CWORK, LWORK, RWORK, LRWORK, INFO )
+           CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
-           SCALEM  = RWORK(1)
+ *
-           NUMRANK = NINT(RWORK(2))
+       ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
-           CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
+ *
- *
+ *        .. Singular Values and Left Singular Vectors                 ..
-       ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
+ *
- *
+ *        .. second preconditioning step to avoid need to accumulate
- *        .. Singular Values and Left Singular Vectors                 ..
+ *        Jacobi rotations in the Jacobi iterations.
- *
+          DO 1965 p = 1, NR
- *        .. second preconditioning step to avoid need to accumulate
+             CALL ZCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
- *        Jacobi rotations in the Jacobi iterations.
+             CALL ZLACGV( N-p+1, U(p,p), 1 )
-          DO 1965 p = 1, NR
+   CONTINUE
-             CALL ZCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
+          CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
-             CALL ZLACGV( N-p+1, U(p,p), 1 )
+ *
-   CONTINUE
+          CALL ZGEQRF( N, NR, U, LDU, CWORK(N+1), CWORK(2*N+1),
-          CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
+      $              LWORK-2*N, IERR )
  *
-          CALL ZGEQRF( N, NR, U, LDU, CWORK(N+1), CWORK(2*N+1),
+          DO 1967 p = 1, NR - 1
-      $              LWORK-2*N, IERR )
+             CALL ZCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
- *
+             CALL ZLACGV( N-p+1, U(p,p), 1 )
-          DO 1967 p = 1, NR - 1
+   CONTINUE
-             CALL ZCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
+          CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
-             CALL ZLACGV( N-p+1, U(p,p), 1 )
+ *
-   CONTINUE
+          CALL ZGESVJ( 'L', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
-          CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
+      $        LDA, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
- *
+          SCALEM  = RWORK(1)
-          CALL ZGESVJ( 'L', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
+          NUMRANK = NINT(RWORK(2))
-      $        LDA, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
+ *
-          SCALEM  = RWORK(1)
+          IF ( NR .LT. M ) THEN
-          NUMRANK = NINT(RWORK(2))
+             CALL ZLASET( 'A',  M-NR, NR,CZERO, CZERO, U(NR+1,1), LDU )
- *
+             IF ( NR .LT. N1 ) THEN
-          IF ( NR .LT. M ) THEN
+                CALL ZLASET( 'A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU )
-             CALL ZLASET( 'A',  M-NR, NR,CZERO, CZERO, U(NR+1,1), LDU )
+                CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,U(NR+1,NR+1),LDU )
-             IF ( NR .LT. N1 ) THEN
+             END IF
-                CALL ZLASET( 'A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU )
+          END IF
-                CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,U(NR+1,NR+1),LDU )
+ *
-             END IF
+          CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
-          END IF
+      $               LDU, CWORK(N+1), LWORK-N, IERR )
  *
-          CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
+          IF ( ROWPIV )
-      $               LDU, CWORK(N+1), LWORK-N, IERR )
+      $       CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
  *
-          IF ( ROWPIV )
+          DO 1974 p = 1, N1
-      $       CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
+             XSC = ONE / DZNRM2( M, U(1,p), 1 )
- *
+             CALL ZDSCAL( M, XSC, U(1,p), 1 )
-          DO 1974 p = 1, N1
+   CONTINUE
-             XSC = ONE / DZNRM2( M, U(1,p), 1 )
+ *
-             CALL ZDSCAL( M, XSC, U(1,p), 1 )
+          IF ( TRANSP ) THEN
-   CONTINUE
+             CALL ZLACPY( 'A', N, N, U, LDU, V, LDV )
- *
+          END IF
-          IF ( TRANSP ) THEN
+ *
-             CALL ZLACPY( 'A', N, N, U, LDU, V, LDV )
+       ELSE
-          END IF
+ *
- *
+ *        .. Full SVD ..
-       ELSE
+ *
- *
+          IF ( .NOT. JRACC ) THEN
- *        .. Full SVD ..
+ *
- *
+          IF ( .NOT. ALMORT ) THEN
-          IF ( .NOT. JRACC ) THEN
+ *
- *
+ *           Second Preconditioning Step (QRF [with pivoting])
-          IF ( .NOT. ALMORT ) THEN
+ *           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
- *
+ *           equivalent to an LQF CALL. Since in many libraries the QRF
- *           Second Preconditioning Step (QRF [with pivoting])
+ *           seems to be better optimized than the LQF, we do explicit
- *           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
+ *           transpose and use the QRF. This is subject to changes in an
- *           equivalent to an LQF CALL. Since in many libraries the QRF
+ *           optimized implementation of ZGEJSV.
- *           seems to be better optimized than the LQF, we do explicit
+ *
- *           transpose and use the QRF. This is subject to changes in an
+             DO 1968 p = 1, NR
- *           optimized implementation of ZGEJSV.
+                CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
- *
+                CALL ZLACGV( N-p+1, V(p,p), 1 )
-             DO 1968 p = 1, NR
+      CONTINUE
-                CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
+ *
-                CALL ZLACGV( N-p+1, V(p,p), 1 )
+ *           .. the following two loops perturb small entries to avoid
-      CONTINUE
+ *           denormals in the second QR factorization, where they are
- *
+ *           as good as zeros. This is done to avoid painfully slow
- *           .. the following two loops perturb small entries to avoid
+ *           computation with denormals. The relative size of the perturbation
- *           denormals in the second QR factorization, where they are
+ *           is a parameter that can be changed by the implementer.
- *           as good as zeros. This is done to avoid painfully slow
+ *           This perturbation device will be obsolete on machines with
- *           computation with denormals. The relative size of the perturbation
+ *           properly implemented arithmetic.
- *           is a parameter that can be changed by the implementer.
+ *           To switch it off, set L2PERT=.FALSE. To remove it from  the
- *           This perturbation device will be obsolete on machines with
+ *           code, remove the action under L2PERT=.TRUE., leave the ELSE part.
- *           properly implemented arithmetic.
+ *           The following two loops should be blocked and fused with the
- *           To switch it off, set L2PERT=.FALSE. To remove it from  the
+ *           transposed copy above.
- *           code, remove the action under L2PERT=.TRUE., leave the ELSE part.
+ *
- *           The following two loops should be blocked and fused with the
+             IF ( L2PERT ) THEN
- *           transposed copy above.
+                XSC = SQRT(SMALL)
- *
+                DO 2969 q = 1, NR
-             IF ( L2PERT ) THEN
+                   CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
-                XSC = SQRT(SMALL)
+                   DO 2968 p = 1, N
-                DO 2969 q = 1, NR
+                      IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
-                   CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
+      $                   .OR. ( p .LT. q ) )
-                   DO 2968 p = 1, N
+ *     $                   V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
-                      IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
+      $                   V(p,q) = CTEMP
-      $                   .OR. ( p .LT. q ) )
+                      IF ( p .LT. q ) V(p,q) = - V(p,q)
- *     $                   V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
+            CONTINUE
-      $                   V(p,q) = CTEMP
+         CONTINUE
-                      IF ( p .LT. q ) V(p,q) = - V(p,q)
+             ELSE
-            CONTINUE
+                CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
-         CONTINUE
+             END IF
-             ELSE
+ *
-                CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
+ *           Estimate the row scaled condition number of R1
-             END IF
+ *           (If R1 is rectangular, N > NR, then the condition number
- *
+ *           of the leading NR x NR submatrix is estimated.)
- *           Estimate the row scaled condition number of R1
+ *
- *           (If R1 is rectangular, N > NR, then the condition number
+             CALL ZLACPY( 'L', NR, NR, V, LDV, CWORK(2*N+1), NR )
- *           of the leading NR x NR submatrix is estimated.)
+             DO 3950 p = 1, NR
- *
+                TEMP1 = DZNRM2(NR-p+1,CWORK(2*N+(p-1)*NR+p),1)
-             CALL ZLACPY( 'L', NR, NR, V, LDV, CWORK(2*N+1), NR )
+                CALL ZDSCAL(NR-p+1,ONE/TEMP1,CWORK(2*N+(p-1)*NR+p),1)
-             DO 3950 p = 1, NR
+      CONTINUE
-                TEMP1 = DZNRM2(NR-p+1,CWORK(2*N+(p-1)*NR+p),1)
+             CALL ZPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1,
-                CALL ZDSCAL(NR-p+1,ONE/TEMP1,CWORK(2*N+(p-1)*NR+p),1)
+      $                   CWORK(2*N+NR*NR+1),RWORK,IERR)
-      CONTINUE
+             CONDR1 = ONE / SQRT(TEMP1)
-             CALL ZPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1,
+ *           .. here need a second opinion on the condition number
-      $                   CWORK(2*N+NR*NR+1),RWORK,IERR)
+ *           .. then assume worst case scenario
-             CONDR1 = ONE / SQRT(TEMP1)
+ *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
- *           .. here need a second oppinion on the condition number
+ *           more conservative    <=> CONDR1 .LT. SQRT(DBLE(N))
- *           .. then assume worst case scenario
+ *
- *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
+             COND_OK = SQRT(SQRT(DBLE(NR)))
- *           more conservative    <=> CONDR1 .LT. SQRT(DBLE(N))
+ *[TP]       COND_OK is a tuning parameter.
  *
-             COND_OK = SQRT(SQRT(DBLE(NR)))
+             IF ( CONDR1 .LT. COND_OK ) THEN
- *[TP]       COND_OK is a tuning parameter.
+ *              .. the second QRF without pivoting. Note: in an optimized
- *
+ *              implementation, this QRF should be implemented as the QRF
-             IF ( CONDR1 .LT. COND_OK ) THEN
+ *              of a lower triangular matrix.
- *              .. the second QRF without pivoting. Note: in an optimized
+ *              R1^* = Q2 * R2
- *              implementation, this QRF should be implemented as the QRF
+                CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
- *              of a lower triangular matrix.
+      $              LWORK-2*N, IERR )
- *              R1^* = Q2 * R2
+ *
-                CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
+                IF ( L2PERT ) THEN
-      $              LWORK-2*N, IERR )
+                   XSC = SQRT(SMALL)/EPSLN
- *
+                   DO 3959 p = 2, NR
-                IF ( L2PERT ) THEN
+                      DO 3958 q = 1, p - 1
-                   XSC = SQRT(SMALL)/EPSLN
+                         CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
-                   DO 3959 p = 2, NR
+      $                              ZERO)
-                      DO 3958 q = 1, p - 1
+                         IF ( ABS(V(q,p)) .LE. TEMP1 )
-                         CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
+ *     $                     V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
-      $                              ZERO)
+      $                     V(q,p) = CTEMP
-                         IF ( ABS(V(q,p)) .LE. TEMP1 )
+               CONTINUE
- *     $                     V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
+            CONTINUE
-      $                     V(q,p) = CTEMP
+                END IF
-               CONTINUE
+ *
-            CONTINUE
+                IF ( NR .NE. N )
-                END IF
+      $         CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
- *
+ *              .. save ...
-                IF ( NR .NE. N )
+ *
-      $         CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
+ *           .. this transposed copy should be better than naive
- *              .. save ...
+                DO 1969 p = 1, NR - 1
- *
+                   CALL ZCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
- *           .. this transposed copy should be better than naive
+                   CALL ZLACGV(NR-p+1, V(p,p), 1 )
-                DO 1969 p = 1, NR - 1
+         CONTINUE
-                   CALL ZCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
+                V(NR,NR)=CONJG(V(NR,NR))
-                   CALL ZLACGV(NR-p+1, V(p,p), 1 )
+ *
-         CONTINUE
+                CONDR2 = CONDR1
-                V(NR,NR)=CONJG(V(NR,NR))
+ *
- *
+             ELSE
-                CONDR2 = CONDR1
+ *
- *
+ *              .. ill-conditioned case: second QRF with pivoting
-             ELSE
+ *              Note that windowed pivoting would be equally good
- *
+ *              numerically, and more run-time efficient. So, in
- *              .. ill-conditioned case: second QRF with pivoting
+ *              an optimal implementation, the next call to ZGEQP3
- *              Note that windowed pivoting would be equaly good
+ *              should be replaced with eg. CALL ZGEQPX (ACM TOMS #782)
- *              numerically, and more run-time efficient. So, in
+ *              with properly (carefully) chosen parameters.
- *              an optimal implementation, the next call to ZGEQP3
+ *
- *              should be replaced with eg. CALL ZGEQPX (ACM TOMS #782)
+ *              R1^* * P2 = Q2 * R2
- *              with properly (carefully) chosen parameters.
+                DO 3003 p = 1, NR
- *
+                   IWORK(N+p) = 0
- *              R1^* * P2 = Q2 * R2
+         CONTINUE
-                DO 3003 p = 1, NR
+                CALL ZGEQP3( N, NR, V, LDV, IWORK(N+1), CWORK(N+1),
-                   IWORK(N+p) = 0
+      $                  CWORK(2*N+1), LWORK-2*N, RWORK, IERR )
-         CONTINUE
+ **               CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
-                CALL ZGEQP3( N, NR, V, LDV, IWORK(N+1), CWORK(N+1),
+ **     $              LWORK-2*N, IERR )
-      $                  CWORK(2*N+1), LWORK-2*N, RWORK, IERR )
+                IF ( L2PERT ) THEN
- **               CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
+                   XSC = SQRT(SMALL)
- **     $              LWORK-2*N, IERR )
+                   DO 3969 p = 2, NR
-                IF ( L2PERT ) THEN
+                      DO 3968 q = 1, p - 1
-                   XSC = SQRT(SMALL)
+                         CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
-                   DO 3969 p = 2, NR
+      $                                ZERO)
-                      DO 3968 q = 1, p - 1
+                         IF ( ABS(V(q,p)) .LE. TEMP1 )
-                         CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
+ *     $                     V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
-      $                                ZERO)
+      $                     V(q,p) = CTEMP
-                         IF ( ABS(V(q,p)) .LE. TEMP1 )
+               CONTINUE
- *     $                     V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
+            CONTINUE
-      $                     V(q,p) = CTEMP
+                END IF
-               CONTINUE
+ *
-            CONTINUE
+                CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
-                END IF
+ *
- *
+                IF ( L2PERT ) THEN
-                CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
+                   XSC = SQRT(SMALL)
- *
+                   DO 8970 p = 2, NR
-                IF ( L2PERT ) THEN
+                      DO 8971 q = 1, p - 1
-                   XSC = SQRT(SMALL)
+                         CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
-                   DO 8970 p = 2, NR
+      $                               ZERO)
-                      DO 8971 q = 1, p - 1
+ *                        V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) )
-                         CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
+                         V(p,q) = - CTEMP
-      $                               ZERO)
+               CONTINUE
- *                        V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) )
+            CONTINUE
-                         V(p,q) = - CTEMP
+                ELSE
-               CONTINUE
+                   CALL ZLASET( 'L',NR-1,NR-1,CZERO,CZERO,V(2,1),LDV )
-            CONTINUE
+                END IF
-                ELSE
+ *              Now, compute R2 = L3 * Q3, the LQ factorization.
-                   CALL ZLASET( 'L',NR-1,NR-1,CZERO,CZERO,V(2,1),LDV )
+                CALL ZGELQF( NR, NR, V, LDV, CWORK(2*N+N*NR+1),
-                END IF
+      $               CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
- *              Now, compute R2 = L3 * Q3, the LQ factorization.
+ *              .. and estimate the condition number
-                CALL ZGELQF( NR, NR, V, LDV, CWORK(2*N+N*NR+1),
+                CALL ZLACPY( 'L',NR,NR,V,LDV,CWORK(2*N+N*NR+NR+1),NR )
-      $               CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
+                DO 4950 p = 1, NR
- *              .. and estimate the condition number
+                   TEMP1 = DZNRM2( p, CWORK(2*N+N*NR+NR+p), NR )
-                CALL ZLACPY( 'L',NR,NR,V,LDV,CWORK(2*N+N*NR+NR+1),NR )
+                   CALL ZDSCAL( p, ONE/TEMP1, CWORK(2*N+N*NR+NR+p), NR )
-                DO 4950 p = 1, NR
+         CONTINUE
-                   TEMP1 = DZNRM2( p, CWORK(2*N+N*NR+NR+p), NR )
+                CALL ZPOCON( 'L',NR,CWORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
-                   CALL ZDSCAL( p, ONE/TEMP1, CWORK(2*N+N*NR+NR+p), NR )
+      $              CWORK(2*N+N*NR+NR+NR*NR+1),RWORK,IERR )
-         CONTINUE
+                CONDR2 = ONE / SQRT(TEMP1)
-                CALL ZPOCON( 'L',NR,CWORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
+ *
-      $              CWORK(2*N+N*NR+NR+NR*NR+1),RWORK,IERR )
+ *
-                CONDR2 = ONE / SQRT(TEMP1)
+                IF ( CONDR2 .GE. COND_OK ) THEN
- *
+ *                 .. save the Householder vectors used for Q3
- *
+ *                 (this overwrites the copy of R2, as it will not be
-                IF ( CONDR2 .GE. COND_OK ) THEN
+ *                 needed in this branch, but it does not overwritte the
- *                 .. save the Householder vectors used for Q3
+ *                 Huseholder vectors of Q2.).
- *                 (this overwrittes the copy of R2, as it will not be
+                   CALL ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N )
- *                 needed in this branch, but it does not overwritte the
+ *                 .. and the rest of the information on Q3 is in
- *                 Huseholder vectors of Q2.).
+ *                 WORK(2*N+N*NR+1:2*N+N*NR+N)
-                   CALL ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N )
+                END IF
- *                 .. and the rest of the information on Q3 is in
+ *
- *                 WORK(2*N+N*NR+1:2*N+N*NR+N)
+             END IF
-                END IF
+ *
- *
+             IF ( L2PERT ) THEN
-             END IF
+                XSC = SQRT(SMALL)
- *
+                DO 4968 q = 2, NR
-             IF ( L2PERT ) THEN
+                   CTEMP = XSC * V(q,q)
-                XSC = SQRT(SMALL)
+                   DO 4969 p = 1, q - 1
-                DO 4968 q = 2, NR
+ *                     V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) )
-                   CTEMP = XSC * V(q,q)
+                      V(p,q) = - CTEMP
-                   DO 4969 p = 1, q - 1
+            CONTINUE
- *                     V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) )
+         CONTINUE
-                      V(p,q) = - CTEMP
+             ELSE
-            CONTINUE
+                CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
-         CONTINUE
+             END IF
-             ELSE
+ *
-                CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
+ *        Second preconditioning finished; continue with Jacobi SVD
-             END IF
+ *        The input matrix is lower trinagular.
  *
- *        Second preconditioning finished; continue with Jacobi SVD
+ *        Recover the right singular vectors as solution of a well
- *        The input matrix is lower trinagular.
+ *        conditioned triangular matrix equation.
  *
- *        Recover the right singular vectors as solution of a well
+             IF ( CONDR1 .LT. COND_OK ) THEN
- *        conditioned triangular matrix equation.
+ *
- *
+                CALL ZGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, LDU,
-             IF ( CONDR1 .LT. COND_OK ) THEN
+      $              CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,RWORK,
- *
+      $              LRWORK, INFO )
-                CALL ZGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, LDU,
+                SCALEM  = RWORK(1)
-      $              CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,RWORK,
+                NUMRANK = NINT(RWORK(2))
-      $              LRWORK, INFO )
+                DO 3970 p = 1, NR
-                SCALEM  = RWORK(1)
+                   CALL ZCOPY(  NR, V(1,p), 1, U(1,p), 1 )
-                NUMRANK = NINT(RWORK(2))
+                   CALL ZDSCAL( NR, SVA(p),    V(1,p), 1 )
-                DO 3970 p = 1, NR
+         CONTINUE
-                   CALL ZCOPY(  NR, V(1,p), 1, U(1,p), 1 )
-                   CALL ZDSCAL( NR, SVA(p),    V(1,p), 1 )
+ *        .. pick the right matrix equation and solve it
-         CONTINUE
+ *
+                IF ( NR .EQ. N ) THEN
- *        .. pick the right matrix equation and solve it
+ * :))             .. best case, R1 is inverted. The solution of this matrix
- *
+ *                 equation is Q2*V2 = the product of the Jacobi rotations
-                IF ( NR .EQ. N ) THEN
+ *                 used in ZGESVJ, premultiplied with the orthogonal matrix
- * :))             .. best case, R1 is inverted. The solution of this matrix
+ *                 from the second QR factorization.
- *                 equation is Q2*V2 = the product of the Jacobi rotations
+                   CALL ZTRSM('L','U','N','N', NR,NR,CONE, A,LDA, V,LDV)
- *                 used in ZGESVJ, premultiplied with the orthogonal matrix
+                ELSE
- *                 from the second QR factorization.
+ *                 .. R1 is well conditioned, but non-square. Adjoint of R2
-                   CALL ZTRSM('L','U','N','N', NR,NR,CONE, A,LDA, V,LDV)
+ *                 is inverted to get the product of the Jacobi rotations
-                ELSE
+ *                 used in ZGESVJ. The Q-factor from the second QR
- *                 .. R1 is well conditioned, but non-square. Adjoint of R2
+ *                 factorization is then built in explicitly.
- *                 is inverted to get the product of the Jacobi rotations
+                   CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),
- *                 used in ZGESVJ. The Q-factor from the second QR
+      $                 N,V,LDV)
- *                 factorization is then built in explicitly.
+                   IF ( NR .LT. N ) THEN
-                   CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),
+                   CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
-      $                 N,V,LDV)
+                   CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
-                   IF ( NR .LT. N ) THEN
+                   CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
-                   CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
+                   END IF
-                   CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
+                   CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
-                   CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
+      $                V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
-                   END IF
+                END IF
-                   CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
+ *
-      $                V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
+             ELSE IF ( CONDR2 .LT. COND_OK ) THEN
-                END IF
+ *
- *
+ *              The matrix R2 is inverted. The solution of the matrix equation
-             ELSE IF ( CONDR2 .LT. COND_OK ) THEN
+ *              is Q3^* * V3 = the product of the Jacobi rotations (appplied to
- *
+ *              the lower triangular L3 from the LQ factorization of
- *              The matrix R2 is inverted. The solution of the matrix equation
+ *              R2=L3*Q3), pre-multiplied with the transposed Q3.
- *              is Q3^* * V3 = the product of the Jacobi rotations (appplied to
+                CALL ZGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
- *              the lower triangular L3 from the LQ factorization of
+      $              LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
- *              R2=L3*Q3), pre-multiplied with the transposed Q3.
+      $              RWORK, LRWORK, INFO )
-                CALL ZGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
+                SCALEM  = RWORK(1)
-      $              LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
+                NUMRANK = NINT(RWORK(2))
-      $              RWORK, LRWORK, INFO )
+                DO 3870 p = 1, NR
-                SCALEM  = RWORK(1)
+                   CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 )
-                NUMRANK = NINT(RWORK(2))
+                   CALL ZDSCAL( NR, SVA(p),    U(1,p), 1 )
-                DO 3870 p = 1, NR
+         CONTINUE
-                   CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 )
+                CALL ZTRSM('L','U','N','N',NR,NR,CONE,CWORK(2*N+1),N,
-                   CALL ZDSCAL( NR, SVA(p),    U(1,p), 1 )
+      $                    U,LDU)
-         CONTINUE
+ *              .. apply the permutation from the second QR factorization
-                CALL ZTRSM('L','U','N','N',NR,NR,CONE,CWORK(2*N+1),N,
+                DO 873 q = 1, NR
-      $                    U,LDU)
+                   DO 872 p = 1, NR
- *              .. apply the permutation from the second QR factorization
+                      CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
-                DO 873 q = 1, NR
+             CONTINUE
-                   DO 872 p = 1, NR
+                   DO 874 p = 1, NR
-                      CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
+                      U(p,q) = CWORK(2*N+N*NR+NR+p)
              CONTINUE
-                   DO 874 p = 1, NR
+          CONTINUE
-                      U(p,q) = CWORK(2*N+N*NR+NR+p)
+                IF ( NR .LT. N ) THEN
-             CONTINUE
+                   CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
-          CONTINUE
+                   CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
-                IF ( NR .LT. N ) THEN
+                   CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
-                   CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
+                END IF
-                   CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
+                CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
-                   CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
+      $              V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
-                END IF
+             ELSE
-                CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
+ *              Last line of defense.
-      $              V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
+ * #:(          This is a rather pathological case: no scaled condition
-             ELSE
+ *              improvement after two pivoted QR factorizations. Other
- *              Last line of defense.
+ *              possibility is that the rank revealing QR factorization
- * #:(          This is a rather pathological case: no scaled condition
+ *              or the condition estimator has failed, or the COND_OK
- *              improvement after two pivoted QR factorizations. Other
+ *              is set very close to ONE (which is unnecessary). Normally,
- *              possibility is that the rank revealing QR factorization
+ *              this branch should never be executed, but in rare cases of
- *              or the condition estimator has failed, or the COND_OK
+ *              failure of the RRQR or condition estimator, the last line of
- *              is set very close to ONE (which is unnecessary). Normally,
+ *              defense ensures that ZGEJSV completes the task.
- *              this branch should never be executed, but in rare cases of
+ *              Compute the full SVD of L3 using ZGESVJ with explicit
- *              failure of the RRQR or condition estimator, the last line of
+ *              accumulation of Jacobi rotations.
- *              defense ensures that ZGEJSV completes the task.
+                CALL ZGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
- *              Compute the full SVD of L3 using ZGESVJ with explicit
+      $              LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
- *              accumulation of Jacobi rotations.
+      $                         RWORK, LRWORK, INFO )
-                CALL ZGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
+                SCALEM  = RWORK(1)
-      $              LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
+                NUMRANK = NINT(RWORK(2))
-      $                         RWORK, LRWORK, INFO )
+                IF ( NR .LT. N ) THEN
-                SCALEM  = RWORK(1)
+                   CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
-                NUMRANK = NINT(RWORK(2))
+                   CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
-                IF ( NR .LT. N ) THEN
+                   CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
-                   CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
+                END IF
-                   CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
+                CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
-                   CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
+      $              V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
-                END IF
+ *
-                CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
+                CALL ZUNMLQ( 'L', 'C', NR, NR, NR, CWORK(2*N+1), N,
-      $              V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
+      $              CWORK(2*N+N*NR+1), U, LDU, CWORK(2*N+N*NR+NR+1),
- *
+      $              LWORK-2*N-N*NR-NR, IERR )
-                CALL ZUNMLQ( 'L', 'C', NR, NR, NR, CWORK(2*N+1), N,
+                DO 773 q = 1, NR
-      $              CWORK(2*N+N*NR+1), U, LDU, CWORK(2*N+N*NR+NR+1),
+                   DO 772 p = 1, NR
-      $              LWORK-2*N-N*NR-NR, IERR )
+                      CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
-                DO 773 q = 1, NR
+             CONTINUE
-                   DO 772 p = 1, NR
+                   DO 774 p = 1, NR
-                      CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
+                      U(p,q) = CWORK(2*N+N*NR+NR+p)
              CONTINUE
-                   DO 774 p = 1, NR
+          CONTINUE
-                      U(p,q) = CWORK(2*N+N*NR+NR+p)
+ *
-             CONTINUE
+             END IF
-          CONTINUE
+ *
- *
+ *           Permute the rows of V using the (column) permutation from the
-             END IF
+ *           first QRF. Also, scale the columns to make them unit in
- *
+ *           Euclidean norm. This applies to all cases.
- *           Permute the rows of V using the (column) permutation from the
+ *
- *           first QRF. Also, scale the columns to make them unit in
+             TEMP1 = SQRT(DBLE(N)) * EPSLN
- *           Euclidean norm. This applies to all cases.
+             DO 1972 q = 1, N
- *
+                DO 972 p = 1, N
-             TEMP1 = SQRT(DBLE(N)) * EPSLN
+                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
-             DO 1972 q = 1, N
+         CONTINUE
-                DO 972 p = 1, N
+                DO 973 p = 1, N
-                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
+                   V(p,q) = CWORK(2*N+N*NR+NR+p)
          CONTINUE
-                DO 973 p = 1, N
+                XSC = ONE / DZNRM2( N, V(1,q), 1 )
-                   V(p,q) = CWORK(2*N+N*NR+NR+p)
+                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
-         CONTINUE
+      $           CALL ZDSCAL( N, XSC, V(1,q), 1 )
-                XSC = ONE / DZNRM2( N, V(1,q), 1 )
+      CONTINUE
-                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
+ *           At this moment, V contains the right singular vectors of A.
-      $           CALL ZDSCAL( N, XSC, V(1,q), 1 )
+ *           Next, assemble the left singular vector matrix U (M x N).
-      CONTINUE
+             IF ( NR .LT. M ) THEN
- *           At this moment, V contains the right singular vectors of A.
+                CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU)
- *           Next, assemble the left singular vector matrix U (M x N).
+                IF ( NR .LT. N1 ) THEN
-             IF ( NR .LT. M ) THEN
+                   CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
-                CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU)
+                   CALL ZLASET('A',M-NR,N1-NR,CZERO,CONE,
-                IF ( NR .LT. N1 ) THEN
+      $                        U(NR+1,NR+1),LDU)
-                   CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
+                END IF
-                   CALL ZLASET('A',M-NR,N1-NR,CZERO,CONE,
+             END IF
-      $                        U(NR+1,NR+1),LDU)
+ *
-                END IF
+ *           The Q matrix from the first QRF is built into the left singular
-             END IF
+ *           matrix U. This applies to all cases.
  *
- *           The Q matrix from the first QRF is built into the left singular
+             CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
- *           matrix U. This applies to all cases.
+      $           LDU, CWORK(N+1), LWORK-N, IERR )
- *
-             CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
+ *           The columns of U are normalized. The cost is O(M*N) flops.
-      $           LDU, CWORK(N+1), LWORK-N, IERR )
+             TEMP1 = SQRT(DBLE(M)) * EPSLN
+             DO 1973 p = 1, NR
- *           The columns of U are normalized. The cost is O(M*N) flops.
+                XSC = ONE / DZNRM2( M, U(1,p), 1 )
-             TEMP1 = SQRT(DBLE(M)) * EPSLN
+                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
-             DO 1973 p = 1, NR
+      $          CALL ZDSCAL( M, XSC, U(1,p), 1 )
-                XSC = ONE / DZNRM2( M, U(1,p), 1 )
+      CONTINUE
-                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
+ *
-      $          CALL ZDSCAL( M, XSC, U(1,p), 1 )
+ *           If the initial QRF is computed with row pivoting, the left
-      CONTINUE
+ *           singular vectors must be adjusted.
  *
- *           If the initial QRF is computed with row pivoting, the left
+             IF ( ROWPIV )
- *           singular vectors must be adjusted.
+      $          CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
  *
-             IF ( ROWPIV )
+          ELSE
-      $          CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
+ *
- *
+ *        .. the initial matrix A has almost orthogonal columns and
-          ELSE
+ *        the second QRF is not needed
  *
- *        .. the initial matrix A has almost orthogonal columns and
+             CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
- *        the second QRF is not needed
+             IF ( L2PERT ) THEN
- *
+                XSC = SQRT(SMALL)
-             CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
+                DO 5970 p = 2, N
-             IF ( L2PERT ) THEN
+                   CTEMP = XSC * CWORK( N + (p-1)*N + p )
-                XSC = SQRT(SMALL)
+                   DO 5971 q = 1, p - 1
-                DO 5970 p = 2, N
+ *                     CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) /
-                   CTEMP = XSC * CWORK( N + (p-1)*N + p )
+ *     $                                        ABS(CWORK(N+(p-1)*N+q)) )
-                   DO 5971 q = 1, p - 1
+                      CWORK(N+(q-1)*N+p)=-CTEMP
- *                     CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) /
+            CONTINUE
- *     $                                        ABS(CWORK(N+(p-1)*N+q)) )
+         CONTINUE
-                      CWORK(N+(q-1)*N+p)=-CTEMP
+             ELSE
-            CONTINUE
+                CALL ZLASET( 'L',N-1,N-1,CZERO,CZERO,CWORK(N+2),N )
-         CONTINUE
+             END IF
-             ELSE
+ *
-                CALL ZLASET( 'L',N-1,N-1,CZERO,CZERO,CWORK(N+2),N )
+             CALL ZGESVJ( 'U', 'U', 'N', N, N, CWORK(N+1), N, SVA,
-             END IF
+      $           N, U, LDU, CWORK(N+N*N+1), LWORK-N-N*N, RWORK, LRWORK,
- *
+      $       INFO )
-             CALL ZGESVJ( 'U', 'U', 'N', N, N, CWORK(N+1), N, SVA,
+ *
-      $           N, U, LDU, CWORK(N+N*N+1), LWORK-N-N*N, RWORK, LRWORK,
+             SCALEM  = RWORK(1)
-      $       INFO )
+             NUMRANK = NINT(RWORK(2))
- *
+             DO 6970 p = 1, N
-             SCALEM  = RWORK(1)
+                CALL ZCOPY( N, CWORK(N+(p-1)*N+1), 1, U(1,p), 1 )
-             NUMRANK = NINT(RWORK(2))
+                CALL ZDSCAL( N, SVA(p), CWORK(N+(p-1)*N+1), 1 )
-             DO 6970 p = 1, N
+      CONTINUE
-                CALL ZCOPY( N, CWORK(N+(p-1)*N+1), 1, U(1,p), 1 )
+ *
-                CALL ZDSCAL( N, SVA(p), CWORK(N+(p-1)*N+1), 1 )
+             CALL ZTRSM( 'L', 'U', 'N', 'N', N, N,
-      CONTINUE
+      $           CONE, A, LDA, CWORK(N+1), N )
- *
+             DO 6972 p = 1, N
-             CALL ZTRSM( 'L', 'U', 'N', 'N', N, N,
+                CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV )
-      $           CONE, A, LDA, CWORK(N+1), N )
+      CONTINUE
-             DO 6972 p = 1, N
+             TEMP1 = SQRT(DBLE(N))*EPSLN
-                CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV )
+             DO 6971 p = 1, N
-      CONTINUE
+                XSC = ONE / DZNRM2( N, V(1,p), 1 )
-             TEMP1 = SQRT(DBLE(N))*EPSLN
+                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
-             DO 6971 p = 1, N
+      $            CALL ZDSCAL( N, XSC, V(1,p), 1 )
-                XSC = ONE / DZNRM2( N, V(1,p), 1 )
+      CONTINUE
-                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
+ *
-      $            CALL ZDSCAL( N, XSC, V(1,p), 1 )
+ *           Assemble the left singular vector matrix U (M x N).
-      CONTINUE
+ *
- *
+             IF ( N .LT. M ) THEN
- *           Assemble the left singular vector matrix U (M x N).
+                CALL ZLASET( 'A',  M-N, N, CZERO, CZERO, U(N+1,1), LDU )
- *
+                IF ( N .LT. N1 ) THEN
-             IF ( N .LT. M ) THEN
+                   CALL ZLASET('A',N,  N1-N, CZERO, CZERO,  U(1,N+1),LDU)
-                CALL ZLASET( 'A',  M-N, N, CZERO, CZERO, U(N+1,1), LDU )
+                   CALL ZLASET( 'A',M-N,N1-N, CZERO, CONE,U(N+1,N+1),LDU)
-                IF ( N .LT. N1 ) THEN
+                END IF
-                   CALL ZLASET('A',N,  N1-N, CZERO, CZERO,  U(1,N+1),LDU)
+             END IF
-                   CALL ZLASET( 'A',M-N,N1-N, CZERO, CONE,U(N+1,N+1),LDU)
+             CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
-                END IF
+      $           LDU, CWORK(N+1), LWORK-N, IERR )
-             END IF
+             TEMP1 = SQRT(DBLE(M))*EPSLN
-             CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
+             DO 6973 p = 1, N1
-      $           LDU, CWORK(N+1), LWORK-N, IERR )
+                XSC = ONE / DZNRM2( M, U(1,p), 1 )
-             TEMP1 = SQRT(DBLE(M))*EPSLN
+                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
-             DO 6973 p = 1, N1
+      $            CALL ZDSCAL( M, XSC, U(1,p), 1 )
-                XSC = ONE / DZNRM2( M, U(1,p), 1 )
+      CONTINUE
-                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
+ *
-      $            CALL ZDSCAL( M, XSC, U(1,p), 1 )
+             IF ( ROWPIV )
-      CONTINUE
+      $         CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
  *
-             IF ( ROWPIV )
+          END IF
-      $         CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
+ *
- *
+ *        end of the  >> almost orthogonal case <<  in the full SVD
-          END IF
+ *
- *
+          ELSE
- *        end of the  >> almost orthogonal case <<  in the full SVD
+ *
- *
+ *        This branch deploys a preconditioned Jacobi SVD with explicitly
-          ELSE
+ *        accumulated rotations. It is included as optional, mainly for
- *
+ *        experimental purposes. It does perform well, and can also be used.
- *        This branch deploys a preconditioned Jacobi SVD with explicitly
+ *        In this implementation, this branch will be automatically activated
- *        accumulated rotations. It is included as optional, mainly for
+ *        if the  condition number sigma_max(A) / sigma_min(A) is predicted
- *        experimental purposes. It does perfom well, and can also be used.
+ *        to be greater than the overflow threshold. This is because the
- *        In this implementation, this branch will be automatically activated
+ *        a posteriori computation of the singular vectors assumes robust
- *        if the  condition number sigma_max(A) / sigma_min(A) is predicted
+ *        implementation of BLAS and some LAPACK procedures, capable of working
- *        to be greater than the overflow threshold. This is because the
+ *        in presence of extreme values, e.g. when the singular values spread from
- *        a posteriori computation of the singular vectors assumes robust
+ *        the underflow to the overflow threshold.
- *        implementation of BLAS and some LAPACK procedures, capable of working
+ *
- *        in presence of extreme values, e.g. when the singular values spread from
+          DO 7968 p = 1, NR
- *        the underflow to the overflow threshold.
+             CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
- *
+             CALL ZLACGV( N-p+1, V(p,p), 1 )
-          DO 7968 p = 1, NR
+   CONTINUE
-             CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
+ *
-             CALL ZLACGV( N-p+1, V(p,p), 1 )
+          IF ( L2PERT ) THEN
-   CONTINUE
+             XSC = SQRT(SMALL/EPSLN)
- *
+             DO 5969 q = 1, NR
-          IF ( L2PERT ) THEN
+                CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
-             XSC = SQRT(SMALL/EPSLN)
+                DO 5968 p = 1, N
-             DO 5969 q = 1, NR
+                   IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
-                CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
+      $                .OR. ( p .LT. q ) )
-                DO 5968 p = 1, N
+ *     $                V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
-                   IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
+      $                V(p,q) = CTEMP
-      $                .OR. ( p .LT. q ) )
+                   IF ( p .LT. q ) V(p,q) = - V(p,q)
- *     $                V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
+         CONTINUE
-      $                V(p,q) = CTEMP
+      CONTINUE
-                   IF ( p .LT. q ) V(p,q) = - V(p,q)
+          ELSE
-         CONTINUE
+             CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
-      CONTINUE
+          END IF
-          ELSE
-             CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
+          CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
-          END IF
+      $        LWORK-2*N, IERR )
+          CALL ZLACPY( 'L', N, NR, V, LDV, CWORK(2*N+1), N )
-          CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
+ *
-      $        LWORK-2*N, IERR )
+          DO 7969 p = 1, NR
-          CALL ZLACPY( 'L', N, NR, V, LDV, CWORK(2*N+1), N )
+             CALL ZCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
- *
+             CALL ZLACGV( NR-p+1, U(p,p), 1 )
-          DO 7969 p = 1, NR
+   CONTINUE
-             CALL ZCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
-             CALL ZLACGV( NR-p+1, U(p,p), 1 )
+          IF ( L2PERT ) THEN
-   CONTINUE
+             XSC = SQRT(SMALL/EPSLN)
+             DO 9970 q = 2, NR
-          IF ( L2PERT ) THEN
+                DO 9971 p = 1, q - 1
-             XSC = SQRT(SMALL/EPSLN)
+                   CTEMP = DCMPLX(XSC * MIN(ABS(U(p,p)),ABS(U(q,q))),
-             DO 9970 q = 2, NR
+      $                            ZERO)
-                DO 9971 p = 1, q - 1
+ *                  U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) )
-                   CTEMP = DCMPLX(XSC * MIN(ABS(U(p,p)),ABS(U(q,q))),
+                   U(p,q) = - CTEMP
-      $                            ZERO)
+         CONTINUE
- *                  U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) )
+      CONTINUE
-                   U(p,q) = - CTEMP
+          ELSE
-         CONTINUE
+             CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
-      CONTINUE
+          END IF
-          ELSE
-             CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
+          CALL ZGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA,
-          END IF
+      $        N, V, LDV, CWORK(2*N+N*NR+1), LWORK-2*N-N*NR,
+      $         RWORK, LRWORK, INFO )
-          CALL ZGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA,
+          SCALEM  = RWORK(1)
-      $        N, V, LDV, CWORK(2*N+N*NR+1), LWORK-2*N-N*NR,
+          NUMRANK = NINT(RWORK(2))
-      $         RWORK, LRWORK, INFO )
-          SCALEM  = RWORK(1)
+          IF ( NR .LT. N ) THEN
-          NUMRANK = NINT(RWORK(2))
+             CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
+             CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
-          IF ( NR .LT. N ) THEN
+             CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV )
-             CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
+          END IF
-             CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
-             CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV )
+          CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
-          END IF
+      $        V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
+ *
-          CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
+ *           Permute the rows of V using the (column) permutation from the
-      $        V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
+ *           first QRF. Also, scale the columns to make them unit in
- *
+ *           Euclidean norm. This applies to all cases.
- *           Permute the rows of V using the (column) permutation from the
+ *
- *           first QRF. Also, scale the columns to make them unit in
+             TEMP1 = SQRT(DBLE(N)) * EPSLN
- *           Euclidean norm. This applies to all cases.
+             DO 7972 q = 1, N
- *
+                DO 8972 p = 1, N
-             TEMP1 = SQRT(DBLE(N)) * EPSLN
+                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
-             DO 7972 q = 1, N
+         CONTINUE
-                DO 8972 p = 1, N
+                DO 8973 p = 1, N
-                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
+                   V(p,q) = CWORK(2*N+N*NR+NR+p)
          CONTINUE
-                DO 8973 p = 1, N
+                XSC = ONE / DZNRM2( N, V(1,q), 1 )
-                   V(p,q) = CWORK(2*N+N*NR+NR+p)
+                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
-         CONTINUE
+      $           CALL ZDSCAL( N, XSC, V(1,q), 1 )
-                XSC = ONE / DZNRM2( N, V(1,q), 1 )
+      CONTINUE
-                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
+ *
-      $           CALL ZDSCAL( N, XSC, V(1,q), 1 )
+ *           At this moment, V contains the right singular vectors of A.
-      CONTINUE
+ *           Next, assemble the left singular vector matrix U (M x N).
  *
- *           At this moment, V contains the right singular vectors of A.
+          IF ( NR .LT. M ) THEN
- *           Next, assemble the left singular vector matrix U (M x N).
+             CALL ZLASET( 'A',  M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU )
- *
+             IF ( NR .LT. N1 ) THEN
-          IF ( NR .LT. M ) THEN
+                CALL ZLASET('A',NR,  N1-NR, CZERO, CZERO,  U(1,NR+1),LDU)
-             CALL ZLASET( 'A',  M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU )
+                CALL ZLASET('A',M-NR,N1-NR, CZERO, CONE,U(NR+1,NR+1),LDU)
-             IF ( NR .LT. N1 ) THEN
+             END IF
-                CALL ZLASET('A',NR,  N1-NR, CZERO, CZERO,  U(1,NR+1),LDU)
+          END IF
-                CALL ZLASET('A',M-NR,N1-NR, CZERO, CONE,U(NR+1,NR+1),LDU)
+ *
-             END IF
+          CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
-          END IF
+      $        LDU, CWORK(N+1), LWORK-N, IERR )
  *
-          CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
+             IF ( ROWPIV )
-      $        LDU, CWORK(N+1), LWORK-N, IERR )
+      $         CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
  *
-             IF ( ROWPIV )
+ *
-      $         CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
+          END IF
- *
+          IF ( TRANSP ) THEN
- *
+ *           .. swap U and V because the procedure worked on A^*
-          END IF
+             DO 6974 p = 1, N
-          IF ( TRANSP ) THEN
+                CALL ZSWAP( N, U(1,p), 1, V(1,p), 1 )
- *           .. swap U and V because the procedure worked on A^*
+      CONTINUE
-             DO 6974 p = 1, N
+          END IF
-                CALL ZSWAP( N, U(1,p), 1, V(1,p), 1 )
+ *
-      CONTINUE
+       END IF
-          END IF
+ *     end of the full SVD
  *
-       END IF
+ *     Undo scaling, if necessary (and possible)
- *     end of the full SVD
+ *
- *
+       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
- *     Undo scaling, if necessary (and possible)
+          CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
- *
+          USCAL1 = ONE
-       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
+          USCAL2 = ONE
-          CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
+       END IF
-          USCAL1 = ONE
+ *
-          USCAL2 = ONE
+       IF ( NR .LT. N ) THEN
-       END IF
+          DO 3004 p = NR+1, N
- *
+             SVA(p) = ZERO
-       IF ( NR .LT. N ) THEN
+   CONTINUE
-          DO 3004 p = NR+1, N
+       END IF
-             SVA(p) = ZERO
+ *
-   CONTINUE
+       RWORK(1) = USCAL2 * SCALEM
-       END IF
+       RWORK(2) = USCAL1
- *
+       IF ( ERREST ) RWORK(3) = SCONDA
-       RWORK(1) = USCAL2 * SCALEM
+       IF ( LSVEC .AND. RSVEC ) THEN
-       RWORK(2) = USCAL1
+          RWORK(4) = CONDR1
-       IF ( ERREST ) RWORK(3) = SCONDA
+          RWORK(5) = CONDR2
-       IF ( LSVEC .AND. RSVEC ) THEN
+       END IF
-          RWORK(4) = CONDR1
+       IF ( L2TRAN ) THEN
-          RWORK(5) = CONDR2
+          RWORK(6) = ENTRA
-       END IF
+          RWORK(7) = ENTRAT
-       IF ( L2TRAN ) THEN
+       END IF
-          RWORK(6) = ENTRA
+ *
-          RWORK(7) = ENTRAT
+       IWORK(1) = NR
-       END IF
+       IWORK(2) = NUMRANK
- *
+       IWORK(3) = WARNING
-       IWORK(1) = NR
+       IF ( TRANSP ) THEN
-       IWORK(2) = NUMRANK
+           IWORK(4) =  1
-       IWORK(3) = WARNING
+       ELSE
-       IF ( TRANSP ) THEN
+           IWORK(4) = -1
-           IWORK(4) =  1
+       END IF
-       ELSE
-           IWORK(4) = -1
+ *
-       END IF
+       RETURN
+ *     ..
- *
+ *     .. END OF ZGEJSV
-       RETURN
+ *     ..
- *     ..
+       END
- *     .. END OF ZGEJSV
+ *
- *     ..
-       END
- *

CVSweb interface <joel.bertrand@systella.fr>

Removed from v.1.4
changed lines
	Added in v.1.9