dgesvdq.f

Go to the documentation of this file.

*> \brief <b> DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices</b>
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download DGESVDQ + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvdq.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvdq.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvdq.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*      SUBROUTINE DGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
*                          S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
*                          WORK, LWORK, RWORK, LRWORK, INFO )
*
*     .. Scalar Arguments ..
*      IMPLICIT    NONE
*      CHARACTER   JOBA, JOBP, JOBR, JOBU, JOBV
*      INTEGER     M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
*                  INFO
*     ..
*     .. Array Arguments ..
*      DOUBLE PRECISION  A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
*      DOUBLE PRECISION  S( * ), RWORK( * )
*      INTEGER     IWORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGESVDQ computes the singular value decomposition (SVD) of a real
*> M-by-N matrix A, where M >= N. The SVD of A is written as
*>                                    [++]   [xx]   [x0]   [xx]
*>              A = U * SIGMA * V^*,  [++] = [xx] * [ox] * [xx]
*>                                    [++]   [xx]
*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
*> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
*> of SIGMA are the singular values of A. The columns of U and V are the
*> left and the right singular vectors of A, respectively.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOBA
*> \verbatim
*>  JOBA is CHARACTER*1
*>  Specifies the level of accuracy in the computed SVD
*>  = 'A' The requested accuracy corresponds to having the backward
*>        error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F,
*>        where EPS = DLAMCH('Epsilon'). This authorises DGESVDQ to
*>        truncate the computed triangular factor in a rank revealing
*>        QR factorization whenever the truncated part is below the
*>        threshold of the order of EPS * ||A||_F. This is aggressive
*>        truncation level.
*>  = 'M' Similarly as with 'A', but the truncation is more gentle: it
*>        is allowed only when there is a drop on the diagonal of the
*>        triangular factor in the QR factorization. This is medium
*>        truncation level.
*>  = 'H' High accuracy requested. No numerical rank determination based
*>        on the rank revealing QR factorization is attempted.
*>  = 'E' Same as 'H', and in addition the condition number of column
*>        scaled A is estimated and returned in  RWORK(1).
*>        N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
*> \endverbatim
*>
*> \param[in] JOBP
*> \verbatim
*>  JOBP is CHARACTER*1
*>  = 'P' The rows of A are ordered in decreasing order with respect to
*>        ||A(i,:)||_\infty. This enhances numerical accuracy at the cost
*>        of extra data movement. Recommended for numerical robustness.
*>  = 'N' No row pivoting.
*> \endverbatim
*>
*> \param[in] JOBR
*> \verbatim
*>          JOBR is CHARACTER*1
*>          = 'T' After the initial pivoted QR factorization, DGESVD is applied to
*>          the transposed R**T of the computed triangular factor R. This involves
*>          some extra data movement (matrix transpositions). Useful for
*>          experiments, research and development.
*>          = 'N' The triangular factor R is given as input to DGESVD. This may be
*>          preferred as it involves less data movement.
*> \endverbatim
*>
*> \param[in] JOBU
*> \verbatim
*>          JOBU is CHARACTER*1
*>          = 'A' All M left singular vectors are computed and returned in the
*>          matrix U. See the description of U.
*>          = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned
*>          in the matrix U. See the description of U.
*>          = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular
*>          vectors are computed and returned in the matrix U.
*>          = 'F' The N left singular vectors are returned in factored form as the
*>          product of the Q factor from the initial QR factorization and the
*>          N left singular vectors of (R**T , 0)**T. If row pivoting is used,
*>          then the necessary information on the row pivoting is stored in
*>          IWORK(N+1:N+M-1).
*>          = 'N' The left singular vectors are not computed.
*> \endverbatim
*>
*> \param[in] JOBV
*> \verbatim
*>          JOBV is CHARACTER*1
*>          = 'A', 'V' All N right singular vectors are computed and returned in
*>          the matrix V.
*>          = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular
*>          vectors are computed and returned in the matrix V. This option is
*>          allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal.
*>          = 'N' The right singular vectors are not computed.
*> \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 of dimensions LDA x N
*>          On entry, the input matrix A.
*>          On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains
*>          the Householder vectors as stored by DGEQP3. If JOBU = 'F', these Householder
*>          vectors together with WORK(1:N) can be used to restore the Q factors from
*>          the initial pivoted QR factorization of A. See the description of U.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER.
*>          The leading dimension of the array A.  LDA >= max(1,M).
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          S is DOUBLE PRECISION array of dimension N.
*>          The singular values of A, ordered so that S(i) >= S(i+1).
*> \endverbatim
*>
*> \param[out] U
*> \verbatim
*>          U is DOUBLE PRECISION array, dimension
*>          LDU x M if JOBU = 'A'; see the description of LDU. In this case,
*>          on exit, U contains the M left singular vectors.
*>          LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this
*>          case, U contains the leading N or the leading NUMRANK left singular vectors.
*>          LDU x N if JOBU = 'F' ; see the description of LDU. In this case U
*>          contains N x N orthogonal matrix that can be used to form the left
*>          singular vectors.
*>          If JOBU = 'N', U is not referenced.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER.
*>          The leading dimension of the array U.
*>          If JOBU = 'A', 'S', 'U', 'R',  LDU >= max(1,M).
*>          If JOBU = 'F',                 LDU >= max(1,N).
*>          Otherwise,                     LDU >= 1.
*> \endverbatim
*>
*> \param[out] V
*> \verbatim
*>          V is DOUBLE PRECISION array, dimension
*>          LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' .
*>          If JOBV = 'A', or 'V',  V contains the N-by-N orthogonal matrix  V**T;
*>          If JOBV = 'R', V contains the first NUMRANK rows of V**T (the right
*>          singular vectors, stored rowwise, of the NUMRANK largest singular values).
*>          If JOBV = 'N' and JOBA = 'E', V is used as a workspace.
*>          If JOBV = 'N', and JOBA.NE.'E', V is not referenced.
*> \endverbatim
*>
*> \param[in] LDV
*> \verbatim
*>          LDV is INTEGER
*>          The leading dimension of the array V.
*>          If JOBV = 'A', 'V', 'R',  or JOBA = 'E', LDV >= max(1,N).
*>          Otherwise,                               LDV >= 1.
*> \endverbatim
*>
*> \param[out] NUMRANK
*> \verbatim
*>          NUMRANK is INTEGER
*>          NUMRANK is the numerical rank first determined after the rank
*>          revealing QR factorization, following the strategy specified by the
*>          value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK
*>          leading singular values and vectors are then requested in the call
*>          of DGESVD. The final value of NUMRANK might be further reduced if
*>          some singular values are computed as zeros.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (max(1, LIWORK)).
*>          On exit, IWORK(1:N) contains column pivoting permutation of the
*>          rank revealing QR factorization.
*>          If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence
*>          of row swaps used in row pivoting. These can be used to restore the
*>          left singular vectors in the case JOBU = 'F'.
*>
*>          If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
*>          IWORK(1) returns the minimal LIWORK.
*> \endverbatim
*>
*> \param[in] LIWORK
*> \verbatim
*>          LIWORK is INTEGER
*>          The dimension of the array IWORK.
*>          LIWORK >= N + M - 1,     if JOBP = 'P' and JOBA .NE. 'E';
*>          LIWORK >= N              if JOBP = 'N' and JOBA .NE. 'E';
*>          LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E';
*>          LIWORK >= N + N          if JOBP = 'N' and JOBA = 'E'.
*>
*>          If LIWORK = -1, then a workspace query is assumed; the routine
*>          only calculates and returns the optimal and minimal sizes
*>          for the WORK, IWORK, and RWORK arrays, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (max(2, LWORK)), used as a workspace.
*>          On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters
*>          needed to recover the Q factor from the QR factorization computed by
*>          DGEQP3.
*>
*>          If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
*>          WORK(1) returns the optimal LWORK, and
*>          WORK(2) returns the minimal LWORK.
*> \endverbatim
*>
*> \param[in,out] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK. It is determined as follows:
*>          Let  LWQP3 = 3*N+1,  LWCON = 3*N, and let
*>          LWORQ = { MAX( N, 1 ),  if JOBU = 'R', 'S', or 'U'
*>                  { MAX( M, 1 ),  if JOBU = 'A'
*>          LWSVD = MAX( 5*N, 1 )
*>          LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ),
*>          LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
*>          Then the minimal value of LWORK is:
*>          = MAX( N + LWQP3, LWSVD )        if only the singular values are needed;
*>          = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
*>                                   and a scaled condition estimate requested;
*>
*>          = N + MAX( LWQP3, LWSVD, LWORQ ) if the singular values and the left
*>                                   singular vectors are requested;
*>          = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular values and the left
*>                                   singular vectors are requested, and also
*>                                   a scaled condition estimate requested;
*>
*>          = N + MAX( LWQP3, LWSVD )        if the singular values and the right
*>                                   singular vectors are requested;
*>          = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
*>                                   singular vectors are requested, and also
*>                                   a scaled condition etimate requested;
*>
*>          = N + MAX( LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV = 'R';
*>                                   independent of JOBR;
*>          = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested,
*>                                   JOBV = 'R' and, also a scaled condition
*>                                   estimate requested; independent of JOBR;
*>          = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
*>         N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the
*>                         full SVD is requested with JOBV = 'A' or 'V', and
*>                         JOBR ='N'
*>          = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
*>         N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) )
*>                         if the full SVD is requested with JOBV = 'A' or 'V', and
*>                         JOBR ='N', and also a scaled condition number estimate
*>                         requested.
*>          = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
*>         N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the
*>                         full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
*>          = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
*>         N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
*>                         if the full SVD is requested with JOBV = 'A' or 'V', and
*>                         JOBR ='T', and also a scaled condition number estimate
*>                         requested.
*>          Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ).
*>
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates and returns the optimal and minimal sizes
*>          for the WORK, IWORK, and RWORK arrays, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension (max(1, LRWORK)).
*>          On exit,
*>          1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
*>          number of column scaled A. If A = C * D where D is diagonal and C
*>          has unit columns in the Euclidean norm, then, assuming full column rank,
*>          N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
*>          Otherwise, RWORK(1) = -1.
*>          2. RWORK(2) contains the number of singular values computed as
*>          exact zeros in DGESVD applied to the upper triangular or trapezoidal
*>          R (from the initial QR factorization). In case of early exit (no call to
*>          DGESVD, such as in the case of zero matrix) RWORK(2) = -1.
*>
*>          If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
*>          RWORK(1) returns the minimal LRWORK.
*> \endverbatim
*>
*> \param[in] LRWORK
*> \verbatim
*>          LRWORK is INTEGER.
*>          The dimension of the array RWORK.
*>          If JOBP ='P', then LRWORK >= MAX(2, M).
*>          Otherwise, LRWORK >= 2
*>
*>          If LRWORK = -1, then a workspace query is assumed; the routine
*>          only calculates and returns the optimal and minimal sizes
*>          for the WORK, IWORK, and RWORK arrays, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*>          > 0:  if DBDSQR did not converge, INFO specifies how many superdiagonals
*>          of an intermediate bidiagonal form B (computed in DGESVD) did not
*>          converge to zero.
*> \endverbatim
*
*> \par Further Details:
*  ========================
*>
*> \verbatim
*>
*>   1. The data movement (matrix transpose) is coded using simple nested
*>   DO-loops because BLAS and LAPACK do not provide corresponding subroutines.
*>   Those DO-loops are easily identified in this source code - by the CONTINUE
*>   statements labeled with 11**. In an optimized version of this code, the
*>   nested DO loops should be replaced with calls to an optimized subroutine.
*>   2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause
*>   column norm overflow. This is the minial precaution and it is left to the
*>   SVD routine (CGESVD) to do its own preemptive scaling if potential over-
*>   or underflows are detected. To avoid repeated scanning of the array A,
*>   an optimal implementation would do all necessary scaling before calling
*>   CGESVD and the scaling in CGESVD can be switched off.
*>   3. Other comments related to code optimization are given in comments in the
*>   code, enclosed in [[double brackets]].
*> \endverbatim
*
*> \par Bugs, examples and comments
*  ===========================
*
*> \verbatim
*>  Please report all bugs and send interesting examples and/or comments to
*>  drmac@math.hr. Thank you.
*> \endverbatim
*
*> \par References
*  ===============
*
*> \verbatim
*>  [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
*>      Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
*>      44(1): 11:1-11:30 (2017)
*>
*>  SIGMA library, xGESVDQ section updated February 2016.
*>  Developed and coded by Zlatko Drmac, Department of Mathematics
*>  University of Zagreb, Croatia, drmac@math.hr
*> \endverbatim
*
*
*> \par Contributors:
*  ==================
*>
*> \verbatim
*> Developed and coded by Zlatko Drmac, Department of Mathematics
*>  University of Zagreb, Croatia, drmac@math.hr
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup gesvdq
*
*  =====================================================================
      SUBROUTINE dgesvdq( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
     $                    S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
     $                    WORK, LWORK, RWORK, LRWORK, INFO )
*     .. Scalar Arguments ..
      IMPLICIT    NONE
      CHARACTER   JOBA, JOBP, JOBR, JOBU, JOBV
      INTEGER     M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
     $            info
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
      DOUBLE PRECISION S( * ), RWORK( * )
      INTEGER          IWORK( * )
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO,         ONE
      PARAMETER      ( ZERO = 0.0d0, one = 1.0d0 )
*     .. Local Scalars ..
      INTEGER     IERR, IWOFF, NR, N1, OPTRATIO, p, q
      INTEGER     LWCON, LWQP3, LWRK_DGELQF, LWRK_DGESVD, LWRK_DGESVD2,
     $            lwrk_dgeqp3,  lwrk_dgeqrf, lwrk_dormlq, lwrk_dormqr,
     $            lwrk_dormqr2, lwlqf, lwqrf, lwsvd, lwsvd2, lworq,
     $            lworq2, lworlq, minwrk, minwrk2, optwrk, optwrk2,
     $            iminwrk, rminwrk
      LOGICAL     ACCLA,  ACCLM, ACCLH, ASCALED, CONDA, DNTWU,  DNTWV,
     $            LQUERY, LSVC0, LSVEC, ROWPRM,  RSVEC, RTRANS, WNTUA,
     $            wntuf,  wntur, wntus, wntva,   wntvr
      DOUBLE PRECISION BIG, EPSLN, RTMP, SCONDA, SFMIN
*     .. Local Arrays
      DOUBLE PRECISION RDUMMY(1)
*     ..
*     .. External Subroutines (BLAS, LAPACK)
      EXTERNAL    dgelqf, dgeqp3, dgeqrf, dgesvd, dlacpy,
     $                   dlapmt,
     $            dlascl, dlaset, dlaswp, dscal,  dpocon, dormlq,
     $            dormqr, xerbla
*     ..
*     .. External Functions (BLAS, LAPACK)
      LOGICAL    LSAME
      INTEGER    IDAMAX
      DOUBLE PRECISION  DLANGE, DNRM2, DLAMCH
      EXTERNAL    dlange, lsame, idamax, dnrm2, dlamch
*     ..
*     .. Intrinsic Functions ..
*
      INTRINSIC   abs, max, min, dble, sqrt
*
*     Test the input arguments
*
      wntus  = lsame( jobu, 'S' ) .OR. lsame( jobu, 'U' )
      wntur  = lsame( jobu, 'R' )
      wntua  = lsame( jobu, 'A' )
      wntuf  = lsame( jobu, 'F' )
      lsvc0  = wntus .OR. wntur .OR. wntua
      lsvec  = lsvc0 .OR. wntuf
      dntwu  = lsame( jobu, 'N' )
*
      wntvr  = lsame( jobv, 'R' )
      wntva  = lsame( jobv, 'A' ) .OR. lsame( jobv, 'V' )
      rsvec  = wntvr .OR. wntva
      dntwv  = lsame( jobv, 'N' )
*
      accla  = lsame( joba, 'A' )
      acclm  = lsame( joba, 'M' )
      conda  = lsame( joba, 'E' )
      acclh  = lsame( joba, 'H' ) .OR. conda
*
      rowprm = lsame( jobp, 'P' )
      rtrans = lsame( jobr, 'T' )
*
      IF ( rowprm ) THEN
         IF ( conda ) THEN
            iminwrk = max( 1, n + m - 1 + n )
         ELSE
            iminwrk = max( 1, n + m - 1 )
         END IF
         rminwrk = max( 2, m )
      ELSE
         IF ( conda ) THEN
            iminwrk = max( 1, n + n )
         ELSE
            iminwrk = max( 1, n )
         END IF
         rminwrk = 2
      END IF
      lquery = (liwork .EQ. -1 .OR. lwork .EQ. -1 .OR. lrwork .EQ. -1)
      info  = 0
      IF ( .NOT. ( accla .OR. acclm .OR. acclh ) ) THEN
         info = -1
      ELSE IF ( .NOT.( rowprm .OR. lsame( jobp, 'N' ) ) ) THEN
          info = -2
      ELSE IF ( .NOT.( rtrans .OR. lsame( jobr, 'N' ) ) ) THEN
          info = -3
      ELSE IF ( .NOT.( lsvec .OR. dntwu ) ) THEN
         info = -4
      ELSE IF ( wntur .AND. wntva ) THEN
         info = -5
      ELSE IF ( .NOT.( rsvec .OR. dntwv )) THEN
         info = -5
      ELSE IF ( m.LT.0 ) THEN
         info = -6
      ELSE IF ( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
         info = -7
      ELSE IF ( lda.LT.max( 1, m ) ) THEN
         info = -9
      ELSE IF ( ldu.LT.1 .OR. ( lsvc0 .AND. ldu.LT.m ) .OR.
     $       ( wntuf .AND. ldu.LT.n ) ) THEN
         info = -12
      ELSE IF ( ldv.LT.1 .OR. ( rsvec .AND. ldv.LT.n ) .OR.
     $          ( conda .AND. ldv.LT.n ) ) THEN
         info = -14
      ELSE IF ( liwork .LT. iminwrk .AND. .NOT. lquery ) THEN
         info = -17
      END IF
*
*
      IF ( info .EQ. 0 ) THEN
*        .. compute the minimal and the optimal workspace lengths
*        [[The expressions for computing the minimal and the optimal
*        values of LWORK are written with a lot of redundancy and
*        can be simplified. However, this detailed form is easier for
*        maintenance and modifications of the code.]]
*
*        .. minimal workspace length for DGEQP3 of an M x N matrix
         lwqp3 = 3 * n + 1
*        .. minimal workspace length for DORMQR to build left singular vectors
         IF ( wntus .OR. wntur ) THEN
             lworq  = max( n  , 1 )
         ELSE IF ( wntua ) THEN
             lworq = max( m , 1 )
         END IF
*        .. minimal workspace length for DPOCON of an N x N matrix
         lwcon = 3 * n
*        .. DGESVD of an N x N matrix
         lwsvd = max( 5 * n, 1 )
         IF ( lquery ) THEN
             CALL dgeqp3( m, n, a, lda, iwork, rdummy, rdummy, -1,
     $           ierr )
             lwrk_dgeqp3 = int( rdummy(1) )
             IF ( wntus .OR. wntur ) THEN
                 CALL dormqr( 'L', 'N', m, n, n, a, lda, rdummy, u,
     $                ldu, rdummy, -1, ierr )
                 lwrk_dormqr = int( rdummy(1) )
             ELSE IF ( wntua ) THEN
                 CALL dormqr( 'L', 'N', m, m, n, a, lda, rdummy, u,
     $                ldu, rdummy, -1, ierr )
                 lwrk_dormqr = int( rdummy(1) )
             ELSE
                 lwrk_dormqr = 0
             END IF
         END IF
         minwrk = 2
         optwrk = 2
         IF ( .NOT. (lsvec .OR. rsvec )) THEN
*            .. minimal and optimal sizes of the workspace if
*            only the singular values are requested
             IF ( conda ) THEN
                minwrk = max( n+lwqp3, lwcon, lwsvd )
             ELSE
                minwrk = max( n+lwqp3, lwsvd )
             END IF
             IF ( lquery ) THEN
                 CALL dgesvd( 'N', 'N', n, n, a, lda, s, u, ldu,
     $                v, ldv, rdummy, -1, ierr )
                 lwrk_dgesvd = int( rdummy(1) )
                 IF ( conda ) THEN
                    optwrk = max( n+lwrk_dgeqp3, n+lwcon, lwrk_dgesvd )
                 ELSE
                    optwrk = max( n+lwrk_dgeqp3, lwrk_dgesvd )
                 END IF
             END IF
         ELSE IF ( lsvec .AND. (.NOT.rsvec) ) THEN
*            .. minimal and optimal sizes of the workspace if the
*            singular values and the left singular vectors are requested
             IF ( conda ) THEN
                 minwrk = n + max( lwqp3, lwcon, lwsvd, lworq )
             ELSE
                 minwrk = n + max( lwqp3, lwsvd, lworq )
             END IF
             IF ( lquery ) THEN
                IF ( rtrans ) THEN
                   CALL dgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
     $                  v, ldv, rdummy, -1, ierr )
                ELSE
                   CALL dgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
     $                  v, ldv, rdummy, -1, ierr )
                END IF
                lwrk_dgesvd = int( rdummy(1) )
                IF ( conda ) THEN
                    optwrk = n + max( lwrk_dgeqp3, lwcon, lwrk_dgesvd,
     $                               lwrk_dormqr )
                ELSE
                    optwrk = n + max( lwrk_dgeqp3, lwrk_dgesvd,
     $                               lwrk_dormqr )
                END IF
             END IF
         ELSE IF ( rsvec .AND. (.NOT.lsvec) ) THEN
