dgesvd.f

Go to the documentation of this file.

*> \brief <b> DGESVD computes the singular value decomposition (SVD) for GE matrices</b>
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download DGESVD + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
*                          WORK, LWORK, INFO )
* 
*       .. Scalar Arguments ..
*       CHARACTER          JOBU, JOBVT
*       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
*      $                   VT( LDVT, * ), WORK( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGESVD computes the singular value decomposition (SVD) of a real
*> M-by-N matrix A, optionally computing the left and/or right singular
*> vectors. The SVD is written
*>
*>      A = U * SIGMA * transpose(V)
*>
*> where SIGMA is an M-by-N matrix which is zero except for its
*> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
*> V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
*> are the singular values of A; they are real and non-negative, and
*> are returned in descending order.  The first min(m,n) columns of
*> U and V are the left and right singular vectors of A.
*>
*> Note that the routine returns V**T, not V.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOBU
*> \verbatim
*>          JOBU is CHARACTER*1
*>          Specifies options for computing all or part of the matrix U:
*>          = 'A':  all M columns of U are returned in array U:
*>          = 'S':  the first min(m,n) columns of U (the left singular
*>                  vectors) are returned in the array U;
*>          = 'O':  the first min(m,n) columns of U (the left singular
*>                  vectors) are overwritten on the array A;
*>          = 'N':  no columns of U (no left singular vectors) are
*>                  computed.
*> \endverbatim
*>
*> \param[in] JOBVT
*> \verbatim
*>          JOBVT is CHARACTER*1
*>          Specifies options for computing all or part of the matrix
*>          V**T:
*>          = 'A':  all N rows of V**T are returned in the array VT;
*>          = 'S':  the first min(m,n) rows of V**T (the right singular
*>                  vectors) are returned in the array VT;
*>          = 'O':  the first min(m,n) rows of V**T (the right singular
*>                  vectors) are overwritten on the array A;
*>          = 'N':  no rows of V**T (no right singular vectors) are
*>                  computed.
*>
*>          JOBVT and JOBU cannot both be 'O'.
*> \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.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA,N)
*>          On entry, the M-by-N matrix A.
*>          On exit,
*>          if JOBU = 'O',  A is overwritten with the first min(m,n)
*>                          columns of U (the left singular vectors,
*>                          stored columnwise);
*>          if JOBVT = 'O', A is overwritten with the first min(m,n)
*>                          rows of V**T (the right singular vectors,
*>                          stored rowwise);
*>          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
*>                          are destroyed.
*> \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, dimension (min(M,N))
*>          The singular values of A, sorted so that S(i) >= S(i+1).
*> \endverbatim
*>
*> \param[out] U
*> \verbatim
*>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)
*>          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
*>          If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
*>          if JOBU = 'S', U contains the first min(m,n) columns of U
*>          (the left singular vectors, stored columnwise);
*>          if JOBU = 'N' or 'O', U is not referenced.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of the array U.  LDU >= 1; if
*>          JOBU = 'S' or 'A', LDU >= M.
*> \endverbatim
*>
*> \param[out] VT
*> \verbatim
*>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
*>          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
*>          V**T;
*>          if JOBVT = 'S', VT contains the first min(m,n) rows of
*>          V**T (the right singular vectors, stored rowwise);
*>          if JOBVT = 'N' or 'O', VT is not referenced.
*> \endverbatim
*>
*> \param[in] LDVT
*> \verbatim
*>          LDVT is INTEGER
*>          The leading dimension of the array VT.  LDVT >= 1; if
*>          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
*>          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
*>          superdiagonal elements of an upper bidiagonal matrix B
*>          whose diagonal is in S (not necessarily sorted). B
*>          satisfies A = U * B * VT, so it has the same singular values
*>          as A, and singular vectors related by U and VT.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK.
*>          LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code):
*>             - PATH 1  (M much larger than N, JOBU='N') 
*>             - PATH 1t (N much larger than M, JOBVT='N')
*>          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) for the other paths
*>          For good performance, LWORK should generally be larger.
*>
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates the optimal size of the WORK array, returns
*>          this value as the first entry of the WORK array, 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
*>                did not converge to zero. See the description of WORK
*>                above for details.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date April 2012
*
*> \ingroup doubleGEsing
*
*  =====================================================================
      SUBROUTINE dgesvd( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
     $                   vt, ldvt, work, lwork, info )
*
*  -- LAPACK driver routine (version 3.4.1) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     April 2012
*
*     .. Scalar Arguments ..
      CHARACTER          jobu, jobvt
      INTEGER            info, lda, ldu, ldvt, lwork, m, n
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   a( lda, * ), s( * ), u( ldu, * ),
     $                   vt( ldvt, * ), work( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   zero, one
      parameter( zero = 0.0d0, one = 1.0d0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            lquery, wntua, wntuas, wntun, wntuo, wntus,
     $                   wntva, wntvas, wntvn, wntvo, wntvs
      INTEGER            bdspac, blk, chunk, i, ie, ierr, ir, iscl,
     $                   itau, itaup, itauq, iu, iwork, ldwrkr, ldwrku,
     $                   maxwrk, minmn, minwrk, mnthr, ncu, ncvt, nru,
     $                   nrvt, wrkbl
      INTEGER            lwork_dgeqrf, lwork_dorgqr_n, lwork_dorgqr_m,
     $                   lwork_dgebrd, lwork_dorgbr_p, lwork_dorgbr_q,
     $                   lwork_dgelqf, lwork_dorglq_n, lwork_dorglq_m
      DOUBLE PRECISION   anrm, bignum, eps, smlnum
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   dum( 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           dbdsqr, dgebrd, dgelqf, dgemm, dgeqrf, dlacpy,
     $                   dlascl, dlaset, dorgbr, dorglq, dorgqr, dormbr,
     $                   xerbla
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      INTEGER            ilaenv
      DOUBLE PRECISION   dlamch, dlange
      EXTERNAL           lsame, ilaenv, dlamch, dlange
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      info = 0
      minmn = min( m, n )
      wntua = lsame( jobu, 'A' )
      wntus = lsame( jobu, 'S' )
      wntuas = wntua .OR. wntus
      wntuo = lsame( jobu, 'O' )
      wntun = lsame( jobu, 'N' )
      wntva = lsame( jobvt, 'A' )
      wntvs = lsame( jobvt, 'S' )
      wntvas = wntva .OR. wntvs
      wntvo = lsame( jobvt, 'O' )
      wntvn = lsame( jobvt, 'N' )
      lquery = ( lwork.EQ.-1 )
*
      IF( .NOT.( wntua .OR. wntus .OR. wntuo .OR. wntun ) ) THEN
         info = -1
      ELSE IF( .NOT.( wntva .OR. wntvs .OR. wntvo .OR. wntvn ) .OR.
     $         ( wntvo .AND. wntuo ) ) THEN
         info = -2
      ELSE IF( m.LT.0 ) THEN
         info = -3
      ELSE IF( n.LT.0 ) THEN
         info = -4
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -6
      ELSE IF( ldu.LT.1 .OR. ( wntuas .AND. ldu.LT.m ) ) THEN
         info = -9
      ELSE IF( ldvt.LT.1 .OR. ( wntva .AND. ldvt.LT.n ) .OR.
     $         ( wntvs .AND. ldvt.LT.minmn ) ) THEN
         info = -11
      END IF
*
*     Compute workspace
*      (Note: Comments in the code beginning "Workspace:" describe the
*       minimal amount of workspace needed at that point in the code,
*       as well as the preferred amount for good performance.
*       NB refers to the optimal block size for the immediately
*       following subroutine, as returned by ILAENV.)
*
      IF( info.EQ.0 ) THEN
         minwrk = 1
         maxwrk = 1
         IF( m.GE.n .AND. minmn.GT.0 ) THEN
*
*           Compute space needed for DBDSQR
*
            mnthr = ilaenv( 6, 'DGESVD', jobu // jobvt, m, n, 0, 0 )
            bdspac = 5*n
*           Compute space needed for DGEQRF
            CALL dgeqrf( m, n, a, lda, dum(1), dum(1), -1, ierr )
            lwork_dgeqrf=dum(1)
*           Compute space needed for DORGQR
            CALL dorgqr( m, n, n, a, lda, dum(1), dum(1), -1, ierr )
            lwork_dorgqr_n=dum(1)
            CALL dorgqr( m, m, n, a, lda, dum(1), dum(1), -1, ierr )
            lwork_dorgqr_m=dum(1)
*           Compute space needed for DGEBRD
            CALL dgebrd( n, n, a, lda, s, dum(1), dum(1),
     $                   dum(1), dum(1), -1, ierr )
            lwork_dgebrd=dum(1)
*           Compute space needed for DORGBR P
            CALL dorgbr( 'P', n, n, n, a, lda, dum(1),
     $                   dum(1), -1, ierr )
            lwork_dorgbr_p=dum(1)
*           Compute space needed for DORGBR Q
            CALL dorgbr( 'Q', n, n, n, a, lda, dum(1),
     $                   dum(1), -1, ierr )
            lwork_dorgbr_q=dum(1)
*
            IF( m.GE.mnthr ) THEN
               IF( wntun ) THEN
*
*                 Path 1 (M much larger than N, JOBU='N')
*
                  maxwrk = n + lwork_dgeqrf
                  maxwrk = max( maxwrk, 3*n+lwork_dgebrd )
                  IF( wntvo .OR. wntvas )
     $               maxwrk = max( maxwrk, 3*n+lwork_dorgbr_p )
                  maxwrk = max( maxwrk, bdspac )
                  minwrk = max( 4*n, bdspac )
               ELSE IF( wntuo .AND. wntvn ) THEN
*
*                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_n )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = max( n*n+wrkbl, n*n+m*n+n )
                  minwrk = max( 3*n+m, bdspac )
               ELSE IF( wntuo .AND. wntvas ) THEN
*
*                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
*                 'A')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_n )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = max( n*n+wrkbl, n*n+m*n+n )
                  minwrk = max( 3*n+m, bdspac )
               ELSE IF( wntus .AND. wntvn ) THEN
*
*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_n )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = n*n + wrkbl
                  minwrk = max( 3*n+m, bdspac )
               ELSE IF( wntus .AND. wntvo ) THEN
