dgejsv.f

Go to the documentation of this file.

*> \brief \b DGEJSV
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DGEJSV + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgejsv.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgejsv.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgejsv.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
*                          M, N, A, LDA, SVA, U, LDU, V, LDV,
*                          WORK, LWORK, IWORK, INFO )
*
*       .. Scalar Arguments ..
*       IMPLICIT    NONE
*       INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
*      $            WORK( LWORK )
*       INTEGER     IWORK( * )
*       CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGEJSV computes the singular value decomposition (SVD) of a real M-by-N
*> matrix [A], where M >= N. The SVD of [A] is written as
*>
*>              [A] = [U] * [SIGMA] * [V]^t,
*>
*> 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) orthonormal matrix, and
*> [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are
*> the singular values of [A]. The columns of [U] and [V] are the left and
*> the right singular vectors of [A], respectively. 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.
*> DGEJSV can sometimes compute tiny singular values and their singular vectors much
*> more accurately than other SVD routines, see below under Further Details.
*> \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
*>             condition number of B, where A=D*B. If A has heavily weighted
*>             rows, then using this condition number gives too pessimistic
*>             error bound.
*>       = 'A': Small singular values are the noise and the matrix is treated
*>             as numerically rank deficient. The error in the computed
*>             singular values is bounded by f(m,n)*epsilon*||A||.
*>             The computed SVD A = U * S * V^t 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. This option is
*>             allowed only if JOBU .NE. 'N', i.e. in computing the full SVD.
*>       = '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 DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
*>        the licence to kill columns of A whose norm in c*A is less than
*>        DSQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
*>        where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
*>       = 'N': Do not kill small columns of c*A. This option assumes that
*>             BLAS and QR factorizations and triangular solvers are
*>             implemented to work in that range. If the condition of A
*>             is greater than BIG, use DGESVJ.
*>       = 'R': RESTRICTED range for sigma(c*A) is [DSQRT(SFMIN), DSQRT(BIG)]
*>             (roughly, as described above). This option is recommended.
*>                                            ~~~~~~~~~~~~~~~~~~~~~~~~~~~
*>        For computing the singular values in the FULL range [SFMIN,BIG]
*>        use DGESVJ.
*> \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^t seems to be better with respect to convergence.
*>        If the matrix is not square, JOBT is ignored. This is subject to
*>        changes in the future.
*>        The decision is based on two values of entropy over the adjoint
*>        orbit of A^t * A. See the descriptions of WORK(6) and WORK(7).
*>       = 'T': transpose if entropy test indicates possibly faster
*>        convergence of Jacobi process if A^t is taken as input. If A is
*>        replaced with A^t, then the row pivoting is included automatically.
*>       = 'N': do not speculate.
*>        This option can be used to compute only the singular values, or the
*>        full SVD (U, SIGMA and V). For only one set of singular vectors
*>        (U or V), the caller should provide both U and V, as one of the
*>        matrices is used as workspace if the matrix A is transposed.
*>        The implementer can easily remove this constraint and make the
*>        code more complicated. See the descriptions of U and V.
*> \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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( LDU, N ) or ( LDU, M )
*>          If JOBU = 'U', then U contains on exit the M-by-N matrix of
*>                         the left singular vectors.
*>          If JOBU = 'F', then U contains on exit the M-by-M matrix of
*>                         the left singular vectors, including an ONB
*>                         of the orthogonal complement of the Range(A).
*>          If JOBU = 'W'  .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
*>                         then U is used as workspace if the procedure
*>                         replaces A with A^t. In that case, [V] is computed
*>                         in U as left singular vectors of A^t 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 DOUBLE PRECISION array, dimension ( LDV, N )
*>          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
*>                         the right singular vectors;
*>          If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
*>                         then V is used as workspace if the procedure
*>                         replaces A with A^t. In that case, [U] is computed
*>                         in V as right singular vectors of A^t 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] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (MAX(7,LWORK))
*>          On exit, if N > 0 .AND. M > 0 (else not referenced),
*>          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
*>                    that SCALE*SVA(1:N) are the computed singular values
*>                    of A. (See the description of SVA().)
*>          WORK(2) = See the description of WORK(1).
*>          WORK(3) = SCONDA is an estimate for the condition number of
*>                    column equilibrated A. (If JOBA = 'E' or 'G')
*>                    SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
*>                    It is computed using DPOCON. It holds
*>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
*>                    where R is the triangular factor from the QRF of A.
*>                    However, if R is truncated and the numerical rank is
*>                    determined to be strictly smaller than N, SCONDA is
*>                    returned as -1, thus indicating that the smallest
*>                    singular values might be lost.
*>
*>          If full SVD is needed, the following two condition numbers are
*>          useful for the analysis of the algorithm. They are provided for
*>          a developer/implementer who is familiar with the details of
*>          the method.
*>
*>          WORK(4) = an estimate of the scaled condition number of the
*>                    triangular factor in the first QR factorization.
*>          WORK(5) = an estimate of the scaled condition number of the
*>                    triangular factor in the second QR factorization.
*>          The following two parameters are computed if JOBT = 'T'.
*>          They are provided for a developer/implementer who is familiar
*>          with the details of the method.
*>
*>          WORK(6) = the entropy of A^t*A :: this is the Shannon entropy
*>                    of diag(A^t*A) / Trace(A^t*A) taken as point in the
*>                    probability simplex.
*>          WORK(7) = the entropy of A*A^t.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          Length of WORK to confirm proper allocation of work space.
*>          LWORK depends on the job:
*>
*>          If only SIGMA is needed (JOBU = 'N', JOBV = 'N') and
*>            -> .. no scaled condition estimate required (JOBE = 'N'):
*>               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
*>               ->> For optimal performance (blocked code) the optimal value
*>               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
*>               block size for DGEQP3 and DGEQRF.
*>               In general, optimal LWORK is computed as
*>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7).
*>            -> .. an estimate of the scaled condition number of A is
*>               required (JOBA='E', 'G'). In this case, LWORK is the maximum
*>               of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7).
*>               ->> For optimal performance (blocked code) the optimal value
*>               is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7).
*>               In general, the optimal length LWORK is computed as
*>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
*>                                                     N+N*N+LWORK(DPOCON),7).
*>
*>          If SIGMA and the right singular vectors are needed (JOBV = 'V'),
*>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
*>            -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
*>               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF,
*>               DORMLQ. In general, the optimal length LWORK is computed as
*>               LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON),
*>                       N+LWORK(DGELQF), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
*>
*>          If SIGMA and the left singular vectors are needed
*>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
*>            -> For optimal performance:
*>               if JOBU = 'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
*>               if JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
*>               where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
*>               In general, the optimal length LWORK is computed as
*>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
*>                        2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
*>               Here LWORK(DORMQR) equals N*NB (for JOBU = 'U') or
*>               M*NB (for JOBU = 'F').
*>
*>          If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
*>            -> if JOBV = 'V'
*>               the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
*>            -> if JOBV = 'J' the minimal requirement is
*>               LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
*>            -> For optimal performance, LWORK should be additionally
*>               larger than N+M*NB, where NB is the optimal block size
*>               for DORMQR.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (MAX(3,M+3*N)).
*>          On exit,
*>          IWORK(1) = the numerical rank determined after the initial
*>                     QR factorization with pivoting. See the descriptions
*>                     of JOBA and JOBR.
*>          IWORK(2) = the number of the computed nonzero singular values
*>          IWORK(3) = if nonzero, a warning message:
*>                     If IWORK(3) = 1 then some of the column norms of A
*>                     were denormalized floats. The requested high accuracy
*>                     is not warranted by the data.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>           < 0:  if INFO = -i, then the i-th argument had an illegal value.
*>           = 0:  successful exit;
*>           > 0:  DGEJSV  did not converge in the maximal allowed number
*>                 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.
*
*> \ingroup gejsv
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses DGEQP3,
*>  DGEQRF, and DGELQF as preprocessors and preconditioners. Optionally, an
*>  additional row pivoting can be used as a preprocessor, which in some
*>  cases results in much higher accuracy. An example is matrix A with the
*>  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
*>  diagonal matrices and C is well-conditioned matrix. In that case, complete
*>  pivoting in the first QR factorizations provides accuracy dependent on the
*>  condition number of C, and independent of D1, D2. Such higher accuracy is
*>  not completely understood theoretically, but it works well in practice.
*>  Further, if A can be written as A = B*D, with well-conditioned B and some
*>  diagonal D, then the high accuracy is guaranteed, both theoretically and
*>  in software, independent of D. For more details see [1], [2].
*>     The computational range for the singular values can be the full range
*>  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
*>  & LAPACK routines called by DGEJSV are implemented to work in that range.
*>  If that is not the case, then the restriction for safe computation with
*>  the singular values in the range of normalized IEEE numbers is that the
*>  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
*>  overflow. This code (DGEJSV) is best used in this restricted range,
*>  meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
*>  returned as zeros. See JOBR for details on this.
*>     Further, this implementation is somewhat slower than the one described
*>  in [1,2] due to replacement of some non-LAPACK components, and because
*>  the choice of some tuning parameters in the iterative part (DGESVJ) is
*>  left to the implementer on a particular machine.
*>     The rank revealing QR factorization (in this code: DGEQP3) should be
*>  implemented as in [3]. We have a new version of DGEQP3 under development
*>  that is more robust than the current one in LAPACK, with a cleaner cut in
*>  rank deficient cases. It will be available in the SIGMA library [4].
*>  If M is much larger than N, it is obvious that the initial QRF with
*>  column pivoting can be preprocessed by the QRF without pivoting. That
*>  well known trick is not used in DGEJSV because in some cases heavy row
*>  weighting can be treated with complete pivoting. The overhead in cases
*>  M much larger than N is then only due to pivoting, but the benefits in
*>  terms of accuracy have prevailed. The implementer/user can incorporate
*>  this extra QRF step easily. The implementer can also improve data movement
*>  (matrix transpose, matrix copy, matrix transposed copy) - this
*>  implementation of DGEJSV uses only the simplest, naive data movement.
*> \endverbatim
*
*> \par Contributors:
*  ==================
*>
*>  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
*
*> \par References:
*  ================
*>
*> \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.
*>     LAPACK Working note 169.
*> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
*>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
*>     LAPACK Working note 170.
*> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
*>     factorization software - a case study.
*>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
*>     LAPACK Working note 176.
*> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
*>     QSVD, (H,K)-SVD computations.
*>     Department of Mathematics, University of Zagreb, 2008.
*> \endverbatim
*
*>  \par Bugs, examples and comments:
*   =================================
*>
*>  Please report all bugs and send interesting examples and/or comments to
*>  drmac@math.hr. Thank you.
*>
*  =====================================================================
      SUBROUTINE dgejsv( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
     $                   M, N, A, LDA, SVA, U, LDU, V, LDV,
     $                   WORK, LWORK, IWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      IMPLICIT    NONE
      INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
     $            WORK( LWORK )
      INTEGER     IWORK( * )
      CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
*     ..
*
*  ===========================================================================
*
*     .. Local Parameters ..
      DOUBLE PRECISION   ZERO,  ONE
      PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
     $        CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,  MAXPRJ, SCALEM,
     $        sconda, sfmin,  small,  temp1,  uscal1, uscal2, xsc
      INTEGER IERR,   N1,     NR,     NUMRANK,        p, q,   WARNING
      LOGICAL ALMORT, DEFR,   ERREST, GOSCAL, JRACC,  KILL,   LSVEC,
     $        l2aber, l2kill, l2pert, l2rank, l2tran,
     $        noscal, rowpiv, rsvec,  transp
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC dabs, dlog, max, min, dble, idnint, dsign, dsqrt
*     ..
*     .. External Functions ..
      DOUBLE PRECISION  DLAMCH, DNRM2
      INTEGER   IDAMAX
      LOGICAL   LSAME
      EXTERNAL  idamax, lsame, dlamch, dnrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL  dcopy,  dgelqf, dgeqp3, dgeqrf, dlacpy, dlascl,
     $          dlaset, dlassq, dlaswp, dorgqr, dormlq,
     $          dormqr, dpocon, dscal,  dswap,  dtrsm,  xerbla
*
      EXTERNAL  dgesvj
*     ..
*
*     Test the input arguments
*
      lsvec  = lsame( jobu, 'U' ) .OR. lsame( jobu, 'F' )
      jracc  = lsame( jobv, 'J' )
      rsvec  = lsame( jobv, 'V' ) .OR. jracc
      rowpiv = lsame( joba, 'F' ) .OR. lsame( joba, 'G' )
      l2rank = lsame( joba, 'R' )
      l2aber = lsame( joba, 'A' )
      errest = lsame( joba, 'E' ) .OR. lsame( joba, 'G' )
      l2tran = lsame( jobt, 'T' )
      l2kill = lsame( jobr, 'R' )
      defr   = lsame( jobr, 'N' )
      l2pert = lsame( jobp, 'P' )
*
      IF ( .NOT.(rowpiv .OR. l2rank .OR. l2aber .OR.
     $     errest .OR. lsame( joba, 'C' ) )) THEN
         info = - 1
      ELSE IF ( .NOT.( lsvec  .OR. lsame( jobu, 'N' ) .OR.
     $                             lsame( jobu, 'W' )) ) THEN
         info = - 2
      ELSE IF ( .NOT.( rsvec .OR. lsame( jobv, 'N' ) .OR.
     $   lsame( jobv, 'W' )) .OR. ( jracc .AND. (.NOT.lsvec) ) ) THEN
         info = - 3
      ELSE IF ( .NOT. ( l2kill .OR. defr ) )    THEN
         info = - 4
      ELSE IF ( .NOT. ( l2tran .OR. lsame( jobt, 'N' ) ) ) THEN
         info = - 5
      ELSE IF ( .NOT. ( l2pert .OR. lsame( jobp, 'N' ) ) ) THEN
         info = - 6
      ELSE IF ( m .LT. 0 ) THEN
         info = - 7
      ELSE IF ( ( n .LT. 0 ) .OR. ( n .GT. m ) ) THEN
         info = - 8
      ELSE IF ( lda .LT. m ) THEN
         info = - 10
      ELSE IF ( lsvec .AND. ( ldu .LT. m ) ) THEN
         info = - 13
      ELSE IF ( rsvec .AND. ( ldv .LT. n ) ) THEN
         info = - 15
      ELSE IF ( (.NOT.(lsvec .OR. rsvec .OR. errest).AND.
     &                           (lwork .LT. max(7,4*n+1,2*m+n))) .OR.
     & (.NOT.(lsvec .OR. rsvec) .AND. errest .AND.
     &                         (lwork .LT. max(7,4*n+n*n,2*m+n))) .OR.
     & (lsvec .AND. (.NOT.rsvec) .AND. (lwork .LT. max(7,2*m+n,4*n+1)))
     & .OR.
     & (rsvec .AND. (.NOT.lsvec) .AND. (lwork .LT. max(7,2*m+n,4*n+1)))
     & .OR.
     & (lsvec .AND. rsvec .AND. (.NOT.jracc) .AND.
     &                          (lwork.LT.max(2*m+n,6*n+2*n*n)))
     & .OR. (lsvec .AND. rsvec .AND. jracc .AND.
     &                          lwork.LT.max(2*m+n,4*n+n*n,2*n+n*n+6)))
     &   THEN
         info = - 17
      ELSE
*        #:)
         info = 0
      END IF
*
      IF ( info .NE. 0 ) THEN
*       #:(
         CALL xerbla( 'DGEJSV', - info )
         RETURN
      END IF
*
*     Quick return for void matrix (Y3K safe)
* #:)
      IF ( ( m .EQ. 0 ) .OR. ( n .EQ. 0 ) ) THEN
         iwork(1:3) = 0
         work(1:7) = 0
         RETURN
      ENDIF
*
*     Determine whether the matrix U should be M x N or M x M
*
      IF ( lsvec ) THEN
         n1 = n
         IF ( lsame( jobu, 'F' ) ) n1 = m
      END IF
*
*     Set numerical parameters
*
*!    NOTE: Make sure DLAMCH() does not fail on the target architecture.
*
      epsln = dlamch('Epsilon')
      sfmin = dlamch('SafeMinimum')
      small = sfmin / epsln
      big   = dlamch('O')
*     BIG   = ONE / SFMIN
*
*     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
*
*(!)  If necessary, scale SVA() to protect the largest norm from
*     overflow. It is possible that this scaling pushes the smallest
*     column norm left from the underflow threshold (extreme case).
*
      scalem  = one / dsqrt(dble(m)*dble(n))
      noscal  = .true.
      goscal  = .true.
      DO 1874 p = 1, n
         aapp = zero
         aaqq = one
         CALL dlassq( m, a(1,p), 1, aapp, aaqq )
         IF ( aapp .GT. big ) THEN
            info = - 9
            CALL xerbla( 'DGEJSV', -info )
            RETURN
         END IF
         aaqq = dsqrt(aaqq)
         IF ( ( aapp .LT. (big / aaqq) ) .AND. noscal  ) THEN
            sva(p)  = aapp * aaqq
         ELSE
            noscal  = .false.
            sva(p)  = aapp * ( aaqq * scalem )
            IF ( goscal ) THEN
               goscal = .false.
               CALL dscal( p-1, scalem, sva, 1 )
            END IF
         END IF
 1874 CONTINUE
*
      IF ( noscal ) scalem = one
*
      aapp = zero
      aaqq = big
      DO 4781 p = 1, n
         aapp = max( aapp, sva(p) )
         IF ( sva(p) .NE. zero ) aaqq = min( aaqq, sva(p) )
 4781 CONTINUE
*
*     Quick return for zero M x N matrix
* #:)
      IF ( aapp .EQ. zero ) THEN
         IF ( lsvec ) CALL dlaset( 'G', m, n1, zero, one, u, ldu )
         IF ( rsvec ) CALL dlaset( 'G', n, n,  zero, one, v, ldv )
         work(1) = one
         work(2) = one
         IF ( errest ) work(3) = one
         IF ( lsvec .AND. rsvec ) THEN
            work(4) = one
            work(5) = one
         END IF
         IF ( l2tran ) THEN
            work(6) = zero
            work(7) = zero
         END IF
         iwork(1) = 0
         iwork(2) = 0
         iwork(3) = 0
         RETURN
      END IF
*
*     Issue warning if denormalized column norms detected. Override the
*     high relative accuracy request. Issue licence to kill columns
*     (set them to zero) whose norm is less than sigma_max / BIG (roughly).
* #:(
      warning = 0
      IF ( aaqq .LE. sfmin ) THEN
         l2rank = .true.
         l2kill = .true.
         warning = 1
      END IF
*
*     Quick return for one-column matrix
* #:)
      IF ( n .EQ. 1 ) THEN
*
         IF ( lsvec ) THEN
            CALL dlascl( 'G',0,0,sva(1),scalem, m,1,a(1,1),lda,ierr )
            CALL dlacpy( 'A', m, 1, a, lda, u, ldu )
*           computing all M left singular vectors of the M x 1 matrix
            IF ( n1 .NE. n  ) THEN
               CALL dgeqrf( m, n, u,ldu, work, work(n+1),lwork-n,ierr )
               CALL dorgqr( m,n1,1, u,ldu,work,work(n+1),lwork-n,ierr )
               CALL dcopy( m, a(1,1), 1, u(1,1), 1 )
            END IF
         END IF
         IF ( rsvec ) THEN
             v(1,1) = one
         END IF
         IF ( sva(1) .LT. (big*scalem) ) THEN
            sva(1)  = sva(1) / scalem
            scalem  = one
         END IF
         work(1) = one / scalem
         work(2) = one
         IF ( sva(1) .NE. zero ) THEN
            iwork(1) = 1
            IF ( ( sva(1) / scalem) .GE. sfmin ) THEN
               iwork(2) = 1
            ELSE
               iwork(2) = 0
            END IF
         ELSE
            iwork(1) = 0
            iwork(2) = 0
         END IF
         iwork(3) = 0
         IF ( errest ) work(3) = one
         IF ( lsvec .AND. rsvec ) THEN
            work(4) = one
            work(5) = one
         END IF
         IF ( l2tran ) THEN
            work(6) = zero
            work(7) = zero
         END IF
         RETURN
*
      END IF
*
      transp = .false.
      l2tran = l2tran .AND. ( m .EQ. n )
*
      aatmax = -one
      aatmin =  big
      IF ( rowpiv .OR. l2tran ) THEN
*
*     Compute the row norms, needed to determine row pivoting sequence
*     (in the case of heavily row weighted A, row pivoting is strongly
*     advised) and to collect information needed to compare the
*     structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.).
*
         IF ( l2tran ) THEN
            DO 1950 p = 1, m
               xsc   = zero
               temp1 = one
               CALL dlassq( n, a(p,1), lda, xsc, temp1 )