*            .. minimal and optimal sizes of the workspace if the
*            singular values and the right singular vectors are requested
             IF ( conda ) THEN
                 minwrk = n + max( lwqp3, lwcon, lwsvd )
             ELSE
                 minwrk = n + max( lwqp3, lwsvd )
             END IF
             IF ( lquery ) THEN
                 IF ( rtrans ) THEN
                     CALL dgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
     $                    v, ldv, rdummy, -1, ierr )
                 ELSE
                     CALL dgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
     $                    v, ldv, rdummy, -1, ierr )
                 END IF
                 lwrk_dgesvd = int( rdummy(1) )
                 IF ( conda ) THEN
                     optwrk = n + max( lwrk_dgeqp3, lwcon, lwrk_dgesvd )
                 ELSE
                     optwrk = n + max( lwrk_dgeqp3, lwrk_dgesvd )
                 END IF
             END IF
         ELSE
*            .. minimal and optimal sizes of the workspace if the
*            full SVD is requested
             IF ( rtrans ) THEN
                 minwrk = max( lwqp3, lwsvd, lworq )
                 IF ( conda ) minwrk = max( minwrk, lwcon )
                 minwrk = minwrk + n
                 IF ( wntva ) THEN
*                   .. minimal workspace length for N x N/2 DGEQRF
                    lwqrf  = max( n/2, 1 )
