SRC/dgesdd.f

   1 *> \brief \b DGESDD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGESDD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
  22 *                          WORK, LWORK, IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ
  26 *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
  31 *      $                   VT( LDVT, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DGESDD computes the singular value decomposition (SVD) of a real
  41 *> M-by-N matrix A, optionally computing the left and right singular
  42 *> vectors.  If singular vectors are desired, it uses a
  43 *> divide-and-conquer algorithm.
  44 *>
  45 *> The SVD is written
  46 *>
  47 *>      A = U * SIGMA * transpose(V)
  48 *>
  49 *> where SIGMA is an M-by-N matrix which is zero except for its
  50 *> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
  51 *> V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
  52 *> are the singular values of A; they are real and non-negative, and
  53 *> are returned in descending order.  The first min(m,n) columns of
  54 *> U and V are the left and right singular vectors of A.
  55 *>
  56 *> Note that the routine returns VT = V**T, not V.
  57 *>
  58 *> The divide and conquer algorithm makes very mild assumptions about
  59 *> floating point arithmetic. It will work on machines with a guard
  60 *> digit in add/subtract, or on those binary machines without guard
  61 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
  62 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
  63 *> without guard digits, but we know of none.
  64 *> \endverbatim
  65 *
  66 *  Arguments:
  67 *  ==========
  68 *
  69 *> \param[in] JOBZ
  70 *> \verbatim
  71 *>          JOBZ is CHARACTER*1
  72 *>          Specifies options for computing all or part of the matrix U:
  73 *>          = 'A':  all M columns of U and all N rows of V**T are
  74 *>                  returned in the arrays U and VT;
  75 *>          = 'S':  the first min(M,N) columns of U and the first
  76 *>                  min(M,N) rows of V**T are returned in the arrays U
  77 *>                  and VT;
  78 *>          = 'O':  If M >= N, the first N columns of U are overwritten
  79 *>                  on the array A and all rows of V**T are returned in
  80 *>                  the array VT;
  81 *>                  otherwise, all columns of U are returned in the
  82 *>                  array U and the first M rows of V**T are overwritten
  83 *>                  in the array A;
  84 *>          = 'N':  no columns of U or rows of V**T are computed.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] M
  88 *> \verbatim
  89 *>          M is INTEGER
  90 *>          The number of rows of the input matrix A.  M >= 0.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] N
  94 *> \verbatim
  95 *>          N is INTEGER
  96 *>          The number of columns of the input matrix A.  N >= 0.
  97 *> \endverbatim
  98 *>
  99 *> \param[in,out] A
 100 *> \verbatim
 101 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 102 *>          On entry, the M-by-N matrix A.
 103 *>          On exit,
 104 *>          if JOBZ = 'O',  A is overwritten with the first N columns
 105 *>                          of U (the left singular vectors, stored
 106 *>                          columnwise) if M >= N;
 107 *>                          A is overwritten with the first M rows
 108 *>                          of V**T (the right singular vectors, stored
 109 *>                          rowwise) otherwise.
 110 *>          if JOBZ .ne. 'O', the contents of A are destroyed.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] LDA
 114 *> \verbatim
 115 *>          LDA is INTEGER
 116 *>          The leading dimension of the array A.  LDA >= max(1,M).
 117 *> \endverbatim
 118 *>
 119 *> \param[out] S
 120 *> \verbatim
 121 *>          S is DOUBLE PRECISION array, dimension (min(M,N))
 122 *>          The singular values of A, sorted so that S(i) >= S(i+1).
 123 *> \endverbatim
 124 *>
 125 *> \param[out] U
 126 *> \verbatim
 127 *>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)
 128 *>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
 129 *>          UCOL = min(M,N) if JOBZ = 'S'.
 130 *>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
 131 *>          orthogonal matrix U;
 132 *>          if JOBZ = 'S', U contains the first min(M,N) columns of U
 133 *>          (the left singular vectors, stored columnwise);
 134 *>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
 135 *> \endverbatim
 136 *>
 137 *> \param[in] LDU
 138 *> \verbatim
 139 *>          LDU is INTEGER
 140 *>          The leading dimension of the array U.  LDU >= 1; if
 141 *>          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
 142 *> \endverbatim
 143 *>
 144 *> \param[out] VT
 145 *> \verbatim
 146 *>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
 147 *>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
 148 *>          N-by-N orthogonal matrix V**T;
 149 *>          if JOBZ = 'S', VT contains the first min(M,N) rows of
 150 *>          V**T (the right singular vectors, stored rowwise);
 151 *>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
 152 *> \endverbatim
 153 *>
 154 *> \param[in] LDVT
 155 *> \verbatim
 156 *>          LDVT is INTEGER
 157 *>          The leading dimension of the array VT.  LDVT >= 1;
 158 *>          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
 159 *>          if JOBZ = 'S', LDVT >= min(M,N).
 160 *> \endverbatim
 161 *>
 162 *> \param[out] WORK
 163 *> \verbatim
 164 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 165 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
 166 *> \endverbatim
 167 *>
 168 *> \param[in] LWORK
 169 *> \verbatim
 170 *>          LWORK is INTEGER
 171 *>          The dimension of the array WORK. LWORK >= 1.
 172 *>          If LWORK = -1, a workspace query is assumed.  The optimal
 173 *>          size for the WORK array is calculated and stored in WORK(1),
 174 *>          and no other work except argument checking is performed.
 175 *>
 176 *>          Let mx = max(M,N) and mn = min(M,N).
 177 *>          If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ).
 178 *>          If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ).
 179 *>          If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn.
 180 *>          If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx.
 181 *>          These are not tight minimums in all cases; see comments inside code.
 182 *>          For good performance, LWORK should generally be larger;
 183 *>          a query is recommended.
 184 *> \endverbatim
 185 *>
 186 *> \param[out] IWORK
 187 *> \verbatim
 188 *>          IWORK is INTEGER array, dimension (8*min(M,N))
 189 *> \endverbatim
 190 *>
 191 *> \param[out] INFO
 192 *> \verbatim
 193 *>          INFO is INTEGER
 194 *>          = 0:  successful exit.
 195 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 196 *>          > 0:  DBDSDC did not converge, updating process failed.
 197 *> \endverbatim
 198 *
 199 *  Authors:
 200 *  ========
 201 *
 202 *> \author Univ. of Tennessee
 203 *> \author Univ. of California Berkeley
 204 *> \author Univ. of Colorado Denver
 205 *> \author NAG Ltd.
 206 *
 207 *> \date June 2016
 208 *
 209 *> \ingroup doubleGEsing
 210 *
 211 *> \par Contributors:
 212 *  ==================
 213 *>
 214 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 215 *>     California at Berkeley, USA
 216 *>
 217 *> @precisions fortran d -> s
 218 *  =====================================================================
 219       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 220      $                   WORK, LWORK, IWORK, INFO )
 221       implicit none
 222 *
 223 *  -- LAPACK driver routine (version 3.6.1) --
 224 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 225 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 226 *     June 2016
 227 *
 228 *     .. Scalar Arguments ..
 229       CHARACTER          JOBZ
 230       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
 231 *     ..
 232 *     .. Array Arguments ..
 233       INTEGER            IWORK( * )
 234       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
 235      $                   VT( LDVT, * ), WORK( * )
 236 *     ..
 237 *
 238 *  =====================================================================
 239 *
 240 *     .. Parameters ..
 241       DOUBLE PRECISION   ZERO, ONE
 242       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 243 *     ..
 244 *     .. Local Scalars ..
 245       LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
 246       INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
 247      $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
 248      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
 249      $                   MNTHR, NWORK, WRKBL
 250       INTEGER            LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
 251      $                   LWORK_DGEBRD_NN, LWORK_DGELQF_MN,
 252      $                   LWORK_DGEQRF_MN,
 253      $                   LWORK_DORGBR_P_MM, LWORK_DORGBR_Q_NN,
 254      $                   LWORK_DORGLQ_MN, LWORK_DORGLQ_NN,
 255      $                   LWORK_DORGQR_MM, LWORK_DORGQR_MN,
 256      $                   LWORK_DORMBR_PRT_MM, LWORK_DORMBR_QLN_MM,
 257      $                   LWORK_DORMBR_PRT_MN, LWORK_DORMBR_QLN_MN,
 258      $                   LWORK_DORMBR_PRT_NN, LWORK_DORMBR_QLN_NN
 259       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
 260 *     ..
 261 *     .. Local Arrays ..
 262       INTEGER            IDUM( 1 )
 263       DOUBLE PRECISION   DUM( 1 )
 264 *     ..
 265 *     .. External Subroutines ..
 266       EXTERNAL           DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
 267      $                   DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
 268      $                   XERBLA
 269 *     ..
 270 *     .. External Functions ..
 271       LOGICAL            LSAME
 272       DOUBLE PRECISION   DLAMCH, DLANGE
 273       EXTERNAL           DLAMCH, DLANGE, LSAME
 274 *     ..
 275 *     .. Intrinsic Functions ..
 276       INTRINSIC          INT, MAX, MIN, SQRT
 277 *     ..
 278 *     .. Executable Statements ..
 279 *
 280 *     Test the input arguments
 281 *
 282       INFO   = 0
 283       MINMN  = MIN( M, N )
 284       WNTQA  = LSAME( JOBZ, 'A' )
 285       WNTQS  = LSAME( JOBZ, 'S' )
 286       WNTQAS = WNTQA .OR. WNTQS
 287       WNTQO  = LSAME( JOBZ, 'O' )
 288       WNTQN  = LSAME( JOBZ, 'N' )
 289       LQUERY = ( LWORK.EQ.-1 )
 290 *
 291       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
 292          INFO = -1
 293       ELSE IF( M.LT.0 ) THEN
 294          INFO = -2
 295       ELSE IF( N.LT.0 ) THEN
 296          INFO = -3
 297       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 298          INFO = -5
 299       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
 300      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
 301          INFO = -8
 302       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
 303      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
 304      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
 305          INFO = -10
 306       END IF
 307 *
 308 *     Compute workspace
 309 *       Note: Comments in the code beginning "Workspace:" describe the
 310 *       minimal amount of workspace allocated at that point in the code,
 311 *       as well as the preferred amount for good performance.
 312 *       NB refers to the optimal block size for the immediately
 313 *       following subroutine, as returned by ILAENV.
 314 *
 315       IF( INFO.EQ.0 ) THEN
 316          MINWRK = 1
 317          MAXWRK = 1
 318          BDSPAC = 0
 319          MNTHR  = INT( MINMN*11.0D0 / 6.0D0 )
 320          IF( M.GE.N .AND. MINMN.GT.0 ) THEN
 321 *
 322 *           Compute space needed for DBDSDC
 323 *
 324             IF( WNTQN ) THEN
 325 *              dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
 326 *              keep 7*N for backwards compatability.
 327                BDSPAC = 7*N
 328             ELSE
 329                BDSPAC = 3*N*N + 4*N
 330             END IF
 331 *
 332 *           Compute space preferred for each routine
 333             CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
 334      $                   DUM(1), DUM(1), -1, IERR )
 335             LWORK_DGEBRD_MN = INT( DUM(1) )
 336 *
 337             CALL DGEBRD( N, N, DUM(1), N, DUM(1), DUM(1), DUM(1),
 338      $                   DUM(1), DUM(1), -1, IERR )
 339             LWORK_DGEBRD_NN = INT( DUM(1) )
 340 *
 341             CALL DGEQRF( M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
 342             LWORK_DGEQRF_MN = INT( DUM(1) )
 343 *
 344             CALL DORGBR( 'Q', N, N, N, DUM(1), N, DUM(1), DUM(1), -1,
 345      $                   IERR )
 346             LWORK_DORGBR_Q_NN = INT( DUM(1) )
 347 *
 348             CALL DORGQR( M, M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
 349             LWORK_DORGQR_MM = INT( DUM(1) )
 350 *
 351             CALL DORGQR( M, N, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
 352             LWORK_DORGQR_MN = INT( DUM(1) )
 353 *
 354             CALL DORMBR( 'P', 'R', 'T', N, N, N, DUM(1), N,
 355      $                   DUM(1), DUM(1), N, DUM(1), -1, IERR )
 356             LWORK_DORMBR_PRT_NN = INT( DUM(1) )
 357 *
 358             CALL DORMBR( 'Q', 'L', 'N', N, N, N, DUM(1), N,
 359      $                   DUM(1), DUM(1), N, DUM(1), -1, IERR )
 360             LWORK_DORMBR_QLN_NN = INT( DUM(1) )
 361 *
 362             CALL DORMBR( 'Q', 'L', 'N', M, N, N, DUM(1), M,
 363      $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
 364             LWORK_DORMBR_QLN_MN = INT( DUM(1) )
 365 *
 366             CALL DORMBR( 'Q', 'L', 'N', M, M, N, DUM(1), M,
 367      $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
 368             LWORK_DORMBR_QLN_MM = INT( DUM(1) )
 369 *
 370             IF( M.GE.MNTHR ) THEN
 371                IF( WNTQN ) THEN
 372 *
 373 *                 Path 1 (M >> N, JOBZ='N')
 374 *
 375                   WRKBL = N + LWORK_DGEQRF_MN
 376                   WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
 377                   MAXWRK = MAX( WRKBL, BDSPAC + N )
 378                   MINWRK = BDSPAC + N
 379                ELSE IF( WNTQO ) THEN
 380 *
 381 *                 Path 2 (M >> N, JOBZ='O')
 382 *
 383                   WRKBL = N + LWORK_DGEQRF_MN
 384                   WRKBL = MAX( WRKBL,   N + LWORK_DORGQR_MN )
 385                   WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
 386                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
 387                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
 388                   WRKBL = MAX( WRKBL, 3*N + BDSPAC )
 389                   MAXWRK = WRKBL + 2*N*N
 390                   MINWRK = BDSPAC + 2*N*N + 3*N
 391                ELSE IF( WNTQS ) THEN
 392 *
 393 *                 Path 3 (M >> N, JOBZ='S')
 394 *
 395                   WRKBL = N + LWORK_DGEQRF_MN
 396                   WRKBL = MAX( WRKBL,   N + LWORK_DORGQR_MN )
 397                   WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
 398                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
 399                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
 400                   WRKBL = MAX( WRKBL, 3*N + BDSPAC )
 401                   MAXWRK = WRKBL + N*N
 402                   MINWRK = BDSPAC + N*N + 3*N
 403                ELSE IF( WNTQA ) THEN
 404 *
 405 *                 Path 4 (M >> N, JOBZ='A')
 406 *
 407                   WRKBL = N + LWORK_DGEQRF_MN
 408                   WRKBL = MAX( WRKBL,   N + LWORK_DORGQR_MM )
 409                   WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
 410                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
 411                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
 412                   WRKBL = MAX( WRKBL, 3*N + BDSPAC )
 413                   MAXWRK = WRKBL + N*N
 414                   MINWRK = N*N + MAX( 3*N + BDSPAC, N + M )
 415                END IF
 416             ELSE
 417 *
 418 *              Path 5 (M >= N, but not much larger)
 419 *
 420                WRKBL = 3*N + LWORK_DGEBRD_MN
 421                IF( WNTQN ) THEN
 422 *                 Path 5n (M >= N, jobz='N')
 423                   MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
 424                   MINWRK = 3*N + MAX( M, BDSPAC )
 425                ELSE IF( WNTQO ) THEN
 426 *                 Path 5o (M >= N, jobz='O')
 427                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
 428                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
 429                   WRKBL = MAX( WRKBL, 3*N + BDSPAC )
 430                   MAXWRK = WRKBL + M*N
 431                   MINWRK = 3*N + MAX( M, N*N + BDSPAC )
 432                ELSE IF( WNTQS ) THEN
 433 *                 Path 5s (M >= N, jobz='S')
 434                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
 435                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
 436                   MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
 437                   MINWRK = 3*N + MAX( M, BDSPAC )
 438                ELSE IF( WNTQA ) THEN
 439 *                 Path 5a (M >= N, jobz='A')
 440                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MM )
 441                   WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
 442                   MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
 443                   MINWRK = 3*N + MAX( M, BDSPAC )
 444                END IF
 445             END IF
 446          ELSE IF( MINMN.GT.0 ) THEN
 447 *
 448 *           Compute space needed for DBDSDC
 449 *
 450             IF( WNTQN ) THEN
 451 *              dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
 452 *              keep 7*N for backwards compatability.
 453                BDSPAC = 7*M
 454             ELSE
 455                BDSPAC = 3*M*M + 4*M
 456             END IF
 457 *
 458 *           Compute space preferred for each routine
 459             CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
 460      $                   DUM(1), DUM(1), -1, IERR )
 461             LWORK_DGEBRD_MN = INT( DUM(1) )
 462 *
 463             CALL DGEBRD( M, M, A, M, S, DUM(1), DUM(1),
 464      $                   DUM(1), DUM(1), -1, IERR )
 465             LWORK_DGEBRD_MM = INT( DUM(1) )
 466 *
 467             CALL DGELQF( M, N, A, M, DUM(1), DUM(1), -1, IERR )
 468             LWORK_DGELQF_MN = INT( DUM(1) )
 469 *
 470             CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
 471             LWORK_DORGLQ_NN = INT( DUM(1) )
 472 *
 473             CALL DORGLQ( M, N, M, A, M, DUM(1), DUM(1), -1, IERR )
 474             LWORK_DORGLQ_MN = INT( DUM(1) )
 475 *
 476             CALL DORGBR( 'P', M, M, M, A, N, DUM(1), DUM(1), -1, IERR )
 477             LWORK_DORGBR_P_MM = INT( DUM(1) )
 478 *
 479             CALL DORMBR( 'P', 'R', 'T', M, M, M, DUM(1), M,
 480      $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
 481             LWORK_DORMBR_PRT_MM = INT( DUM(1) )
 482 *
 483             CALL DORMBR( 'P', 'R', 'T', M, N, M, DUM(1), M,
 484      $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
 485             LWORK_DORMBR_PRT_MN = INT( DUM(1) )
 486 *
 487             CALL DORMBR( 'P', 'R', 'T', N, N, M, DUM(1), N,
 488      $                   DUM(1), DUM(1), N, DUM(1), -1, IERR )
 489             LWORK_DORMBR_PRT_NN = INT( DUM(1) )
 490 *
 491             CALL DORMBR( 'Q', 'L', 'N', M, M, M, DUM(1), M,
 492      $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
 493             LWORK_DORMBR_QLN_MM = INT( DUM(1) )
 494 *
 495             IF( N.GE.MNTHR ) THEN
 496                IF( WNTQN ) THEN
 497 *
 498 *                 Path 1t (N >> M, JOBZ='N')
 499 *
 500                   WRKBL = M + LWORK_DGELQF_MN
 501                   WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
 502                   MAXWRK = MAX( WRKBL, BDSPAC + M )
 503                   MINWRK = BDSPAC + M
 504                ELSE IF( WNTQO ) THEN
 505 *
 506 *                 Path 2t (N >> M, JOBZ='O')
 507 *
 508                   WRKBL = M + LWORK_DGELQF_MN
 509                   WRKBL = MAX( WRKBL,   M + LWORK_DORGLQ_MN )
 510                   WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
 511                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
 512                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
 513                   WRKBL = MAX( WRKBL, 3*M + BDSPAC )
 514                   MAXWRK = WRKBL + 2*M*M
 515                   MINWRK = BDSPAC + 2*M*M + 3*M
 516                ELSE IF( WNTQS ) THEN
 517 *
 518 *                 Path 3t (N >> M, JOBZ='S')
 519 *
 520                   WRKBL = M + LWORK_DGELQF_MN
 521                   WRKBL = MAX( WRKBL,   M + LWORK_DORGLQ_MN )
 522                   WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
 523                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
 524                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
 525                   WRKBL = MAX( WRKBL, 3*M + BDSPAC )
 526                   MAXWRK = WRKBL + M*M
 527                   MINWRK = BDSPAC + M*M + 3*M
 528                ELSE IF( WNTQA ) THEN
 529 *
 530 *                 Path 4t (N >> M, JOBZ='A')
 531 *
 532                   WRKBL = M + LWORK_DGELQF_MN
 533                   WRKBL = MAX( WRKBL,   M + LWORK_DORGLQ_NN )
 534                   WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
 535                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
 536                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
 537                   WRKBL = MAX( WRKBL, 3*M + BDSPAC )
 538                   MAXWRK = WRKBL + M*M
 539                   MINWRK = M*M + MAX( 3*M + BDSPAC, M + N )
 540                END IF
 541             ELSE
 542 *
 543 *              Path 5t (N > M, but not much larger)
 544 *
 545                WRKBL = 3*M + LWORK_DGEBRD_MN
 546                IF( WNTQN ) THEN
 547 *                 Path 5tn (N > M, jobz='N')
 548                   MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
 549                   MINWRK = 3*M + MAX( N, BDSPAC )
 550                ELSE IF( WNTQO ) THEN
 551 *                 Path 5to (N > M, jobz='O')
 552                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
 553                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
 554                   WRKBL = MAX( WRKBL, 3*M + BDSPAC )
 555                   MAXWRK = WRKBL + M*N
 556                   MINWRK = 3*M + MAX( N, M*M + BDSPAC )
 557                ELSE IF( WNTQS ) THEN
 558 *                 Path 5ts (N > M, jobz='S')
 559                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
 560                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
 561                   MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
 562                   MINWRK = 3*M + MAX( N, BDSPAC )
 563                ELSE IF( WNTQA ) THEN
 564 *                 Path 5ta (N > M, jobz='A')
 565                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
 566                   WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_NN )
 567                   MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
 568                   MINWRK = 3*M + MAX( N, BDSPAC )
 569                END IF
 570             END IF
 571          END IF
 572
 573          MAXWRK = MAX( MAXWRK, MINWRK )
 574          WORK( 1 ) = MAXWRK
 575 *
 576          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 577             INFO = -12
 578          END IF
 579       END IF
 580 *
 581       IF( INFO.NE.0 ) THEN
 582          CALL XERBLA( 'DGESDD', -INFO )
 583          RETURN
 584       ELSE IF( LQUERY ) THEN
 585          RETURN
 586       END IF
 587 *
 588 *     Quick return if possible
 589 *
 590       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 591          RETURN
 592       END IF
 593 *
 594 *     Get machine constants
 595 *
 596       EPS = DLAMCH( 'P' )
 597       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
 598       BIGNUM = ONE / SMLNUM
 599 *
 600 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 601 *
 602       ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
 603       ISCL = 0
 604       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 605          ISCL = 1
 606          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
 607       ELSE IF( ANRM.GT.BIGNUM ) THEN
 608          ISCL = 1
 609          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
 610       END IF
 611 *
 612       IF( M.GE.N ) THEN
 613 *
 614 *        A has at least as many rows as columns. If A has sufficiently
 615 *        more rows than columns, first reduce using the QR
 616 *        decomposition (if sufficient workspace available)
 617 *
 618          IF( M.GE.MNTHR ) THEN
 619 *
 620             IF( WNTQN ) THEN
 621 *
 622 *              Path 1 (M >> N, JOBZ='N')
 623 *              No singular vectors to be computed
 624 *
 625                ITAU = 1
 626                NWORK = ITAU + N
 627 *
 628 *              Compute A=Q*R
 629 *              Workspace: need   N [tau] + N    [work]
 630 *              Workspace: prefer N [tau] + N*NB [work]
 631 *
 632                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 633      $                      LWORK - NWORK + 1, IERR )
 634 *
 635 *              Zero out below R
 636 *
 637                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
 638                IE = 1
 639                ITAUQ = IE + N
 640                ITAUP = ITAUQ + N
 641                NWORK = ITAUP + N
 642 *
 643 *              Bidiagonalize R in A
 644 *              Workspace: need   3*N [e, tauq, taup] + N      [work]
 645 *              Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work]
 646 *
 647                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 648      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 649      $                      IERR )
 650                NWORK = IE + N
 651 *
 652 *              Perform bidiagonal SVD, computing singular values only
 653 *              Workspace: need   N [e] + BDSPAC
 654 *
 655                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
 656      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
 657 *
 658             ELSE IF( WNTQO ) THEN
 659 *
 660 *              Path 2 (M >> N, JOBZ = 'O')
 661 *              N left singular vectors to be overwritten on A and
 662 *              N right singular vectors to be computed in VT
 663 *
 664                IR = 1
 665 *
 666 *              WORK(IR) is LDWRKR by N
 667 *
 668                IF( LWORK .GE. LDA*N + N*N + 3*N + BDSPAC ) THEN
 669                   LDWRKR = LDA
 670                ELSE
 671                   LDWRKR = ( LWORK - N*N - 3*N - BDSPAC ) / N
 672                END IF
 673                ITAU = IR + LDWRKR*N
 674                NWORK = ITAU + N
 675 *
 676 *              Compute A=Q*R
 677 *              Workspace: need   N*N [R] + N [tau] + N    [work]
 678 *              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
 679 *
 680                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 681      $                      LWORK - NWORK + 1, IERR )
 682 *
 683 *              Copy R to WORK(IR), zeroing out below it
 684 *
 685                CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
 686                CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
 687      $                      LDWRKR )
 688 *
 689 *              Generate Q in A
 690 *              Workspace: need   N*N [R] + N [tau] + N    [work]
 691 *              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
 692 *
 693                CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
 694      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
 695                IE = ITAU
 696                ITAUQ = IE + N
 697                ITAUP = ITAUQ + N
 698                NWORK = ITAUP + N
 699 *
 700 *              Bidiagonalize R in WORK(IR)
 701 *              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N      [work]
 702 *              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
 703 *
 704                CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
 705      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
 706      $                      LWORK - NWORK + 1, IERR )
 707 *
 708 *              WORK(IU) is N by N
 709 *
 710                IU = NWORK
 711                NWORK = IU + N*N
 712 *
 713 *              Perform bidiagonal SVD, computing left singular vectors
 714 *              of bidiagonal matrix in WORK(IU) and computing right
 715 *              singular vectors of bidiagonal matrix in VT
 716 *              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC
 717 *
 718                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
 719      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 720      $                      INFO )
 721 *
 722 *              Overwrite WORK(IU) by left singular vectors of R
 723 *              and VT by right singular vectors of R
 724 *              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N    [work]
 725 *              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
 726 *
 727                CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
 728      $                      WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
 729      $                      LWORK - NWORK + 1, IERR )
 730                CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
 731      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 732      $                      LWORK - NWORK + 1, IERR )
 733 *
 734 *              Multiply Q in A by left singular vectors of R in
 735 *              WORK(IU), storing result in WORK(IR) and copying to A
 736 *              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N*N [U]
 737 *              Workspace: prefer M*N [R] + 3*N [e, tauq, taup] + N*N [U]
 738 *
 739                DO 10 I = 1, M, LDWRKR
 740                   CHUNK = MIN( M - I + 1, LDWRKR )
 741                   CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
 742      $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
 743      $                        LDWRKR )
 744                   CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
 745      $                         A( I, 1 ), LDA )
 746    10          CONTINUE
 747 *
 748             ELSE IF( WNTQS ) THEN
 749 *
 750 *              Path 3 (M >> N, JOBZ='S')
 751 *              N left singular vectors to be computed in U and
 752 *              N right singular vectors to be computed in VT
 753 *
 754                IR = 1
 755 *
 756 *              WORK(IR) is N by N
 757 *
 758                LDWRKR = N
 759                ITAU = IR + LDWRKR*N
 760                NWORK = ITAU + N
 761 *
 762 *              Compute A=Q*R
 763 *              Workspace: need   N*N [R] + N [tau] + N    [work]
 764 *              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
 765 *
 766                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 767      $                      LWORK - NWORK + 1, IERR )
 768 *
 769 *              Copy R to WORK(IR), zeroing out below it
 770 *
 771                CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
 772                CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
 773      $                      LDWRKR )
 774 *
 775 *              Generate Q in A
 776 *              Workspace: need   N*N [R] + N [tau] + N    [work]
 777 *              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
 778 *
 779                CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
 780      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
 781                IE = ITAU
 782                ITAUQ = IE + N
 783                ITAUP = ITAUQ + N
 784                NWORK = ITAUP + N
 785 *
 786 *              Bidiagonalize R in WORK(IR)
 787 *              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N      [work]
 788 *              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
 789 *
 790                CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
 791      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
 792      $                      LWORK - NWORK + 1, IERR )
 793 *
 794 *              Perform bidiagonal SVD, computing left singular vectors
 795 *              of bidiagoal matrix in U and computing right singular
 796 *              vectors of bidiagonal matrix in VT
 797 *              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + BDSPAC
 798 *
 799                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
 800      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 801      $                      INFO )
 802 *
 803 *              Overwrite U by left singular vectors of R and VT
 804 *              by right singular vectors of R
 805 *              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N    [work]
 806 *              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*NB [work]
 807 *
 808                CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
 809      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
 810      $                      LWORK - NWORK + 1, IERR )
 811 *
 812                CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
 813      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 814      $                      LWORK - NWORK + 1, IERR )
 815 *
 816 *              Multiply Q in A by left singular vectors of R in
 817 *              WORK(IR), storing result in U
 818 *              Workspace: need   N*N [R]
 819 *
 820                CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
 821                CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
 822      $                     LDWRKR, ZERO, U, LDU )
 823 *
 824             ELSE IF( WNTQA ) THEN
 825 *
 826 *              Path 4 (M >> N, JOBZ='A')
 827 *              M left singular vectors to be computed in U and
 828 *              N right singular vectors to be computed in VT
 829 *
 830                IU = 1
 831 *
 832 *              WORK(IU) is N by N
 833 *
 834                LDWRKU = N
 835                ITAU = IU + LDWRKU*N
 836                NWORK = ITAU + N
 837 *
 838 *              Compute A=Q*R, copying result to U
 839 *              Workspace: need   N*N [U] + N [tau] + N    [work]
 840 *              Workspace: prefer N*N [U] + N [tau] + N*NB [work]
 841 *
 842                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 843      $                      LWORK - NWORK + 1, IERR )
 844                CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 845 *
 846 *              Generate Q in U
 847 *              Workspace: need   N*N [U] + N [tau] + M    [work]
 848 *              Workspace: prefer N*N [U] + N [tau] + M*NB [work]
 849                CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
 850      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
 851 *
 852 *              Produce R in A, zeroing out other entries
 853 *
 854                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
 855                IE = ITAU
 856                ITAUQ = IE + N
 857                ITAUP = ITAUQ + N
 858                NWORK = ITAUP + N
 859 *
 860 *              Bidiagonalize R in A
 861 *              Workspace: need   N*N [U] + 3*N [e, tauq, taup] + N      [work]
 862 *              Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + 2*N*NB [work]
 863 *
 864                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 865      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 866      $                      IERR )
 867 *
 868 *              Perform bidiagonal SVD, computing left singular vectors
 869 *              of bidiagonal matrix in WORK(IU) and computing right
 870 *              singular vectors of bidiagonal matrix in VT
 871 *              Workspace: need   N*N [U] + 3*N [e, tauq, taup] + BDSPAC
 872 *
 873                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
 874      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 875      $                      INFO )
 876 *
 877 *              Overwrite WORK(IU) by left singular vectors of R and VT
 878 *              by right singular vectors of R
 879 *              Workspace: need   N*N [U] + 3*N [e, tauq, taup] + N    [work]
 880 *              Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + N*NB [work]
 881 *
 882                CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
 883      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
 884      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
 885                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
 886      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 887      $                      LWORK - NWORK + 1, IERR )
 888 *
 889 *              Multiply Q in U by left singular vectors of R in
 890 *              WORK(IU), storing result in A
 891 *              Workspace: need   N*N [U]
 892 *
 893                CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
 894      $                     LDWRKU, ZERO, A, LDA )
 895 *
 896 *              Copy left singular vectors of A from A to U
 897 *
 898                CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
 899 *
 900             END IF
 901 *
 902          ELSE
 903 *
 904 *           M .LT. MNTHR
 905 *
 906 *           Path 5 (M >= N, but not much larger)
 907 *           Reduce to bidiagonal form without QR decomposition
 908 *
 909             IE = 1
 910             ITAUQ = IE + N
 911             ITAUP = ITAUQ + N
 912             NWORK = ITAUP + N
 913 *
 914 *           Bidiagonalize A
 915 *           Workspace: need   3*N [e, tauq, taup] + M        [work]
 916 *           Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work]
 917 *
 918             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 919      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 920      $                   IERR )
 921             IF( WNTQN ) THEN
 922 *
 923 *              Path 5n (M >= N, JOBZ='N')
 924 *              Perform bidiagonal SVD, only computing singular values
 925 *              Workspace: need   3*N [e, tauq, taup] + BDSPAC
 926 *
 927                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
 928      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
 929             ELSE IF( WNTQO ) THEN
 930 *              Path 5o (M >= N, JOBZ='O')
 931                IU = NWORK
 932                IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
 933 *
 934 *                 WORK( IU ) is M by N
 935 *
 936                   LDWRKU = M
 937                   NWORK = IU + LDWRKU*N
 938                   CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
 939      $                         LDWRKU )
 940 *                 IR is unused; silence compile warnings
 941                   IR = -1
 942                ELSE
 943 *
 944 *                 WORK( IU ) is N by N
 945 *
 946                   LDWRKU = N
 947                   NWORK = IU + LDWRKU*N
 948 *
 949 *                 WORK(IR) is LDWRKR by N
 950 *
 951                   IR = NWORK
 952                   LDWRKR = ( LWORK - N*N - 3*N ) / N
 953                END IF
 954                NWORK = IU + LDWRKU*N
 955 *
 956 *              Perform bidiagonal SVD, computing left singular vectors
 957 *              of bidiagonal matrix in WORK(IU) and computing right
 958 *              singular vectors of bidiagonal matrix in VT
 959 *              Workspace: need   3*N [e, tauq, taup] + N*N [U] + BDSPAC
 960 *
 961                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
 962      $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
 963      $                      IWORK, INFO )
 964 *
 965 *              Overwrite VT by right singular vectors of A
 966 *              Workspace: need   3*N [e, tauq, taup] + N*N [U] + N    [work]
 967 *              Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
 968 *
 969                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
 970      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 971      $                      LWORK - NWORK + 1, IERR )
 972 *
 973                IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
 974 *
 975 *                 Path 5o-fast
 976 *                 Overwrite WORK(IU) by left singular vectors of A
 977 *                 Workspace: need   3*N [e, tauq, taup] + M*N [U] + N    [work]
 978 *                 Workspace: prefer 3*N [e, tauq, taup] + M*N [U] + N*NB [work]
 979 *
 980                   CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
 981      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
 982      $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
 983 *
 984 *                 Copy left singular vectors of A from WORK(IU) to A
 985 *
 986                   CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
 987                ELSE
 988 *
 989 *                 Path 5o-slow
 990 *                 Generate Q in A
 991 *                 Workspace: need   3*N [e, tauq, taup] + N*N [U] + N    [work]
 992 *                 Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
 993 *
 994                   CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
 995      $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
 996 *
 997 *                 Multiply Q in A by left singular vectors of
 998 *                 bidiagonal matrix in WORK(IU), storing result in
 999 *                 WORK(IR) and copying to A