*              DLASSQ gets both the ell_2 and the ell_infinity norm
*              in one pass through the vector
               work(m+n+p)  = xsc * scalem
               work(n+p)    = xsc * (scalem*dsqrt(temp1))
               aatmax = max( aatmax, work(n+p) )
               IF (work(n+p) .NE. zero) aatmin = min(aatmin,work(n+p))
 1950       CONTINUE
         ELSE
            DO 1904 p = 1, m
               work(m+n+p) = scalem*dabs( a(p,idamax(n,a(p,1),lda)) )
               aatmax = max( aatmax, work(m+n+p) )
               aatmin = min( aatmin, work(m+n+p) )
 1904       CONTINUE
         END IF
*
      END IF
*
*     For square matrix A try to determine whether A^t  would be  better
*     input for the preconditioned Jacobi SVD, with faster convergence.
*     The decision is based on an O(N) function of the vector of column
*     and row norms of A, based on the Shannon entropy. This should give
*     the right choice in most cases when the difference actually matters.
*     It may fail and pick the slower converging side.
*
      entra  = zero
      entrat = zero
      IF ( l2tran ) THEN
*
         xsc   = zero
         temp1 = one
         CALL dlassq( n, sva, 1, xsc, temp1 )
         temp1 = one / temp1
*
         entra = zero
         DO 1113 p = 1, n
            big1  = ( ( sva(p) / xsc )**2 ) * temp1
            IF ( big1 .NE. zero ) entra = entra + big1 * dlog(big1)
 1113    CONTINUE
         entra = - entra / dlog(dble(n))