*                   .. minimal workspace length for N/2 x N/2 DGESVD
                    lwsvd2 = max( 5 * (n/2), 1 )
                    lworq2 = max( n, 1 )
                    minwrk2 = max( lwqp3, n/2+lwqrf, n/2+lwsvd2,
     $                        n/2+lworq2, lworq )
                    IF ( conda ) minwrk2 = max( minwrk2, lwcon )
                    minwrk2 = n + minwrk2
                    minwrk = max( minwrk, minwrk2 )
                 END IF
             ELSE
                 minwrk = max( lwqp3, lwsvd, lworq )
                 IF ( conda ) minwrk = max( minwrk, lwcon )
                 minwrk = minwrk + n
                 IF ( wntva ) THEN
*                   .. minimal workspace length for N/2 x N DGELQF
                    lwlqf  = max( n/2, 1 )
                    lwsvd2 = max( 5 * (n/2), 1 )
                    lworlq = max( n , 1 )
                    minwrk2 = max( lwqp3, n/2+lwlqf, n/2+lwsvd2,
     $                        n/2+lworlq, lworq )
                    IF ( conda ) minwrk2 = max( minwrk2, lwcon )
                    minwrk2 = n + minwrk2
                    minwrk = max( minwrk, minwrk2 )
                 END IF
             END IF
             IF ( lquery ) THEN
                IF ( rtrans ) THEN
                   CALL dgesvd( 'O', 'A', n, n, a, lda, s, u, ldu,
     $                  v, ldv, rdummy, -1, ierr )
                   lwrk_dgesvd = int( rdummy(1) )
                   optwrk = max(lwrk_dgeqp3,lwrk_dgesvd,lwrk_dormqr)
                   IF ( conda ) optwrk = max( optwrk, lwcon )
                   optwrk = n + optwrk
                   IF ( wntva ) THEN
                       CALL dgeqrf(n,n/2,u,ldu,rdummy,rdummy,-1,ierr)
                       lwrk_dgeqrf = int( rdummy(1) )
                       CALL dgesvd( 'S', 'O', n/2,n/2, v,ldv, s, u,
     $                              ldu,
     $                      v, ldv, rdummy, -1, ierr )
                       lwrk_dgesvd2 = int( rdummy(1) )
                       CALL dormqr( 'R', 'C', n, n, n/2, u, ldu,
     $                              rdummy,
     $                      v, ldv, rdummy, -1, ierr )
                       lwrk_dormqr2 = int( rdummy(1) )
                       optwrk2 = max( lwrk_dgeqp3, n/2+lwrk_dgeqrf,
     $                           n/2+lwrk_dgesvd2, n/2+lwrk_dormqr2 )
                       IF ( conda ) optwrk2 = max( optwrk2, lwcon )
                       optwrk2 = n + optwrk2
                       optwrk = max( optwrk, optwrk2 )
                   END IF
                ELSE
                   CALL dgesvd( 'S', 'O', n, n, a, lda, s, u, ldu,
     $                  v, ldv, rdummy, -1, ierr )
                   lwrk_dgesvd = int( rdummy(1) )
                   optwrk = max(lwrk_dgeqp3,lwrk_dgesvd,lwrk_dormqr)
                   IF ( conda ) optwrk = max( optwrk, lwcon )
                   optwrk = n + optwrk
                   IF ( wntva ) THEN
                      CALL dgelqf(n/2,n,u,ldu,rdummy,rdummy,-1,ierr)
                      lwrk_dgelqf = int( rdummy(1) )
                      CALL dgesvd( 'S','O', n/2,n/2, v, ldv, s, u,
     $                             ldu,
     $                     v, ldv, rdummy, -1, ierr )
                      lwrk_dgesvd2 = int( rdummy(1) )
                      CALL dormlq( 'R', 'N', n, n, n/2, u, ldu,
     $                             rdummy,
     $                     v, ldv, rdummy,-1,ierr )
                      lwrk_dormlq = int( rdummy(1) )
                      optwrk2 = max( lwrk_dgeqp3, n/2+lwrk_dgelqf,
     $                           n/2+lwrk_dgesvd2, n/2+lwrk_dormlq )
                       IF ( conda ) optwrk2 = max( optwrk2, lwcon )
                       optwrk2 = n + optwrk2
                       optwrk = max( optwrk, optwrk2 )
                   END IF
                END IF
             END IF
         END IF