1000 *                 Workspace: need   3*N [e, tauq, taup] + N*N [U] + NB*N [R]
1001 *                 Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + M*N  [R]
1002 *
1003                   DO 20 I = 1, M, LDWRKR
1004                      CHUNK = MIN( M - I + 1, LDWRKR )
1005                      CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
1006      $                           LDA, WORK( IU ), LDWRKU, ZERO,
1007      $                           WORK( IR ), LDWRKR )
1008                      CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
1009      $                            A( I, 1 ), LDA )
1010    20             CONTINUE
1011                END IF
1012 *
1013             ELSE IF( WNTQS ) THEN
1014 *
1015 *              Path 5s (M >= N, JOBZ='S')
1016 *              Perform bidiagonal SVD, computing left singular vectors
1017 *              of bidiagonal matrix in U and computing right singular
1018 *              vectors of bidiagonal matrix in VT
1019 *              Workspace: need   3*N [e, tauq, taup] + BDSPAC
1020 *
1021                CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
1022                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
1023      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1024      $                      INFO )
1025 *
1026 *              Overwrite U by left singular vectors of A and VT
1027 *              by right singular vectors of A
1028 *              Workspace: need   3*N [e, tauq, taup] + N    [work]
1029 *              Workspace: prefer 3*N [e, tauq, taup] + N*NB [work]
1030 *
1031                CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1032      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1033      $                      LWORK - NWORK + 1, IERR )
1034                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
1035      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1036      $                      LWORK - NWORK + 1, IERR )
1037             ELSE IF( WNTQA ) THEN
1038 *
1039 *              Path 5a (M >= N, JOBZ='A')
1040 *              Perform bidiagonal SVD, computing left singular vectors
1041 *              of bidiagonal matrix in U and computing right singular
1042 *              vectors of bidiagonal matrix in VT
1043 *              Workspace: need   3*N [e, tauq, taup] + BDSPAC
1044 *
1045                CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
1046                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
1047      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1048      $                      INFO )
1049 *
1050 *              Set the right corner of U to identity matrix
1051 *
1052                IF( M.GT.N ) THEN
1053                   CALL DLASET( 'F', M - N, M - N, ZERO, ONE, U(N+1,N+1),
1054      $                         LDU )
1055                END IF
1056 *
1057 *              Overwrite U by left singular vectors of A and VT
1058 *              by right singular vectors of A
1059 *              Workspace: need   3*N [e, tauq, taup] + M    [work]
1060 *              Workspace: prefer 3*N [e, tauq, taup] + M*NB [work]
1061 *
1062                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1063      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1064      $                      LWORK - NWORK + 1, IERR )
1065                CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1066      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1067      $                      LWORK - NWORK + 1, IERR )
1068             END IF
1069 *
1070          END IF
1071 *
1072       ELSE
1073 *
1074 *        A has more columns than rows. If A has sufficiently more
1075 *        columns than rows, first reduce using the LQ decomposition (if
1076 *        sufficient workspace available)
1077 *
1078          IF( N.GE.MNTHR ) THEN
1079 *
1080             IF( WNTQN ) THEN
1081 *
1082 *              Path 1t (N >> M, JOBZ='N')
1083 *              No singular vectors to be computed
1084 *
1085                ITAU = 1
1086                NWORK = ITAU + M
1087 *
1088 *              Compute A=L*Q
1089 *              Workspace: need   M [tau] + M [work]
1090 *              Workspace: prefer M [tau] + M*NB [work]
1091 *
1092                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1093      $                      LWORK - NWORK + 1, IERR )
1094 *
1095 *              Zero out above L
1096 *
1097                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1098                IE = 1
1099                ITAUQ = IE + M
1100                ITAUP = ITAUQ + M
1101                NWORK = ITAUP + M
1102 *
1103 *              Bidiagonalize L in A
1104 *              Workspace: need   3*M [e, tauq, taup] + M      [work]
1105 *              Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work]
1106 *
1107                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1108      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1109      $                      IERR )
1110                NWORK = IE + M
1111 *
1112 *              Perform bidiagonal SVD, computing singular values only
1113 *              Workspace: need   M [e] + BDSPAC
1114 *
1115                CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1116      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1117 *
1118             ELSE IF( WNTQO ) THEN
1119 *
1120 *              Path 2t (N >> M, JOBZ='O')
1121 *              M right singular vectors to be overwritten on A and
1122 *              M left singular vectors to be computed in U
1123 *
1124                IVT = 1
1125 *
1126 *              WORK(IVT) is M by M
1127 *              WORK(IL)  is M by M; it is later resized to M by chunk for gemm
1128 *
1129                IL = IVT + M*M
1130                IF( LWORK .GE. M*N + M*M + 3*M + BDSPAC ) THEN
1131                   LDWRKL = M
1132                   CHUNK = N
1133                ELSE
1134                   LDWRKL = M
1135                   CHUNK = ( LWORK - M*M ) / M
1136                END IF
1137                ITAU = IL + LDWRKL*M
1138                NWORK = ITAU + M
1139 *
1140 *              Compute A=L*Q
1141 *              Workspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
1142 *              Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1143 *
1144                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1145      $                      LWORK - NWORK + 1, IERR )
1146 *
1147 *              Copy L to WORK(IL), zeroing about above it
1148 *
1149                CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1150                CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1151      $                      WORK( IL + LDWRKL ), LDWRKL )
1152 *
1153 *              Generate Q in A
1154 *              Workspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
1155 *              Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1156 *
1157                CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1158      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1159                IE = ITAU
1160                ITAUQ = IE + M
1161                ITAUP = ITAUQ + M
1162                NWORK = ITAUP + M
1163 *
1164 *              Bidiagonalize L in WORK(IL)
1165 *              Workspace: need   M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M      [work]
1166 *              Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1167 *
1168                CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1169      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1170      $                      LWORK - NWORK + 1, IERR )
1171 *
1172 *              Perform bidiagonal SVD, computing left singular vectors
1173 *              of bidiagonal matrix in U, and computing right singular
1174 *              vectors of bidiagonal matrix in WORK(IVT)
1175 *              Workspace: need   M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1176 *
1177                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1178      $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
1179      $                      IWORK, INFO )
1180 *
1181 *              Overwrite U by left singular vectors of L and WORK(IVT)
1182 *              by right singular vectors of L
1183 *              Workspace: need   M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M    [work]
1184 *              Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1185 *
1186                CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1187      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1188      $                      LWORK - NWORK + 1, IERR )
1189                CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1190      $                      WORK( ITAUP ), WORK( IVT ), M,
1191      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1192 *
1193 *              Multiply right singular vectors of L in WORK(IVT) by Q
1194 *              in A, storing result in WORK(IL) and copying to A
1195 *              Workspace: need   M*M [VT] + M*M [L]
1196 *              Workspace: prefer M*M [VT] + M*N [L]
1197 *              At this point, L is resized as M by chunk.
1198 *
1199                DO 30 I = 1, N, CHUNK
1200                   BLK = MIN( N - I + 1, CHUNK )
1201                   CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
1202      $                        A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1203                   CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1204      $                         A( 1, I ), LDA )
1205    30          CONTINUE
1206 *
1207             ELSE IF( WNTQS ) THEN
1208 *
1209 *              Path 3t (N >> M, JOBZ='S')
1210 *              M right singular vectors to be computed in VT and
1211 *              M left singular vectors to be computed in U
1212 *
1213                IL = 1
1214 *
1215 *              WORK(IL) is M by M
1216 *
1217                LDWRKL = M
1218                ITAU = IL + LDWRKL*M
1219                NWORK = ITAU + M
1220 *
1221 *              Compute A=L*Q
1222 *              Workspace: need   M*M [L] + M [tau] + M    [work]
1223 *              Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1224 *
1225                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1226      $                      LWORK - NWORK + 1, IERR )
1227 *
1228 *              Copy L to WORK(IL), zeroing out above it
1229 *
1230                CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1231                CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1232      $                      WORK( IL + LDWRKL ), LDWRKL )
1233 *
1234 *              Generate Q in A
1235 *              Workspace: need   M*M [L] + M [tau] + M    [work]
1236 *              Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1237 *
1238                CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1239      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1240                IE = ITAU
1241                ITAUQ = IE + M
1242                ITAUP = ITAUQ + M
1243                NWORK = ITAUP + M
1244 *
1245 *              Bidiagonalize L in WORK(IU).
1246 *              Workspace: need   M*M [L] + 3*M [e, tauq, taup] + M      [work]
1247 *              Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1248 *
1249                CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1250      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1251      $                      LWORK - NWORK + 1, IERR )
1252 *
1253 *              Perform bidiagonal SVD, computing left singular vectors
1254 *              of bidiagonal matrix in U and computing right singular
1255 *              vectors of bidiagonal matrix in VT
1256 *              Workspace: need   M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1257 *
1258                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
1259      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1260      $                      INFO )
1261 *
1262 *              Overwrite U by left singular vectors of L and VT
1263 *              by right singular vectors of L
1264 *              Workspace: need   M*M [L] + 3*M [e, tauq, taup] + M    [work]
1265 *              Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1266 *
1267                CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1268      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1269      $                      LWORK - NWORK + 1, IERR )
1270                CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1271      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1272      $                      LWORK - NWORK + 1, IERR )
1273 *
1274 *              Multiply right singular vectors of L in WORK(IL) by
1275 *              Q in A, storing result in VT
1276 *              Workspace: need   M*M [L]
1277 *
1278                CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1279                CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
1280      $                     A, LDA, ZERO, VT, LDVT )
1281 *
1282             ELSE IF( WNTQA ) THEN
1283 *
1284 *              Path 4t (N >> M, JOBZ='A')
1285 *              N right singular vectors to be computed in VT and
1286 *              M left singular vectors to be computed in U
1287 *
1288                IVT = 1
1289 *
1290 *              WORK(IVT) is M by M
1291 *
1292                LDWKVT = M
1293                ITAU = IVT + LDWKVT*M
1294                NWORK = ITAU + M
1295 *
1296 *              Compute A=L*Q, copying result to VT
1297 *              Workspace: need   M*M [VT] + M [tau] + M    [work]
1298 *              Workspace: prefer M*M [VT] + M [tau] + M*NB [work]
1299 *
1300                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1301      $                      LWORK - NWORK + 1, IERR )
1302                CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
1303 *
1304 *              Generate Q in VT
1305 *              Workspace: need   M*M [VT] + M [tau] + N    [work]
1306 *              Workspace: prefer M*M [VT] + M [tau] + N*NB [work]
1307 *
1308                CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1309      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1310 *
1311 *              Produce L in A, zeroing out other entries
1312 *
1313                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1314                IE = ITAU
1315                ITAUQ = IE + M
1316                ITAUP = ITAUQ + M
1317                NWORK = ITAUP + M
1318 *
1319 *              Bidiagonalize L in A
1320 *              Workspace: need   M*M [VT] + 3*M [e, tauq, taup] + M      [work]
1321 *              Workspace: prefer M*M [VT] + 3*M [e, tauq, taup] + 2*M*NB [work]
1322 *
1323                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1324      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1325      $                      IERR )
1326 *
1327 *              Perform bidiagonal SVD, computing left singular vectors
1328 *              of bidiagonal matrix in U and computing right singular
1329 *              vectors of bidiagonal matrix in WORK(IVT)
1330 *              Workspace: need   M*M [VT] + 3*M [e, tauq, taup] + BDSPAC
1331 *
1332                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1333      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1334      $                      WORK( NWORK ), IWORK, INFO )
1335 *
1336 *              Overwrite U by left singular vectors of L and WORK(IVT)
1337 *              by right singular vectors of L
1338 *              Workspace: need   M*M [VT] + 3*M [e, tauq, taup]+ M    [work]
1339 *              Workspace: prefer M*M [VT] + 3*M [e, tauq, taup]+ M*NB [work]
1340 *
1341                CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1342      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1343      $                      LWORK - NWORK + 1, IERR )
1344                CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1345      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1346      $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1347 *
1348 *              Multiply right singular vectors of L in WORK(IVT) by
1349 *              Q in VT, storing result in A
1350 *              Workspace: need   M*M [VT]
1351 *
1352                CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
1353      $                     VT, LDVT, ZERO, A, LDA )
1354 *
1355 *              Copy right singular vectors of A from A to VT
1356 *
1357                CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
1358 *
1359             END IF
1360 *
1361          ELSE
1362 *
1363 *           N .LT. MNTHR
1364 *
1365 *           Path 5t (N > M, but not much larger)
1366 *           Reduce to bidiagonal form without LQ decomposition
1367 *
1368             IE = 1
1369             ITAUQ = IE + M
1370             ITAUP = ITAUQ + M
1371             NWORK = ITAUP + M
1372 *
1373 *           Bidiagonalize A
1374 *           Workspace: need   3*M [e, tauq, taup] + N        [work]
1375 *           Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work]
1376 *
1377             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1378      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1379      $                   IERR )
1380             IF( WNTQN ) THEN
1381 *
1382 *              Path 5tn (N > M, JOBZ='N')
1383 *              Perform bidiagonal SVD, only computing singular values
1384 *              Workspace: need   3*M [e, tauq, taup] + BDSPAC
1385 *
1386                CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1387      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1388             ELSE IF( WNTQO ) THEN
1389 *              Path 5to (N > M, JOBZ='O')
1390                LDWKVT = M
1391                IVT = NWORK
1392                IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1393 *
1394 *                 WORK( IVT ) is M by N
1395 *
1396                   CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
1397      $                         LDWKVT )
1398                   NWORK = IVT + LDWKVT*N
1399 *                 IL is unused; silence compile warnings
1400                   IL = -1
1401                ELSE
1402 *
1403 *                 WORK( IVT ) is M by M
1404 *
1405                   NWORK = IVT + LDWKVT*M
1406                   IL = NWORK
1407 *
1408 *                 WORK(IL) is M by CHUNK
1409 *
1410                   CHUNK = ( LWORK - M*M - 3*M ) / M
1411                END IF
1412 *
1413 *              Perform bidiagonal SVD, computing left singular vectors
1414 *              of bidiagonal matrix in U and computing right singular
1415 *              vectors of bidiagonal matrix in WORK(IVT)
1416 *              Workspace: need   3*M [e, tauq, taup] + M*M [VT] + BDSPAC
1417 *
1418                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
1419      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1420      $                      WORK( NWORK ), IWORK, INFO )
1421 *
1422 *              Overwrite U by left singular vectors of A
1423 *              Workspace: need   3*M [e, tauq, taup] + M*M [VT] + M    [work]
1424 *              Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1425 *
1426                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1427      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1428      $                      LWORK - NWORK + 1, IERR )
1429 *
1430                IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1431 *
1432 *                 Path 5to-fast
1433 *                 Overwrite WORK(IVT) by left singular vectors of A
1434 *                 Workspace: need   3*M [e, tauq, taup] + M*N [VT] + M    [work]
1435 *                 Workspace: prefer 3*M [e, tauq, taup] + M*N [VT] + M*NB [work]
1436 *
1437                   CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1438      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
1439      $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
1440 *
1441 *                 Copy right singular vectors of A from WORK(IVT) to A
1442 *
1443                   CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1444                ELSE
1445 *
1446 *                 Path 5to-slow
1447 *                 Generate P**T in A
1448 *                 Workspace: need   3*M [e, tauq, taup] + M*M [VT] + M    [work]
1449 *                 Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1450 *
1451                   CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1452      $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
1453 *
1454 *                 Multiply Q in A by right singular vectors of
1455 *                 bidiagonal matrix in WORK(IVT), storing result in
1456 *                 WORK(IL) and copying to A
1457 *                 Workspace: need   3*M [e, tauq, taup] + M*M [VT] + M*NB [L]
1458 *                 Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*N  [L]
1459 *
1460                   DO 40 I = 1, N, CHUNK
1461                      BLK = MIN( N - I + 1, CHUNK )
1462                      CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
1463      $                           LDWKVT, A( 1, I ), LDA, ZERO,
1464      $                           WORK( IL ), M )
1465                      CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
1466      $                            LDA )
1467    40             CONTINUE
1468                END IF
1469             ELSE IF( WNTQS ) THEN
1470 *
1471 *              Path 5ts (N > M, JOBZ='S')
1472 *              Perform bidiagonal SVD, computing left singular vectors
1473 *              of bidiagonal matrix in U and computing right singular
1474 *              vectors of bidiagonal matrix in VT
1475 *              Workspace: need   3*M [e, tauq, taup] + BDSPAC
1476 *
1477                CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
1478                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1479      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1480      $                      INFO )
1481 *
1482 *              Overwrite U by left singular vectors of A and VT
1483 *              by right singular vectors of A
1484 *              Workspace: need   3*M [e, tauq, taup] + M    [work]
1485 *              Workspace: prefer 3*M [e, tauq, taup] + M*NB [work]
1486 *
1487                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1488      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1489      $                      LWORK - NWORK + 1, IERR )
1490                CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1491      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1492      $                      LWORK - NWORK + 1, IERR )
1493             ELSE IF( WNTQA ) THEN
1494 *
1495 *              Path 5ta (N > M, JOBZ='A')
1496 *              Perform bidiagonal SVD, computing left singular vectors
1497 *              of bidiagonal matrix in U and computing right singular
1498 *              vectors of bidiagonal matrix in VT
1499 *              Workspace: need   3*M [e, tauq, taup] + BDSPAC
1500 *
1501                CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
1502                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1503      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1504      $                      INFO )
1505 *
1506 *              Set the right corner of VT to identity matrix
1507 *
1508                IF( N.GT.M ) THEN
1509                   CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT(M+1,M+1),
1510      $                         LDVT )
1511                END IF
1512 *
1513 *              Overwrite U by left singular vectors of A and VT
1514 *              by right singular vectors of A
1515 *              Workspace: need   3*M [e, tauq, taup] + N    [work]
1516 *              Workspace: prefer 3*M [e, tauq, taup] + N*NB [work]
1517 *
1518                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1519      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1520      $                      LWORK - NWORK + 1, IERR )
1521                CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1522      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1523      $                      LWORK - NWORK + 1, IERR )
1524             END IF
1525 *
1526          END IF
1527 *
1528       END IF
1529 *
1530 *     Undo scaling if necessary
1531 *
1532       IF( ISCL.EQ.1 ) THEN
1533          IF( ANRM.GT.BIGNUM )
1534      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1535      $                   IERR )
1536          IF( ANRM.LT.SMLNUM )
1537      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1538      $                   IERR )
1539       END IF
1540 *
1541 *     Return optimal workspace in WORK(1)
1542 *
1543       WORK( 1 ) = MAXWRK
1544 *
1545       RETURN
1546 *
1547 *     End of DGESDD
1548 *
1549       END