*
*        Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex.
*        It is derived from the diagonal of  A^t * A.  Do the same with the
*        diagonal of A * A^t, compute the entropy of the corresponding
*        probability distribution. Note that A * A^t and A^t * A have the
*        same trace.
*
         entrat = zero
         DO 1114 p = n+1, n+m
            big1 = ( ( work(p) / xsc )**2 ) * temp1
            IF ( big1 .NE. zero ) entrat = entrat + big1 * dlog(big1)
 1114    CONTINUE
         entrat = - entrat / dlog(dble(m))
*
*        Analyze the entropies and decide A or A^t. Smaller entropy
*        usually means better input for the algorithm.
*
         transp = ( entrat .LT. entra )
*
*        If A^t is better than A, transpose A.
*
         IF ( transp ) THEN
*           In an optimal implementation, this trivial transpose
*           should be replaced with faster transpose.
            DO 1115 p = 1, n - 1
               DO 1116 q = p + 1, n
                   temp1 = a(q,p)
                  a(q,p) = a(p,q)
                  a(p,q) = temp1
 1116          CONTINUE
 1115       CONTINUE
            DO 1117 p = 1, n
               work(m+n+p) = sva(p)
               sva(p)      = work(n+p)
 1117       CONTINUE
            temp1  = aapp
            aapp   = aatmax
            aatmax = temp1
            temp1  = aaqq
            aaqq   = aatmin
            aatmin = temp1
            kill   = lsvec
            lsvec  = rsvec
            rsvec  = kill
            IF ( lsvec ) n1 = n