*
*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_n )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = 2*n*n + wrkbl
                  minwrk = max( 3*n+m, bdspac )
               ELSE IF( wntus .AND. wntvas ) THEN
*
*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
*                 'A')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_n )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = n*n + wrkbl
                  minwrk = max( 3*n+m, bdspac )
               ELSE IF( wntua .AND. wntvn ) THEN
*
*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_m )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = n*n + wrkbl
                  minwrk = max( 3*n+m, bdspac )
               ELSE IF( wntua .AND. wntvo ) THEN
*
*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_m )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = 2*n*n + wrkbl
                  minwrk = max( 3*n+m, bdspac )
               ELSE IF( wntua .AND. wntvas ) THEN
*
*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
*                 'A')
*
                  wrkbl = n + lwork_dgeqrf
                  wrkbl = max( wrkbl, n+lwork_dorgqr_m )
                  wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = n*n + wrkbl
                  minwrk = max( 3*n+m, bdspac )
               END IF
            ELSE
*
*              Path 10 (M at least N, but not much larger)
*
               CALL dgebrd( m, n, a, lda, s, dum(1), dum(1),
     $                   dum(1), dum(1), -1, ierr )
               lwork_dgebrd=dum(1)
               maxwrk = 3*n + lwork_dgebrd
               IF( wntus .OR. wntuo ) THEN
                  CALL dorgbr( 'Q', m, n, n, a, lda, dum(1),
     $                   dum(1), -1, ierr )
                  lwork_dorgbr_q=dum(1)
                  maxwrk = max( maxwrk, 3*n+lwork_dorgbr_q )
               END IF
               IF( wntua ) THEN
                  CALL dorgbr( 'Q', m, m, n, a, lda, dum(1),
     $                   dum(1), -1, ierr )
                  lwork_dorgbr_q=dum(1)
                  maxwrk = max( maxwrk, 3*n+lwork_dorgbr_q )
               END IF
               IF( .NOT.wntvn ) THEN
                 maxwrk = max( maxwrk, 3*n+lwork_dorgbr_p )
               END IF
               maxwrk = max( maxwrk, bdspac )
               minwrk = max( 3*n+m, bdspac )
            END IF
         ELSE IF( minmn.GT.0 ) THEN
*
*           Compute space needed for DBDSQR
*
            mnthr = ilaenv( 6, 'DGESVD', jobu // jobvt, m, n, 0, 0 )
            bdspac = 5*m
*           Compute space needed for DGELQF
            CALL dgelqf( m, n, a, lda, dum(1), dum(1), -1, ierr )
            lwork_dgelqf=dum(1)
*           Compute space needed for DORGLQ
            CALL dorglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
            lwork_dorglq_n=dum(1)
            CALL dorglq( m, n, m, a, lda, dum(1), dum(1), -1, ierr )
            lwork_dorglq_m=dum(1)
*           Compute space needed for DGEBRD
            CALL dgebrd( m, m, a, lda, s, dum(1), dum(1),
     $                   dum(1), dum(1), -1, ierr )
            lwork_dgebrd=dum(1)
*            Compute space needed for DORGBR P
            CALL dorgbr( 'P', m, m, m, a, n, dum(1),
     $                   dum(1), -1, ierr )
            lwork_dorgbr_p=dum(1)
*           Compute space needed for DORGBR Q
            CALL dorgbr( 'Q', m, m, m, a, n, dum(1),
     $                   dum(1), -1, ierr )
            lwork_dorgbr_q=dum(1)
            IF( n.GE.mnthr ) THEN
               IF( wntvn ) THEN
*
*                 Path 1t(N much larger than M, JOBVT='N')
*
                  maxwrk = m + lwork_dgelqf
                  maxwrk = max( maxwrk, 3*m+lwork_dgebrd )
                  IF( wntuo .OR. wntuas )
     $               maxwrk = max( maxwrk, 3*m+lwork_dorgbr_q )
                  maxwrk = max( maxwrk, bdspac )
                  minwrk = max( 4*m, bdspac )
               ELSE IF( wntvo .AND. wntun ) THEN
*
*                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_m )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = max( m*m+wrkbl, m*m+m*n+m )
                  minwrk = max( 3*m+n, bdspac )
               ELSE IF( wntvo .AND. wntuas ) THEN
*
*                 Path 3t(N much larger than M, JOBU='S' or 'A',
*                 JOBVT='O')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_m )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = max( m*m+wrkbl, m*m+m*n+m )
                  minwrk = max( 3*m+n, bdspac )
               ELSE IF( wntvs .AND. wntun ) THEN
*
*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_m )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = m*m + wrkbl
                  minwrk = max( 3*m+n, bdspac )
               ELSE IF( wntvs .AND. wntuo ) THEN
*
*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_m )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = 2*m*m + wrkbl
                  minwrk = max( 3*m+n, bdspac )
               ELSE IF( wntvs .AND. wntuas ) THEN
*
*                 Path 6t(N much larger than M, JOBU='S' or 'A',
*                 JOBVT='S')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_m )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = m*m + wrkbl
                  minwrk = max( 3*m+n, bdspac )
               ELSE IF( wntva .AND. wntun ) THEN
*
*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_n )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = m*m + wrkbl
                  minwrk = max( 3*m+n, bdspac )
               ELSE IF( wntva .AND. wntuo ) THEN
*
*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_n )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = 2*m*m + wrkbl
                  minwrk = max( 3*m+n, bdspac )
               ELSE IF( wntva .AND. wntuas ) THEN
*
*                 Path 9t(N much larger than M, JOBU='S' or 'A',
*                 JOBVT='A')
*
                  wrkbl = m + lwork_dgelqf
                  wrkbl = max( wrkbl, m+lwork_dorglq_n )
                  wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
                  wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
                  wrkbl = max( wrkbl, bdspac )
                  maxwrk = m*m + wrkbl
                  minwrk = max( 3*m+n, bdspac )
               END IF
            ELSE
*
*              Path 10t(N greater than M, but not much larger)
*
               CALL dgebrd( m, n, a, lda, s, dum(1), dum(1),
     $                   dum(1), dum(1), -1, ierr )
               lwork_dgebrd=dum(1)
               maxwrk = 3*m + lwork_dgebrd
               IF( wntvs .OR. wntvo ) THEN
*                Compute space needed for DORGBR P
                 CALL dorgbr( 'P', m, n, m, a, n, dum(1),
     $                   dum(1), -1, ierr )
                 lwork_dorgbr_p=dum(1)
                 maxwrk = max( maxwrk, 3*m+lwork_dorgbr_p )
               END IF
               IF( wntva ) THEN
                 CALL dorgbr( 'P', n, n, m, a, n, dum(1),
     $                   dum(1), -1, ierr )
                 lwork_dorgbr_p=dum(1)
                 maxwrk = max( maxwrk, 3*m+lwork_dorgbr_p )
               END IF
               IF( .NOT.wntun ) THEN
                  maxwrk = max( maxwrk, 3*m+lwork_dorgbr_q )
               END IF
               maxwrk = max( maxwrk, bdspac )
               minwrk = max( 3*m+n, bdspac )
            END IF
         END IF
         maxwrk = max( maxwrk, minwrk )
         work( 1 ) = maxwrk
*
         IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
            info = -13
         END IF
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DGESVD', -info )
         return
      ELSE IF( lquery ) THEN
         return
      END IF
*
*     Quick return if possible
*
      IF( m.EQ.0 .OR. n.EQ.0 ) THEN
         return
      END IF
*
*     Get machine constants
*
      eps = dlamch( 'P' )
      smlnum = sqrt( dlamch( 'S' ) ) / eps
      bignum = one / smlnum
*
*     Scale A if max element outside range [SMLNUM,BIGNUM]
*
      anrm = dlange( 'M', m, n, a, lda, dum )
      iscl = 0
      IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
         iscl = 1
         CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
      ELSE IF( anrm.GT.bignum ) THEN
         iscl = 1
         CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
      END IF
*
      IF( m.GE.n ) THEN