*
         minwrk = max( 2, minwrk )
         optwrk = max( 2, optwrk )
         IF ( lwork .LT. minwrk .AND. (.NOT.lquery) ) info = -19
*
      END IF
*
      IF (info .EQ. 0 .AND. lrwork .LT. rminwrk .AND. .NOT. lquery) THEN
         info = -21
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DGESVDQ', -info )
         RETURN
      ELSE IF ( lquery ) THEN
*
*     Return optimal workspace
*
          iwork(1) = iminwrk
          work(1) = optwrk
          work(2) = minwrk
          rwork(1) = rminwrk
          RETURN
      END IF
*
*     Quick return if the matrix is void.
*
      IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) ) THEN
*     .. all output is void.
         RETURN
      END IF
*
      big = dlamch('O')
      ascaled = .false.
      iwoff = 1
      IF ( rowprm ) THEN
            iwoff = m
*           .. reordering the rows in decreasing sequence in the
*           ell-infinity norm - this enhances numerical robustness in
*           the case of differently scaled rows.
            DO 1904 p = 1, m
*               RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
*               [[DLANGE will return NaN if an entry of the p-th row is Nan]]
                rwork(p) = dlange( 'M', 1, n, a(p,1), lda, rdummy )
*               .. check for NaN's and Inf's
                IF ( ( rwork(p) .NE. rwork(p) ) .OR.
     $               ( (rwork(p)*zero) .NE. zero ) ) THEN
                    info = -8
                    CALL xerbla( 'DGESVDQ', -info )
                    RETURN
                END IF
 1904       CONTINUE
            DO 1952 p = 1, m - 1
            q = idamax( m-p+1, rwork(p), 1 ) + p - 1
            iwork(n+p) = q
            IF ( p .NE. q ) THEN
               rtmp     = rwork(p)
               rwork(p) = rwork(q)
               rwork(q) = rtmp
            END IF
 1952       CONTINUE