*
            rowpiv = .true.
         END IF
*
      END IF
*     END IF L2TRAN
*
*     Scale the matrix so that its maximal singular value remains less
*     than DSQRT(BIG) -- the matrix is scaled so that its maximal column
*     has Euclidean norm equal to DSQRT(BIG/N). The only reason to keep
*     DSQRT(BIG) instead of BIG is the fact that DGEJSV uses LAPACK and
*     BLAS routines that, in some implementations, are not capable of
*     working in the full interval [SFMIN,BIG] and that they may provoke
*     overflows in the intermediate results. If the singular values spread
*     from SFMIN to BIG, then DGESVJ will compute them. So, in that case,
*     one should use DGESVJ instead of DGEJSV.
*
      big1   = dsqrt( big )
      temp1  = dsqrt( big / dble(n) )
*
      CALL dlascl( 'G', 0, 0, aapp, temp1, n, 1, sva, n, ierr )
      IF ( aaqq .GT. (aapp * sfmin) ) THEN
          aaqq = ( aaqq / aapp ) * temp1
      ELSE
          aaqq = ( aaqq * temp1 ) / aapp
      END IF
      temp1 = temp1 * scalem
      CALL dlascl( 'G', 0, 0, aapp, temp1, m, n, a, lda, ierr )
*
*     To undo scaling at the end of this procedure, multiply the
*     computed singular values with USCAL2 / USCAL1.
*
      uscal1 = temp1
      uscal2 = aapp
*
      IF ( l2kill ) THEN
*        L2KILL enforces computation of nonzero singular values in
*        the restricted range of condition number of the initial A,
*        sigma_max(A) / sigma_min(A) approx. DSQRT(BIG)/DSQRT(SFMIN).
         xsc = dsqrt( sfmin )
      ELSE
         xsc = small
*
*        Now, if the condition number of A is too big,
*        sigma_max(A) / sigma_min(A) .GT. DSQRT(BIG/N) * EPSLN / SFMIN,
*        as a precaution measure, the full SVD is computed using DGESVJ
*        with accumulated Jacobi rotations. This provides numerically
*        more robust computation, at the cost of slightly increased run
*        time. Depending on the concrete implementation of BLAS and LAPACK
*        (i.e. how they behave in presence of extreme ill-conditioning) the
*        implementor may decide to remove this switch.
         IF ( ( aaqq.LT.dsqrt(sfmin) ) .AND. lsvec .AND. rsvec ) THEN
            jracc = .true.
         END IF
*
      END IF
      IF ( aaqq .LT. xsc ) THEN
         DO 700 p = 1, n
            IF ( sva(p) .LT. xsc ) THEN
               CALL dlaset( 'A', m, 1, zero, zero, a(1,p), lda )
               sva(p) = zero
            END IF
 700     CONTINUE
      END IF
*
*     Preconditioning using QR factorization with pivoting
*
      IF ( rowpiv ) THEN
*        Optional row permutation (Bjoerck row pivoting):
*        A result by Cox and Higham shows that the Bjoerck's
*        row pivoting combined with standard column pivoting
*        has similar effect as Powell-Reid complete pivoting.
*        The ell-infinity norms of A are made nonincreasing.
         DO 1952 p = 1, m - 1
            q = idamax( m-p+1, work(m+n+p), 1 ) + p - 1
            iwork(2*n+p) = q
            IF ( p .NE. q ) THEN
               temp1       = work(m+n+p)
               work(m+n+p) = work(m+n+q)
               work(m+n+q) = temp1
            END IF
 1952    CONTINUE
         CALL dlaswp( n, a, lda, 1, m-1, iwork(2*n+1), 1 )
      END IF