*
*        A has at least as many rows as columns. If A has sufficiently
*        more rows than columns, first reduce using the QR
*        decomposition (if sufficient workspace available)
*
         IF( m.GE.mnthr ) THEN
*
            IF( wntun ) THEN
*
*              Path 1 (M much larger than N, JOBU='N')
*              No left singular vectors to be computed
*
               itau = 1
               iwork = itau + n
*
*              Compute A=Q*R
*              (Workspace: need 2*N, prefer N+N*NB)
*
               CALL dgeqrf( m, n, a, lda, work( itau ), work( iwork ),
     $                      lwork-iwork+1, ierr )
*
*              Zero out below R
*
               CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
               ie = 1
               itauq = ie + n
               itaup = itauq + n
               iwork = itaup + n
*
*              Bidiagonalize R in A
*              (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
               CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
     $                      work( itaup ), work( iwork ), lwork-iwork+1,
     $                      ierr )
               ncvt = 0
               IF( wntvo .OR. wntvas ) THEN
*
*                 If right singular vectors desired, generate P'.
*                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
                  CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  ncvt = n
               END IF
               iwork = ie + n
*
*              Perform bidiagonal QR iteration, computing right
*              singular vectors of A in A if desired
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'U', n, ncvt, 0, 0, s, work( ie ), a, lda,
     $                      dum, 1, dum, 1, work( iwork ), info )
*
*              If right singular vectors desired in VT, copy them there
*
               IF( wntvas )
     $            CALL dlacpy( 'F', n, n, a, lda, vt, ldvt )
*
            ELSE IF( wntuo .AND. wntvn ) THEN
*
*              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
*              N left singular vectors to be overwritten on A and
*              no right singular vectors to be computed
*
               IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
*
*                 Sufficient workspace for a fast algorithm
*
                  ir = 1
                  IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n ) THEN
*
*                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
*
                     ldwrku = lda
                     ldwrkr = lda
                  ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n ) THEN
*
*                    WORK(IU) is LDA by N, WORK(IR) is N by N
*
                     ldwrku = lda
                     ldwrkr = n
                  ELSE
*
*                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N
*
                     ldwrku = ( lwork-n*n-n ) / n
                     ldwrkr = n
                  END IF
                  itau = ir + ldwrkr*n
                  iwork = itau + n
*
*                 Compute A=Q*R
*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                  CALL dgeqrf( m, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Copy R to WORK(IR) and zero out below it
*
                  CALL dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
                  CALL dlaset( 'L', n-1, n-1, zero, zero, work( ir+1 ),
     $                         ldwrkr )
*
*                 Generate Q in A
*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                  CALL dorgqr( m, n, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  ie = itau
                  itauq = ie + n
                  itaup = itauq + n
                  iwork = itaup + n
*
*                 Bidiagonalize R in WORK(IR)
*                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
                  CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Generate left vectors bidiagonalizing R
*                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
*
                  CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
     $                         work( itauq ), work( iwork ),
     $                         lwork-iwork+1, ierr )
                  iwork = ie + n
*
*                 Perform bidiagonal QR iteration, computing left
*                 singular vectors of R in WORK(IR)
*                 (Workspace: need N*N+BDSPAC)
*
                  CALL dbdsqr( 'U', n, 0, n, 0, s, work( ie ), dum, 1,
     $                         work( ir ), ldwrkr, dum, 1,
     $                         work( iwork ), info )
                  iu = ie + n
*
*                 Multiply Q in A by left singular vectors of R in
*                 WORK(IR), storing result in WORK(IU) and copying to A
*                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
*
                  DO 10 i = 1, m, ldwrku
                     chunk = min( m-i+1, ldwrku )
                     CALL dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
     $                           lda, work( ir ), ldwrkr, zero,
     $                           work( iu ), ldwrku )
                     CALL dlacpy( 'F', chunk, n, work( iu ), ldwrku,
     $                            a( i, 1 ), lda )
   10             continue
*
               ELSE
*
*                 Insufficient workspace for a fast algorithm
*
                  ie = 1
                  itauq = ie + n
                  itaup = itauq + n
                  iwork = itaup + n
*
*                 Bidiagonalize A
*                 (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
*
                  CALL dgebrd( m, n, a, lda, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Generate left vectors bidiagonalizing A
*                 (Workspace: need 4*N, prefer 3*N+N*NB)
*
                  CALL dorgbr( 'Q', m, n, n, a, lda, work( itauq ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  iwork = ie + n
*
*                 Perform bidiagonal QR iteration, computing left
*                 singular vectors of A in A
*                 (Workspace: need BDSPAC)
*
                  CALL dbdsqr( 'U', n, 0, m, 0, s, work( ie ), dum, 1,
     $                         a, lda, dum, 1, work( iwork ), info )
*
               END IF
*
            ELSE IF( wntuo .AND. wntvas ) THEN
*
*              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
*              N left singular vectors to be overwritten on A and
*              N right singular vectors to be computed in VT
*
               IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
*
*                 Sufficient workspace for a fast algorithm
*
                  ir = 1
                  IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n ) THEN
*
*                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
*
                     ldwrku = lda
                     ldwrkr = lda
                  ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n ) THEN
*
*                    WORK(IU) is LDA by N and WORK(IR) is N by N
*
                     ldwrku = lda
                     ldwrkr = n
                  ELSE
*
*                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N
*
                     ldwrku = ( lwork-n*n-n ) / n
                     ldwrkr = n
                  END IF
                  itau = ir + ldwrkr*n
                  iwork = itau + n
*
*                 Compute A=Q*R
*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                  CALL dgeqrf( m, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Copy R to VT, zeroing out below it
*
                  CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
                  IF( n.GT.1 )
     $               CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            vt( 2, 1 ), ldvt )
*
*                 Generate Q in A
*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                  CALL dorgqr( m, n, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  ie = itau
                  itauq = ie + n
                  itaup = itauq + n
                  iwork = itaup + n
*
*                 Bidiagonalize R in VT, copying result to WORK(IR)
*                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
                  CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  CALL dlacpy( 'L', n, n, vt, ldvt, work( ir ), ldwrkr )
*
*                 Generate left vectors bidiagonalizing R in WORK(IR)
*                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
*
                  CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
     $                         work( itauq ), work( iwork ),
     $                         lwork-iwork+1, ierr )
*
*                 Generate right vectors bidiagonalizing R in VT
*                 (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)
*
                  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  iwork = ie + n
*
*                 Perform bidiagonal QR iteration, computing left
*                 singular vectors of R in WORK(IR) and computing right
*                 singular vectors of R in VT
*                 (Workspace: need N*N+BDSPAC)
*
                  CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ), vt, ldvt,
     $                         work( ir ), ldwrkr, dum, 1,
     $                         work( iwork ), info )
                  iu = ie + n
*
*                 Multiply Q in A by left singular vectors of R in
*                 WORK(IR), storing result in WORK(IU) and copying to A
*                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
*
                  DO 20 i = 1, m, ldwrku
                     chunk = min( m-i+1, ldwrku )
                     CALL dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
     $                           lda, work( ir ), ldwrkr, zero,
     $                           work( iu ), ldwrku )
                     CALL dlacpy( 'F', chunk, n, work( iu ), ldwrku,
     $                            a( i, 1 ), lda )
   20             continue
*
               ELSE
*
*                 Insufficient workspace for a fast algorithm
*
                  itau = 1
                  iwork = itau + n
*
*                 Compute A=Q*R
*                 (Workspace: need 2*N, prefer N+N*NB)
*
                  CALL dgeqrf( m, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Copy R to VT, zeroing out below it
*
                  CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
                  IF( n.GT.1 )
     $               CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            vt( 2, 1 ), ldvt )
*
*                 Generate Q in A
*                 (Workspace: need 2*N, prefer N+N*NB)
*
                  CALL dorgqr( m, n, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  ie = itau
                  itauq = ie + n
                  itaup = itauq + n
                  iwork = itaup + n
*
*                 Bidiagonalize R in VT
*                 (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
                  CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Multiply Q in A by left vectors bidiagonalizing R
*                 (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
                  CALL dormbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
     $                         work( itauq ), a, lda, work( iwork ),
     $                         lwork-iwork+1, ierr )
*
*                 Generate right vectors bidiagonalizing R in VT
*                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
                  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  iwork = ie + n
*
*                 Perform bidiagonal QR iteration, computing left
*                 singular vectors of A in A and computing right
*                 singular vectors of A in VT
*                 (Workspace: need BDSPAC)
*
                  CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), vt, ldvt,
     $                         a, lda, dum, 1, work( iwork ), info )
*
               END IF
*
            ELSE IF( wntus ) THEN
*
               IF( wntvn ) THEN
*
*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
*                 N left singular vectors to be computed in U and
*                 no right singular vectors to be computed
*
                  IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     ir = 1
                     IF( lwork.GE.wrkbl+lda*n ) THEN
*
*                       WORK(IR) is LDA by N
*
                        ldwrkr = lda
                     ELSE
*
*                       WORK(IR) is N by N
*
                        ldwrkr = n
                     END IF
                     itau = ir + ldwrkr*n
                     iwork = itau + n