*
            IF ( rwork(1) .EQ. zero ) THEN
*              Quick return: A is the M x N zero matrix.
               numrank = 0
               CALL dlaset( 'G', n, 1, zero, zero, s, n )
               IF ( wntus ) CALL dlaset('G', m, n, zero, one, u, ldu)
               IF ( wntua ) CALL dlaset('G', m, m, zero, one, u, ldu)
               IF ( wntva ) CALL dlaset('G', n, n, zero, one, v, ldv)
               IF ( wntuf ) THEN
                   CALL dlaset( 'G', n, 1, zero, zero, work, n )
                   CALL dlaset( 'G', m, n, zero,  one, u, ldu )
               END IF
               DO 5001 p = 1, n
                   iwork(p) = p
 5001          CONTINUE
               IF ( rowprm ) THEN
                   DO 5002 p = n + 1, n + m - 1
                       iwork(p) = p - n
 5002              CONTINUE
               END IF
               IF ( conda ) rwork(1) = -1
               rwork(2) = -1
               RETURN
            END IF
*
            IF ( rwork(1) .GT. big / sqrt(dble(m)) ) THEN
*               .. to prevent overflow in the QR factorization, scale the
*               matrix by 1/sqrt(M) if too large entry detected
                CALL dlascl('G',0,0,sqrt(dble(m)),one, m,n, a,lda,
     $                       ierr)
                ascaled = .true.
            END IF
            CALL dlaswp( n, a, lda, 1, m-1, iwork(n+1), 1 )
      END IF
*
*    .. At this stage, preemptive scaling is done only to avoid column
*    norms overflows during the QR factorization. The SVD procedure should
*    have its own scaling to save the singular values from overflows and
*    underflows. That depends on the SVD procedure.
*
      IF ( .NOT.rowprm ) THEN
          rtmp = dlange( 'M', m, n, a, lda, rdummy )
          IF ( ( rtmp .NE. rtmp ) .OR.
     $         ( (rtmp*zero) .NE. zero ) ) THEN
               info = -8
               CALL xerbla( 'DGESVDQ', -info )
               RETURN
          END IF
          IF ( rtmp .GT. big / sqrt(dble(m)) ) THEN
*             .. to prevent overflow in the QR factorization, scale the
*             matrix by 1/sqrt(M) if too large entry detected
              CALL dlascl('G',0,0, sqrt(dble(m)),one, m,n, a,lda,
     $                     ierr)
              ascaled = .true.
          END IF
      END IF
*
*     .. QR factorization with column pivoting
*
*     A * P = Q * [ R ]
*                 [ 0 ]
*
      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 )
*
*    If the user requested accuracy level allows truncation in the
*    computed upper triangular factor, the matrix R is examined and,
*    if possible, replaced with its leading upper trapezoidal part.
*
      epsln = dlamch('E')
      sfmin = dlamch('S')
*     SMALL = SFMIN / EPSLN
      nr = n
*
      IF ( accla ) THEN
*
*        Standard absolute error bound suffices. All sigma_i with
*        sigma_i < N*EPS*||A||_F are flushed to zero. This is an
*        aggressive enforcement of lower numerical rank by introducing a
*        backward error of the order of N*EPS*||A||_F.
         nr = 1
         rtmp = sqrt(dble(n))*epsln
         DO 3001 p = 2, n
            IF ( abs(a(p,p)) .LT. (rtmp*abs(a(1,1))) ) GO TO 3002
               nr = nr + 1
 3001    CONTINUE
 3002    CONTINUE
*
      ELSEIF ( acclm ) THEN
