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