*
*                    Compute A=Q*R
*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy R to WORK(IR), zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, work( ir ),
     $                            ldwrkr )
                     CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            work( ir+1 ), ldwrkr )
*
*                    Generate Q in A
*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                     CALL dorgqr( m, n, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in WORK(IR)
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
                     CALL dgebrd( n, n, work( ir ), ldwrkr, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate left vectors bidiagonalizing R in WORK(IR)
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
*
                     CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of R in WORK(IR)
*                    (Workspace: need N*N+BDSPAC)
*
                     CALL dbdsqr( 'U', n, 0, n, 0, s, work( ie ), dum,
     $                            1, work( ir ), ldwrkr, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply Q in A by left singular vectors of R in
*                    WORK(IR), storing result in U
*                    (Workspace: need N*N)
*
                     CALL dgemm( 'N', 'N', m, n, n, one, a, lda,
     $                           work( ir ), ldwrkr, zero, u, ldu )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dorgqr( m, n, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Zero out below R in A
*
                     CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
     $                            lda )
*
*                    Bidiagonalize R in A
*                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
                     CALL dgebrd( n, n, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply Q in U by left vectors bidiagonalizing R
*                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
                     CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
     $                            work( itauq ), u, ldu, work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', n, 0, m, 0, s, work( ie ), dum,
     $                            1, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntvo ) THEN
*
*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
*                 N left singular vectors to be computed in U and
*                 N right singular vectors to be overwritten on A
*
                  IF( lwork.GE.2*n*n+max( 4*n, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+2*lda*n ) THEN
*
*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
*
                        ldwrku = lda
                        ir = iu + ldwrku*n
                        ldwrkr = lda
                     ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
*
*                       WORK(IU) is LDA by N and WORK(IR) is N by N
*
                        ldwrku = lda
                        ir = iu + ldwrku*n
                        ldwrkr = n
                     ELSE
*
*                       WORK(IU) is N by N and WORK(IR) is N by N
*
                        ldwrku = n
                        ir = iu + ldwrku*n
                        ldwrkr = n
                     END IF
                     itau = ir + ldwrkr*n
                     iwork = itau + n
*
*                    Compute A=Q*R
*                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy R to WORK(IU), zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            work( iu+1 ), ldwrku )
*
*                    Generate Q in A
*                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
*
                     CALL dorgqr( m, n, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in WORK(IU), copying result to
*                    WORK(IR)
*                    (Workspace: need 2*N*N+4*N,
*                                prefer 2*N*N+3*N+2*N*NB)
*
                     CALL dgebrd( n, n, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', n, n, work( iu ), ldwrku,
     $                            work( ir ), ldwrkr )
*
*                    Generate left bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
*
                     CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in WORK(IR)
*                    (Workspace: need 2*N*N+4*N-1,
*                                prefer 2*N*N+3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, work( ir ), ldwrkr,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of R in WORK(IU) and computing
*                    right singular vectors of R in WORK(IR)
*                    (Workspace: need 2*N*N+BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ),
     $                            work( ir ), ldwrkr, work( iu ),
     $                            ldwrku, dum, 1, work( iwork ), info )
*
*                    Multiply Q in A by left singular vectors of R in
*                    WORK(IU), storing result in U
*                    (Workspace: need N*N)
*
                     CALL dgemm( 'N', 'N', m, n, n, one, a, lda,
     $                           work( iu ), ldwrku, zero, u, ldu )
*
*                    Copy right singular vectors of R to A
*                    (Workspace: need N*N)
*
                     CALL dlacpy( 'F', n, n, work( ir ), ldwrkr, a,
     $                            lda )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dorgqr( m, n, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Zero out below R in A
*
                     CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
     $                            lda )
*
*                    Bidiagonalize R in A
*                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
                     CALL dgebrd( n, n, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply Q in U by left vectors bidiagonalizing R
*                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
                     CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
     $                            work( itauq ), u, ldu, work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right vectors bidiagonalizing R in A
*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U and computing right
*                    singular vectors of A in A
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), a,
     $                            lda, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntvas ) THEN
*
*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
*                         or 'A')
*                 N left singular vectors to be computed in U and
*                 N right singular vectors to be computed in VT
*
                  IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+lda*n ) THEN
*
*                       WORK(IU) is LDA by N
*
                        ldwrku = lda
                     ELSE
*
*                       WORK(IU) is N by N
*
                        ldwrku = n
                     END IF
                     itau = iu + ldwrku*n
                     iwork = itau + n
*
*                    Compute A=Q*R
*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy R to WORK(IU), zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            work( iu+1 ), ldwrku )
*
*                    Generate Q in A
*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                     CALL dorgqr( m, n, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in WORK(IU), copying result to VT
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
                     CALL dgebrd( n, n, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', n, n, work( iu ), ldwrku, vt,
     $                            ldvt )
*
*                    Generate left bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
*
                     CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in VT
*                    (Workspace: need N*N+4*N-1,
*                                prefer N*N+3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of R in WORK(IU) and computing
*                    right singular vectors of R in VT
*                    (Workspace: need N*N+BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ), vt,
     $                            ldvt, work( iu ), ldwrku, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply Q in A by left singular vectors of R in
*                    WORK(IU), storing result in U
*                    (Workspace: need N*N)
*
                     CALL dgemm( 'N', 'N', m, n, n, one, a, lda,
     $                           work( iu ), ldwrku, zero, u, ldu )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dorgqr( m, n, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy R to VT, zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
                     IF( n.GT.1 )
     $                  CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                               vt( 2, 1 ), ldvt )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in VT
*                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
                     CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply Q in U by left bidiagonalizing vectors
*                    in VT
*                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
                     CALL dormbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
     $                            work( itauq ), u, ldu, work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in VT
*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U and computing right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), vt,
     $                            ldvt, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               END IF
*
            ELSE IF( wntua ) THEN
*
               IF( wntvn ) THEN
*
*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
*                 M left singular vectors to be computed in U and
*                 no right singular vectors to be computed
*
                  IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     ir = 1
                     IF( lwork.GE.wrkbl+lda*n ) THEN
*
*                       WORK(IR) is LDA by N
*
                        ldwrkr = lda
                     ELSE
*
*                       WORK(IR) is N by N
*
                        ldwrkr = n
                     END IF
                     itau = ir + ldwrkr*n
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Copy R to WORK(IR), zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, work( ir ),
     $                            ldwrkr )
                     CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            work( ir+1 ), ldwrkr )
*
*                    Generate Q in U
*                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
*
                     CALL dorgqr( m, m, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in WORK(IR)
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
                     CALL dgebrd( n, n, work( ir ), ldwrkr, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in WORK(IR)
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
*
                     CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of R in WORK(IR)
*                    (Workspace: need N*N+BDSPAC)
*
                     CALL dbdsqr( 'U', n, 0, n, 0, s, work( ie ), dum,
     $                            1, work( ir ), ldwrkr, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply Q in U by left singular vectors of R in
*                    WORK(IR), storing result in A
*                    (Workspace: need N*N)
*
                     CALL dgemm( 'N', 'N', m, n, n, one, u, ldu,
     $                           work( ir ), ldwrkr, zero, a, lda )
*
*                    Copy left singular vectors of A from A to U
*
                     CALL dlacpy( 'F', m, n, a, lda, u, ldu )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need N+M, prefer N+M*NB)
*
                     CALL dorgqr( m, m, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Zero out below R in A
*
                     CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
     $                            lda )
*
*                    Bidiagonalize R in A
*                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
                     CALL dgebrd( n, n, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply Q in U by left bidiagonalizing vectors
*                    in A
*                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
                     CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
     $                            work( itauq ), u, ldu, work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', n, 0, m, 0, s, work( ie ), dum,
     $                            1, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntvo ) THEN
*
*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
*                 M left singular vectors to be computed in U and
*                 N right singular vectors to be overwritten on A
*
                  IF( lwork.GE.2*n*n+max( n+m, 4*n, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+2*lda*n ) THEN
*
*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
*
                        ldwrku = lda
                        ir = iu + ldwrku*n
                        ldwrkr = lda
                     ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
*
*                       WORK(IU) is LDA by N and WORK(IR) is N by N
*
                        ldwrku = lda
                        ir = iu + ldwrku*n
                        ldwrkr = n
                     ELSE
*
*                       WORK(IU) is N by N and WORK(IR) is N by N
*
                        ldwrku = n
                        ir = iu + ldwrku*n
                        ldwrkr = n
                     END IF
                     itau = ir + ldwrkr*n
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
*
                     CALL dorgqr( m, m, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy R to WORK(IU), zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            work( iu+1 ), ldwrku )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in WORK(IU), copying result to
*                    WORK(IR)
*                    (Workspace: need 2*N*N+4*N,
*                                prefer 2*N*N+3*N+2*N*NB)
*
                     CALL dgebrd( n, n, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', n, n, work( iu ), ldwrku,
     $                            work( ir ), ldwrkr )
*
*                    Generate left bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
*
                     CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in WORK(IR)
*                    (Workspace: need 2*N*N+4*N-1,
*                                prefer 2*N*N+3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, work( ir ), ldwrkr,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of R in WORK(IU) and computing
*                    right singular vectors of R in WORK(IR)
*                    (Workspace: need 2*N*N+BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ),
     $                            work( ir ), ldwrkr, work( iu ),
     $                            ldwrku, dum, 1, work( iwork ), info )
*
*                    Multiply Q in U by left singular vectors of R in
*                    WORK(IU), storing result in A
*                    (Workspace: need N*N)
*
                     CALL dgemm( 'N', 'N', m, n, n, one, u, ldu,
     $                           work( iu ), ldwrku, zero, a, lda )
*
*                    Copy left singular vectors of A from A to U
*
                     CALL dlacpy( 'F', m, n, a, lda, u, ldu )
*
*                    Copy right singular vectors of R from WORK(IR) to A
*
                     CALL dlacpy( 'F', n, n, work( ir ), ldwrkr, a,
     $                            lda )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need N+M, prefer N+M*NB)
*
                     CALL dorgqr( m, m, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Zero out below R in A
*
                     CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
     $                            lda )
*
*                    Bidiagonalize R in A
*                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
                     CALL dgebrd( n, n, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply Q in U by left bidiagonalizing vectors
*                    in A
*                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
                     CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
     $                            work( itauq ), u, ldu, work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in A
*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U and computing right
*                    singular vectors of A in A
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), a,
     $                            lda, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntvas ) THEN
*
*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
*                         or 'A')
*                 M left singular vectors to be computed in U and
*                 N right singular vectors to be computed in VT
*
                  IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+lda*n ) THEN