*        .. similarly as above, only slightly more gentle (less aggressive).
*        Sudden drop on the diagonal of R is used as the criterion for being
*        close-to-rank-deficient. The threshold is set to EPSLN=DLAMCH('E').
*        [[This can be made more flexible by replacing this hard-coded value
*        with a user specified threshold.]] Also, the values that underflow
*        will be truncated.
         nr = 1
         DO 3401 p = 2, n
            IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
     $           ( abs(a(p,p)) .LT. sfmin ) ) GO TO 3402
            nr = nr + 1
 3401    CONTINUE
 3402    CONTINUE
*
      ELSE
*        .. RRQR not authorized to determine numerical rank except in the
*        obvious case of zero pivots.
*        .. inspect R for exact zeros on the diagonal;
*        R(i,i)=0 => R(i:N,i:N)=0.
         nr = 1
         DO 3501 p = 2, n
            IF ( abs(a(p,p)) .EQ. zero ) GO TO 3502
            nr = nr + 1
 3501    CONTINUE
 3502    CONTINUE
*
         IF ( conda ) THEN
*           Estimate the scaled condition number of A. Use the fact that it is
*           the same as the scaled condition number of R.
*              .. V is used as workspace
               CALL dlacpy( 'U', n, n, a, lda, v, ldv )
*              Only the leading NR x NR submatrix of the triangular factor
*              is considered. Only if NR=N will this give a reliable error
*              bound. However, even for NR < N, this can be used on an
*              expert level and obtain useful information in the sense of
*              perturbation theory.
               DO 3053 p = 1, nr
                  rtmp = dnrm2( p, v(1,p), 1 )
                  CALL dscal( p, one/rtmp, v(1,p), 1 )
 3053          CONTINUE
               IF ( .NOT. ( lsvec .OR. rsvec ) ) THEN
                   CALL dpocon( 'U', nr, v, ldv, one, rtmp,
     $                  work, iwork(n+iwoff), ierr )
               ELSE
                   CALL dpocon( 'U', nr, v, ldv, one, rtmp,
     $                  work(n+1), iwork(n+iwoff), ierr )
               END IF
               sconda = one / sqrt(rtmp)
*           For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
*           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
*           See the reference [1] for more details.
         END IF
*
      ENDIF
*
      IF ( wntur ) THEN
          n1 = nr
      ELSE IF ( wntus .OR. wntuf) THEN
          n1 = n
      ELSE IF ( wntua ) THEN
          n1 = m
      END IF
*
      IF ( .NOT. ( rsvec .OR. lsvec ) ) THEN
*.......................................................................
*        .. only the singular values are requested
*.......................................................................
         IF ( rtrans ) THEN
*
*         .. compute the singular values of R**T = [A](1:NR,1:N)**T
*           .. set the lower triangle of [A] to [A](1:NR,1:N)**T and
*           the upper triangle of [A] to zero.
            DO 1146 p = 1, min( n, nr )
               DO 1147 q = p + 1, n
                  a(q,p) = a(p,q)
                  IF ( q .LE. nr ) a(p,q) = zero
 1147          CONTINUE
 1146       CONTINUE
*
            CALL dgesvd( 'N', 'N', n, nr, a, lda, s, u, ldu,
     $           v, ldv, work, lwork, info )
*
         ELSE
*
*           .. compute the singular values of R = [A](1:NR,1:N)
*
            IF ( nr .GT. 1 )
     $          CALL dlaset( 'L', nr-1,nr-1, zero,zero, a(2,1), lda )
            CALL dgesvd( 'N', 'N', nr, n, a, lda, s, u, ldu,
     $           v, ldv, work, lwork, info )
*
         END IF
*
      ELSE IF ( lsvec .AND. ( .NOT. rsvec) ) THEN
*.......................................................................
*       .. the singular values and the left singular vectors requested
*.......................................................................""""""""
         IF ( rtrans ) THEN
*            .. apply DGESVD to R**T
*            .. copy R**T into [U] and overwrite [U] with the right singular
*            vectors of R
            DO 1192 p = 1, nr
               DO 1193 q = p, n
                  u(q,p) = a(p,q)
 1193          CONTINUE
 1192       CONTINUE
            IF ( nr .GT. 1 )
     $          CALL dlaset( 'U', nr-1,nr-1, zero,zero, u(1,2), ldu )
*           .. the left singular vectors not computed, the NR right singular
*           vectors overwrite [U](1:NR,1:NR) as transposed. These
*           will be pre-multiplied by Q to build the left singular vectors of A.
               CALL dgesvd( 'N', 'O', n, nr, u, ldu, s, u, ldu,
     $              u, ldu, work(n+1), lwork-n, info )
*
               DO 1119 p = 1, nr
                   DO 1120 q = p + 1, nr
                      rtmp   = u(q,p)
                      u(q,p) = u(p,q)
                      u(p,q) = rtmp
 1120              CONTINUE
 1119          CONTINUE
*
         ELSE
*            .. apply DGESVD to R
*            .. copy R into [U] and overwrite [U] with the left singular vectors
             CALL dlacpy( 'U', nr, n, a, lda, u, ldu )
             IF ( nr .GT. 1 )
     $         CALL dlaset( 'L', nr-1, nr-1, zero, zero, u(2,1),
     $                      ldu )
*            .. the right singular vectors not computed, the NR left singular
*            vectors overwrite [U](1:NR,1:NR)
                CALL dgesvd( 'O', 'N', nr, n, u, ldu, s, u, ldu,
     $               v, ldv, work(n+1), lwork-n, info )
*               .. now [U](1:NR,1:NR) contains the NR left singular vectors of
*               R. These will be pre-multiplied by Q to build the left singular
*               vectors of A.
         END IF
*
*           .. assemble the left singular vector matrix U of dimensions
*              (M x NR) or (M x N) or (M x M).
         IF ( ( nr .LT. m ) .AND. ( .NOT.wntuf ) ) 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
*           vectors matrix U.
*
         IF ( .NOT.wntuf )
     $       CALL dormqr( 'L', 'N', m, n1, n, a, lda, work, u,
     $            ldu, work(n+1), lwork-n, ierr )
         IF ( rowprm .AND. .NOT.wntuf )
     $          CALL dlaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
*
      ELSE IF ( rsvec .AND. ( .NOT. lsvec ) ) THEN
*.......................................................................
*       .. the singular values and the right singular vectors requested
*.......................................................................
          IF ( rtrans ) THEN