*
*     End of the preparation phase (scaling, optional sorting and
*     transposing, optional flushing of small columns).
*
*     Preconditioning
*
*     If the full SVD is needed, the right singular vectors are computed
*     from a matrix equation, and for that we need theoretical analysis
*     of the Businger-Golub pivoting. So we use DGEQP3 as the first RR QRF.
*     In all other cases the first RR QRF can be chosen by other criteria
*     (eg speed by replacing global with restricted window pivoting, such
*     as in SGEQPX from TOMS # 782). Good results will be obtained using
*     SGEQPX with properly (!) chosen numerical parameters.
*     Any improvement of DGEQP3 improves overall performance of DGEJSV.
*
*     A * P1 = Q1 * [ R1^t 0]^t:
      DO 1963 p = 1, n
*        .. all columns are free columns
         iwork(p) = 0
 1963 CONTINUE
      CALL dgeqp3( m,n,a,lda, iwork,work, work(n+1),lwork-n, ierr )
*
*     The upper triangular matrix R1 from the first QRF is inspected for
*     rank deficiency and possibilities for deflation, or possible
*     ill-conditioning. Depending on the user specified flag L2RANK,
*     the procedure explores possibilities to reduce the numerical
*     rank by inspecting the computed upper triangular factor. If
*     L2RANK or L2ABER are up, then DGEJSV will compute the SVD of
*     A + dA, where ||dA|| <= f(M,N)*EPSLN.
*
      nr = 1
      IF ( l2aber ) THEN
*        Standard absolute error bound suffices. All sigma_i with
*        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
*        aggressive enforcement of lower numerical rank by introducing a
*        backward error of the order of N*EPSLN*||A||.
         temp1 = dsqrt(dble(n))*epsln
         DO 3001 p = 2, n
            IF ( dabs(a(p,p)) .GE. (temp1*dabs(a(1,1))) ) THEN
               nr = nr + 1
            ELSE
               GO TO 3002
            END IF
 3001    CONTINUE
 3002    CONTINUE
      ELSE IF ( l2rank ) THEN
*        .. similarly as above, only slightly more gentle (less aggressive).
*        Sudden drop on the diagonal of R1 is used as the criterion for
*        close-to-rank-deficient.
         temp1 = dsqrt(sfmin)
         DO 3401 p = 2, n
            IF ( ( dabs(a(p,p)) .LT. (epsln*dabs(a(p-1,p-1))) ) .OR.
     $           ( dabs(a(p,p)) .LT. small ) .OR.
     $           ( l2kill .AND. (dabs(a(p,p)) .LT. temp1) ) ) GO TO 3402
            nr = nr + 1
 3401    CONTINUE
 3402    CONTINUE
*
      ELSE
*        The goal is high relative accuracy. However, if the matrix
*        has high scaled condition number the relative accuracy is in
*        general not feasible. Later on, a condition number estimator
*        will be deployed to estimate the scaled condition number.
*        Here we just remove the underflowed part of the triangular
*        factor. This prevents the situation in which the code is
*        working hard to get the accuracy not warranted by the data.
         temp1  = dsqrt(sfmin)
         DO 3301 p = 2, n
            IF ( ( dabs(a(p,p)) .LT. small ) .OR.
     $          ( l2kill .AND. (dabs(a(p,p)) .LT. temp1) ) ) GO TO 3302
            nr = nr + 1
 3301    CONTINUE
 3302    CONTINUE
*
      END IF
*
      almort = .false.
      IF ( nr .EQ. n ) THEN
         maxprj = one
         DO 3051 p = 2, n
            temp1  = dabs(a(p,p)) / sva(iwork(p))
            maxprj = min( maxprj, temp1 )
 3051    CONTINUE
         IF ( maxprj**2 .GE. one - dble(n)*epsln ) almort = .true.
      END IF
*
*
      sconda = - one
      condr1 = - one
      condr2 = - one
*
      IF ( errest ) THEN
         IF ( n .EQ. nr ) THEN
            IF ( rsvec ) THEN
*              .. V is available as workspace
               CALL dlacpy( 'U', n, n, a, lda, v, ldv )
               DO 3053 p = 1, n
                  temp1 = sva(iwork(p))
                  CALL dscal( p, one/temp1, v(1,p), 1 )
 3053          CONTINUE
               CALL dpocon( 'U', n, v, ldv, one, temp1,
     $              work(n+1), iwork(2*n+m+1), ierr )
            ELSE IF ( lsvec ) THEN
*              .. U is available as workspace
               CALL dlacpy( 'U', n, n, a, lda, u, ldu )
               DO 3054 p = 1, n
                  temp1 = sva(iwork(p))
                  CALL dscal( p, one/temp1, u(1,p), 1 )
 3054          CONTINUE
               CALL dpocon( 'U', n, u, ldu, one, temp1,
     $              work(n+1), iwork(2*n+m+1), ierr )
            ELSE
               CALL dlacpy( 'U', n, n, a, lda, work(n+1), n )
               DO 3052 p = 1, n
                  temp1 = sva(iwork(p))
                  CALL dscal( p, one/temp1, work(n+(p-1)*n+1), 1 )
 3052          CONTINUE
*           .. the columns of R are scaled to have unit Euclidean lengths.
               CALL dpocon( 'U', n, work(n+1), n, one, temp1,
     $              work(n+n*n+1), iwork(2*n+m+1), ierr )
            END IF
            sconda = one / dsqrt(temp1)
*           SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
*           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
         ELSE
            sconda = - one
         END IF
      END IF
*
      l2pert = l2pert .AND. ( dabs( a(1,1)/a(nr,nr) ) .GT. dsqrt(big1) )
*     If there is no violent scaling, artificial perturbation is not needed.
*
*     Phase 3:
*
      IF ( .NOT. ( rsvec .OR. lsvec ) ) THEN
*
*         Singular Values only
*
*         .. transpose A(1:NR,1:N)
         DO 1946 p = 1, min( n-1, nr )
            CALL dcopy( n-p, a(p,p+1), lda, a(p+1,p), 1 )
 1946    CONTINUE
*
*        The following two DO-loops introduce small relative perturbation
*        into the strict upper triangle of the lower triangular matrix.
*        Small entries below the main diagonal are also changed.
*        This modification is useful if the computing environment does not
*        provide/allow FLUSH TO ZERO underflow, for it prevents many
*        annoying denormalized numbers in case of strongly scaled matrices.
*        The perturbation is structured so that it does not introduce any
*        new perturbation of the singular values, and it does not destroy
*        the job done by the preconditioner.
*        The licence for this perturbation is in the variable L2PERT, which
*        should be .FALSE. if FLUSH TO ZERO underflow is active.
*
         IF ( .NOT. almort ) THEN
*
            IF ( l2pert ) THEN
*              XSC = DSQRT(SMALL)
               xsc = epsln / dble(n)
               DO 4947 q = 1, nr
                  temp1 = xsc*dabs(a(q,q))
                  DO 4949 p = 1, n
                     IF ( ( (p.GT.q) .AND. (dabs(a(p,q)).LE.temp1) )
     $                    .OR. ( p .LT. q ) )
     $                     a(p,q) = dsign( temp1, a(p,q) )
 4949             CONTINUE
 4947          CONTINUE
            ELSE
               CALL dlaset( 'U', nr-1,nr-1, zero,zero, a(1,2),lda )
            END IF
*
*            .. second preconditioning using the QR factorization
*
            CALL dgeqrf( n,nr, a,lda, work, work(n+1),lwork-n, ierr )
*
*           .. and transpose upper to lower triangular
            DO 1948 p = 1, nr - 1
               CALL dcopy( nr-p, a(p,p+1), lda, a(p+1,p), 1 )
 1948       CONTINUE
*
         END IF
*
*           Row-cyclic Jacobi SVD algorithm with column pivoting
*
*           .. again some perturbation (a "background noise") is added
*           to drown denormals
            IF ( l2pert ) THEN
*              XSC = DSQRT(SMALL)
               xsc = epsln / dble(n)
               DO 1947 q = 1, nr
                  temp1 = xsc*dabs(a(q,q))
                  DO 1949 p = 1, nr
                     IF ( ( (p.GT.q) .AND. (dabs(a(p,q)).LE.temp1) )
     $                       .OR. ( p .LT. q ) )
     $                   a(p,q) = dsign( temp1, a(p,q) )
 1949             CONTINUE
 1947          CONTINUE
            ELSE
               CALL dlaset( 'U', nr-1, nr-1, zero, zero, a(1,2), lda )
            END IF
*
*           .. and one-sided Jacobi rotations are started on a lower
*           triangular matrix (plus perturbation which is ignored in
*           the part which destroys triangular form (confusing?!))
*
            CALL dgesvj( 'L', 'NoU', 'NoV', nr, nr, a, lda, sva,
     $                      n, v, ldv, work, lwork, info )
*
            scalem  = work(1)
            numrank = idnint(work(2))
*
*
      ELSE IF ( rsvec .AND. ( .NOT. lsvec ) ) THEN
*
*        -> Singular Values and Right Singular Vectors <-
*
         IF ( almort ) THEN
*
*           .. in this case NR equals N
            DO 1998 p = 1, nr
               CALL dcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
 1998       CONTINUE
            CALL dlaset( 'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
*
            CALL dgesvj( 'L','U','N', n, nr, v,ldv, sva, nr, a,lda,
     $                  work, lwork, info )
            scalem  = work(1)
            numrank = idnint(work(2))
 
         ELSE
*
*        .. two more QR factorizations ( one QRF is not enough, two require
*        accumulated product of Jacobi rotations, three are perfect )
*
            CALL dlaset( 'Lower', nr-1, nr-1, zero, zero, a(2,1), lda )
            CALL dgelqf( nr, n, a, lda, work, work(n+1), lwork-n, ierr)
            CALL dlacpy( 'Lower', nr, nr, a, lda, v, ldv )
            CALL dlaset( 'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
            CALL dgeqrf( nr, nr, v, ldv, work(n+1), work(2*n+1),
     $                   lwork-2*n, ierr )
            DO 8998 p = 1, nr
               CALL dcopy( nr-p+1, v(p,p), ldv, v(p,p), 1 )
 8998       CONTINUE
            CALL dlaset( 'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
*
            CALL dgesvj( 'Lower', 'U','N', nr, nr, v,ldv, sva, nr, u,
     $                  ldu, work(n+1), lwork, info )
            scalem  = work(n+1)
            numrank = idnint(work(n+2))
            IF ( nr .LT. n ) THEN
               CALL dlaset( 'A',n-nr, nr, zero,zero, v(nr+1,1),   ldv )
               CALL dlaset( 'A',nr, n-nr, zero,zero, v(1,nr+1),   ldv )
               CALL dlaset( 'A',n-nr,n-nr,zero,one, v(nr+1,nr+1), ldv )
            END IF
*
         CALL dormlq( 'Left', 'Transpose', n, n, nr, a, lda, work,
     $               v, ldv, work(n+1), lwork-n, ierr )
*
         END IF
*
         DO 8991 p = 1, n
            CALL dcopy( n, v(p,1), ldv, a(iwork(p),1), lda )
 8991    CONTINUE
         CALL dlacpy( 'All', n, n, a, lda, v, ldv )
*
         IF ( transp ) THEN
            CALL dlacpy( 'All', n, n, v, ldv, u, ldu )
         END IF
*
      ELSE IF ( lsvec .AND. ( .NOT. rsvec ) ) THEN
*
*        .. Singular Values and Left Singular Vectors                 ..
*
*        .. second preconditioning step to avoid need to accumulate
*        Jacobi rotations in the Jacobi iterations.
         DO 1965 p = 1, nr
            CALL dcopy( n-p+1, a(p,p), lda, u(p,p), 1 )
 1965    CONTINUE
         CALL dlaset( 'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
*
         CALL dgeqrf( n, nr, u, ldu, work(n+1), work(2*n+1),
     $              lwork-2*n, ierr )
*
         DO 1967 p = 1, nr - 1
            CALL dcopy( nr-p, u(p,p+1), ldu, u(p+1,p), 1 )
 1967    CONTINUE
         CALL dlaset( 'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
*
         CALL dgesvj( 'Lower', 'U', 'N', nr,nr, u, ldu, sva, nr, a,
     $        lda, work(n+1), lwork-n, info )
         scalem  = work(n+1)
         numrank = idnint(work(n+2))
*
         IF ( nr .LT. m ) THEN
            CALL dlaset( 'A',  m-nr, nr,zero, zero, u(nr+1,1), ldu )
            IF ( nr .LT. n1 ) THEN
               CALL dlaset( 'A',nr, n1-nr, zero, zero, u(1,nr+1), ldu )
               CALL dlaset( 'A',m-nr,n1-nr,zero,one,u(nr+1,nr+1), ldu )
            END IF
         END IF
*
         CALL dormqr( 'Left', 'No Tr', m, n1, n, a, lda, work, u,
     $               ldu, work(n+1), lwork-n, ierr )
*
         IF ( rowpiv )
     $       CALL dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
*
         DO 1974 p = 1, n1
            xsc = one / dnrm2( m, u(1,p), 1 )
            CALL dscal( m, xsc, u(1,p), 1 )
 1974    CONTINUE
*
         IF ( transp ) THEN
            CALL dlacpy( 'All', n, n, u, ldu, v, ldv )
         END IF
*
      ELSE
*
*        .. Full SVD ..
*
         IF ( .NOT. jracc ) THEN
*
         IF ( .NOT. almort ) THEN
*
*           Second Preconditioning Step (QRF [with pivoting])
*           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
*           equivalent to an LQF CALL. Since in many libraries the QRF
*           seems to be better optimized than the LQF, we do explicit
*           transpose and use the QRF. This is subject to changes in an
*           optimized implementation of DGEJSV.
*
            DO 1968 p = 1, nr
               CALL dcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
 1968       CONTINUE
*
*           .. the following two loops perturb small entries to avoid
*           denormals in the second QR factorization, where they are
*           as good as zeros. This is done to avoid painfully slow
*           computation with denormals. The relative size of the perturbation
*           is a parameter that can be changed by the implementer.
*           This perturbation device will be obsolete on machines with
*           properly implemented arithmetic.
*           To switch it off, set L2PERT=.FALSE. To remove it from  the
*           code, remove the action under L2PERT=.TRUE., leave the ELSE part.
*           The following two loops should be blocked and fused with the
*           transposed copy above.
*
            IF ( l2pert ) THEN
               xsc = dsqrt(small)
               DO 2969 q = 1, nr
                  temp1 = xsc*dabs( v(q,q) )
                  DO 2968 p = 1, n
                     IF ( ( p .GT. q ) .AND. ( dabs(v(p,q)) .LE. temp1 )
     $                   .OR. ( p .LT. q ) )
     $                   v(p,q) = dsign( temp1, v(p,q) )
                     IF ( p .LT. q ) v(p,q) = - v(p,q)
 2968             CONTINUE
 2969          CONTINUE
            ELSE
               CALL dlaset( 'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
            END IF
*
*           Estimate the row scaled condition number of R1
*           (If R1 is rectangular, N > NR, then the condition number
*           of the leading NR x NR submatrix is estimated.)
*
            CALL dlacpy( 'L', nr, nr, v, ldv, work(2*n+1), nr )
            DO 3950 p = 1, nr
               temp1 = dnrm2(nr-p+1,work(2*n+(p-1)*nr+p),1)
               CALL dscal(nr-p+1,one/temp1,work(2*n+(p-1)*nr+p),1)
 3950       CONTINUE
            CALL dpocon('Lower',nr,work(2*n+1),nr,one,temp1,
     $                   work(2*n+nr*nr+1),iwork(m+2*n+1),ierr)
            condr1 = one / dsqrt(temp1)
*           .. here need a second opinion on the condition number
*           .. then assume worst case scenario
*           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
*           more conservative    <=> CONDR1 .LT. DSQRT(DBLE(N))
*
            cond_ok = dsqrt(dble(nr))
*[TP]       COND_OK is a tuning parameter.
 
            IF ( condr1 .LT. cond_ok ) THEN
*              .. the second QRF without pivoting. Note: in an optimized
*              implementation, this QRF should be implemented as the QRF
*              of a lower triangular matrix.
*              R1^t = Q2 * R2
               CALL dgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
     $              lwork-2*n, ierr )
*
               IF ( l2pert ) THEN
                  xsc = dsqrt(small)/epsln
                  DO 3959 p = 2, nr
                     DO 3958 q = 1, p - 1
                        temp1 = xsc * min(dabs(v(p,p)),dabs(v(q,q)))
                        IF ( dabs(v(q,p)) .LE. temp1 )
     $                     v(q,p) = dsign( temp1, v(q,p) )
 3958                CONTINUE
 3959             CONTINUE
               END IF
*
               IF ( nr .NE. n )
     $         CALL dlacpy( 'A', n, nr, v, ldv, work(2*n+1), n )
*              .. save ...
*
*           .. this transposed copy should be better than naive
               DO 1969 p = 1, nr - 1
                  CALL dcopy( nr-p, v(p,p+1), ldv, v(p+1,p), 1 )
 1969          CONTINUE
*
               condr2 = condr1
*
            ELSE
*
*              .. ill-conditioned case: second QRF with pivoting
*              Note that windowed pivoting would be equally good
*              numerically, and more run-time efficient. So, in
*              an optimal implementation, the next call to DGEQP3
*              should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
*              with properly (carefully) chosen parameters.
*
*              R1^t * P2 = Q2 * R2
               DO 3003 p = 1, nr
                  iwork(n+p) = 0
 3003          CONTINUE
               CALL dgeqp3( n, nr, v, ldv, iwork(n+1), work(n+1),
     $                  work(2*n+1), lwork-2*n, ierr )
**               CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
**     $              LWORK-2*N, IERR )
               IF ( l2pert ) THEN
                  xsc = dsqrt(small)
                  DO 3969 p = 2, nr
                     DO 3968 q = 1, p - 1
                        temp1 = xsc * min(dabs(v(p,p)),dabs(v(q,q)))
                        IF ( dabs(v(q,p)) .LE. temp1 )
     $                     v(q,p) = dsign( temp1, v(q,p) )
 3968                CONTINUE
 3969             CONTINUE
               END IF
*
               CALL dlacpy( 'A', n, nr, v, ldv, work(2*n+1), n )
*
               IF ( l2pert ) THEN
                  xsc = dsqrt(small)
                  DO 8970 p = 2, nr
                     DO 8971 q = 1, p - 1
                        temp1 = xsc * min(dabs(v(p,p)),dabs(v(q,q)))
                        v(p,q) = - dsign( temp1, v(q,p) )
 8971                CONTINUE
 8970             CONTINUE
               ELSE
                  CALL dlaset( 'L',nr-1,nr-1,zero,zero,v(2,1),ldv )
               END IF
*              Now, compute R2 = L3 * Q3, the LQ factorization.
               CALL dgelqf( nr, nr, v, ldv, work(2*n+n*nr+1),
     $               work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, ierr )
*              .. and estimate the condition number
               CALL dlacpy( 'L',nr,nr,v,ldv,work(2*n+n*nr+nr+1),nr )
               DO 4950 p = 1, nr
                  temp1 = dnrm2( p, work(2*n+n*nr+nr+p), nr )
                  CALL dscal( p, one/temp1, work(2*n+n*nr+nr+p), nr )
 4950          CONTINUE
               CALL dpocon( 'L',nr,work(2*n+n*nr+nr+1),nr,one,temp1,
     $              work(2*n+n*nr+nr+nr*nr+1),iwork(m+2*n+1),ierr )
               condr2 = one / dsqrt(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
*                 needed in this branch, but it does not overwrite the
*                 Huseholder vectors of Q2.).
                  CALL dlacpy( 'U', nr, nr, v, ldv, work(2*n+1), n )
*                 .. 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
               xsc = dsqrt(small)
               DO 4968 q = 2, nr
                  temp1 = xsc * v(q,q)
                  DO 4969 p = 1, q - 1
*                    V(p,q) = - DSIGN( TEMP1, V(q,p) )
                     v(p,q) = - dsign( temp1, v(p,q) )
 4969             CONTINUE
 4968          CONTINUE
            ELSE
               CALL dlaset( 'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
            END IF
*
*        Second preconditioning finished; continue with Jacobi SVD
*        The input matrix is lower triangular.
*
*        Recover the right singular vectors as solution of a well
*        conditioned triangular matrix equation.
*
            IF ( condr1 .LT. cond_ok ) THEN
*
               CALL dgesvj( 'L','U','N',nr,nr,v,ldv,sva,nr,u,
     $              ldu,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,info )
               scalem  = work(2*n+n*nr+nr+1)
               numrank = idnint(work(2*n+n*nr+nr+2))
               DO 3970 p = 1, nr
                  CALL dcopy( nr, v(1,p), 1, u(1,p), 1 )
                  CALL dscal( nr, sva(p),    v(1,p), 1 )
 3970          CONTINUE
 
*        .. pick the right matrix equation and solve it
*
               IF ( nr .EQ. n ) THEN
* :))             .. best case, R1 is inverted. The solution of this matrix
*                 equation is Q2*V2 = the product of the Jacobi rotations
*                 used in DGESVJ, premultiplied with the orthogonal matrix
*                 from the second QR factorization.
                  CALL dtrsm( 'L','U','N','N', nr,nr,one, a,lda, v,ldv )
               ELSE
*                 .. R1 is well conditioned, but non-square. Transpose(R2)
*                 is inverted to get the product of the Jacobi rotations
*                 used in DGESVJ. The Q-factor from the second QR
*                 factorization is then built in explicitly.
                  CALL dtrsm('L','U','T','N',nr,nr,one,work(2*n+1),
     $                 n,v,ldv)
                  IF ( nr .LT. n ) THEN
                    CALL dlaset('A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
                    CALL dlaset('A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
                    CALL dlaset('A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv)
                  END IF
                  CALL dormqr('L','N',n,n,nr,work(2*n+1),n,work(n+1),
     $                 v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr)
               END IF
*
            ELSE IF ( condr2 .LT. cond_ok ) THEN
*
* :)           .. the input matrix A is very likely a relative of
*              the Kahan matrix :)
*              The matrix R2 is inverted. The solution of the matrix equation
*              is Q3^T*V3 = the product of the Jacobi rotations (applied to
*              the lower triangular L3 from the LQ factorization of
*              R2=L3*Q3), pre-multiplied with the transposed Q3.
               CALL dgesvj( 'L', 'U', 'N', nr, nr, v, ldv, sva, nr, u,
     $              ldu, work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, info )
               scalem  = work(2*n+n*nr+nr+1)
               numrank = idnint(work(2*n+n*nr+nr+2))
               DO 3870 p = 1, nr
                  CALL dcopy( nr, v(1,p), 1, u(1,p), 1 )
                  CALL dscal( nr, sva(p),    u(1,p), 1 )
 3870          CONTINUE
               CALL dtrsm('L','U','N','N',nr,nr,one,work(2*n+1),n,u,ldu)
*              .. apply the permutation from the second QR factorization
               DO 873 q = 1, nr
                  DO 872 p = 1, nr
                     work(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
 872              CONTINUE
                  DO 874 p = 1, nr
                     u(p,q) = work(2*n+n*nr+nr+p)
 874              CONTINUE
 873           CONTINUE
               IF ( nr .LT. n ) THEN
                  CALL dlaset( 'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
                  CALL dlaset( 'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
                  CALL dlaset( 'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
               END IF
               CALL dormqr( 'L','N',n,n,nr,work(2*n+1),n,work(n+1),
     $              v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
            ELSE
*              Last line of defense.
* #:(          This is a rather pathological case: no scaled condition
*              improvement after two pivoted QR factorizations. Other
*              possibility is that the rank revealing QR factorization
*              or the condition estimator has failed, or the COND_OK
*              is set very close to ONE (which is unnecessary). Normally,
*              this branch should never be executed, but in rare cases of
*              failure of the RRQR or condition estimator, the last line of
*              defense ensures that DGEJSV completes the task.
*              Compute the full SVD of L3 using DGESVJ with explicit
*              accumulation of Jacobi rotations.
               CALL dgesvj( 'L', 'U', 'V', nr, nr, v, ldv, sva, nr, u,
     $              ldu, work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, info )
               scalem  = work(2*n+n*nr+nr+1)
               numrank = idnint(work(2*n+n*nr+nr+2))
               IF ( nr .LT. n ) THEN
                  CALL dlaset( 'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
                  CALL dlaset( 'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
                  CALL dlaset( 'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
               END IF
               CALL dormqr( 'L','N',n,n,nr,work(2*n+1),n,work(n+1),
     $              v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
*
               CALL dormlq( 'L', 'T', nr, nr, nr, work(2*n+1), n,
     $              work(2*n+n*nr+1), u, ldu, work(2*n+n*nr+nr+1),
     $              lwork-2*n-n*nr-nr, ierr )
               DO 773 q = 1, nr
                  DO 772 p = 1, nr
                     work(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
 772              CONTINUE
                  DO 774 p = 1, nr
                     u(p,q) = work(2*n+n*nr+nr+p)
 774              CONTINUE
 773           CONTINUE
*
            END IF
*
*           Permute the rows of V using the (column) permutation from the
*           first QRF. Also, scale the columns to make them unit in
*           Euclidean norm. This applies to all cases.
*
            temp1 = dsqrt(dble(n)) * epsln
            DO 1972 q = 1, n
               DO 972 p = 1, n
                  work(2*n+n*nr+nr+iwork(p)) = v(p,q)
  972          CONTINUE
               DO 973 p = 1, n
                  v(p,q) = work(2*n+n*nr+nr+p)
  973          CONTINUE
               xsc = one / dnrm2( n, v(1,q), 1 )
               IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
     $           CALL dscal( n, xsc, v(1,q), 1 )
 1972       CONTINUE
*           At this moment, V contains the right singular vectors of A.
*           Next, assemble the left singular vector matrix U (M x N).
            IF ( nr .LT. m ) THEN
               CALL dlaset( 'A', m-nr, nr, zero, zero, u(nr+1,1), ldu )
               IF ( nr .LT. n1 ) THEN
                  CALL dlaset('A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
                  CALL dlaset('A',m-nr,n1-nr,zero,one,u(nr+1,nr+1),ldu)
               END IF
            END IF
*
*           The Q matrix from the first QRF is built into the left singular
*           matrix U. This applies to all cases.
*
            CALL dormqr( 'Left', 'No_Tr', m, n1, n, a, lda, work, u,
     $           ldu, work(n+1), lwork-n, ierr )
 
*           The columns of U are normalized. The cost is O(M*N) flops.
            temp1 = dsqrt(dble(m)) * epsln
            DO 1973 p = 1, nr
               xsc = one / dnrm2( m, u(1,p), 1 )
               IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
     $          CALL dscal( m, xsc, u(1,p), 1 )
 1973       CONTINUE
*
*           If the initial QRF is computed with row pivoting, the left
*           singular vectors must be adjusted.
*
            IF ( rowpiv )
     $          CALL dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
*
         ELSE
*
*        .. the initial matrix A has almost orthogonal columns and
*        the second QRF is not needed
*
            CALL dlacpy( 'Upper', n, n, a, lda, work(n+1), n )
            IF ( l2pert ) THEN
               xsc = dsqrt(small)
               DO 5970 p = 2, n
                  temp1 = xsc * work( n + (p-1)*n + p )
                  DO 5971 q = 1, p - 1
                     work(n+(q-1)*n+p)=-dsign(temp1,work(n+(p-1)*n+q))
 5971             CONTINUE
 5970          CONTINUE
            ELSE
               CALL dlaset( 'Lower',n-1,n-1,zero,zero,work(n+2),n )
            END IF
*
            CALL dgesvj( 'Upper', 'U', 'N', n, n, work(n+1), n, sva,
     $           n, u, ldu, work(n+n*n+1), lwork-n-n*n, info )
*
            scalem  = work(n+n*n+1)
            numrank = idnint(work(n+n*n+2))
            DO 6970 p = 1, n
               CALL dcopy( n, work(n+(p-1)*n+1), 1, u(1,p), 1 )
               CALL dscal( n, sva(p), work(n+(p-1)*n+1), 1 )
 6970       CONTINUE
*
            CALL dtrsm( 'Left', 'Upper', 'NoTrans', 'No UD', n, n,
     $           one, a, lda, work(n+1), n )
            DO 6972 p = 1, n
               CALL dcopy( n, work(n+p), n, v(iwork(p),1), ldv )
 6972       CONTINUE
            temp1 = dsqrt(dble(n))*epsln
            DO 6971 p = 1, n
               xsc = one / dnrm2( n, v(1,p), 1 )
               IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
     $            CALL dscal( n, xsc, v(1,p), 1 )
 6971       CONTINUE
*
*           Assemble the left singular vector matrix U (M x N).
*
            IF ( n .LT. m ) THEN
               CALL dlaset( 'A',  m-n, n, zero, zero, u(n+1,1), ldu )
               IF ( n .LT. n1 ) THEN
                  CALL dlaset( 'A',n,  n1-n, zero, zero,  u(1,n+1),ldu )
                  CALL dlaset( 'A',m-n,n1-n, zero, one,u(n+1,n+1),ldu )
               END IF
            END IF
            CALL dormqr( 'Left', 'No Tr', m, n1, n, a, lda, work, u,
     $           ldu, work(n+1), lwork-n, ierr )
            temp1 = dsqrt(dble(m))*epsln
            DO 6973 p = 1, n1
               xsc = one / dnrm2( m, u(1,p), 1 )
               IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
     $            CALL dscal( m, xsc, u(1,p), 1 )
 6973       CONTINUE
*
            IF ( rowpiv )
     $         CALL dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
*
         END IF
*
*        end of the  >> almost orthogonal case <<  in the full SVD
*
         ELSE
*
*        This branch deploys a preconditioned Jacobi SVD with explicitly
*        accumulated rotations. It is included as optional, mainly for
*        experimental purposes. It does perform well, and can also be used.
*        In this implementation, this branch will be automatically activated
*        if the  condition number sigma_max(A) / sigma_min(A) is predicted
*        to be greater than the overflow threshold. This is because the
*        a posteriori computation of the singular vectors assumes robust
*        implementation of BLAS and some LAPACK procedures, capable of working
*        in presence of extreme values. Since that is not always the case, ...
*
         DO 7968 p = 1, nr
            CALL dcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
 7968    CONTINUE
*
         IF ( l2pert ) THEN
            xsc = dsqrt(small/epsln)
            DO 5969 q = 1, nr
               temp1 = xsc*dabs( v(q,q) )
               DO 5968 p = 1, n
                  IF ( ( p .GT. q ) .AND. ( dabs(v(p,q)) .LE. temp1 )
     $                .OR. ( p .LT. q ) )
     $                v(p,q) = dsign( temp1, v(p,q) )
                  IF ( p .LT. q ) v(p,q) = - v(p,q)
 5968          CONTINUE
 5969       CONTINUE
         ELSE
            CALL dlaset( 'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
         END IF
 
         CALL dgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
     $        lwork-2*n, ierr )
         CALL dlacpy( 'L', n, nr, v, ldv, work(2*n+1), n )
*
         DO 7969 p = 1, nr
            CALL dcopy( nr-p+1, v(p,p), ldv, u(p,p), 1 )
 7969    CONTINUE
 
         IF ( l2pert ) THEN
            xsc = dsqrt(small/epsln)
            DO 9970 q = 2, nr
               DO 9971 p = 1, q - 1
                  temp1 = xsc * min(dabs(u(p,p)),dabs(u(q,q)))
                  u(p,q) = - dsign( temp1, u(q,p) )
 9971          CONTINUE
 9970       CONTINUE
         ELSE
            CALL dlaset('U', nr-1, nr-1, zero, zero, u(1,2), ldu )
         END IF
 
         CALL dgesvj( 'G', 'U', 'V', nr, nr, u, ldu, sva,
     $        n, v, ldv, work(2*n+n*nr+1), lwork-2*n-n*nr, info )
         scalem  = work(2*n+n*nr+1)
         numrank = idnint(work(2*n+n*nr+2))
 
         IF ( nr .LT. n ) THEN
            CALL dlaset( 'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
            CALL dlaset( 'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
            CALL dlaset( 'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
         END IF
 
         CALL dormqr( 'L','N',n,n,nr,work(2*n+1),n,work(n+1),
     $        v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
*
*           Permute the rows of V using the (column) permutation from the
*           first QRF. Also, scale the columns to make them unit in
*           Euclidean norm. This applies to all cases.
*
            temp1 = dsqrt(dble(n)) * epsln
            DO 7972 q = 1, n
               DO 8972 p = 1, n
                  work(2*n+n*nr+nr+iwork(p)) = v(p,q)
 8972          CONTINUE
               DO 8973 p = 1, n
                  v(p,q) = work(2*n+n*nr+nr+p)
 8973          CONTINUE
               xsc = one / dnrm2( n, v(1,q), 1 )
               IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
     $           CALL dscal( n, xsc, v(1,q), 1 )
 7972       CONTINUE
*
*           At this moment, V contains the right singular vectors of A.
*           Next, assemble the left singular vector matrix U (M x N).
*
         IF ( nr .LT. m ) THEN
            CALL dlaset( 'A',  m-nr, nr, zero, zero, u(nr+1,1), ldu )
            IF ( nr .LT. n1 ) THEN
               CALL dlaset( 'A',nr,  n1-nr, zero, zero,  u(1,nr+1),ldu )
               CALL dlaset( 'A',m-nr,n1-nr, zero, one,u(nr+1,nr+1),ldu )
            END IF
         END IF
*
         CALL dormqr( 'Left', 'No Tr', m, n1, n, a, lda, work, u,
     $        ldu, work(n+1), lwork-n, ierr )
*
            IF ( rowpiv )
     $         CALL dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
*
*
         END IF
         IF ( transp ) THEN
*           .. swap U and V because the procedure worked on A^t
            DO 6974 p = 1, n
               CALL dswap( n, u(1,p), 1, v(1,p), 1 )
 6974       CONTINUE
         END IF
*
      END IF
*     end of the full SVD
*
*     Undo scaling, if necessary (and possible)
*
      IF ( uscal2 .LE. (big/sva(1))*uscal1 ) THEN
         CALL dlascl( 'G', 0, 0, uscal1, uscal2, nr, 1, sva, n, ierr )
         uscal1 = one
         uscal2 = one
      END IF
*
      IF ( nr .LT. n ) THEN
         DO 3004 p = nr+1, n
            sva(p) = zero
 3004    CONTINUE
      END IF
*
      work(1) = uscal2 * scalem
      work(2) = uscal1
      IF ( errest ) work(3) = sconda
      IF ( lsvec .AND. rsvec ) THEN
         work(4) = condr1
         work(5) = condr2
      END IF
      IF ( l2tran ) THEN
         work(6) = entra
         work(7) = entrat
      END IF
*
      iwork(1) = nr
      iwork(2) = numrank
      iwork(3) = warning
*
      RETURN
*     ..
*     .. END OF DGEJSV
*     ..
      END
*