*
*                       WORK(IU) is LDA by N
*
                        ldwrku = lda
                     ELSE
*
*                       WORK(IU) is N by N
*
                        ldwrku = n
                     END IF
                     itau = iu + ldwrku*n
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
*
                     CALL dorgqr( m, m, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy R to WORK(IU), zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                            work( iu+1 ), ldwrku )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in WORK(IU), copying result to VT
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
                     CALL dgebrd( n, n, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', n, n, work( iu ), ldwrku, vt,
     $                            ldvt )
*
*                    Generate left bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
*
                     CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in VT
*                    (Workspace: need N*N+4*N-1,
*                                prefer N*N+3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of R in WORK(IU) and computing
*                    right singular vectors of R in VT
*                    (Workspace: need N*N+BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ), vt,
     $                            ldvt, work( iu ), ldwrku, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply Q in U by left singular vectors of R in
*                    WORK(IU), storing result in A
*                    (Workspace: need N*N)
*
                     CALL dgemm( 'N', 'N', m, n, n, one, u, ldu,
     $                           work( iu ), ldwrku, zero, a, lda )
*
*                    Copy left singular vectors of A from A to U
*
                     CALL dlacpy( 'F', m, n, a, lda, u, ldu )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + n
*
*                    Compute A=Q*R, copying result to U
*                    (Workspace: need 2*N, prefer N+N*NB)
*
                     CALL dgeqrf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, n, a, lda, u, ldu )
*
*                    Generate Q in U
*                    (Workspace: need N+M, prefer N+M*NB)
*
                     CALL dorgqr( m, m, n, u, ldu, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy R from A to VT, zeroing out below it
*
                     CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
                     IF( n.GT.1 )
     $                  CALL dlaset( 'L', n-1, n-1, zero, zero,
     $                               vt( 2, 1 ), ldvt )
                     ie = itau
                     itauq = ie + n
                     itaup = itauq + n
                     iwork = itaup + n
*
*                    Bidiagonalize R in VT
*                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
                     CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply Q in U by left bidiagonalizing vectors
*                    in VT
*                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
                     CALL dormbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
     $                            work( itauq ), u, ldu, work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in VT
*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
                     CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + n
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U and computing right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), vt,
     $                            ldvt, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               END IF
*
            END IF
*
         ELSE
*
*           M .LT. MNTHR
*
*           Path 10 (M at least N, but not much larger)
*           Reduce to bidiagonal form without QR decomposition
*
            ie = 1
            itauq = ie + n
            itaup = itauq + n
            iwork = itaup + n
*
*           Bidiagonalize A
*           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
*
            CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
     $                   work( itaup ), work( iwork ), lwork-iwork+1,
     $                   ierr )
            IF( wntuas ) THEN
*
*              If left singular vectors desired in U, copy result to U
*              and generate left bidiagonalizing vectors in U
*              (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)
*
               CALL dlacpy( 'L', m, n, a, lda, u, ldu )
               IF( wntus )
     $            ncu = n
               IF( wntua )
     $            ncu = m
               CALL dorgbr( 'Q', m, ncu, n, u, ldu, work( itauq ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            IF( wntvas ) THEN
*
*              If right singular vectors desired in VT, copy result to
*              VT and generate right bidiagonalizing vectors in VT
*              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
               CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
               CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            IF( wntuo ) THEN
*
*              If left singular vectors desired in A, generate left
*              bidiagonalizing vectors in A
*              (Workspace: need 4*N, prefer 3*N+N*NB)
*
               CALL dorgbr( 'Q', m, n, n, a, lda, work( itauq ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            IF( wntvo ) THEN
*
*              If right singular vectors desired in A, generate right
*              bidiagonalizing vectors in A
*              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
               CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            iwork = ie + n
            IF( wntuas .OR. wntuo )
     $         nru = m
            IF( wntun )
     $         nru = 0
            IF( wntvas .OR. wntvo )
     $         ncvt = n
            IF( wntvn )
     $         ncvt = 0
            IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
*
*              Perform bidiagonal QR iteration, if desired, computing
*              left singular vectors in U and computing right singular
*              vectors in VT
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'U', n, ncvt, nru, 0, s, work( ie ), vt,
     $                      ldvt, u, ldu, dum, 1, work( iwork ), info )
            ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
*
*              Perform bidiagonal QR iteration, if desired, computing
*              left singular vectors in U and computing right singular
*              vectors in A
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'U', n, ncvt, nru, 0, s, work( ie ), a, lda,
     $                      u, ldu, dum, 1, work( iwork ), info )
            ELSE
*
*              Perform bidiagonal QR iteration, if desired, computing
*              left singular vectors in A and computing right singular
*              vectors in VT
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'U', n, ncvt, nru, 0, s, work( ie ), vt,
     $                      ldvt, a, lda, dum, 1, work( iwork ), info )
            END IF
*
         END IF
*
      ELSE
*
*        A has more columns than rows. If A has sufficiently more
*        columns than rows, first reduce using the LQ decomposition (if
*        sufficient workspace available)
*
         IF( n.GE.mnthr ) THEN
*
            IF( wntvn ) THEN
*
*              Path 1t(N much larger than M, JOBVT='N')
*              No right singular vectors to be computed
*
               itau = 1
               iwork = itau + m
*
*              Compute A=L*Q
*              (Workspace: need 2*M, prefer M+M*NB)
*
               CALL dgelqf( m, n, a, lda, work( itau ), work( iwork ),
     $                      lwork-iwork+1, ierr )
*
*              Zero out above L
*
               CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
               ie = 1
               itauq = ie + m
               itaup = itauq + m
               iwork = itaup + m
*
*              Bidiagonalize L in A
*              (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
               CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
     $                      work( itaup ), work( iwork ), lwork-iwork+1,
     $                      ierr )
               IF( wntuo .OR. wntuas ) THEN
*
*                 If left singular vectors desired, generate Q
*                 (Workspace: need 4*M, prefer 3*M+M*NB)
*
                  CALL dorgbr( 'Q', m, m, m, a, lda, work( itauq ),
     $                         work( iwork ), lwork-iwork+1, ierr )
               END IF
               iwork = ie + m
               nru = 0
               IF( wntuo .OR. wntuas )
     $            nru = m
*
*              Perform bidiagonal QR iteration, computing left singular
*              vectors of A in A if desired
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'U', m, 0, nru, 0, s, work( ie ), dum, 1, a,
     $                      lda, dum, 1, work( iwork ), info )
*
*              If left singular vectors desired in U, copy them there
*
               IF( wntuas )
     $            CALL dlacpy( 'F', m, m, a, lda, u, ldu )
*
            ELSE IF( wntvo .AND. wntun ) THEN
*
*              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
*              M right singular vectors to be overwritten on A and
*              no left singular vectors to be computed
*
               IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
*
*                 Sufficient workspace for a fast algorithm
*
                  ir = 1
                  IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m ) THEN
*
*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
*
                     ldwrku = lda
                     chunk = n
                     ldwrkr = lda
                  ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m ) THEN
*
*                    WORK(IU) is LDA by N and WORK(IR) is M by M
*
                     ldwrku = lda
                     chunk = n
                     ldwrkr = m
                  ELSE
*
*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
*
                     ldwrku = m
                     chunk = ( lwork-m*m-m ) / m
                     ldwrkr = m
                  END IF
                  itau = ir + ldwrkr*m
                  iwork = itau + m