*            .. apply DGESVD to R**T
*            .. copy R**T into V and overwrite V with the left singular vectors
            DO 1165 p = 1, nr
               DO 1166 q = p, n
                  v(q,p) = (a(p,q))
 1166          CONTINUE
 1165       CONTINUE
            IF ( nr .GT. 1 )
     $          CALL dlaset( 'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
*           .. the left singular vectors of R**T overwrite V, the right singular
*           vectors not computed
            IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
               CALL dgesvd( 'O', 'N', n, nr, v, ldv, s, u, ldu,
     $              u, ldu, work(n+1), lwork-n, info )
*
               DO 1121 p = 1, nr
                   DO 1122 q = p + 1, nr
                      rtmp   = v(q,p)
                      v(q,p) = v(p,q)
                      v(p,q) = rtmp
 1122              CONTINUE
 1121          CONTINUE
*
               IF ( nr .LT. n ) THEN
                   DO 1103 p = 1, nr
                      DO 1104 q = nr + 1, n
                          v(p,q) = v(q,p)
 1104                 CONTINUE
 1103              CONTINUE
               END IF
               CALL dlapmt( .false., nr, n, v, ldv, iwork )
            ELSE
*               .. need all N right singular vectors and NR < N
*               [!] This is simple implementation that augments [V](1:N,1:NR)
*               by padding a zero block. In the case NR << N, a more efficient
*               way is to first use the QR factorization. For more details
*               how to implement this, see the " FULL SVD " branch.
                CALL dlaset('G', n, n-nr, zero, zero, v(1,nr+1), ldv)
                CALL dgesvd( 'O', 'N', n, n, v, ldv, s, u, ldu,
     $               u, ldu, work(n+1), lwork-n, info )
*
                DO 1123 p = 1, n
                   DO 1124 q = p + 1, n
                      rtmp   = v(q,p)
                      v(q,p) = v(p,q)
                      v(p,q) = rtmp
 1124              CONTINUE
 1123           CONTINUE
                CALL dlapmt( .false., n, n, v, ldv, iwork )
            END IF
*
          ELSE
*            .. aply DGESVD to R
*            .. copy R into V and overwrite V with the right singular vectors
             CALL dlacpy( 'U', nr, n, a, lda, v, ldv )
             IF ( nr .GT. 1 )
     $         CALL dlaset( 'L', nr-1, nr-1, zero, zero, v(2,1),
     $                      ldv )
*            .. the right singular vectors overwrite V, the NR left singular
*            vectors stored in U(1:NR,1:NR)
             IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
                CALL dgesvd( 'N', 'O', nr, n, v, ldv, s, u, ldu,
     $               v, ldv, work(n+1), lwork-n, info )
                CALL dlapmt( .false., nr, n, v, ldv, iwork )
*               .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
             ELSE
*               .. need all N right singular vectors and NR < N
*               [!] This is simple implementation that augments [V](1:NR,1:N)
*               by padding a zero block. In the case NR << N, a more efficient
*               way is to first use the LQ factorization. For more details
*               how to implement this, see the " FULL SVD " branch.
                 CALL dlaset('G', n-nr, n, zero,zero, v(nr+1,1), ldv)
                 CALL dgesvd( 'N', 'O', n, n, v, ldv, s, u, ldu,
     $                v, ldv, work(n+1), lwork-n, info )
                 CALL dlapmt( .false., n, n, v, ldv, iwork )
             END IF
*            .. now [V] contains the transposed matrix of the right singular
*            vectors of A.
          END IF
*
      ELSE
*.......................................................................
*       .. FULL SVD requested
*.......................................................................
         IF ( rtrans ) THEN
*
*            .. apply DGESVD to R**T [[this option is left for R&D&T]]
*
            IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
*            .. copy R**T into [V] and overwrite [V] with the left singular
*            vectors of R**T
            DO 1168 p = 1, nr
               DO 1169 q = p, n
                  v(q,p) = a(p,q)
 1169          CONTINUE
 1168       CONTINUE
            IF ( nr .GT. 1 )
     $          CALL dlaset( 'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
*
*           .. the left singular vectors of R**T overwrite [V], the NR right
*           singular vectors of R**T stored in [U](1:NR,1:NR) as transposed
               CALL dgesvd( 'O', 'A', n, nr, v, ldv, s, v, ldv,
     $              u, ldu, work(n+1), lwork-n, info )
*              .. assemble V
               DO 1115 p = 1, nr
                  DO 1116 q = p + 1, nr
                     rtmp   = v(q,p)
                     v(q,p) = v(p,q)
                     v(p,q) = rtmp
 1116             CONTINUE
 1115          CONTINUE
               IF ( nr .LT. n ) THEN
                   DO 1101 p = 1, nr
                      DO 1102 q = nr+1, n
                         v(p,q) = v(q,p)
 1102                 CONTINUE
 1101              CONTINUE
               END IF
               CALL dlapmt( .false., nr, n, v, ldv, iwork )
*
                DO 1117 p = 1, nr
                   DO 1118 q = p + 1, nr
                      rtmp   = u(q,p)
                      u(q,p) = u(p,q)
                      u(p,q) = rtmp
 1118              CONTINUE
 1117           CONTINUE
*
                IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) 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
*
            ELSE
*               .. need all N right singular vectors and NR < N
*            .. copy R**T into [V] and overwrite [V] with the left singular
*            vectors of R**T
*               [[The optimal ratio N/NR for using QRF instead of padding
*                 with zeros. Here hard coded to 2; it must be at least
*                 two due to work space constraints.]]
*               OPTRATIO = ILAENV(6, 'DGESVD', 'S' // 'O', NR,N,0,0)
*               OPTRATIO = MAX( OPTRATIO, 2 )
                optratio = 2
                IF ( optratio*nr .GT. n ) THEN
                   DO 1198 p = 1, nr
                      DO 1199 q = p, n
                         v(q,p) = a(p,q)
 1199                 CONTINUE
 1198              CONTINUE
                   IF ( nr .GT. 1 )
     $             CALL dlaset('U',nr-1,nr-1, zero,zero, v(1,2),ldv)
*
                   CALL dlaset('A',n,n-nr,zero,zero,v(1,nr+1),ldv)
                   CALL dgesvd( 'O', 'A', n, n, v, ldv, s, v, ldv,
     $                  u, ldu, work(n+1), lwork-n, info )
*
                   DO 1113 p = 1, n
                      DO 1114 q = p + 1, n
                         rtmp   = v(q,p)
                         v(q,p) = v(p,q)
                         v(p,q) = rtmp
 1114                 CONTINUE
 1113              CONTINUE
                   CALL dlapmt( .false., n, n, v, ldv, iwork )
*              .. assemble the left singular vector matrix U of dimensions
*              (M x N1), i.e. (M x N) or (M x M).
*
                   DO 1111 p = 1, n
                      DO 1112 q = p + 1, n
                         rtmp   = u(q,p)
                         u(q,p) = u(p,q)
                         u(p,q) = rtmp
 1112                 CONTINUE
 1111              CONTINUE
*
                   IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) 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
                ELSE
*                  .. copy R**T into [U] and overwrite [U] with the right
*                  singular vectors of R
                   DO 1196 p = 1, nr
                      DO 1197 q = p, n
                         u(q,nr+p) = a(p,q)
 1197                 CONTINUE
 1196              CONTINUE
                   IF ( nr .GT. 1 )
     $             CALL dlaset('U',nr-1,nr-1,zero,zero,u(1,nr+2),ldu)
                   CALL dgeqrf( n, nr, u(1,nr+1), ldu, work(n+1),
     $                  work(n+nr+1), lwork-n-nr, ierr )
                   DO 1143 p = 1, nr
                       DO 1144 q = 1, n
                           v(q,p) = u(p,nr+q)
 1144                  CONTINUE
 1143              CONTINUE
                  CALL dlaset('U',nr-1,nr-1,zero,zero,v(1,2),ldv)
                  CALL dgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
     $                 v,ldv, work(n+nr+1),lwork-n-nr, info )
                  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)
                  CALL dormqr('R','C', n, n, nr, u(1,nr+1), ldu,
     $                 work(n+1),v,ldv,work(n+nr+1),lwork-n-nr,ierr)
                  CALL dlapmt( .false., n, n, v, ldv, iwork )
*                 .. assemble the left singular vector matrix U of dimensions
*                 (M x NR) or (M x N) or (M x M).
                  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) 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
                END IF
            END IF
*
         ELSE
*
*            .. apply DGESVD to R [[this is the recommended option]]
*
             IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
*                .. copy R into [V] and overwrite V with the right singular vectors
                 CALL dlacpy( 'U', nr, n, a, lda, v, ldv )
                IF ( nr .GT. 1 )
     $          CALL dlaset( 'L', nr-1,nr-1, zero,zero, v(2,1), ldv )
*               .. the right singular vectors of R overwrite [V], the NR left
*               singular vectors of R stored in [U](1:NR,1:NR)
                CALL dgesvd( 'S', 'O', nr, n, v, ldv, s, u, ldu,
     $               v, ldv, work(n+1), lwork-n, info )
                CALL dlapmt( .false., nr, n, v, ldv, iwork )
*               .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
*               .. assemble the left singular vector matrix U of dimensions
*              (M x NR) or (M x N) or (M x M).
               IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) 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
*
             ELSE
*              .. need all N right singular vectors and NR < N
*              .. the requested number of the left singular vectors
*               is then N1 (N or M)
*               [[The optimal ratio N/NR for using LQ instead of padding
*                 with zeros. Here hard coded to 2; it must be at least
*                 two due to work space constraints.]]
*               OPTRATIO = ILAENV(6, 'DGESVD', 'S' // 'O', NR,N,0,0)
*               OPTRATIO = MAX( OPTRATIO, 2 )
               optratio = 2
               IF ( optratio * nr .GT. n ) THEN
                  CALL dlacpy( 'U', nr, n, a, lda, v, ldv )
                  IF ( nr .GT. 1 )
     $            CALL dlaset('L', nr-1,nr-1, zero,zero, v(2,1),ldv)
*              .. the right singular vectors of R overwrite [V], the NR left
*                 singular vectors of R stored in [U](1:NR,1:NR)
                  CALL dlaset('A', n-nr,n, zero,zero, v(nr+1,1),ldv)
                  CALL dgesvd( 'S', 'O', n, n, v, ldv, s, u, ldu,
     $                 v, ldv, work(n+1), lwork-n, info )
                  CALL dlapmt( .false., n, n, v, ldv, iwork )
*                 .. now [V] contains the transposed matrix of the right
*                 singular vectors of A. The leading N left singular vectors
*                 are in [U](1:N,1:N)
*                 .. assemble the left singular vector matrix U of dimensions
*                 (M x N1), i.e. (M x N) or (M x M).
                  IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) 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
               ELSE
                  CALL dlacpy( 'U', nr, n, a, lda, u(nr+1,1), ldu )
                  IF ( nr .GT. 1 )
     $            CALL dlaset('L',nr-1,nr-1,zero,zero,u(nr+2,1),ldu)
                  CALL dgelqf( nr, n, u(nr+1,1), ldu, work(n+1),
     $                 work(n+nr+1), lwork-n-nr, ierr )
                  CALL dlacpy('L',nr,nr,u(nr+1,1),ldu,v,ldv)
                  IF ( nr .GT. 1 )
     $            CALL dlaset('U',nr-1,nr-1,zero,zero,v(1,2),ldv)
                  CALL dgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
     $                 v, ldv, work(n+nr+1), lwork-n-nr, info )
                  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)
                  CALL dormlq('R','N',n,n,nr,u(nr+1,1),ldu,work(n+1),
     $                 v, ldv, work(n+nr+1),lwork-n-nr,ierr)
                  CALL dlapmt( .false., n, n, v, ldv, iwork )
*               .. assemble the left singular vector matrix U of dimensions
*              (M x NR) or (M x N) or (M x M).
                  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) 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
               END IF
             END IF
*        .. end of the "R**T or R" branch
         END IF
*
*           The Q matrix from the first QRF is built into the left singular
*           vectors matrix U.
*
         IF ( .NOT. wntuf )
     $       CALL dormqr( 'L', 'N', m, n1, n, a, lda, work, u,
     $            ldu, work(n+1), lwork-n, ierr )
         IF ( rowprm .AND. .NOT.wntuf )
     $          CALL dlaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
*
*     ... end of the "full SVD" branch
      END IF
*
*     Check whether some singular values are returned as zeros, e.g.
*     due to underflow, and update the numerical rank.
      p = nr
      DO 4001 q = p, 1, -1
          IF ( s(q) .GT. zero ) GO TO 4002
          nr = nr - 1
 4001 CONTINUE
 4002 CONTINUE
*
*     .. if numerical rank deficiency is detected, the truncated
*     singular values are set to zero.
      IF ( nr .LT. n ) CALL dlaset( 'G', n-nr,1, zero,zero, s(nr+1),
     $     n )
*     .. undo scaling; this may cause overflow in the largest singular
*     values.
      IF ( ascaled )
     $   CALL dlascl( 'G',0,0, one,sqrt(dble(m)), nr,1, s, n, ierr )
      IF ( conda ) rwork(1) = sconda
      rwork(2) = p - nr
*     .. p-NR is the number of singular values that are computed as
*     exact zeros in DGESVD() applied to the (possibly truncated)
*     full row rank triangular (trapezoidal) factor of A.
      numrank = nr
*
      RETURN
*
*     End of DGESVDQ
*
      END