*
*                 Compute A=L*Q
*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                  CALL dgelqf( m, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Copy L to WORK(IR) and zero out above it
*
                  CALL dlacpy( 'L', m, m, a, lda, work( ir ), ldwrkr )
                  CALL dlaset( 'U', m-1, m-1, zero, zero,
     $                         work( ir+ldwrkr ), ldwrkr )
*
*                 Generate Q in A
*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                  CALL dorglq( m, n, m, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  ie = itau
                  itauq = ie + m
                  itaup = itauq + m
                  iwork = itaup + m
*
*                 Bidiagonalize L in WORK(IR)
*                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
                  CALL dgebrd( m, m, work( ir ), ldwrkr, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Generate right vectors bidiagonalizing L
*                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
*
                  CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
     $                         work( itaup ), work( iwork ),
     $                         lwork-iwork+1, ierr )
                  iwork = ie + m
*
*                 Perform bidiagonal QR iteration, computing right
*                 singular vectors of L in WORK(IR)
*                 (Workspace: need M*M+BDSPAC)
*
                  CALL dbdsqr( 'U', m, m, 0, 0, s, work( ie ),
     $                         work( ir ), ldwrkr, dum, 1, dum, 1,
     $                         work( iwork ), info )
                  iu = ie + m
*
*                 Multiply right singular vectors of L in WORK(IR) by Q
*                 in A, storing result in WORK(IU) and copying to A
*                 (Workspace: need M*M+2*M, prefer M*M+M*N+M)
*
                  DO 30 i = 1, n, chunk
                     blk = min( n-i+1, chunk )
                     CALL dgemm( 'N', 'N', m, blk, m, one, work( ir ),
     $                           ldwrkr, a( 1, i ), lda, zero,
     $                           work( iu ), ldwrku )
                     CALL dlacpy( 'F', m, blk, work( iu ), ldwrku,
     $                            a( 1, i ), lda )
   30             continue
*
               ELSE
*
*                 Insufficient workspace for a fast algorithm
*
                  ie = 1
                  itauq = ie + m
                  itaup = itauq + m
                  iwork = itaup + m
*
*                 Bidiagonalize A
*                 (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
*
                  CALL dgebrd( m, n, a, lda, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Generate right vectors bidiagonalizing A
*                 (Workspace: need 4*M, prefer 3*M+M*NB)
*
                  CALL dorgbr( 'P', m, n, m, a, lda, work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  iwork = ie + m
*
*                 Perform bidiagonal QR iteration, computing right
*                 singular vectors of A in A
*                 (Workspace: need BDSPAC)
*
                  CALL dbdsqr( 'L', m, n, 0, 0, s, work( ie ), a, lda,
     $                         dum, 1, dum, 1, work( iwork ), info )
*
               END IF
*
            ELSE IF( wntvo .AND. wntuas ) THEN
*
*              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
*              M right singular vectors to be overwritten on A and
*              M left singular vectors to be computed in U
*
               IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
*
*                 Sufficient workspace for a fast algorithm
*
                  ir = 1
                  IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m ) THEN
*
*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
*
                     ldwrku = lda
                     chunk = n
                     ldwrkr = lda
                  ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m ) THEN
*
*                    WORK(IU) is LDA by N and WORK(IR) is M by M
*
                     ldwrku = lda
                     chunk = n
                     ldwrkr = m
                  ELSE
*
*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
*
                     ldwrku = m
                     chunk = ( lwork-m*m-m ) / m
                     ldwrkr = m
                  END IF
                  itau = ir + ldwrkr*m
                  iwork = itau + m
*
*                 Compute A=L*Q
*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                  CALL dgelqf( m, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Copy L to U, zeroing about above it
*
                  CALL dlacpy( 'L', m, m, a, lda, u, ldu )
                  CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
     $                         ldu )
*
*                 Generate Q in A
*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                  CALL dorglq( m, n, m, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  ie = itau
                  itauq = ie + m
                  itaup = itauq + m
                  iwork = itaup + m
*
*                 Bidiagonalize L in U, copying result to WORK(IR)
*                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
                  CALL dgebrd( m, m, u, ldu, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  CALL dlacpy( 'U', m, m, u, ldu, work( ir ), ldwrkr )
*
*                 Generate right vectors bidiagonalizing L in WORK(IR)
*                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
*
                  CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
     $                         work( itaup ), work( iwork ),
     $                         lwork-iwork+1, ierr )
*
*                 Generate left vectors bidiagonalizing L in U
*                 (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
*
                  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  iwork = ie + m
*
*                 Perform bidiagonal QR iteration, computing left
*                 singular vectors of L in U, and computing right
*                 singular vectors of L in WORK(IR)
*                 (Workspace: need M*M+BDSPAC)
*
                  CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
     $                         work( ir ), ldwrkr, u, ldu, dum, 1,
     $                         work( iwork ), info )
                  iu = ie + m
*
*                 Multiply right singular vectors of L in WORK(IR) by Q
*                 in A, storing result in WORK(IU) and copying to A
*                 (Workspace: need M*M+2*M, prefer M*M+M*N+M))
*
                  DO 40 i = 1, n, chunk
                     blk = min( n-i+1, chunk )
                     CALL dgemm( 'N', 'N', m, blk, m, one, work( ir ),
     $                           ldwrkr, a( 1, i ), lda, zero,
     $                           work( iu ), ldwrku )
                     CALL dlacpy( 'F', m, blk, work( iu ), ldwrku,
     $                            a( 1, i ), lda )
   40             continue
*
               ELSE
*
*                 Insufficient workspace for a fast algorithm
*
                  itau = 1
                  iwork = itau + m
*
*                 Compute A=L*Q
*                 (Workspace: need 2*M, prefer M+M*NB)
*
                  CALL dgelqf( m, n, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Copy L to U, zeroing out above it
*
                  CALL dlacpy( 'L', m, m, a, lda, u, ldu )
                  CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
     $                         ldu )
*
*                 Generate Q in A
*                 (Workspace: need 2*M, prefer M+M*NB)
*
                  CALL dorglq( m, n, m, a, lda, work( itau ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  ie = itau
                  itauq = ie + m
                  itaup = itauq + m
                  iwork = itaup + m
*
*                 Bidiagonalize L in U
*                 (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
                  CALL dgebrd( m, m, u, ldu, s, work( ie ),
     $                         work( itauq ), work( itaup ),
     $                         work( iwork ), lwork-iwork+1, ierr )
*
*                 Multiply right vectors bidiagonalizing L by Q in A
*                 (Workspace: need 3*M+N, prefer 3*M+N*NB)
*
                  CALL dormbr( 'P', 'L', 'T', m, n, m, u, ldu,
     $                         work( itaup ), a, lda, work( iwork ),
     $                         lwork-iwork+1, ierr )
*
*                 Generate left vectors bidiagonalizing L in U
*                 (Workspace: need 4*M, prefer 3*M+M*NB)
*
                  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
     $                         work( iwork ), lwork-iwork+1, ierr )
                  iwork = ie + m
*
*                 Perform bidiagonal QR iteration, computing left
*                 singular vectors of A in U and computing right
*                 singular vectors of A in A
*                 (Workspace: need BDSPAC)
*
                  CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), a, lda,
     $                         u, ldu, dum, 1, work( iwork ), info )
*
               END IF
*
            ELSE IF( wntvs ) THEN
*
               IF( wntun ) THEN
*
*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
*                 M right singular vectors to be computed in VT and
*                 no left singular vectors to be computed
*
                  IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     ir = 1
                     IF( lwork.GE.wrkbl+lda*m ) THEN
*
*                       WORK(IR) is LDA by M
*
                        ldwrkr = lda
                     ELSE
*
*                       WORK(IR) is M by M
*
                        ldwrkr = m
                     END IF
                     itau = ir + ldwrkr*m
                     iwork = itau + m
*
*                    Compute A=L*Q
*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy L to WORK(IR), zeroing out above it
*
                     CALL dlacpy( 'L', m, m, a, lda, work( ir ),
     $                            ldwrkr )
                     CALL dlaset( 'U', m-1, m-1, zero, zero,
     $                            work( ir+ldwrkr ), ldwrkr )
*
*                    Generate Q in A
*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                     CALL dorglq( m, n, m, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in WORK(IR)
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
                     CALL dgebrd( m, m, work( ir ), ldwrkr, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right vectors bidiagonalizing L in
*                    WORK(IR)
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
*
                     CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing right
*                    singular vectors of L in WORK(IR)
*                    (Workspace: need M*M+BDSPAC)
*
                     CALL dbdsqr( 'U', m, m, 0, 0, s, work( ie ),
     $                            work( ir ), ldwrkr, dum, 1, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply right singular vectors of L in WORK(IR) by
*                    Q in A, storing result in VT
*                    (Workspace: need M*M)
*
                     CALL dgemm( 'N', 'N', m, n, m, one, work( ir ),
     $                           ldwrkr, a, lda, zero, vt, ldvt )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + m
*
*                    Compute A=L*Q
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy result to VT
*
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dorglq( m, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Zero out above L in A
*
                     CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
     $                            lda )
*
*                    Bidiagonalize L in A
*                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
                     CALL dgebrd( m, m, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply right vectors bidiagonalizing L by Q in VT
*                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
*
                     CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
     $                            work( itaup ), vt, ldvt,
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', m, n, 0, 0, s, work( ie ), vt,
     $                            ldvt, dum, 1, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntuo ) THEN
*
*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
*                 M right singular vectors to be computed in VT and
*                 M left singular vectors to be overwritten on A
*
                  IF( lwork.GE.2*m*m+max( 4*m, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+2*lda*m ) THEN
*
*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
*
                        ldwrku = lda
                        ir = iu + ldwrku*m
                        ldwrkr = lda
                     ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
*
*                       WORK(IU) is LDA by M and WORK(IR) is M by M
*
                        ldwrku = lda
                        ir = iu + ldwrku*m
                        ldwrkr = m
                     ELSE
*
*                       WORK(IU) is M by M and WORK(IR) is M by M
*
                        ldwrku = m
                        ir = iu + ldwrku*m
                        ldwrkr = m
                     END IF
                     itau = ir + ldwrkr*m
                     iwork = itau + m
*
*                    Compute A=L*Q
*                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy L to WORK(IU), zeroing out below it
*
                     CALL dlacpy( 'L', m, m, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'U', m-1, m-1, zero, zero,
     $                            work( iu+ldwrku ), ldwrku )
*
*                    Generate Q in A
*                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
*
                     CALL dorglq( m, n, m, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in WORK(IU), copying result to
*                    WORK(IR)
*                    (Workspace: need 2*M*M+4*M,
*                                prefer 2*M*M+3*M+2*M*NB)
*
                     CALL dgebrd( m, m, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, m, work( iu ), ldwrku,
     $                            work( ir ), ldwrkr )
*
*                    Generate right bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need 2*M*M+4*M-1,
*                                prefer 2*M*M+3*M+(M-1)*NB)
*
                     CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in WORK(IR)
*                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, work( ir ), ldwrkr,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of L in WORK(IR) and computing
*                    right singular vectors of L in WORK(IU)
*                    (Workspace: need 2*M*M+BDSPAC)
*
                     CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
     $                            work( iu ), ldwrku, work( ir ),
     $                            ldwrkr, dum, 1, work( iwork ), info )
*
*                    Multiply right singular vectors of L in WORK(IU) by
*                    Q in A, storing result in VT
*                    (Workspace: need M*M)
*
                     CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
     $                           ldwrku, a, lda, zero, vt, ldvt )
*
*                    Copy left singular vectors of L to A
*                    (Workspace: need M*M)
*
                     CALL dlacpy( 'F', m, m, work( ir ), ldwrkr, a,
     $                            lda )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dorglq( m, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Zero out above L in A
*
                     CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
     $                            lda )
*
*                    Bidiagonalize L in A
*                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
                     CALL dgebrd( m, m, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply right vectors bidiagonalizing L by Q in VT
*                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
*
                     CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
     $                            work( itaup ), vt, ldvt,
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors of L in A
*                    (Workspace: need 4*M, prefer 3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, a, lda, work( itauq ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, compute left
*                    singular vectors of A in A and compute right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
     $                            ldvt, a, lda, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntuas ) THEN
*
*                 Path 6t(N much larger than M, JOBU='S' or 'A',
*                         JOBVT='S')
*                 M right singular vectors to be computed in VT and
*                 M left singular vectors to be computed in U
*
                  IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+lda*m ) THEN
*
*                       WORK(IU) is LDA by N
*
                        ldwrku = lda
                     ELSE
*
*                       WORK(IU) is LDA by M
*
                        ldwrku = m
                     END IF
                     itau = iu + ldwrku*m
                     iwork = itau + m
*
*                    Compute A=L*Q
*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy L to WORK(IU), zeroing out above it
*
                     CALL dlacpy( 'L', m, m, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'U', m-1, m-1, zero, zero,
     $                            work( iu+ldwrku ), ldwrku )
*
*                    Generate Q in A
*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                     CALL dorglq( m, n, m, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in WORK(IU), copying result to U
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
                     CALL dgebrd( m, m, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, m, work( iu ), ldwrku, u,
     $                            ldu )
*
*                    Generate right bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need M*M+4*M-1,
*                                prefer M*M+3*M+(M-1)*NB)
*
                     CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in U
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of L in U and computing right
*                    singular vectors of L in WORK(IU)
*                    (Workspace: need M*M+BDSPAC)
*
                     CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
     $                            work( iu ), ldwrku, u, ldu, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply right singular vectors of L in WORK(IU) by
*                    Q in A, storing result in VT
*                    (Workspace: need M*M)
*
                     CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
     $                           ldwrku, a, lda, zero, vt, ldvt )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dorglq( m, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy L to U, zeroing out above it
*
                     CALL dlacpy( 'L', m, m, a, lda, u, ldu )
                     CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
     $                            ldu )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in U
*                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
                     CALL dgebrd( m, m, u, ldu, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply right bidiagonalizing vectors in U by Q
*                    in VT
*                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
*
                     CALL dormbr( 'P', 'L', 'T', m, n, m, u, ldu,
     $                            work( itaup ), vt, ldvt,
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in U
*                    (Workspace: need 4*M, prefer 3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U and computing right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
     $                            ldvt, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               END IF
*
            ELSE IF( wntva ) THEN
*
               IF( wntun ) THEN
*
*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
*                 N right singular vectors to be computed in VT and
*                 no left singular vectors to be computed
*
                  IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     ir = 1
                     IF( lwork.GE.wrkbl+lda*m ) THEN
*
*                       WORK(IR) is LDA by M
*
                        ldwrkr = lda
                     ELSE
*
*                       WORK(IR) is M by M
*
                        ldwrkr = m
                     END IF
                     itau = ir + ldwrkr*m
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Copy L to WORK(IR), zeroing out above it
*
                     CALL dlacpy( 'L', m, m, a, lda, work( ir ),
     $                            ldwrkr )
                     CALL dlaset( 'U', m-1, m-1, zero, zero,
     $                            work( ir+ldwrkr ), ldwrkr )
*
*                    Generate Q in VT
*                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
*
                     CALL dorglq( n, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in WORK(IR)
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
                     CALL dgebrd( m, m, work( ir ), ldwrkr, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate right bidiagonalizing vectors in WORK(IR)
*                    (Workspace: need M*M+4*M-1,
*                                prefer M*M+3*M+(M-1)*NB)
*
                     CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing right
*                    singular vectors of L in WORK(IR)
*                    (Workspace: need M*M+BDSPAC)
*
                     CALL dbdsqr( 'U', m, m, 0, 0, s, work( ie ),
     $                            work( ir ), ldwrkr, dum, 1, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply right singular vectors of L in WORK(IR) by
*                    Q in VT, storing result in A
*                    (Workspace: need M*M)
*
                     CALL dgemm( 'N', 'N', m, n, m, one, work( ir ),
     $                           ldwrkr, vt, ldvt, zero, a, lda )
*
*                    Copy right singular vectors of A from A to VT
*
                     CALL dlacpy( 'F', m, n, a, lda, vt, ldvt )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need M+N, prefer M+N*NB)
*
                     CALL dorglq( n, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Zero out above L in A
*
                     CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
     $                            lda )
*
*                    Bidiagonalize L in A
*                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
                     CALL dgebrd( m, m, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply right bidiagonalizing vectors in A by Q
*                    in VT
*                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
*
                     CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
     $                            work( itaup ), vt, ldvt,
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', m, n, 0, 0, s, work( ie ), vt,
     $                            ldvt, dum, 1, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntuo ) THEN
*
*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
*                 N right singular vectors to be computed in VT and
*                 M left singular vectors to be overwritten on A
*
                  IF( lwork.GE.2*m*m+max( n+m, 4*m, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+2*lda*m ) THEN
*
*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
*
                        ldwrku = lda
                        ir = iu + ldwrku*m
                        ldwrkr = lda
                     ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
*
*                       WORK(IU) is LDA by M and WORK(IR) is M by M
*
                        ldwrku = lda
                        ir = iu + ldwrku*m
                        ldwrkr = m
                     ELSE
*
*                       WORK(IU) is M by M and WORK(IR) is M by M
*
                        ldwrku = m
                        ir = iu + ldwrku*m
                        ldwrkr = m
                     END IF
                     itau = ir + ldwrkr*m
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
*
                     CALL dorglq( n, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy L to WORK(IU), zeroing out above it
*
                     CALL dlacpy( 'L', m, m, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'U', m-1, m-1, zero, zero,
     $                            work( iu+ldwrku ), ldwrku )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in WORK(IU), copying result to
*                    WORK(IR)
*                    (Workspace: need 2*M*M+4*M,
*                                prefer 2*M*M+3*M+2*M*NB)
*
                     CALL dgebrd( m, m, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, m, work( iu ), ldwrku,
     $                            work( ir ), ldwrkr )
*
*                    Generate right bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need 2*M*M+4*M-1,
*                                prefer 2*M*M+3*M+(M-1)*NB)
*
                     CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in WORK(IR)
*                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, work( ir ), ldwrkr,
     $                            work( itauq ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of L in WORK(IR) and computing
*                    right singular vectors of L in WORK(IU)
*                    (Workspace: need 2*M*M+BDSPAC)
*
                     CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
     $                            work( iu ), ldwrku, work( ir ),
     $                            ldwrkr, dum, 1, work( iwork ), info )
*
*                    Multiply right singular vectors of L in WORK(IU) by
*                    Q in VT, storing result in A
*                    (Workspace: need M*M)
*
                     CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
     $                           ldwrku, vt, ldvt, zero, a, lda )
*
*                    Copy right singular vectors of A from A to VT
*
                     CALL dlacpy( 'F', m, n, a, lda, vt, ldvt )
*
*                    Copy left singular vectors of A from WORK(IR) to A
*
                     CALL dlacpy( 'F', m, m, work( ir ), ldwrkr, a,
     $                            lda )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need M+N, prefer M+N*NB)
*
                     CALL dorglq( n, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Zero out above L in A
*
                     CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
     $                            lda )
*
*                    Bidiagonalize L in A
*                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
                     CALL dgebrd( m, m, a, lda, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply right bidiagonalizing vectors in A by Q
*                    in VT
*                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
*
                     CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
     $                            work( itaup ), vt, ldvt,
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in A
*                    (Workspace: need 4*M, prefer 3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, a, lda, work( itauq ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in A and computing right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
     $                            ldvt, a, lda, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               ELSE IF( wntuas ) THEN
*
*                 Path 9t(N much larger than M, JOBU='S' or 'A',
*                         JOBVT='A')
*                 N right singular vectors to be computed in VT and
*                 M left singular vectors to be computed in U
*
                  IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) ) THEN
*
*                    Sufficient workspace for a fast algorithm
*
                     iu = 1
                     IF( lwork.GE.wrkbl+lda*m ) THEN
*
*                       WORK(IU) is LDA by M
*
                        ldwrku = lda
                     ELSE
*
*                       WORK(IU) is M by M
*
                        ldwrku = m
                     END IF
                     itau = iu + ldwrku*m
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
*
                     CALL dorglq( n, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy L to WORK(IU), zeroing out above it
*
                     CALL dlacpy( 'L', m, m, a, lda, work( iu ),
     $                            ldwrku )
                     CALL dlaset( 'U', m-1, m-1, zero, zero,
     $                            work( iu+ldwrku ), ldwrku )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in WORK(IU), copying result to U
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
                     CALL dgebrd( m, m, work( iu ), ldwrku, s,
     $                            work( ie ), work( itauq ),
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
                     CALL dlacpy( 'L', m, m, work( iu ), ldwrku, u,
     $                            ldu )
*
*                    Generate right bidiagonalizing vectors in WORK(IU)
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
*
                     CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
     $                            work( itaup ), work( iwork ),
     $                            lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in U
*                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of L in U and computing right
*                    singular vectors of L in WORK(IU)
*                    (Workspace: need M*M+BDSPAC)
*
                     CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
     $                            work( iu ), ldwrku, u, ldu, dum, 1,
     $                            work( iwork ), info )
*
*                    Multiply right singular vectors of L in WORK(IU) by
*                    Q in VT, storing result in A
*                    (Workspace: need M*M)
*
                     CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
     $                           ldwrku, vt, ldvt, zero, a, lda )
*
*                    Copy right singular vectors of A from A to VT
*
                     CALL dlacpy( 'F', m, n, a, lda, vt, ldvt )
*
                  ELSE
*
*                    Insufficient workspace for a fast algorithm
*
                     itau = 1
                     iwork = itau + m
*
*                    Compute A=L*Q, copying result to VT
*                    (Workspace: need 2*M, prefer M+M*NB)
*
                     CALL dgelqf( m, n, a, lda, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
*
*                    Generate Q in VT
*                    (Workspace: need M+N, prefer M+N*NB)
*
                     CALL dorglq( n, n, m, vt, ldvt, work( itau ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Copy L to U, zeroing out above it
*
                     CALL dlacpy( 'L', m, m, a, lda, u, ldu )
                     CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
     $                            ldu )
                     ie = itau
                     itauq = ie + m
                     itaup = itauq + m
                     iwork = itaup + m
*
*                    Bidiagonalize L in U
*                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
                     CALL dgebrd( m, m, u, ldu, s, work( ie ),
     $                            work( itauq ), work( itaup ),
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Multiply right bidiagonalizing vectors in U by Q
*                    in VT
*                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
*
                     CALL dormbr( 'P', 'L', 'T', m, n, m, u, ldu,
     $                            work( itaup ), vt, ldvt,
     $                            work( iwork ), lwork-iwork+1, ierr )
*
*                    Generate left bidiagonalizing vectors in U
*                    (Workspace: need 4*M, prefer 3*M+M*NB)
*
                     CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
     $                            work( iwork ), lwork-iwork+1, ierr )
                     iwork = ie + m
*
*                    Perform bidiagonal QR iteration, computing left
*                    singular vectors of A in U and computing right
*                    singular vectors of A in VT
*                    (Workspace: need BDSPAC)
*
                     CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
     $                            ldvt, u, ldu, dum, 1, work( iwork ),
     $                            info )
*
                  END IF
*
               END IF
*
            END IF
*
         ELSE
*
*           N .LT. MNTHR
*
*           Path 10t(N greater than M, but not much larger)
*           Reduce to bidiagonal form without LQ decomposition
*
            ie = 1
            itauq = ie + m
            itaup = itauq + m
            iwork = itaup + m
*
*           Bidiagonalize A
*           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
*
            CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
     $                   work( itaup ), work( iwork ), lwork-iwork+1,
     $                   ierr )
            IF( wntuas ) THEN
*
*              If left singular vectors desired in U, copy result to U
*              and generate left bidiagonalizing vectors in U
*              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
*
               CALL dlacpy( 'L', m, m, a, lda, u, ldu )
               CALL dorgbr( 'Q', m, m, n, u, ldu, work( itauq ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            IF( wntvas ) THEN
*
*              If right singular vectors desired in VT, copy result to
*              VT and generate right bidiagonalizing vectors in VT
*              (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)
*
               CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
               IF( wntva )
     $            nrvt = n
               IF( wntvs )
     $            nrvt = m
               CALL dorgbr( 'P', nrvt, n, m, vt, ldvt, work( itaup ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            IF( wntuo ) THEN
*
*              If left singular vectors desired in A, generate left
*              bidiagonalizing vectors in A
*              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
*
               CALL dorgbr( 'Q', m, m, n, a, lda, work( itauq ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            IF( wntvo ) THEN
*
*              If right singular vectors desired in A, generate right
*              bidiagonalizing vectors in A
*              (Workspace: need 4*M, prefer 3*M+M*NB)
*
               CALL dorgbr( 'P', m, n, m, a, lda, work( itaup ),
     $                      work( iwork ), lwork-iwork+1, ierr )
            END IF
            iwork = ie + m
            IF( wntuas .OR. wntuo )
     $         nru = m
            IF( wntun )
     $         nru = 0
            IF( wntvas .OR. wntvo )
     $         ncvt = n
            IF( wntvn )
     $         ncvt = 0
            IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
*
*              Perform bidiagonal QR iteration, if desired, computing
*              left singular vectors in U and computing right singular
*              vectors in VT
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'L', m, ncvt, nru, 0, s, work( ie ), vt,
     $                      ldvt, u, ldu, dum, 1, work( iwork ), info )
            ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
*
*              Perform bidiagonal QR iteration, if desired, computing
*              left singular vectors in U and computing right singular
*              vectors in A
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'L', m, ncvt, nru, 0, s, work( ie ), a, lda,
     $                      u, ldu, dum, 1, work( iwork ), info )
            ELSE
*
*              Perform bidiagonal QR iteration, if desired, computing
*              left singular vectors in A and computing right singular
*              vectors in VT
*              (Workspace: need BDSPAC)
*
               CALL dbdsqr( 'L', m, ncvt, nru, 0, s, work( ie ), vt,
     $                      ldvt, a, lda, dum, 1, work( iwork ), info )
            END IF
*
         END IF
*
      END IF
*
*     If DBDSQR failed to converge, copy unconverged superdiagonals
*     to WORK( 2:MINMN )
*
      IF( info.NE.0 ) THEN
         IF( ie.GT.2 ) THEN
            DO 50 i = 1, minmn - 1
               work( i+1 ) = work( i+ie-1 )
   50       continue
         END IF
         IF( ie.LT.2 ) THEN
            DO 60 i = minmn - 1, 1, -1
               work( i+1 ) = work( i+ie-1 )
   60       continue
         END IF
      END IF
*
*     Undo scaling if necessary
*
      IF( iscl.EQ.1 ) THEN
         IF( anrm.GT.bignum )
     $      CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
     $                   ierr )
         IF( info.NE.0 .AND. anrm.GT.bignum )
     $      CALL dlascl( 'G', 0, 0, bignum, anrm, minmn-1, 1, work( 2 ),
     $                   minmn, ierr )
         IF( anrm.LT.smlnum )
     $      CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
     $                   ierr )
         IF( info.NE.0 .AND. anrm.LT.smlnum )
     $      CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn-1, 1, work( 2 ),
     $                   minmn, ierr )
      END IF
*
*     Return optimal workspace in WORK(1)
*
      work( 1 ) = maxwrk
*
      return
*
*     End of DGESVD
*
      END