SRC/zgesvd.f

   1 *> \brief <b> ZGESVD computes the singular value decomposition (SVD) for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGESVD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
  22 *                          WORK, LWORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBU, JOBVT
  26 *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   RWORK( * ), S( * )
  30 *       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
  31 *      $                   WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZGESVD computes the singular value decomposition (SVD) of a complex
  41 *> M-by-N matrix A, optionally computing the left and/or right singular
  42 *> vectors. The SVD is written
  43 *>
  44 *>      A = U * SIGMA * conjugate-transpose(V)
  45 *>
  46 *> where SIGMA is an M-by-N matrix which is zero except for its
  47 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
  48 *> V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
  49 *> are the singular values of A; they are real and non-negative, and
  50 *> are returned in descending order.  The first min(m,n) columns of
  51 *> U and V are the left and right singular vectors of A.
  52 *>
  53 *> Note that the routine returns V**H, not V.
  54 *> \endverbatim
  55 *
  56 *  Arguments:
  57 *  ==========
  58 *
  59 *> \param[in] JOBU
  60 *> \verbatim
  61 *>          JOBU is CHARACTER*1
  62 *>          Specifies options for computing all or part of the matrix U:
  63 *>          = 'A':  all M columns of U are returned in array U:
  64 *>          = 'S':  the first min(m,n) columns of U (the left singular
  65 *>                  vectors) are returned in the array U;
  66 *>          = 'O':  the first min(m,n) columns of U (the left singular
  67 *>                  vectors) are overwritten on the array A;
  68 *>          = 'N':  no columns of U (no left singular vectors) are
  69 *>                  computed.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] JOBVT
  73 *> \verbatim
  74 *>          JOBVT is CHARACTER*1
  75 *>          Specifies options for computing all or part of the matrix
  76 *>          V**H:
  77 *>          = 'A':  all N rows of V**H are returned in the array VT;
  78 *>          = 'S':  the first min(m,n) rows of V**H (the right singular
  79 *>                  vectors) are returned in the array VT;
  80 *>          = 'O':  the first min(m,n) rows of V**H (the right singular
  81 *>                  vectors) are overwritten on the array A;
  82 *>          = 'N':  no rows of V**H (no right singular vectors) are
  83 *>                  computed.
  84 *>
  85 *>          JOBVT and JOBU cannot both be 'O'.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] M
  89 *> \verbatim
  90 *>          M is INTEGER
  91 *>          The number of rows of the input matrix A.  M >= 0.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] N
  95 *> \verbatim
  96 *>          N is INTEGER
  97 *>          The number of columns of the input matrix A.  N >= 0.
  98 *> \endverbatim
  99 *>
 100 *> \param[in,out] A
 101 *> \verbatim
 102 *>          A is COMPLEX*16 array, dimension (LDA,N)
 103 *>          On entry, the M-by-N matrix A.
 104 *>          On exit,
 105 *>          if JOBU = 'O',  A is overwritten with the first min(m,n)
 106 *>                          columns of U (the left singular vectors,
 107 *>                          stored columnwise);
 108 *>          if JOBVT = 'O', A is overwritten with the first min(m,n)
 109 *>                          rows of V**H (the right singular vectors,
 110 *>                          stored rowwise);
 111 *>          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
 112 *>                          are destroyed.
 113 *> \endverbatim
 114 *>
 115 *> \param[in] LDA
 116 *> \verbatim
 117 *>          LDA is INTEGER
 118 *>          The leading dimension of the array A.  LDA >= max(1,M).
 119 *> \endverbatim
 120 *>
 121 *> \param[out] S
 122 *> \verbatim
 123 *>          S is DOUBLE PRECISION array, dimension (min(M,N))
 124 *>          The singular values of A, sorted so that S(i) >= S(i+1).
 125 *> \endverbatim
 126 *>
 127 *> \param[out] U
 128 *> \verbatim
 129 *>          U is COMPLEX*16 array, dimension (LDU,UCOL)
 130 *>          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
 131 *>          If JOBU = 'A', U contains the M-by-M unitary matrix U;
 132 *>          if JOBU = 'S', U contains the first min(m,n) columns of U
 133 *>          (the left singular vectors, stored columnwise);
 134 *>          if JOBU = 'N' or 'O', 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 *>          JOBU = 'S' or 'A', LDU >= M.
 142 *> \endverbatim
 143 *>
 144 *> \param[out] VT
 145 *> \verbatim
 146 *>          VT is COMPLEX*16 array, dimension (LDVT,N)
 147 *>          If JOBVT = 'A', VT contains the N-by-N unitary matrix
 148 *>          V**H;
 149 *>          if JOBVT = 'S', VT contains the first min(m,n) rows of
 150 *>          V**H (the right singular vectors, stored rowwise);
 151 *>          if JOBVT = 'N' or 'O', 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; if
 158 *>          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
 159 *> \endverbatim
 160 *>
 161 *> \param[out] WORK
 162 *> \verbatim
 163 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 164 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 165 *> \endverbatim
 166 *>
 167 *> \param[in] LWORK
 168 *> \verbatim
 169 *>          LWORK is INTEGER
 170 *>          The dimension of the array WORK.
 171 *>          LWORK >=  MAX(1,2*MIN(M,N)+MAX(M,N)).
 172 *>          For good performance, LWORK should generally be larger.
 173 *>
 174 *>          If LWORK = -1, then a workspace query is assumed; the routine
 175 *>          only calculates the optimal size of the WORK array, returns
 176 *>          this value as the first entry of the WORK array, and no error
 177 *>          message related to LWORK is issued by XERBLA.
 178 *> \endverbatim
 179 *>
 180 *> \param[out] RWORK
 181 *> \verbatim
 182 *>          RWORK is DOUBLE PRECISION array, dimension (5*min(M,N))
 183 *>          On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
 184 *>          unconverged superdiagonal elements of an upper bidiagonal
 185 *>          matrix B whose diagonal is in S (not necessarily sorted).
 186 *>          B satisfies A = U * B * VT, so it has the same singular
 187 *>          values as A, and singular vectors related by U and VT.
 188 *> \endverbatim
 189 *>
 190 *> \param[out] INFO
 191 *> \verbatim
 192 *>          INFO is INTEGER
 193 *>          = 0:  successful exit.
 194 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 195 *>          > 0:  if ZBDSQR did not converge, INFO specifies how many
 196 *>                superdiagonals of an intermediate bidiagonal form B
 197 *>                did not converge to zero. See the description of RWORK
 198 *>                above for details.
 199 *> \endverbatim
 200 *
 201 *  Authors:
 202 *  ========
 203 *
 204 *> \author Univ. of Tennessee
 205 *> \author Univ. of California Berkeley
 206 *> \author Univ. of Colorado Denver
 207 *> \author NAG Ltd.
 208 *
 209 *> \date April 2012
 210 *
 211 *> \ingroup complex16GEsing
 212 *
 213 *  =====================================================================
 214       SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
 215      $                   VT, LDVT, WORK, LWORK, RWORK, INFO )
 216 *
 217 *  -- LAPACK driver routine (version 3.6.1) --
 218 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 219 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 220 *     April 2012
 221 *
 222 *     .. Scalar Arguments ..
 223       CHARACTER          JOBU, JOBVT
 224       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
 225 *     ..
 226 *     .. Array Arguments ..
 227       DOUBLE PRECISION   RWORK( * ), S( * )
 228       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
 229      $                   WORK( * )
 230 *     ..
 231 *
 232 *  =====================================================================
 233 *
 234 *     .. Parameters ..
 235       COMPLEX*16         CZERO, CONE
 236       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
 237      $                   CONE = ( 1.0D0, 0.0D0 ) )
 238       DOUBLE PRECISION   ZERO, ONE
 239       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 240 *     ..
 241 *     .. Local Scalars ..
 242       LOGICAL            LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
 243      $                   WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS
 244       INTEGER            BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL,
 245      $                   ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
 246      $                   MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
 247      $                   NRVT, WRKBL
 248       INTEGER            LWORK_ZGEQRF, LWORK_ZUNGQR_N, LWORK_ZUNGQR_M,
 249      $                   LWORK_ZGEBRD, LWORK_ZUNGBR_P, LWORK_ZUNGBR_Q,
 250      $                   LWORK_ZGELQF, LWORK_ZUNGLQ_N, LWORK_ZUNGLQ_M
 251       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
 252 *     ..
 253 *     .. Local Arrays ..
 254       DOUBLE PRECISION   DUM( 1 )
 255       COMPLEX*16         CDUM( 1 )
 256 *     ..
 257 *     .. External Subroutines ..
 258       EXTERNAL           DLASCL, XERBLA, ZBDSQR, ZGEBRD, ZGELQF, ZGEMM,
 259      $                   ZGEQRF, ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNGLQ,
 260      $                   ZUNGQR, ZUNMBR
 261 *     ..
 262 *     .. External Functions ..
 263       LOGICAL            LSAME
 264       INTEGER            ILAENV
 265       DOUBLE PRECISION   DLAMCH, ZLANGE
 266       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
 267 *     ..
 268 *     .. Intrinsic Functions ..
 269       INTRINSIC          MAX, MIN, SQRT
 270 *     ..
 271 *     .. Executable Statements ..
 272 *
 273 *     Test the input arguments
 274 *
 275       INFO = 0
 276       MINMN = MIN( M, N )
 277       WNTUA = LSAME( JOBU, 'A' )
 278       WNTUS = LSAME( JOBU, 'S' )
 279       WNTUAS = WNTUA .OR. WNTUS
 280       WNTUO = LSAME( JOBU, 'O' )
 281       WNTUN = LSAME( JOBU, 'N' )
 282       WNTVA = LSAME( JOBVT, 'A' )
 283       WNTVS = LSAME( JOBVT, 'S' )
 284       WNTVAS = WNTVA .OR. WNTVS
 285       WNTVO = LSAME( JOBVT, 'O' )
 286       WNTVN = LSAME( JOBVT, 'N' )
 287       LQUERY = ( LWORK.EQ.-1 )
 288 *
 289       IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN
 290          INFO = -1
 291       ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR.
 292      $         ( WNTVO .AND. WNTUO ) ) THEN
 293          INFO = -2
 294       ELSE IF( M.LT.0 ) THEN
 295          INFO = -3
 296       ELSE IF( N.LT.0 ) THEN
 297          INFO = -4
 298       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 299          INFO = -6
 300       ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN
 301          INFO = -9
 302       ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR.
 303      $         ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN
 304          INFO = -11
 305       END IF
 306 *
 307 *     Compute workspace
 308 *      (Note: Comments in the code beginning "Workspace:" describe the
 309 *       minimal amount of workspace needed at that point in the code,
 310 *       as well as the preferred amount for good performance.
 311 *       CWorkspace refers to complex workspace, and RWorkspace to
 312 *       real workspace. NB refers to the optimal block size for the
 313 *       immediately following subroutine, as returned by ILAENV.)
 314 *
 315       IF( INFO.EQ.0 ) THEN
 316          MINWRK = 1
 317          MAXWRK = 1
 318          IF( M.GE.N .AND. MINMN.GT.0 ) THEN
 319 *
 320 *           Space needed for ZBDSQR is BDSPAC = 5*N
 321 *
 322             MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
 323 *           Compute space needed for ZGEQRF
 324             CALL ZGEQRF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
 325             LWORK_ZGEQRF = INT( CDUM(1) )
 326 *           Compute space needed for ZUNGQR
 327             CALL ZUNGQR( M, N, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
 328             LWORK_ZUNGQR_N = INT( CDUM(1) )
 329             CALL ZUNGQR( M, M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
 330             LWORK_ZUNGQR_M = INT( CDUM(1) )
 331 *           Compute space needed for ZGEBRD
 332             CALL ZGEBRD( N, N, A, LDA, S, DUM(1), CDUM(1),
 333      $                   CDUM(1), CDUM(1), -1, IERR )
 334             LWORK_ZGEBRD = INT( CDUM(1) )
 335 *           Compute space needed for ZUNGBR
 336             CALL ZUNGBR( 'P', N, N, N, A, LDA, CDUM(1),
 337      $                   CDUM(1), -1, IERR )
 338             LWORK_ZUNGBR_P = INT( CDUM(1) )
 339             CALL ZUNGBR( 'Q', N, N, N, A, LDA, CDUM(1),
 340      $                   CDUM(1), -1, IERR )
 341             LWORK_ZUNGBR_Q = INT( CDUM(1) )
 342 *
 343             IF( M.GE.MNTHR ) THEN
 344                IF( WNTUN ) THEN
 345 *
 346 *                 Path 1 (M much larger than N, JOBU='N')
 347 *
 348                   MAXWRK = N + LWORK_ZGEQRF
 349                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZGEBRD )
 350                   IF( WNTVO .OR. WNTVAS )
 351      $               MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P )
 352                   MINWRK = 3*N
 353                ELSE IF( WNTUO .AND. WNTVN ) THEN
 354 *
 355 *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
 356 *
 357                   WRKBL = N + LWORK_ZGEQRF
 358                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N )
 359                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 360                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 361                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
 362                   MINWRK = 2*N + M
 363                ELSE IF( WNTUO .AND. WNTVAS ) THEN
 364 *
 365 *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
 366 *                 'A')
 367 *
 368                   WRKBL = N + LWORK_ZGEQRF
 369                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N )
 370                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 371                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 372                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P )
 373                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
 374                   MINWRK = 2*N + M
 375                ELSE IF( WNTUS .AND. WNTVN ) THEN
 376 *
 377 *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
 378 *
 379                   WRKBL = N + LWORK_ZGEQRF
 380                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N )
 381                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 382                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 383                   MAXWRK = N*N + WRKBL
 384                   MINWRK = 2*N + M
 385                ELSE IF( WNTUS .AND. WNTVO ) THEN
 386 *
 387 *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
 388 *
 389                   WRKBL = N + LWORK_ZGEQRF
 390                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N )
 391                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 392                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 393                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P )
 394                   MAXWRK = 2*N*N + WRKBL
 395                   MINWRK = 2*N + M
 396                ELSE IF( WNTUS .AND. WNTVAS ) THEN
 397 *
 398 *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
 399 *                 'A')
 400 *
 401                   WRKBL = N + LWORK_ZGEQRF
 402                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N )
 403                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 404                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 405                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P )
 406                   MAXWRK = N*N + WRKBL
 407                   MINWRK = 2*N + M
 408                ELSE IF( WNTUA .AND. WNTVN ) THEN
 409 *
 410 *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
 411 *
 412                   WRKBL = N + LWORK_ZGEQRF
 413                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M )
 414                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 415                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 416                   MAXWRK = N*N + WRKBL
 417                   MINWRK = 2*N + M
 418                ELSE IF( WNTUA .AND. WNTVO ) THEN
 419 *
 420 *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
 421 *
 422                   WRKBL = N + LWORK_ZGEQRF
 423                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M )
 424                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 425                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 426                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P )
 427                   MAXWRK = 2*N*N + WRKBL
 428                   MINWRK = 2*N + M
 429                ELSE IF( WNTUA .AND. WNTVAS ) THEN
 430 *
 431 *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
 432 *                 'A')
 433 *
 434                   WRKBL = N + LWORK_ZGEQRF
 435                   WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M )
 436                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD )
 437                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q )
 438                   WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P )
 439                   MAXWRK = N*N + WRKBL
 440                   MINWRK = 2*N + M
 441                END IF
 442             ELSE
 443 *
 444 *              Path 10 (M at least N, but not much larger)
 445 *
 446                CALL ZGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1),
 447      $                   CDUM(1), CDUM(1), -1, IERR )
 448                LWORK_ZGEBRD = INT( CDUM(1) )
 449                MAXWRK = 2*N + LWORK_ZGEBRD
 450                IF( WNTUS .OR. WNTUO ) THEN
 451                   CALL ZUNGBR( 'Q', M, N, N, A, LDA, CDUM(1),
 452      $                   CDUM(1), -1, IERR )
 453                   LWORK_ZUNGBR_Q = INT( CDUM(1) )
 454                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q )
 455                END IF
 456                IF( WNTUA ) THEN
 457                   CALL ZUNGBR( 'Q', M, M, N, A, LDA, CDUM(1),
 458      $                   CDUM(1), -1, IERR )
 459                   LWORK_ZUNGBR_Q = INT( CDUM(1) )
 460                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q )
 461                END IF
 462                IF( .NOT.WNTVN ) THEN
 463                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P )
 464                END IF
 465                MINWRK = 2*N + M
 466             END IF
 467          ELSE IF( MINMN.GT.0 ) THEN
 468 *
 469 *           Space needed for ZBDSQR is BDSPAC = 5*M
 470 *
 471             MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
 472 *           Compute space needed for ZGELQF
 473             CALL ZGELQF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
 474             LWORK_ZGELQF = INT( CDUM(1) )
 475 *           Compute space needed for ZUNGLQ
 476             CALL ZUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1), -1,
 477      $                   IERR )
 478             LWORK_ZUNGLQ_N = INT( CDUM(1) )
 479             CALL ZUNGLQ( M, N, M, A, LDA, CDUM(1), CDUM(1), -1, IERR )
 480             LWORK_ZUNGLQ_M = INT( CDUM(1) )
 481 *           Compute space needed for ZGEBRD
 482             CALL ZGEBRD( M, M, A, LDA, S, DUM(1), CDUM(1),
 483      $                   CDUM(1), CDUM(1), -1, IERR )
 484             LWORK_ZGEBRD = INT( CDUM(1) )
 485 *            Compute space needed for ZUNGBR P
 486             CALL ZUNGBR( 'P', M, M, M, A, N, CDUM(1),
 487      $                   CDUM(1), -1, IERR )
 488             LWORK_ZUNGBR_P = INT( CDUM(1) )
 489 *           Compute space needed for ZUNGBR Q
 490             CALL ZUNGBR( 'Q', M, M, M, A, N, CDUM(1),
 491      $                   CDUM(1), -1, IERR )
 492             LWORK_ZUNGBR_Q = INT( CDUM(1) )
 493             IF( N.GE.MNTHR ) THEN
 494                IF( WNTVN ) THEN
 495 *
 496 *                 Path 1t(N much larger than M, JOBVT='N')
 497 *
 498                   MAXWRK = M + LWORK_ZGELQF
 499                   MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZGEBRD )
 500                   IF( WNTUO .OR. WNTUAS )
 501      $               MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q )
 502                   MINWRK = 3*M
 503                ELSE IF( WNTVO .AND. WNTUN ) THEN
 504 *
 505 *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
 506 *
 507                   WRKBL = M + LWORK_ZGELQF
 508                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M )
 509                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 510                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 511                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
 512                   MINWRK = 2*M + N
 513                ELSE IF( WNTVO .AND. WNTUAS ) THEN
 514 *
 515 *                 Path 3t(N much larger than M, JOBU='S' or 'A',
 516 *                 JOBVT='O')
 517 *
 518                   WRKBL = M + LWORK_ZGELQF
 519                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M )
 520                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 521                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 522                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q )
 523                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
 524                   MINWRK = 2*M + N
 525                ELSE IF( WNTVS .AND. WNTUN ) THEN
 526 *
 527 *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
 528 *
 529                   WRKBL = M + LWORK_ZGELQF
 530                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M )
 531                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 532                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 533                   MAXWRK = M*M + WRKBL
 534                   MINWRK = 2*M + N
 535                ELSE IF( WNTVS .AND. WNTUO ) THEN
 536 *
 537 *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
 538 *
 539                   WRKBL = M + LWORK_ZGELQF
 540                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M )
 541                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 542                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 543                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q )
 544                   MAXWRK = 2*M*M + WRKBL
 545                   MINWRK = 2*M + N
 546                ELSE IF( WNTVS .AND. WNTUAS ) THEN
 547 *
 548 *                 Path 6t(N much larger than M, JOBU='S' or 'A',
 549 *                 JOBVT='S')
 550 *
 551                   WRKBL = M + LWORK_ZGELQF
 552                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M )
 553                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 554                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 555                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q )
 556                   MAXWRK = M*M + WRKBL
 557                   MINWRK = 2*M + N
 558                ELSE IF( WNTVA .AND. WNTUN ) THEN
 559 *
 560 *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
 561 *
 562                   WRKBL = M + LWORK_ZGELQF
 563                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N )
 564                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 565                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 566                   MAXWRK = M*M + WRKBL
 567                   MINWRK = 2*M + N
 568                ELSE IF( WNTVA .AND. WNTUO ) THEN
 569 *
 570 *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
 571 *
 572                   WRKBL = M + LWORK_ZGELQF
 573                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N )
 574                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 575                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 576                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q )
 577                   MAXWRK = 2*M*M + WRKBL
 578                   MINWRK = 2*M + N
 579                ELSE IF( WNTVA .AND. WNTUAS ) THEN
 580 *
 581 *                 Path 9t(N much larger than M, JOBU='S' or 'A',
 582 *                 JOBVT='A')
 583 *
 584                   WRKBL = M + LWORK_ZGELQF
 585                   WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N )
 586                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD )
 587                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P )
 588                   WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q )
 589                   MAXWRK = M*M + WRKBL
 590                   MINWRK = 2*M + N
 591                END IF
 592             ELSE
 593 *
 594 *              Path 10t(N greater than M, but not much larger)
 595 *
 596                CALL ZGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1),
 597      $                   CDUM(1), CDUM(1), -1, IERR )
 598                LWORK_ZGEBRD = INT( CDUM(1) )
 599                MAXWRK = 2*M + LWORK_ZGEBRD
 600                IF( WNTVS .OR. WNTVO ) THEN
 601 *                Compute space needed for ZUNGBR P
 602                  CALL ZUNGBR( 'P', M, N, M, A, N, CDUM(1),
 603      $                   CDUM(1), -1, IERR )
 604                  LWORK_ZUNGBR_P = INT( CDUM(1) )
 605                  MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P )
 606                END IF
 607                IF( WNTVA ) THEN
 608                  CALL ZUNGBR( 'P', N,  N, M, A, N, CDUM(1),
 609      $                   CDUM(1), -1, IERR )
 610                  LWORK_ZUNGBR_P = INT( CDUM(1) )
 611                  MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P )
 612                END IF
 613                IF( .NOT.WNTUN ) THEN
 614                   MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q )
 615                END IF
 616                MINWRK = 2*M + N
 617             END IF
 618          END IF
 619          MAXWRK = MAX( MAXWRK, MINWRK )
 620          WORK( 1 ) = MAXWRK
 621 *
 622          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 623             INFO = -13
 624          END IF
 625       END IF
 626 *
 627       IF( INFO.NE.0 ) THEN
 628          CALL XERBLA( 'ZGESVD', -INFO )
 629          RETURN
 630       ELSE IF( LQUERY ) THEN
 631          RETURN
 632       END IF
 633 *
 634 *     Quick return if possible
 635 *
 636       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 637          RETURN
 638       END IF
 639 *
 640 *     Get machine constants
 641 *
 642       EPS = DLAMCH( 'P' )
 643       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
 644       BIGNUM = ONE / SMLNUM
 645 *
 646 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 647 *
 648       ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
 649       ISCL = 0
 650       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 651          ISCL = 1
 652          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
 653       ELSE IF( ANRM.GT.BIGNUM ) THEN
 654          ISCL = 1
 655          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
 656       END IF
 657 *
 658       IF( M.GE.N ) THEN
 659 *
 660 *        A has at least as many rows as columns. If A has sufficiently
 661 *        more rows than columns, first reduce using the QR
 662 *        decomposition (if sufficient workspace available)
 663 *
 664          IF( M.GE.MNTHR ) THEN
 665 *
 666             IF( WNTUN ) THEN
 667 *
 668 *              Path 1 (M much larger than N, JOBU='N')
 669 *              No left singular vectors to be computed
 670 *
 671                ITAU = 1
 672                IWORK = ITAU + N
 673 *
 674 *              Compute A=Q*R
 675 *              (CWorkspace: need 2*N, prefer N+N*NB)
 676 *              (RWorkspace: need 0)
 677 *
 678                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
 679      $                      LWORK-IWORK+1, IERR )
 680 *
 681 *              Zero out below R
 682 *
 683                IF( N .GT. 1 ) THEN
 684                   CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
 685      $                         LDA )
 686                END IF
 687                IE = 1
 688                ITAUQ = 1
 689                ITAUP = ITAUQ + N
 690                IWORK = ITAUP + N
 691 *
 692 *              Bidiagonalize R in A
 693 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
 694 *              (RWorkspace: need N)
 695 *
 696                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
 697      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
 698      $                      IERR )
 699                NCVT = 0
 700                IF( WNTVO .OR. WNTVAS ) THEN
 701 *
 702 *                 If right singular vectors desired, generate P'.
 703 *                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 704 *                 (RWorkspace: 0)
 705 *
 706                   CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
 707      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 708                   NCVT = N
 709                END IF
 710                IRWORK = IE + N
 711 *
 712 *              Perform bidiagonal QR iteration, computing right
 713 *              singular vectors of A in A if desired
 714 *              (CWorkspace: 0)
 715 *              (RWorkspace: need BDSPAC)
 716 *
 717                CALL ZBDSQR( 'U', N, NCVT, 0, 0, S, RWORK( IE ), A, LDA,
 718      $                      CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
 719 *
 720 *              If right singular vectors desired in VT, copy them there
 721 *
 722                IF( WNTVAS )
 723      $            CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
 724 *
 725             ELSE IF( WNTUO .AND. WNTVN ) THEN
 726 *
 727 *              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
 728 *              N left singular vectors to be overwritten on A and
 729 *              no right singular vectors to be computed
 730 *
 731                IF( LWORK.GE.N*N+3*N ) THEN
 732 *
 733 *                 Sufficient workspace for a fast algorithm
 734 *
 735                   IR = 1
 736                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
 737 *
 738 *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
 739 *
 740                      LDWRKU = LDA
 741                      LDWRKR = LDA
 742                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
 743 *
 744 *                    WORK(IU) is LDA by N, WORK(IR) is N by N
 745 *
 746                      LDWRKU = LDA
 747                      LDWRKR = N
 748                   ELSE
 749 *
 750 *                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N
 751 *
 752                      LDWRKU = ( LWORK-N*N ) / N
 753                      LDWRKR = N
 754                   END IF
 755                   ITAU = IR + LDWRKR*N
 756                   IWORK = ITAU + N
 757 *
 758 *                 Compute A=Q*R
 759 *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
 760 *                 (RWorkspace: 0)
 761 *
 762                   CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
 763      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 764 *
 765 *                 Copy R to WORK(IR) and zero out below it
 766 *
 767                   CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
 768                   CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
 769      $                         WORK( IR+1 ), LDWRKR )
 770 *
 771 *                 Generate Q in A
 772 *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
 773 *                 (RWorkspace: 0)
 774 *
 775                   CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
 776      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 777                   IE = 1
 778                   ITAUQ = ITAU
 779                   ITAUP = ITAUQ + N
 780                   IWORK = ITAUP + N
 781 *
 782 *                 Bidiagonalize R in WORK(IR)
 783 *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
 784 *                 (RWorkspace: need N)
 785 *
 786                   CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
 787      $                         WORK( ITAUQ ), WORK( ITAUP ),
 788      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 789 *
 790 *                 Generate left vectors bidiagonalizing R
 791 *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
 792 *                 (RWorkspace: need 0)
 793 *
 794                   CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
 795      $                         WORK( ITAUQ ), WORK( IWORK ),
 796      $                         LWORK-IWORK+1, IERR )
 797                   IRWORK = IE + N
 798 *
 799 *                 Perform bidiagonal QR iteration, computing left
 800 *                 singular vectors of R in WORK(IR)
 801 *                 (CWorkspace: need N*N)
 802 *                 (RWorkspace: need BDSPAC)
 803 *
 804                   CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1,
 805      $                         WORK( IR ), LDWRKR, CDUM, 1,
 806      $                         RWORK( IRWORK ), INFO )
 807                   IU = ITAUQ
 808 *
 809 *                 Multiply Q in A by left singular vectors of R in
 810 *                 WORK(IR), storing result in WORK(IU) and copying to A
 811 *                 (CWorkspace: need N*N+N, prefer N*N+M*N)
 812 *                 (RWorkspace: 0)
 813 *
 814                   DO 10 I = 1, M, LDWRKU
 815                      CHUNK = MIN( M-I+1, LDWRKU )
 816                      CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
 817      $                           LDA, WORK( IR ), LDWRKR, CZERO,
 818      $                           WORK( IU ), LDWRKU )
 819                      CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
 820      $                            A( I, 1 ), LDA )
 821    10             CONTINUE
 822 *
 823                ELSE
 824 *
 825 *                 Insufficient workspace for a fast algorithm
 826 *
 827                   IE = 1
 828                   ITAUQ = 1
 829                   ITAUP = ITAUQ + N
 830                   IWORK = ITAUP + N
 831 *
 832 *                 Bidiagonalize A
 833 *                 (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
 834 *                 (RWorkspace: N)
 835 *
 836                   CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ),
 837      $                         WORK( ITAUQ ), WORK( ITAUP ),
 838      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 839 *
 840 *                 Generate left vectors bidiagonalizing A
 841 *                 (CWorkspace: need 3*N, prefer 2*N+N*NB)
 842 *                 (RWorkspace: 0)
 843 *
 844                   CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
 845      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 846                   IRWORK = IE + N
 847 *
 848 *                 Perform bidiagonal QR iteration, computing left
 849 *                 singular vectors of A in A
 850 *                 (CWorkspace: need 0)
 851 *                 (RWorkspace: need BDSPAC)
 852 *
 853                   CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1,
 854      $                         A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
 855 *
 856                END IF
 857 *
 858             ELSE IF( WNTUO .AND. WNTVAS ) THEN
 859 *
 860 *              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
 861 *              N left singular vectors to be overwritten on A and
 862 *              N right singular vectors to be computed in VT
 863 *
 864                IF( LWORK.GE.N*N+3*N ) THEN
 865 *
 866 *                 Sufficient workspace for a fast algorithm
 867 *
 868                   IR = 1
 869                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
 870 *
 871 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
 872 *
 873                      LDWRKU = LDA
 874                      LDWRKR = LDA
 875                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
 876 *
 877 *                    WORK(IU) is LDA by N and WORK(IR) is N by N
 878 *
 879                      LDWRKU = LDA
 880                      LDWRKR = N
 881                   ELSE
 882 *
 883 *                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N
 884 *
 885                      LDWRKU = ( LWORK-N*N ) / N
 886                      LDWRKR = N
 887                   END IF
 888                   ITAU = IR + LDWRKR*N
 889                   IWORK = ITAU + N
 890 *
 891 *                 Compute A=Q*R
 892 *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
 893 *                 (RWorkspace: 0)
 894 *
 895                   CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
 896      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 897 *
 898 *                 Copy R to VT, zeroing out below it
 899 *
 900                   CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
 901                   IF( N.GT.1 )
 902      $               CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
 903      $                            VT( 2, 1 ), LDVT )
 904 *
 905 *                 Generate Q in A
 906 *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
 907 *                 (RWorkspace: 0)
 908 *
 909                   CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
 910      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 911                   IE = 1
 912                   ITAUQ = ITAU
 913                   ITAUP = ITAUQ + N
 914                   IWORK = ITAUP + N
 915 *
 916 *                 Bidiagonalize R in VT, copying result to WORK(IR)
 917 *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
 918 *                 (RWorkspace: need N)
 919 *
 920                   CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
 921      $                         WORK( ITAUQ ), WORK( ITAUP ),
 922      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 923                   CALL ZLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
 924 *
 925 *                 Generate left vectors bidiagonalizing R in WORK(IR)
 926 *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
 927 *                 (RWorkspace: 0)
 928 *
 929                   CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
 930      $                         WORK( ITAUQ ), WORK( IWORK ),
 931      $                         LWORK-IWORK+1, IERR )
 932 *
 933 *                 Generate right vectors bidiagonalizing R in VT
 934 *                 (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
 935 *                 (RWorkspace: 0)
 936 *
 937                   CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
 938      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 939                   IRWORK = IE + N
 940 *
 941 *                 Perform bidiagonal QR iteration, computing left
 942 *                 singular vectors of R in WORK(IR) and computing right
 943 *                 singular vectors of R in VT
 944 *                 (CWorkspace: need N*N)
 945 *                 (RWorkspace: need BDSPAC)
 946 *
 947                   CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
 948      $                         LDVT, WORK( IR ), LDWRKR, CDUM, 1,
 949      $                         RWORK( IRWORK ), INFO )
 950                   IU = ITAUQ
 951 *
 952 *                 Multiply Q in A by left singular vectors of R in
 953 *                 WORK(IR), storing result in WORK(IU) and copying to A
 954 *                 (CWorkspace: need N*N+N, prefer N*N+M*N)
 955 *                 (RWorkspace: 0)
 956 *
 957                   DO 20 I = 1, M, LDWRKU
 958                      CHUNK = MIN( M-I+1, LDWRKU )
 959                      CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
 960      $                           LDA, WORK( IR ), LDWRKR, CZERO,
 961      $                           WORK( IU ), LDWRKU )
 962                      CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
 963      $                            A( I, 1 ), LDA )
 964    20             CONTINUE
 965 *
 966                ELSE
 967 *
 968 *                 Insufficient workspace for a fast algorithm
 969 *
 970                   ITAU = 1
 971                   IWORK = ITAU + N
 972 *
 973 *                 Compute A=Q*R
 974 *                 (CWorkspace: need 2*N, prefer N+N*NB)
 975 *                 (RWorkspace: 0)
 976 *
 977                   CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
 978      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 979 *
 980 *                 Copy R to VT, zeroing out below it
 981 *
 982                   CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
 983                   IF( N.GT.1 )
 984      $               CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
 985      $                            VT( 2, 1 ), LDVT )
 986 *
 987 *                 Generate Q in A
 988 *                 (CWorkspace: need 2*N, prefer N+N*NB)
 989 *                 (RWorkspace: 0)
 990 *
 991                   CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
 992      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 993                   IE = 1
 994                   ITAUQ = ITAU
 995                   ITAUP = ITAUQ + N
 996                   IWORK = ITAUP + N
 997 *
 998 *                 Bidiagonalize R in VT
 999 *                 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1000 *                 (RWorkspace: N)
1001 *
1002                   CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
1003      $                         WORK( ITAUQ ), WORK( ITAUP ),
1004      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
1005 *
1006 *                 Multiply Q in A by left vectors bidiagonalizing R
1007 *                 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1008 *                 (RWorkspace: 0)
1009 *
1010                   CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
1011      $                         WORK( ITAUQ ), A, LDA, WORK( IWORK ),
1012      $                         LWORK-IWORK+1, IERR )
1013 *
1014 *                 Generate right vectors bidiagonalizing R in VT
1015 *                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1016 *                 (RWorkspace: 0)
1017 *
1018                   CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1019      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
1020                   IRWORK = IE + N
1021 *
1022 *                 Perform bidiagonal QR iteration, computing left
1023 *                 singular vectors of A in A and computing right
1024 *                 singular vectors of A in VT
1025 *                 (CWorkspace: 0)
1026 *                 (RWorkspace: need BDSPAC)
1027 *
1028                   CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
1029      $                         LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
1030      $                         INFO )
1031 *
1032                END IF
1033 *
1034             ELSE IF( WNTUS ) THEN
1035 *
1036                IF( WNTVN ) THEN
1037 *
1038 *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
1039 *                 N left singular vectors to be computed in U and
1040 *                 no right singular vectors to be computed
1041 *
1042                   IF( LWORK.GE.N*N+3*N ) THEN
1043 *
1044 *                    Sufficient workspace for a fast algorithm
1045 *
1046                      IR = 1
1047                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
1048 *
1049 *                       WORK(IR) is LDA by N
1050 *
1051                         LDWRKR = LDA
1052                      ELSE
1053 *
1054 *                       WORK(IR) is N by N
1055 *
1056                         LDWRKR = N
1057                      END IF
1058                      ITAU = IR + LDWRKR*N
1059                      IWORK = ITAU + N
1060 *
1061 *                    Compute A=Q*R
1062 *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1063 *                    (RWorkspace: 0)
1064 *
1065                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1066      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1067 *
1068 *                    Copy R to WORK(IR), zeroing out below it
1069 *
1070                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ),
1071      $                            LDWRKR )
1072                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1073      $                            WORK( IR+1 ), LDWRKR )
1074 *
1075 *                    Generate Q in A
1076 *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1077 *                    (RWorkspace: 0)
1078 *
1079                      CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
1080      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1081                      IE = 1
1082                      ITAUQ = ITAU
1083                      ITAUP = ITAUQ + N
1084                      IWORK = ITAUP + N
1085 *
1086 *                    Bidiagonalize R in WORK(IR)
1087 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1088 *                    (RWorkspace: need N)
1089 *
1090                      CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S,
1091      $                            RWORK( IE ), WORK( ITAUQ ),
1092      $                            WORK( ITAUP ), WORK( IWORK ),
1093      $                            LWORK-IWORK+1, IERR )
1094 *
1095 *                    Generate left vectors bidiagonalizing R in WORK(IR)
1096 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1097 *                    (RWorkspace: 0)
1098 *
1099                      CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
1100      $                            WORK( ITAUQ ), WORK( IWORK ),
1101      $                            LWORK-IWORK+1, IERR )
1102                      IRWORK = IE + N
1103 *
1104 *                    Perform bidiagonal QR iteration, computing left
1105 *                    singular vectors of R in WORK(IR)
1106 *                    (CWorkspace: need N*N)
1107 *                    (RWorkspace: need BDSPAC)
1108 *
1109                      CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
1110      $                            1, WORK( IR ), LDWRKR, CDUM, 1,
1111      $                            RWORK( IRWORK ), INFO )
1112 *
1113 *                    Multiply Q in A by left singular vectors of R in
1114 *                    WORK(IR), storing result in U
1115 *                    (CWorkspace: need N*N)
1116 *                    (RWorkspace: 0)
1117 *
1118                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
1119      $                           WORK( IR ), LDWRKR, CZERO, U, LDU )
1120 *
1121                   ELSE
1122 *
1123 *                    Insufficient workspace for a fast algorithm
1124 *
1125                      ITAU = 1
1126                      IWORK = ITAU + N
1127 *
1128 *                    Compute A=Q*R, copying result to U
1129 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1130 *                    (RWorkspace: 0)
1131 *
1132                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1133      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1134                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1135 *
1136 *                    Generate Q in U
1137 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1138 *                    (RWorkspace: 0)
1139 *
1140                      CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
1141      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1142                      IE = 1
1143                      ITAUQ = ITAU
1144                      ITAUP = ITAUQ + N
1145                      IWORK = ITAUP + N
1146 *
1147 *                    Zero out below R in A
1148 *
1149                      IF( N .GT. 1 ) THEN
1150                         CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1151      $                               A( 2, 1 ), LDA )
1152                      END IF
1153 *
1154 *                    Bidiagonalize R in A
1155 *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1156 *                    (RWorkspace: need N)
1157 *
1158                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
1159      $                            WORK( ITAUQ ), WORK( ITAUP ),
1160      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1161 *
1162 *                    Multiply Q in U by left vectors bidiagonalizing R
1163 *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1164 *                    (RWorkspace: 0)
1165 *
1166                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1167      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1168      $                            LWORK-IWORK+1, IERR )
1169                      IRWORK = IE + N
1170 *
1171 *                    Perform bidiagonal QR iteration, computing left
1172 *                    singular vectors of A in U
1173 *                    (CWorkspace: 0)
1174 *                    (RWorkspace: need BDSPAC)
1175 *
1176                      CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
1177      $                            1, U, LDU, CDUM, 1, RWORK( IRWORK ),
1178      $                            INFO )
1179 *
1180                   END IF
1181 *
1182                ELSE IF( WNTVO ) THEN
1183 *
1184 *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1185 *                 N left singular vectors to be computed in U and
1186 *                 N right singular vectors to be overwritten on A
1187 *
1188                   IF( LWORK.GE.2*N*N+3*N ) THEN
1189 *
1190 *                    Sufficient workspace for a fast algorithm
1191 *
1192                      IU = 1
1193                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
1194 *
1195 *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
1196 *
1197                         LDWRKU = LDA
1198                         IR = IU + LDWRKU*N
1199                         LDWRKR = LDA
1200                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
1201 *
1202 *                       WORK(IU) is LDA by N and WORK(IR) is N by N
1203 *
1204                         LDWRKU = LDA
1205                         IR = IU + LDWRKU*N
1206                         LDWRKR = N
1207                      ELSE
1208 *
1209 *                       WORK(IU) is N by N and WORK(IR) is N by N
1210 *
1211                         LDWRKU = N
1212                         IR = IU + LDWRKU*N
1213                         LDWRKR = N
1214                      END IF
1215                      ITAU = IR + LDWRKR*N
1216                      IWORK = ITAU + N
1217 *
1218 *                    Compute A=Q*R
1219 *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1220 *                    (RWorkspace: 0)
1221 *
1222                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1223      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1224 *
1225 *                    Copy R to WORK(IU), zeroing out below it
1226 *
1227                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
1228      $                            LDWRKU )
1229                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1230      $                            WORK( IU+1 ), LDWRKU )
1231 *
1232 *                    Generate Q in A
1233 *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1234 *                    (RWorkspace: 0)
1235 *
1236                      CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
1237      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1238                      IE = 1
1239                      ITAUQ = ITAU
1240                      ITAUP = ITAUQ + N
1241                      IWORK = ITAUP + N
1242 *
1243 *                    Bidiagonalize R in WORK(IU), copying result to
1244 *                    WORK(IR)
1245 *                    (CWorkspace: need   2*N*N+3*N,
1246 *                                 prefer 2*N*N+2*N+2*N*NB)
1247 *                    (RWorkspace: need   N)
1248 *
1249                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
1250      $                            RWORK( IE ), WORK( ITAUQ ),
1251      $                            WORK( ITAUP ), WORK( IWORK ),
1252      $                            LWORK-IWORK+1, IERR )
1253                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU,
1254      $                            WORK( IR ), LDWRKR )
1255 *
1256 *                    Generate left bidiagonalizing vectors in WORK(IU)
1257 *                    (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1258 *                    (RWorkspace: 0)
1259 *
1260                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1261      $                            WORK( ITAUQ ), WORK( IWORK ),
1262      $                            LWORK-IWORK+1, IERR )
1263 *
1264 *                    Generate right bidiagonalizing vectors in WORK(IR)
1265 *                    (CWorkspace: need   2*N*N+3*N-1,
1266 *                                 prefer 2*N*N+2*N+(N-1)*NB)
1267 *                    (RWorkspace: 0)
1268 *
1269                      CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
1270      $                            WORK( ITAUP ), WORK( IWORK ),
1271      $                            LWORK-IWORK+1, IERR )
1272                      IRWORK = IE + N
1273 *
1274 *                    Perform bidiagonal QR iteration, computing left
1275 *                    singular vectors of R in WORK(IU) and computing
1276 *                    right singular vectors of R in WORK(IR)
1277 *                    (CWorkspace: need 2*N*N)
1278 *                    (RWorkspace: need BDSPAC)
1279 *
1280                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
1281      $                            WORK( IR ), LDWRKR, WORK( IU ),
1282      $                            LDWRKU, CDUM, 1, RWORK( IRWORK ),
1283      $                            INFO )
1284 *
1285 *                    Multiply Q in A by left singular vectors of R in
1286 *                    WORK(IU), storing result in U
1287 *                    (CWorkspace: need N*N)
1288 *                    (RWorkspace: 0)
1289 *
1290                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
1291      $                           WORK( IU ), LDWRKU, CZERO, U, LDU )
1292 *
1293 *                    Copy right singular vectors of R to A
1294 *                    (CWorkspace: need N*N)
1295 *                    (RWorkspace: 0)
1296 *
1297                      CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
1298      $                            LDA )
1299 *
1300                   ELSE
1301 *
1302 *                    Insufficient workspace for a fast algorithm
1303 *
1304                      ITAU = 1
1305                      IWORK = ITAU + N
1306 *
1307 *                    Compute A=Q*R, copying result to U
1308 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1309 *                    (RWorkspace: 0)
1310 *
1311                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1312      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1313                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1314 *
1315 *                    Generate Q in U
1316 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1317 *                    (RWorkspace: 0)
1318 *
1319                      CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
1320      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1321                      IE = 1
1322                      ITAUQ = ITAU
1323                      ITAUP = ITAUQ + N
1324                      IWORK = ITAUP + N
1325 *
1326 *                    Zero out below R in A
1327 *
1328                      IF( N .GT. 1 ) THEN
1329                         CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1330      $                               A( 2, 1 ), LDA )
1331                      END IF
1332 *
1333 *                    Bidiagonalize R in A
1334 *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1335 *                    (RWorkspace: need N)
1336 *
1337                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
1338      $                            WORK( ITAUQ ), WORK( ITAUP ),
1339      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1340 *
1341 *                    Multiply Q in U by left vectors bidiagonalizing R
1342 *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1343 *                    (RWorkspace: 0)
1344 *
1345                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1346      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1347      $                            LWORK-IWORK+1, IERR )
1348 *
1349 *                    Generate right vectors bidiagonalizing R in A
1350 *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1351 *                    (RWorkspace: 0)
1352 *
1353                      CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
1354      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1355                      IRWORK = IE + N
1356 *
1357 *                    Perform bidiagonal QR iteration, computing left
1358 *                    singular vectors of A in U and computing right
1359 *                    singular vectors of A in A
1360 *                    (CWorkspace: 0)
1361 *                    (RWorkspace: need BDSPAC)
1362 *
1363                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
1364      $                            LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
1365      $                            INFO )
1366 *
1367                   END IF
1368 *
1369                ELSE IF( WNTVAS ) THEN
1370 *
1371 *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1372 *                         or 'A')
1373 *                 N left singular vectors to be computed in U and
1374 *                 N right singular vectors to be computed in VT
1375 *
1376                   IF( LWORK.GE.N*N+3*N ) THEN
1377 *
1378 *                    Sufficient workspace for a fast algorithm
1379 *
1380                      IU = 1
1381                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
1382 *
1383 *                       WORK(IU) is LDA by N
1384 *
1385                         LDWRKU = LDA
1386                      ELSE
1387 *
1388 *                       WORK(IU) is N by N
1389 *
1390                         LDWRKU = N
1391                      END IF
1392                      ITAU = IU + LDWRKU*N
1393                      IWORK = ITAU + N
1394 *
1395 *                    Compute A=Q*R
1396 *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1397 *                    (RWorkspace: 0)
1398 *
1399                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1400      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1401 *
1402 *                    Copy R to WORK(IU), zeroing out below it
1403 *
1404                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
1405      $                            LDWRKU )
1406                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1407      $                            WORK( IU+1 ), LDWRKU )
1408 *
1409 *                    Generate Q in A
1410 *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1411 *                    (RWorkspace: 0)
1412 *
1413                      CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
1414      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1415                      IE = 1
1416                      ITAUQ = ITAU
1417                      ITAUP = ITAUQ + N
1418                      IWORK = ITAUP + N
1419 *
1420 *                    Bidiagonalize R in WORK(IU), copying result to VT
1421 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1422 *                    (RWorkspace: need N)
1423 *
1424                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
1425      $                            RWORK( IE ), WORK( ITAUQ ),
1426      $                            WORK( ITAUP ), WORK( IWORK ),
1427      $                            LWORK-IWORK+1, IERR )
1428                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
1429      $                            LDVT )
1430 *
1431 *                    Generate left bidiagonalizing vectors in WORK(IU)
1432 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1433 *                    (RWorkspace: 0)
1434 *
1435                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1436      $                            WORK( ITAUQ ), WORK( IWORK ),
1437      $                            LWORK-IWORK+1, IERR )
1438 *
1439 *                    Generate right bidiagonalizing vectors in VT
1440 *                    (CWorkspace: need   N*N+3*N-1,
1441 *                                 prefer N*N+2*N+(N-1)*NB)
1442 *                    (RWorkspace: 0)
1443 *
1444                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1445      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1446                      IRWORK = IE + N
1447 *
1448 *                    Perform bidiagonal QR iteration, computing left
1449 *                    singular vectors of R in WORK(IU) and computing
1450 *                    right singular vectors of R in VT
1451 *                    (CWorkspace: need N*N)
1452 *                    (RWorkspace: need BDSPAC)
1453 *
1454                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
1455      $                            LDVT, WORK( IU ), LDWRKU, CDUM, 1,
1456      $                            RWORK( IRWORK ), INFO )
1457 *
1458 *                    Multiply Q in A by left singular vectors of R in
1459 *                    WORK(IU), storing result in U
1460 *                    (CWorkspace: need N*N)
1461 *                    (RWorkspace: 0)
1462 *
1463                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
1464      $                           WORK( IU ), LDWRKU, CZERO, U, LDU )
1465 *
1466                   ELSE
1467 *
1468 *                    Insufficient workspace for a fast algorithm
1469 *
1470                      ITAU = 1
1471                      IWORK = ITAU + N
1472 *
1473 *                    Compute A=Q*R, copying result to U
1474 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1475 *                    (RWorkspace: 0)
1476 *
1477                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1478      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1479                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1480 *
1481 *                    Generate Q in U
1482 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1483 *                    (RWorkspace: 0)
1484 *
1485                      CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
1486      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1487 *
1488 *                    Copy R to VT, zeroing out below it
1489 *
1490                      CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
1491                      IF( N.GT.1 )
1492      $                  CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1493      $                               VT( 2, 1 ), LDVT )
1494                      IE = 1
1495                      ITAUQ = ITAU
1496                      ITAUP = ITAUQ + N
1497                      IWORK = ITAUP + N
1498 *
1499 *                    Bidiagonalize R in VT
1500 *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1501 *                    (RWorkspace: need N)
1502 *
1503                      CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
1504      $                            WORK( ITAUQ ), WORK( ITAUP ),
1505      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1506 *
1507 *                    Multiply Q in U by left bidiagonalizing vectors
1508 *                    in VT
1509 *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1510 *                    (RWorkspace: 0)
1511 *
1512                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
1513      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1514      $                            LWORK-IWORK+1, IERR )
1515 *
1516 *                    Generate right bidiagonalizing vectors in VT
1517 *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1518 *                    (RWorkspace: 0)
1519 *
1520                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1521      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1522                      IRWORK = IE + N
1523 *
1524 *                    Perform bidiagonal QR iteration, computing left
1525 *                    singular vectors of A in U and computing right
1526 *                    singular vectors of A in VT
1527 *                    (CWorkspace: 0)
1528 *                    (RWorkspace: need BDSPAC)
1529 *
1530                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
1531      $                            LDVT, U, LDU, CDUM, 1,
1532      $                            RWORK( IRWORK ), INFO )
1533 *
1534                   END IF
1535 *
1536                END IF
1537 *
1538             ELSE IF( WNTUA ) THEN
1539 *
1540                IF( WNTVN ) THEN
1541 *
1542 *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1543 *                 M left singular vectors to be computed in U and
1544 *                 no right singular vectors to be computed
1545 *
1546                   IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
1547 *
1548 *                    Sufficient workspace for a fast algorithm
1549 *
1550                      IR = 1
1551                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
1552 *
1553 *                       WORK(IR) is LDA by N
1554 *
1555                         LDWRKR = LDA
1556                      ELSE
1557 *
1558 *                       WORK(IR) is N by N
1559 *
1560                         LDWRKR = N
1561                      END IF
1562                      ITAU = IR + LDWRKR*N
1563                      IWORK = ITAU + N
1564 *
1565 *                    Compute A=Q*R, copying result to U
1566 *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1567 *                    (RWorkspace: 0)
1568 *
1569                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1570      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1571                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1572 *
1573 *                    Copy R to WORK(IR), zeroing out below it
1574 *
1575                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ),
1576      $                            LDWRKR )
1577                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1578      $                            WORK( IR+1 ), LDWRKR )
1579 *
1580 *                    Generate Q in U
1581 *                    (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1582 *                    (RWorkspace: 0)
1583 *
1584                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1585      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1586                      IE = 1
1587                      ITAUQ = ITAU
1588                      ITAUP = ITAUQ + N
1589                      IWORK = ITAUP + N
1590 *
1591 *                    Bidiagonalize R in WORK(IR)
1592 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1593 *                    (RWorkspace: need N)
1594 *
1595                      CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S,
1596      $                            RWORK( IE ), WORK( ITAUQ ),
1597      $                            WORK( ITAUP ), WORK( IWORK ),
1598      $                            LWORK-IWORK+1, IERR )
1599 *
1600 *                    Generate left bidiagonalizing vectors in WORK(IR)
1601 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1602 *                    (RWorkspace: 0)
1603 *
1604                      CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
1605      $                            WORK( ITAUQ ), WORK( IWORK ),
1606      $                            LWORK-IWORK+1, IERR )
1607                      IRWORK = IE + N
1608 *
1609 *                    Perform bidiagonal QR iteration, computing left
1610 *                    singular vectors of R in WORK(IR)
1611 *                    (CWorkspace: need N*N)
1612 *                    (RWorkspace: need BDSPAC)
1613 *
1614                      CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
1615      $                            1, WORK( IR ), LDWRKR, CDUM, 1,
1616      $                            RWORK( IRWORK ), INFO )
1617 *
1618 *                    Multiply Q in U by left singular vectors of R in
1619 *                    WORK(IR), storing result in A
1620 *                    (CWorkspace: need N*N)
1621 *                    (RWorkspace: 0)
1622 *
1623                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
1624      $                           WORK( IR ), LDWRKR, CZERO, A, LDA )
1625 *
1626 *                    Copy left singular vectors of A from A to U
1627 *
1628                      CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
1629 *
1630                   ELSE
1631 *
1632 *                    Insufficient workspace for a fast algorithm
1633 *
1634                      ITAU = 1
1635                      IWORK = ITAU + N
1636 *
1637 *                    Compute A=Q*R, copying result to U
1638 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1639 *                    (RWorkspace: 0)
1640 *
1641                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1642      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1643                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1644 *
1645 *                    Generate Q in U
1646 *                    (CWorkspace: need N+M, prefer N+M*NB)
1647 *                    (RWorkspace: 0)
1648 *
1649                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1650      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1651                      IE = 1
1652                      ITAUQ = ITAU
1653                      ITAUP = ITAUQ + N
1654                      IWORK = ITAUP + N
1655 *
1656 *                    Zero out below R in A
1657 *
1658                      IF( N .GT. 1 ) THEN
1659                         CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1660      $                               A( 2, 1 ), LDA )
1661                      END IF
1662 *
1663 *                    Bidiagonalize R in A
1664 *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1665 *                    (RWorkspace: need N)
1666 *
1667                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
1668      $                            WORK( ITAUQ ), WORK( ITAUP ),
1669      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1670 *
1671 *                    Multiply Q in U by left bidiagonalizing vectors
1672 *                    in A
1673 *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1674 *                    (RWorkspace: 0)
1675 *
1676                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1677      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1678      $                            LWORK-IWORK+1, IERR )
1679                      IRWORK = IE + N
1680 *
1681 *                    Perform bidiagonal QR iteration, computing left
1682 *                    singular vectors of A in U
1683 *                    (CWorkspace: 0)
1684 *                    (RWorkspace: need BDSPAC)
1685 *
1686                      CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
1687      $                            1, U, LDU, CDUM, 1, RWORK( IRWORK ),
1688      $                            INFO )
1689 *
1690                   END IF
1691 *
1692                ELSE IF( WNTVO ) THEN
1693 *
1694 *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
1695 *                 M left singular vectors to be computed in U and
1696 *                 N right singular vectors to be overwritten on A
1697 *
1698                   IF( LWORK.GE.2*N*N+MAX( N+M, 3*N ) ) THEN
1699 *
1700 *                    Sufficient workspace for a fast algorithm
1701 *
1702                      IU = 1
1703                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
1704 *
1705 *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
1706 *
1707                         LDWRKU = LDA
1708                         IR = IU + LDWRKU*N
1709                         LDWRKR = LDA
1710                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
1711 *
1712 *                       WORK(IU) is LDA by N and WORK(IR) is N by N
1713 *
1714                         LDWRKU = LDA
1715                         IR = IU + LDWRKU*N
1716                         LDWRKR = N
1717                      ELSE
1718 *
1719 *                       WORK(IU) is N by N and WORK(IR) is N by N
1720 *
1721                         LDWRKU = N
1722                         IR = IU + LDWRKU*N
1723                         LDWRKR = N
1724                      END IF
1725                      ITAU = IR + LDWRKR*N
1726                      IWORK = ITAU + N
1727 *
1728 *                    Compute A=Q*R, copying result to U
1729 *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1730 *                    (RWorkspace: 0)
1731 *
1732                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1733      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1734                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1735 *
1736 *                    Generate Q in U
1737 *                    (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
1738 *                    (RWorkspace: 0)
1739 *
1740                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1741      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1742 *
1743 *                    Copy R to WORK(IU), zeroing out below it
1744 *
1745                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
1746      $                            LDWRKU )
1747                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1748      $                            WORK( IU+1 ), LDWRKU )
1749                      IE = 1
1750                      ITAUQ = ITAU
1751                      ITAUP = ITAUQ + N
1752                      IWORK = ITAUP + N
1753 *
1754 *                    Bidiagonalize R in WORK(IU), copying result to
1755 *                    WORK(IR)
1756 *                    (CWorkspace: need   2*N*N+3*N,
1757 *                                 prefer 2*N*N+2*N+2*N*NB)
1758 *                    (RWorkspace: need   N)
1759 *
1760                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
1761      $                            RWORK( IE ), WORK( ITAUQ ),
1762      $                            WORK( ITAUP ), WORK( IWORK ),
1763      $                            LWORK-IWORK+1, IERR )
1764                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU,
1765      $                            WORK( IR ), LDWRKR )
1766 *
1767 *                    Generate left bidiagonalizing vectors in WORK(IU)
1768 *                    (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1769 *                    (RWorkspace: 0)
1770 *
1771                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1772      $                            WORK( ITAUQ ), WORK( IWORK ),
1773      $                            LWORK-IWORK+1, IERR )
1774 *
1775 *                    Generate right bidiagonalizing vectors in WORK(IR)
1776 *                    (CWorkspace: need   2*N*N+3*N-1,
1777 *                                 prefer 2*N*N+2*N+(N-1)*NB)
1778 *                    (RWorkspace: 0)
1779 *
1780                      CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
1781      $                            WORK( ITAUP ), WORK( IWORK ),
1782      $                            LWORK-IWORK+1, IERR )
1783                      IRWORK = IE + N
1784 *
1785 *                    Perform bidiagonal QR iteration, computing left
1786 *                    singular vectors of R in WORK(IU) and computing
1787 *                    right singular vectors of R in WORK(IR)
1788 *                    (CWorkspace: need 2*N*N)
1789 *                    (RWorkspace: need BDSPAC)
1790 *
1791                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
1792      $                            WORK( IR ), LDWRKR, WORK( IU ),
1793      $                            LDWRKU, CDUM, 1, RWORK( IRWORK ),
1794      $                            INFO )
1795 *
1796 *                    Multiply Q in U by left singular vectors of R in
1797 *                    WORK(IU), storing result in A
1798 *                    (CWorkspace: need N*N)
1799 *                    (RWorkspace: 0)
1800 *
1801                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
1802      $                           WORK( IU ), LDWRKU, CZERO, A, LDA )
1803 *
1804 *                    Copy left singular vectors of A from A to U
1805 *
1806                      CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
1807 *
1808 *                    Copy right singular vectors of R from WORK(IR) to A
1809 *
1810                      CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
1811      $                            LDA )
1812 *
1813                   ELSE
1814 *
1815 *                    Insufficient workspace for a fast algorithm
1816 *
1817                      ITAU = 1
1818                      IWORK = ITAU + N
1819 *
1820 *                    Compute A=Q*R, copying result to U
1821 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1822 *                    (RWorkspace: 0)
1823 *
1824                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1825      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1826                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1827 *
1828 *                    Generate Q in U
1829 *                    (CWorkspace: need N+M, prefer N+M*NB)
1830 *                    (RWorkspace: 0)
1831 *
1832                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1833      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1834                      IE = 1
1835                      ITAUQ = ITAU
1836                      ITAUP = ITAUQ + N
1837                      IWORK = ITAUP + N
1838 *
1839 *                    Zero out below R in A
1840 *
1841                      IF( N .GT. 1 ) THEN
1842                         CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1843      $                               A( 2, 1 ), LDA )
1844                      END IF
1845 *
1846 *                    Bidiagonalize R in A
1847 *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1848 *                    (RWorkspace: need N)
1849 *
1850                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
1851      $                            WORK( ITAUQ ), WORK( ITAUP ),
1852      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1853 *
1854 *                    Multiply Q in U by left bidiagonalizing vectors
1855 *                    in A
1856 *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1857 *                    (RWorkspace: 0)
1858 *
1859                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1860      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1861      $                            LWORK-IWORK+1, IERR )
1862 *
1863 *                    Generate right bidiagonalizing vectors in A
1864 *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1865 *                    (RWorkspace: 0)
1866 *
1867                      CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
1868      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1869                      IRWORK = IE + N
1870 *
1871 *                    Perform bidiagonal QR iteration, computing left
1872 *                    singular vectors of A in U and computing right
1873 *                    singular vectors of A in A
1874 *                    (CWorkspace: 0)
1875 *                    (RWorkspace: need BDSPAC)
1876 *
1877                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
1878      $                            LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
1879      $                            INFO )
1880 *
1881                   END IF
1882 *
1883                ELSE IF( WNTVAS ) THEN
1884 *
1885 *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
1886 *                         or 'A')
1887 *                 M left singular vectors to be computed in U and
1888 *                 N right singular vectors to be computed in VT
1889 *
1890                   IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
1891 *
1892 *                    Sufficient workspace for a fast algorithm
1893 *
1894                      IU = 1
1895                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
1896 *
1897 *                       WORK(IU) is LDA by N
1898 *
1899                         LDWRKU = LDA
1900                      ELSE
1901 *
1902 *                       WORK(IU) is N by N
1903 *
1904                         LDWRKU = N
1905                      END IF
1906                      ITAU = IU + LDWRKU*N
1907                      IWORK = ITAU + N
1908 *
1909 *                    Compute A=Q*R, copying result to U
1910 *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1911 *                    (RWorkspace: 0)
1912 *
1913                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1914      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1915                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1916 *
1917 *                    Generate Q in U
1918 *                    (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1919 *                    (RWorkspace: 0)
1920 *
1921                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1922      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1923 *
1924 *                    Copy R to WORK(IU), zeroing out below it
1925 *
1926                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
1927      $                            LDWRKU )
1928                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
1929      $                            WORK( IU+1 ), LDWRKU )
1930                      IE = 1
1931                      ITAUQ = ITAU
1932                      ITAUP = ITAUQ + N
1933                      IWORK = ITAUP + N
1934 *
1935 *                    Bidiagonalize R in WORK(IU), copying result to VT
1936 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1937 *                    (RWorkspace: need N)
1938 *
1939                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
1940      $                            RWORK( IE ), WORK( ITAUQ ),
1941      $                            WORK( ITAUP ), WORK( IWORK ),
1942      $                            LWORK-IWORK+1, IERR )
1943                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
1944      $                            LDVT )
1945 *
1946 *                    Generate left bidiagonalizing vectors in WORK(IU)
1947 *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1948 *                    (RWorkspace: 0)
1949 *
1950                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1951      $                            WORK( ITAUQ ), WORK( IWORK ),
1952      $                            LWORK-IWORK+1, IERR )
1953 *
1954 *                    Generate right bidiagonalizing vectors in VT
1955 *                    (CWorkspace: need   N*N+3*N-1,
1956 *                                 prefer N*N+2*N+(N-1)*NB)
1957 *                    (RWorkspace: need   0)
1958 *
1959                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1960      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1961                      IRWORK = IE + N
1962 *
1963 *                    Perform bidiagonal QR iteration, computing left
1964 *                    singular vectors of R in WORK(IU) and computing
1965 *                    right singular vectors of R in VT
1966 *                    (CWorkspace: need N*N)
1967 *                    (RWorkspace: need BDSPAC)
1968 *
1969                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
1970      $                            LDVT, WORK( IU ), LDWRKU, CDUM, 1,
1971      $                            RWORK( IRWORK ), INFO )
1972 *
1973 *                    Multiply Q in U by left singular vectors of R in
1974 *                    WORK(IU), storing result in A
1975 *                    (CWorkspace: need N*N)
1976 *                    (RWorkspace: 0)
1977 *
1978                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
1979      $                           WORK( IU ), LDWRKU, CZERO, A, LDA )
1980 *
1981 *                    Copy left singular vectors of A from A to U
1982 *
1983                      CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
1984 *
1985                   ELSE
1986 *
1987 *                    Insufficient workspace for a fast algorithm
1988 *
1989                      ITAU = 1
1990                      IWORK = ITAU + N
1991 *
1992 *                    Compute A=Q*R, copying result to U
1993 *                    (CWorkspace: need 2*N, prefer N+N*NB)
1994 *                    (RWorkspace: 0)
1995 *
1996                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
1997      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1998                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1999 *
2000 *                    Generate Q in U
2001 *                    (CWorkspace: need N+M, prefer N+M*NB)
2002 *                    (RWorkspace: 0)
2003 *
2004                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
2005      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2006 *
2007 *                    Copy R from A to VT, zeroing out below it
2008 *
2009                      CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
2010                      IF( N.GT.1 )
2011      $                  CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
2012      $                               VT( 2, 1 ), LDVT )
2013                      IE = 1
2014                      ITAUQ = ITAU
2015                      ITAUP = ITAUQ + N
2016                      IWORK = ITAUP + N
2017 *
2018 *                    Bidiagonalize R in VT
2019 *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
2020 *                    (RWorkspace: need N)
2021 *
2022                      CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
2023      $                            WORK( ITAUQ ), WORK( ITAUP ),
2024      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2025 *
2026 *                    Multiply Q in U by left bidiagonalizing vectors
2027 *                    in VT
2028 *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
2029 *                    (RWorkspace: 0)
2030 *
2031                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
2032      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
2033      $                            LWORK-IWORK+1, IERR )
2034 *
2035 *                    Generate right bidiagonalizing vectors in VT
2036 *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2037 *                    (RWorkspace: 0)
2038 *
2039                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
2040      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2041                      IRWORK = IE + N
2042 *
2043 *                    Perform bidiagonal QR iteration, computing left
2044 *                    singular vectors of A in U and computing right
2045 *                    singular vectors of A in VT
2046 *                    (CWorkspace: 0)
2047 *                    (RWorkspace: need BDSPAC)
2048 *
2049                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
2050      $                            LDVT, U, LDU, CDUM, 1,
2051      $                            RWORK( IRWORK ), INFO )
2052 *
2053                   END IF
2054 *
2055                END IF
2056 *
2057             END IF
2058 *
2059          ELSE
2060 *
2061 *           M .LT. MNTHR
2062 *
2063 *           Path 10 (M at least N, but not much larger)
2064 *           Reduce to bidiagonal form without QR decomposition
2065 *
2066             IE = 1
2067             ITAUQ = 1
2068             ITAUP = ITAUQ + N
2069             IWORK = ITAUP + N
2070 *
2071 *           Bidiagonalize A
2072 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
2073 *           (RWorkspace: need N)
2074 *
2075             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
2076      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
2077      $                   IERR )
2078             IF( WNTUAS ) THEN
2079 *
2080 *              If left singular vectors desired in U, copy result to U
2081 *              and generate left bidiagonalizing vectors in U
2082 *              (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
2083 *              (RWorkspace: 0)
2084 *
2085                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
2086                IF( WNTUS )
2087      $            NCU = N
2088                IF( WNTUA )
2089      $            NCU = M
2090                CALL ZUNGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ),
2091      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
2092             END IF
2093             IF( WNTVAS ) THEN
2094 *
2095 *              If right singular vectors desired in VT, copy result to
2096 *              VT and generate right bidiagonalizing vectors in VT
2097 *              (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2098 *              (RWorkspace: 0)
2099 *
2100                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
2101                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
2102      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
2103             END IF
2104             IF( WNTUO ) THEN
2105 *
2106 *              If left singular vectors desired in A, generate left
2107 *              bidiagonalizing vectors in A
2108 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
2109 *              (RWorkspace: 0)
2110 *
2111                CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
2112      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
2113             END IF
2114             IF( WNTVO ) THEN
2115 *
2116 *              If right singular vectors desired in A, generate right
2117 *              bidiagonalizing vectors in A
2118 *              (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2119 *              (RWorkspace: 0)
2120 *
2121                CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
2122      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
2123             END IF
2124             IRWORK = IE + N
2125             IF( WNTUAS .OR. WNTUO )
2126      $         NRU = M
2127             IF( WNTUN )
2128      $         NRU = 0
2129             IF( WNTVAS .OR. WNTVO )
2130      $         NCVT = N
2131             IF( WNTVN )
2132      $         NCVT = 0
2133             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
2134 *
2135 *              Perform bidiagonal QR iteration, if desired, computing
2136 *              left singular vectors in U and computing right singular
2137 *              vectors in VT
2138 *              (CWorkspace: 0)
2139 *              (RWorkspace: need BDSPAC)
2140 *
2141                CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
2142      $                      LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
2143      $                      INFO )
2144             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
2145 *
2146 *              Perform bidiagonal QR iteration, if desired, computing
2147 *              left singular vectors in U and computing right singular
2148 *              vectors in A
2149 *              (CWorkspace: 0)
2150 *              (RWorkspace: need BDSPAC)
2151 *
2152                CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), A,
2153      $                      LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
2154      $                      INFO )
2155             ELSE
2156 *
2157 *              Perform bidiagonal QR iteration, if desired, computing
2158 *              left singular vectors in A and computing right singular
2159 *              vectors in VT
2160 *              (CWorkspace: 0)
2161 *              (RWorkspace: need BDSPAC)
2162 *
2163                CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
2164      $                      LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
2165      $                      INFO )
2166             END IF
2167 *
2168          END IF
2169 *
2170       ELSE
2171 *
2172 *        A has more columns than rows. If A has sufficiently more
2173 *        columns than rows, first reduce using the LQ decomposition (if
2174 *        sufficient workspace available)
2175 *
2176          IF( N.GE.MNTHR ) THEN
2177 *
2178             IF( WNTVN ) THEN
2179 *
2180 *              Path 1t(N much larger than M, JOBVT='N')
2181 *              No right singular vectors to be computed
2182 *
2183                ITAU = 1
2184                IWORK = ITAU + M
2185 *
2186 *              Compute A=L*Q
2187 *              (CWorkspace: need 2*M, prefer M+M*NB)
2188 *              (RWorkspace: 0)
2189 *
2190                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
2191      $                      LWORK-IWORK+1, IERR )
2192 *
2193 *              Zero out above L
2194 *
2195                CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
2196      $                      LDA )
2197                IE = 1
2198                ITAUQ = 1
2199                ITAUP = ITAUQ + M
2200                IWORK = ITAUP + M
2201 *
2202 *              Bidiagonalize L in A
2203 *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2204 *              (RWorkspace: need M)
2205 *
2206                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
2207      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
2208      $                      IERR )
2209                IF( WNTUO .OR. WNTUAS ) THEN
2210 *
2211 *                 If left singular vectors desired, generate Q
2212 *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2213 *                 (RWorkspace: 0)
2214 *
2215                   CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
2216      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2217                END IF
2218                IRWORK = IE + M
2219                NRU = 0
2220                IF( WNTUO .OR. WNTUAS )
2221      $            NRU = M
2222 *
2223 *              Perform bidiagonal QR iteration, computing left singular
2224 *              vectors of A in A if desired
2225 *              (CWorkspace: 0)
2226 *              (RWorkspace: need BDSPAC)
2227 *
2228                CALL ZBDSQR( 'U', M, 0, NRU, 0, S, RWORK( IE ), CDUM, 1,
2229      $                      A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
2230 *
2231 *              If left singular vectors desired in U, copy them there
2232 *
2233                IF( WNTUAS )
2234      $            CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )
2235 *
2236             ELSE IF( WNTVO .AND. WNTUN ) THEN
2237 *
2238 *              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2239 *              M right singular vectors to be overwritten on A and
2240 *              no left singular vectors to be computed
2241 *
2242                IF( LWORK.GE.M*M+3*M ) THEN
2243 *
2244 *                 Sufficient workspace for a fast algorithm
2245 *
2246                   IR = 1
2247                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
2248 *
2249 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
2250 *
2251                      LDWRKU = LDA
2252                      CHUNK = N
2253                      LDWRKR = LDA
2254                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
2255 *
2256 *                    WORK(IU) is LDA by N and WORK(IR) is M by M
2257 *
2258                      LDWRKU = LDA
2259                      CHUNK = N
2260                      LDWRKR = M
2261                   ELSE
2262 *
2263 *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
2264 *
2265                      LDWRKU = M
2266                      CHUNK = ( LWORK-M*M ) / M
2267                      LDWRKR = M
2268                   END IF
2269                   ITAU = IR + LDWRKR*M
2270                   IWORK = ITAU + M
2271 *
2272 *                 Compute A=L*Q
2273 *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2274 *                 (RWorkspace: 0)
2275 *
2276                   CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2277      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2278 *
2279 *                 Copy L to WORK(IR) and zero out above it
2280 *
2281                   CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR )
2282                   CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
2283      $                         WORK( IR+LDWRKR ), LDWRKR )
2284 *
2285 *                 Generate Q in A
2286 *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2287 *                 (RWorkspace: 0)
2288 *
2289                   CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2290      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2291                   IE = 1
2292                   ITAUQ = ITAU
2293                   ITAUP = ITAUQ + M
2294                   IWORK = ITAUP + M
2295 *
2296 *                 Bidiagonalize L in WORK(IR)
2297 *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2298 *                 (RWorkspace: need M)
2299 *
2300                   CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, RWORK( IE ),
2301      $                         WORK( ITAUQ ), WORK( ITAUP ),
2302      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2303 *
2304 *                 Generate right vectors bidiagonalizing L
2305 *                 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2306 *                 (RWorkspace: 0)
2307 *
2308                   CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
2309      $                         WORK( ITAUP ), WORK( IWORK ),
2310      $                         LWORK-IWORK+1, IERR )
2311                   IRWORK = IE + M
2312 *
2313 *                 Perform bidiagonal QR iteration, computing right
2314 *                 singular vectors of L in WORK(IR)
2315 *                 (CWorkspace: need M*M)
2316 *                 (RWorkspace: need BDSPAC)
2317 *
2318                   CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
2319      $                         WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
2320      $                         RWORK( IRWORK ), INFO )
2321                   IU = ITAUQ
2322 *
2323 *                 Multiply right singular vectors of L in WORK(IR) by Q
2324 *                 in A, storing result in WORK(IU) and copying to A
2325 *                 (CWorkspace: need M*M+M, prefer M*M+M*N)
2326 *                 (RWorkspace: 0)
2327 *
2328                   DO 30 I = 1, N, CHUNK
2329                      BLK = MIN( N-I+1, CHUNK )
2330                      CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
2331      $                           LDWRKR, A( 1, I ), LDA, CZERO,
2332      $                           WORK( IU ), LDWRKU )
2333                      CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
2334      $                            A( 1, I ), LDA )
2335    30             CONTINUE
2336 *
2337                ELSE
2338 *
2339 *                 Insufficient workspace for a fast algorithm
2340 *
2341                   IE = 1
2342                   ITAUQ = 1
2343                   ITAUP = ITAUQ + M
2344                   IWORK = ITAUP + M
2345 *
2346 *                 Bidiagonalize A
2347 *                 (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
2348 *                 (RWorkspace: need M)
2349 *
2350                   CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ),
2351      $                         WORK( ITAUQ ), WORK( ITAUP ),
2352      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2353 *
2354 *                 Generate right vectors bidiagonalizing A
2355 *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2356 *                 (RWorkspace: 0)
2357 *
2358                   CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
2359      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2360                   IRWORK = IE + M
2361 *
2362 *                 Perform bidiagonal QR iteration, computing right
2363 *                 singular vectors of A in A
2364 *                 (CWorkspace: 0)
2365 *                 (RWorkspace: need BDSPAC)
2366 *
2367                   CALL ZBDSQR( 'L', M, N, 0, 0, S, RWORK( IE ), A, LDA,
2368      $                         CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
2369 *
2370                END IF
2371 *
2372             ELSE IF( WNTVO .AND. WNTUAS ) THEN
2373 *
2374 *              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2375 *              M right singular vectors to be overwritten on A and
2376 *              M left singular vectors to be computed in U
2377 *
2378                IF( LWORK.GE.M*M+3*M ) THEN
2379 *
2380 *                 Sufficient workspace for a fast algorithm
2381 *
2382                   IR = 1
2383                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
2384 *
2385 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
2386 *
2387                      LDWRKU = LDA
2388                      CHUNK = N
2389                      LDWRKR = LDA
2390                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
2391 *
2392 *                    WORK(IU) is LDA by N and WORK(IR) is M by M
2393 *
2394                      LDWRKU = LDA
2395                      CHUNK = N
2396                      LDWRKR = M
2397                   ELSE
2398 *
2399 *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
2400 *
2401                      LDWRKU = M
2402                      CHUNK = ( LWORK-M*M ) / M
2403                      LDWRKR = M
2404                   END IF
2405                   ITAU = IR + LDWRKR*M
2406                   IWORK = ITAU + M
2407 *
2408 *                 Compute A=L*Q
2409 *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2410 *                 (RWorkspace: 0)
2411 *
2412                   CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2413      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2414 *
2415 *                 Copy L to U, zeroing about above it
2416 *
2417                   CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
2418                   CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
2419      $                         LDU )
2420 *
2421 *                 Generate Q in A
2422 *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2423 *                 (RWorkspace: 0)
2424 *
2425                   CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2426      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2427                   IE = 1
2428                   ITAUQ = ITAU
2429                   ITAUP = ITAUQ + M
2430                   IWORK = ITAUP + M
2431 *
2432 *                 Bidiagonalize L in U, copying result to WORK(IR)
2433 *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2434 *                 (RWorkspace: need M)
2435 *
2436                   CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
2437      $                         WORK( ITAUQ ), WORK( ITAUP ),
2438      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2439                   CALL ZLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
2440 *
2441 *                 Generate right vectors bidiagonalizing L in WORK(IR)
2442 *                 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2443 *                 (RWorkspace: 0)
2444 *
2445                   CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
2446      $                         WORK( ITAUP ), WORK( IWORK ),
2447      $                         LWORK-IWORK+1, IERR )
2448 *
2449 *                 Generate left vectors bidiagonalizing L in U
2450 *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2451 *                 (RWorkspace: 0)
2452 *
2453                   CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
2454      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2455                   IRWORK = IE + M
2456 *
2457 *                 Perform bidiagonal QR iteration, computing left
2458 *                 singular vectors of L in U, and computing right
2459 *                 singular vectors of L in WORK(IR)
2460 *                 (CWorkspace: need M*M)
2461 *                 (RWorkspace: need BDSPAC)
2462 *
2463                   CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
2464      $                         WORK( IR ), LDWRKR, U, LDU, CDUM, 1,
2465      $                         RWORK( IRWORK ), INFO )
2466                   IU = ITAUQ
2467 *
2468 *                 Multiply right singular vectors of L in WORK(IR) by Q
2469 *                 in A, storing result in WORK(IU) and copying to A
2470 *                 (CWorkspace: need M*M+M, prefer M*M+M*N))
2471 *                 (RWorkspace: 0)
2472 *
2473                   DO 40 I = 1, N, CHUNK
2474                      BLK = MIN( N-I+1, CHUNK )
2475                      CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
2476      $                           LDWRKR, A( 1, I ), LDA, CZERO,
2477      $                           WORK( IU ), LDWRKU )
2478                      CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
2479      $                            A( 1, I ), LDA )
2480    40             CONTINUE
2481 *
2482                ELSE
2483 *
2484 *                 Insufficient workspace for a fast algorithm
2485 *
2486                   ITAU = 1
2487                   IWORK = ITAU + M
2488 *
2489 *                 Compute A=L*Q
2490 *                 (CWorkspace: need 2*M, prefer M+M*NB)
2491 *                 (RWorkspace: 0)
2492 *
2493                   CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2494      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2495 *
2496 *                 Copy L to U, zeroing out above it
2497 *
2498                   CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
2499                   CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
2500      $                         LDU )
2501 *
2502 *                 Generate Q in A
2503 *                 (CWorkspace: need 2*M, prefer M+M*NB)
2504 *                 (RWorkspace: 0)
2505 *
2506                   CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2507      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2508                   IE = 1
2509                   ITAUQ = ITAU
2510                   ITAUP = ITAUQ + M
2511                   IWORK = ITAUP + M
2512 *
2513 *                 Bidiagonalize L in U
2514 *                 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2515 *                 (RWorkspace: need M)
2516 *
2517                   CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
2518      $                         WORK( ITAUQ ), WORK( ITAUP ),
2519      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2520 *
2521 *                 Multiply right vectors bidiagonalizing L by Q in A
2522 *                 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2523 *                 (RWorkspace: 0)
2524 *
2525                   CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
2526      $                         WORK( ITAUP ), A, LDA, WORK( IWORK ),
2527      $                         LWORK-IWORK+1, IERR )
2528 *
2529 *                 Generate left vectors bidiagonalizing L in U
2530 *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2531 *                 (RWorkspace: 0)
2532 *
2533                   CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
2534      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2535                   IRWORK = IE + M
2536 *
2537 *                 Perform bidiagonal QR iteration, computing left
2538 *                 singular vectors of A in U and computing right
2539 *                 singular vectors of A in A
2540 *                 (CWorkspace: 0)
2541 *                 (RWorkspace: need BDSPAC)
2542 *
2543                   CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), A, LDA,
2544      $                         U, LDU, CDUM, 1, RWORK( IRWORK ), INFO )
2545 *
2546                END IF
2547 *
2548             ELSE IF( WNTVS ) THEN
2549 *
2550                IF( WNTUN ) THEN
2551 *
2552 *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2553 *                 M right singular vectors to be computed in VT and
2554 *                 no left singular vectors to be computed
2555 *
2556                   IF( LWORK.GE.M*M+3*M ) THEN
2557 *
2558 *                    Sufficient workspace for a fast algorithm
2559 *
2560                      IR = 1
2561                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
2562 *
2563 *                       WORK(IR) is LDA by M
2564 *
2565                         LDWRKR = LDA
2566                      ELSE
2567 *
2568 *                       WORK(IR) is M by M
2569 *
2570                         LDWRKR = M
2571                      END IF
2572                      ITAU = IR + LDWRKR*M
2573                      IWORK = ITAU + M
2574 *
2575 *                    Compute A=L*Q
2576 *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2577 *                    (RWorkspace: 0)
2578 *
2579                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2580      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2581 *
2582 *                    Copy L to WORK(IR), zeroing out above it
2583 *
2584                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ),
2585      $                            LDWRKR )
2586                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
2587      $                            WORK( IR+LDWRKR ), LDWRKR )
2588 *
2589 *                    Generate Q in A
2590 *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2591 *                    (RWorkspace: 0)
2592 *
2593                      CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2594      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2595                      IE = 1
2596                      ITAUQ = ITAU
2597                      ITAUP = ITAUQ + M
2598                      IWORK = ITAUP + M
2599 *
2600 *                    Bidiagonalize L in WORK(IR)
2601 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2602 *                    (RWorkspace: need M)
2603 *
2604                      CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S,
2605      $                            RWORK( IE ), WORK( ITAUQ ),
2606      $                            WORK( ITAUP ), WORK( IWORK ),
2607      $                            LWORK-IWORK+1, IERR )
2608 *
2609 *                    Generate right vectors bidiagonalizing L in
2610 *                    WORK(IR)
2611 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
2612 *                    (RWorkspace: 0)
2613 *
2614                      CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
2615      $                            WORK( ITAUP ), WORK( IWORK ),
2616      $                            LWORK-IWORK+1, IERR )
2617                      IRWORK = IE + M
2618 *
2619 *                    Perform bidiagonal QR iteration, computing right
2620 *                    singular vectors of L in WORK(IR)
2621 *                    (CWorkspace: need M*M)
2622 *                    (RWorkspace: need BDSPAC)
2623 *
2624                      CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
2625      $                            WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
2626      $                            RWORK( IRWORK ), INFO )
2627 *
2628 *                    Multiply right singular vectors of L in WORK(IR) by
2629 *                    Q in A, storing result in VT
2630 *                    (CWorkspace: need M*M)
2631 *                    (RWorkspace: 0)
2632 *
2633                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
2634      $                           LDWRKR, A, LDA, CZERO, VT, LDVT )
2635 *
2636                   ELSE
2637 *
2638 *                    Insufficient workspace for a fast algorithm
2639 *
2640                      ITAU = 1
2641                      IWORK = ITAU + M
2642 *
2643 *                    Compute A=L*Q
2644 *                    (CWorkspace: need 2*M, prefer M+M*NB)
2645 *                    (RWorkspace: 0)
2646 *
2647                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2648      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2649 *
2650 *                    Copy result to VT
2651 *
2652                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
2653 *
2654 *                    Generate Q in VT
2655 *                    (CWorkspace: need 2*M, prefer M+M*NB)
2656 *                    (RWorkspace: 0)
2657 *
2658                      CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
2659      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2660                      IE = 1
2661                      ITAUQ = ITAU
2662                      ITAUP = ITAUQ + M
2663                      IWORK = ITAUP + M
2664 *
2665 *                    Zero out above L in A
2666 *
2667                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
2668      $                            A( 1, 2 ), LDA )
2669 *
2670 *                    Bidiagonalize L in A
2671 *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2672 *                    (RWorkspace: need M)
2673 *
2674                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
2675      $                            WORK( ITAUQ ), WORK( ITAUP ),
2676      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2677 *
2678 *                    Multiply right vectors bidiagonalizing L by Q in VT
2679 *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2680 *                    (RWorkspace: 0)
2681 *
2682                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
2683      $                            WORK( ITAUP ), VT, LDVT,
2684      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2685                      IRWORK = IE + M
2686 *
2687 *                    Perform bidiagonal QR iteration, computing right
2688 *                    singular vectors of A in VT
2689 *                    (CWorkspace: 0)
2690 *                    (RWorkspace: need BDSPAC)
2691 *
2692                      CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
2693      $                            LDVT, CDUM, 1, CDUM, 1,
2694      $                            RWORK( IRWORK ), INFO )
2695 *
2696                   END IF
2697 *
2698                ELSE IF( WNTUO ) THEN
2699 *
2700 *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
2701 *                 M right singular vectors to be computed in VT and
2702 *                 M left singular vectors to be overwritten on A
2703 *
2704                   IF( LWORK.GE.2*M*M+3*M ) THEN
2705 *
2706 *                    Sufficient workspace for a fast algorithm
2707 *
2708                      IU = 1
2709                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
2710 *
2711 *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
2712 *
2713                         LDWRKU = LDA
2714                         IR = IU + LDWRKU*M
2715                         LDWRKR = LDA
2716                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
2717 *
2718 *                       WORK(IU) is LDA by M and WORK(IR) is M by M
2719 *
2720                         LDWRKU = LDA
2721                         IR = IU + LDWRKU*M
2722                         LDWRKR = M
2723                      ELSE
2724 *
2725 *                       WORK(IU) is M by M and WORK(IR) is M by M
2726 *
2727                         LDWRKU = M
2728                         IR = IU + LDWRKU*M
2729                         LDWRKR = M
2730                      END IF
2731                      ITAU = IR + LDWRKR*M
2732                      IWORK = ITAU + M
2733 *
2734 *                    Compute A=L*Q
2735 *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2736 *                    (RWorkspace: 0)
2737 *
2738                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2739      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2740 *
2741 *                    Copy L to WORK(IU), zeroing out below it
2742 *
2743                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
2744      $                            LDWRKU )
2745                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
2746      $                            WORK( IU+LDWRKU ), LDWRKU )
2747 *
2748 *                    Generate Q in A
2749 *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2750 *                    (RWorkspace: 0)
2751 *
2752                      CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2753      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2754                      IE = 1
2755                      ITAUQ = ITAU
2756                      ITAUP = ITAUQ + M
2757                      IWORK = ITAUP + M
2758 *
2759 *                    Bidiagonalize L in WORK(IU), copying result to
2760 *                    WORK(IR)
2761 *                    (CWorkspace: need   2*M*M+3*M,
2762 *                                 prefer 2*M*M+2*M+2*M*NB)
2763 *                    (RWorkspace: need   M)
2764 *
2765                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
2766      $                            RWORK( IE ), WORK( ITAUQ ),
2767      $                            WORK( ITAUP ), WORK( IWORK ),
2768      $                            LWORK-IWORK+1, IERR )
2769                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU,
2770      $                            WORK( IR ), LDWRKR )
2771 *
2772 *                    Generate right bidiagonalizing vectors in WORK(IU)
2773 *                    (CWorkspace: need   2*M*M+3*M-1,
2774 *                                 prefer 2*M*M+2*M+(M-1)*NB)
2775 *                    (RWorkspace: 0)
2776 *
2777                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
2778      $                            WORK( ITAUP ), WORK( IWORK ),
2779      $                            LWORK-IWORK+1, IERR )
2780 *
2781 *                    Generate left bidiagonalizing vectors in WORK(IR)
2782 *                    (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
2783 *                    (RWorkspace: 0)
2784 *
2785                      CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
2786      $                            WORK( ITAUQ ), WORK( IWORK ),
2787      $                            LWORK-IWORK+1, IERR )
2788                      IRWORK = IE + M
2789 *
2790 *                    Perform bidiagonal QR iteration, computing left
2791 *                    singular vectors of L in WORK(IR) and computing
2792 *                    right singular vectors of L in WORK(IU)
2793 *                    (CWorkspace: need 2*M*M)
2794 *                    (RWorkspace: need BDSPAC)
2795 *
2796                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
2797      $                            WORK( IU ), LDWRKU, WORK( IR ),
2798      $                            LDWRKR, CDUM, 1, RWORK( IRWORK ),
2799      $                            INFO )
2800 *
2801 *                    Multiply right singular vectors of L in WORK(IU) by
2802 *                    Q in A, storing result in VT
2803 *                    (CWorkspace: need M*M)
2804 *                    (RWorkspace: 0)
2805 *
2806                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
2807      $                           LDWRKU, A, LDA, CZERO, VT, LDVT )
2808 *
2809 *                    Copy left singular vectors of L to A
2810 *                    (CWorkspace: need M*M)
2811 *                    (RWorkspace: 0)
2812 *
2813                      CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
2814      $                            LDA )
2815 *
2816                   ELSE
2817 *
2818 *                    Insufficient workspace for a fast algorithm
2819 *
2820                      ITAU = 1
2821                      IWORK = ITAU + M
2822 *
2823 *                    Compute A=L*Q, copying result to VT
2824 *                    (CWorkspace: need 2*M, prefer M+M*NB)
2825 *                    (RWorkspace: 0)
2826 *
2827                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2828      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2829                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
2830 *
2831 *                    Generate Q in VT
2832 *                    (CWorkspace: need 2*M, prefer M+M*NB)
2833 *                    (RWorkspace: 0)
2834 *
2835                      CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
2836      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2837                      IE = 1
2838                      ITAUQ = ITAU
2839                      ITAUP = ITAUQ + M
2840                      IWORK = ITAUP + M
2841 *
2842 *                    Zero out above L in A
2843 *
2844                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
2845      $                            A( 1, 2 ), LDA )
2846 *
2847 *                    Bidiagonalize L in A
2848 *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2849 *                    (RWorkspace: need M)
2850 *
2851                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
2852      $                            WORK( ITAUQ ), WORK( ITAUP ),
2853      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2854 *
2855 *                    Multiply right vectors bidiagonalizing L by Q in VT
2856 *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2857 *                    (RWorkspace: 0)
2858 *
2859                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
2860      $                            WORK( ITAUP ), VT, LDVT,
2861      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2862 *
2863 *                    Generate left bidiagonalizing vectors of L in A
2864 *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
2865 *                    (RWorkspace: 0)
2866 *
2867                      CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
2868      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2869                      IRWORK = IE + M
2870 *
2871 *                    Perform bidiagonal QR iteration, computing left
2872 *                    singular vectors of A in A and computing right
2873 *                    singular vectors of A in VT
2874 *                    (CWorkspace: 0)
2875 *                    (RWorkspace: need BDSPAC)
2876 *
2877                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
2878      $                            LDVT, A, LDA, CDUM, 1,
2879      $                            RWORK( IRWORK ), INFO )
2880 *
2881                   END IF
2882 *
2883                ELSE IF( WNTUAS ) THEN
2884 *
2885 *                 Path 6t(N much larger than M, JOBU='S' or 'A',
2886 *                         JOBVT='S')
2887 *                 M right singular vectors to be computed in VT and
2888 *                 M left singular vectors to be computed in U
2889 *
2890                   IF( LWORK.GE.M*M+3*M ) THEN
2891 *
2892 *                    Sufficient workspace for a fast algorithm
2893 *
2894                      IU = 1
2895                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
2896 *
2897 *                       WORK(IU) is LDA by N
2898 *
2899                         LDWRKU = LDA
2900                      ELSE
2901 *
2902 *                       WORK(IU) is LDA by M
2903 *
2904                         LDWRKU = M
2905                      END IF
2906                      ITAU = IU + LDWRKU*M
2907                      IWORK = ITAU + M
2908 *
2909 *                    Compute A=L*Q
2910 *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2911 *                    (RWorkspace: 0)
2912 *
2913                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2914      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2915 *
2916 *                    Copy L to WORK(IU), zeroing out above it
2917 *
2918                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
2919      $                            LDWRKU )
2920                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
2921      $                            WORK( IU+LDWRKU ), LDWRKU )
2922 *
2923 *                    Generate Q in A
2924 *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2925 *                    (RWorkspace: 0)
2926 *
2927                      CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2928      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2929                      IE = 1
2930                      ITAUQ = ITAU
2931                      ITAUP = ITAUQ + M
2932                      IWORK = ITAUP + M
2933 *
2934 *                    Bidiagonalize L in WORK(IU), copying result to U
2935 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2936 *                    (RWorkspace: need M)
2937 *
2938                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
2939      $                            RWORK( IE ), WORK( ITAUQ ),
2940      $                            WORK( ITAUP ), WORK( IWORK ),
2941      $                            LWORK-IWORK+1, IERR )
2942                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
2943      $                            LDU )
2944 *
2945 *                    Generate right bidiagonalizing vectors in WORK(IU)
2946 *                    (CWorkspace: need   M*M+3*M-1,
2947 *                                 prefer M*M+2*M+(M-1)*NB)
2948 *                    (RWorkspace: 0)
2949 *
2950                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
2951      $                            WORK( ITAUP ), WORK( IWORK ),
2952      $                            LWORK-IWORK+1, IERR )
2953 *
2954 *                    Generate left bidiagonalizing vectors in U
2955 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2956 *                    (RWorkspace: 0)
2957 *
2958                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
2959      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2960                      IRWORK = IE + M
2961 *
2962 *                    Perform bidiagonal QR iteration, computing left
2963 *                    singular vectors of L in U and computing right
2964 *                    singular vectors of L in WORK(IU)
2965 *                    (CWorkspace: need M*M)
2966 *                    (RWorkspace: need BDSPAC)
2967 *
2968                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
2969      $                            WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
2970      $                            RWORK( IRWORK ), INFO )
2971 *
2972 *                    Multiply right singular vectors of L in WORK(IU) by
2973 *                    Q in A, storing result in VT
2974 *                    (CWorkspace: need M*M)
2975 *                    (RWorkspace: 0)
2976 *
2977                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
2978      $                           LDWRKU, A, LDA, CZERO, VT, LDVT )
2979 *
2980                   ELSE
2981 *
2982 *                    Insufficient workspace for a fast algorithm
2983 *
2984                      ITAU = 1
2985                      IWORK = ITAU + M
2986 *
2987 *                    Compute A=L*Q, copying result to VT
2988 *                    (CWorkspace: need 2*M, prefer M+M*NB)
2989 *                    (RWorkspace: 0)
2990 *
2991                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
2992      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2993                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
2994 *
2995 *                    Generate Q in VT
2996 *                    (CWorkspace: need 2*M, prefer M+M*NB)
2997 *                    (RWorkspace: 0)
2998 *
2999                      CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
3000      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3001 *
3002 *                    Copy L to U, zeroing out above it
3003 *
3004                      CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
3005                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
3006      $                            U( 1, 2 ), LDU )
3007                      IE = 1
3008                      ITAUQ = ITAU
3009                      ITAUP = ITAUQ + M
3010                      IWORK = ITAUP + M
3011 *
3012 *                    Bidiagonalize L in U
3013 *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3014 *                    (RWorkspace: need M)
3015 *
3016                      CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
3017      $                            WORK( ITAUQ ), WORK( ITAUP ),
3018      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3019 *
3020 *                    Multiply right bidiagonalizing vectors in U by Q
3021 *                    in VT
3022 *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3023 *                    (RWorkspace: 0)
3024 *
3025                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
3026      $                            WORK( ITAUP ), VT, LDVT,
3027      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3028 *
3029 *                    Generate left bidiagonalizing vectors in U
3030 *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
3031 *                    (RWorkspace: 0)
3032 *
3033                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
3034      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3035                      IRWORK = IE + M
3036 *
3037 *                    Perform bidiagonal QR iteration, computing left
3038 *                    singular vectors of A in U and computing right
3039 *                    singular vectors of A in VT
3040 *                    (CWorkspace: 0)
3041 *                    (RWorkspace: need BDSPAC)
3042 *
3043                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
3044      $                            LDVT, U, LDU, CDUM, 1,
3045      $                            RWORK( IRWORK ), INFO )
3046 *
3047                   END IF
3048 *
3049                END IF
3050 *
3051             ELSE IF( WNTVA ) THEN
3052 *
3053                IF( WNTUN ) THEN
3054 *
3055 *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
3056 *                 N right singular vectors to be computed in VT and
3057 *                 no left singular vectors to be computed
3058 *
3059                   IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
3060 *
3061 *                    Sufficient workspace for a fast algorithm
3062 *
3063                      IR = 1
3064                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
3065 *
3066 *                       WORK(IR) is LDA by M
3067 *
3068                         LDWRKR = LDA
3069                      ELSE
3070 *
3071 *                       WORK(IR) is M by M
3072 *
3073                         LDWRKR = M
3074                      END IF
3075                      ITAU = IR + LDWRKR*M
3076                      IWORK = ITAU + M
3077 *
3078 *                    Compute A=L*Q, copying result to VT
3079 *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3080 *                    (RWorkspace: 0)
3081 *
3082                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
3083      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3084                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
3085 *
3086 *                    Copy L to WORK(IR), zeroing out above it
3087 *
3088                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ),
3089      $                            LDWRKR )
3090                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
3091      $                            WORK( IR+LDWRKR ), LDWRKR )
3092 *
3093 *                    Generate Q in VT
3094 *                    (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3095 *                    (RWorkspace: 0)
3096 *
3097                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3098      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3099                      IE = 1
3100                      ITAUQ = ITAU
3101                      ITAUP = ITAUQ + M
3102                      IWORK = ITAUP + M
3103 *
3104 *                    Bidiagonalize L in WORK(IR)
3105 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3106 *                    (RWorkspace: need M)
3107 *
3108                      CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S,
3109      $                            RWORK( IE ), WORK( ITAUQ ),
3110      $                            WORK( ITAUP ), WORK( IWORK ),
3111      $                            LWORK-IWORK+1, IERR )
3112 *
3113 *                    Generate right bidiagonalizing vectors in WORK(IR)
3114 *                    (CWorkspace: need   M*M+3*M-1,
3115 *                                 prefer M*M+2*M+(M-1)*NB)
3116 *                    (RWorkspace: 0)
3117 *
3118                      CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
3119      $                            WORK( ITAUP ), WORK( IWORK ),
3120      $                            LWORK-IWORK+1, IERR )
3121                      IRWORK = IE + M
3122 *
3123 *                    Perform bidiagonal QR iteration, computing right
3124 *                    singular vectors of L in WORK(IR)
3125 *                    (CWorkspace: need M*M)
3126 *                    (RWorkspace: need BDSPAC)
3127 *
3128                      CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
3129      $                            WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
3130      $                            RWORK( IRWORK ), INFO )
3131 *
3132 *                    Multiply right singular vectors of L in WORK(IR) by
3133 *                    Q in VT, storing result in A
3134 *                    (CWorkspace: need M*M)
3135 *                    (RWorkspace: 0)
3136 *
3137                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
3138      $                           LDWRKR, VT, LDVT, CZERO, A, LDA )
3139 *
3140 *                    Copy right singular vectors of A from A to VT
3141 *
3142                      CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
3143 *
3144                   ELSE
3145 *
3146 *                    Insufficient workspace for a fast algorithm
3147 *
3148                      ITAU = 1
3149                      IWORK = ITAU + M
3150 *
3151 *                    Compute A=L*Q, copying result to VT
3152 *                    (CWorkspace: need 2*M, prefer M+M*NB)
3153 *                    (RWorkspace: 0)
3154 *
3155                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
3156      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3157                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
3158 *
3159 *                    Generate Q in VT
3160 *                    (CWorkspace: need M+N, prefer M+N*NB)
3161 *                    (RWorkspace: 0)
3162 *
3163                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3164      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3165                      IE = 1
3166                      ITAUQ = ITAU
3167                      ITAUP = ITAUQ + M
3168                      IWORK = ITAUP + M
3169 *
3170 *                    Zero out above L in A
3171 *
3172                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
3173      $                            A( 1, 2 ), LDA )
3174 *
3175 *                    Bidiagonalize L in A
3176 *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3177 *                    (RWorkspace: need M)
3178 *
3179                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
3180      $                            WORK( ITAUQ ), WORK( ITAUP ),
3181      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3182 *
3183 *                    Multiply right bidiagonalizing vectors in A by Q
3184 *                    in VT
3185 *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3186 *                    (RWorkspace: 0)
3187 *
3188                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
3189      $                            WORK( ITAUP ), VT, LDVT,
3190      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3191                      IRWORK = IE + M
3192 *
3193 *                    Perform bidiagonal QR iteration, computing right
3194 *                    singular vectors of A in VT
3195 *                    (CWorkspace: 0)
3196 *                    (RWorkspace: need BDSPAC)
3197 *
3198                      CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
3199      $                            LDVT, CDUM, 1, CDUM, 1,
3200      $                            RWORK( IRWORK ), INFO )
3201 *
3202                   END IF
3203 *
3204                ELSE IF( WNTUO ) THEN
3205 *
3206 *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3207 *                 N right singular vectors to be computed in VT and
3208 *                 M left singular vectors to be overwritten on A
3209 *
3210                   IF( LWORK.GE.2*M*M+MAX( N+M, 3*M ) ) THEN
3211 *
3212 *                    Sufficient workspace for a fast algorithm
3213 *
3214                      IU = 1
3215                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
3216 *
3217 *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
3218 *
3219                         LDWRKU = LDA
3220                         IR = IU + LDWRKU*M
3221                         LDWRKR = LDA
3222                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
3223 *
3224 *                       WORK(IU) is LDA by M and WORK(IR) is M by M
3225 *
3226                         LDWRKU = LDA
3227                         IR = IU + LDWRKU*M
3228                         LDWRKR = M
3229                      ELSE
3230 *
3231 *                       WORK(IU) is M by M and WORK(IR) is M by M
3232 *
3233                         LDWRKU = M
3234                         IR = IU + LDWRKU*M
3235                         LDWRKR = M
3236                      END IF
3237                      ITAU = IR + LDWRKR*M
3238                      IWORK = ITAU + M
3239 *
3240 *                    Compute A=L*Q, copying result to VT
3241 *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3242 *                    (RWorkspace: 0)
3243 *
3244                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
3245      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3246                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
3247 *
3248 *                    Generate Q in VT
3249 *                    (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3250 *                    (RWorkspace: 0)
3251 *
3252                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3253      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3254 *
3255 *                    Copy L to WORK(IU), zeroing out above it
3256 *
3257                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
3258      $                            LDWRKU )
3259                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
3260      $                            WORK( IU+LDWRKU ), LDWRKU )
3261                      IE = 1
3262                      ITAUQ = ITAU
3263                      ITAUP = ITAUQ + M
3264                      IWORK = ITAUP + M
3265 *
3266 *                    Bidiagonalize L in WORK(IU), copying result to
3267 *                    WORK(IR)
3268 *                    (CWorkspace: need   2*M*M+3*M,
3269 *                                 prefer 2*M*M+2*M+2*M*NB)
3270 *                    (RWorkspace: need   M)
3271 *
3272                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
3273      $                            RWORK( IE ), WORK( ITAUQ ),
3274      $                            WORK( ITAUP ), WORK( IWORK ),
3275      $                            LWORK-IWORK+1, IERR )
3276                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU,
3277      $                            WORK( IR ), LDWRKR )
3278 *
3279 *                    Generate right bidiagonalizing vectors in WORK(IU)
3280 *                    (CWorkspace: need   2*M*M+3*M-1,
3281 *                                 prefer 2*M*M+2*M+(M-1)*NB)
3282 *                    (RWorkspace: 0)
3283 *
3284                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
3285      $                            WORK( ITAUP ), WORK( IWORK ),
3286      $                            LWORK-IWORK+1, IERR )
3287 *
3288 *                    Generate left bidiagonalizing vectors in WORK(IR)
3289 *                    (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3290 *                    (RWorkspace: 0)
3291 *
3292                      CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
3293      $                            WORK( ITAUQ ), WORK( IWORK ),
3294      $                            LWORK-IWORK+1, IERR )
3295                      IRWORK = IE + M
3296 *
3297 *                    Perform bidiagonal QR iteration, computing left
3298 *                    singular vectors of L in WORK(IR) and computing
3299 *                    right singular vectors of L in WORK(IU)
3300 *                    (CWorkspace: need 2*M*M)
3301 *                    (RWorkspace: need BDSPAC)
3302 *
3303                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
3304      $                            WORK( IU ), LDWRKU, WORK( IR ),
3305      $                            LDWRKR, CDUM, 1, RWORK( IRWORK ),
3306      $                            INFO )
3307 *
3308 *                    Multiply right singular vectors of L in WORK(IU) by
3309 *                    Q in VT, storing result in A
3310 *                    (CWorkspace: need M*M)
3311 *                    (RWorkspace: 0)
3312 *
3313                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
3314      $                           LDWRKU, VT, LDVT, CZERO, A, LDA )
3315 *
3316 *                    Copy right singular vectors of A from A to VT
3317 *
3318                      CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
3319 *
3320 *                    Copy left singular vectors of A from WORK(IR) to A
3321 *
3322                      CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
3323      $                            LDA )
3324 *
3325                   ELSE
3326 *
3327 *                    Insufficient workspace for a fast algorithm
3328 *
3329                      ITAU = 1
3330                      IWORK = ITAU + M
3331 *
3332 *                    Compute A=L*Q, copying result to VT
3333 *                    (CWorkspace: need 2*M, prefer M+M*NB)
3334 *                    (RWorkspace: 0)
3335 *
3336                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
3337      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3338                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
3339 *
3340 *                    Generate Q in VT
3341 *                    (CWorkspace: need M+N, prefer M+N*NB)
3342 *                    (RWorkspace: 0)
3343 *
3344                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3345      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3346                      IE = 1
3347                      ITAUQ = ITAU
3348                      ITAUP = ITAUQ + M
3349                      IWORK = ITAUP + M
3350 *
3351 *                    Zero out above L in A
3352 *
3353                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
3354      $                            A( 1, 2 ), LDA )
3355 *
3356 *                    Bidiagonalize L in A
3357 *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3358 *                    (RWorkspace: need M)
3359 *
3360                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
3361      $                            WORK( ITAUQ ), WORK( ITAUP ),
3362      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3363 *
3364 *                    Multiply right bidiagonalizing vectors in A by Q
3365 *                    in VT
3366 *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3367 *                    (RWorkspace: 0)
3368 *
3369                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
3370      $                            WORK( ITAUP ), VT, LDVT,
3371      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3372 *
3373 *                    Generate left bidiagonalizing vectors in A
3374 *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
3375 *                    (RWorkspace: 0)
3376 *
3377                      CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
3378      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3379                      IRWORK = IE + M
3380 *
3381 *                    Perform bidiagonal QR iteration, computing left
3382 *                    singular vectors of A in A and computing right
3383 *                    singular vectors of A in VT
3384 *                    (CWorkspace: 0)
3385 *                    (RWorkspace: need BDSPAC)
3386 *
3387                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
3388      $                            LDVT, A, LDA, CDUM, 1,
3389      $                            RWORK( IRWORK ), INFO )
3390 *
3391                   END IF
3392 *
3393                ELSE IF( WNTUAS ) THEN
3394 *
3395 *                 Path 9t(N much larger than M, JOBU='S' or 'A',
3396 *                         JOBVT='A')
3397 *                 N right singular vectors to be computed in VT and
3398 *                 M left singular vectors to be computed in U
3399 *
3400                   IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
3401 *
3402 *                    Sufficient workspace for a fast algorithm
3403 *
3404                      IU = 1
3405                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
3406 *
3407 *                       WORK(IU) is LDA by M
3408 *
3409                         LDWRKU = LDA
3410                      ELSE
3411 *
3412 *                       WORK(IU) is M by M
3413 *
3414                         LDWRKU = M
3415                      END IF
3416                      ITAU = IU + LDWRKU*M
3417                      IWORK = ITAU + M
3418 *
3419 *                    Compute A=L*Q, copying result to VT
3420 *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3421 *                    (RWorkspace: 0)
3422 *
3423                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
3424      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3425                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
3426 *
3427 *                    Generate Q in VT
3428 *                    (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3429 *                    (RWorkspace: 0)
3430 *
3431                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3432      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3433 *
3434 *                    Copy L to WORK(IU), zeroing out above it
3435 *
3436                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
3437      $                            LDWRKU )
3438                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
3439      $                            WORK( IU+LDWRKU ), LDWRKU )
3440                      IE = 1
3441                      ITAUQ = ITAU
3442                      ITAUP = ITAUQ + M
3443                      IWORK = ITAUP + M
3444 *
3445 *                    Bidiagonalize L in WORK(IU), copying result to U
3446 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3447 *                    (RWorkspace: need M)
3448 *
3449                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
3450      $                            RWORK( IE ), WORK( ITAUQ ),
3451      $                            WORK( ITAUP ), WORK( IWORK ),
3452      $                            LWORK-IWORK+1, IERR )
3453                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
3454      $                            LDU )
3455 *
3456 *                    Generate right bidiagonalizing vectors in WORK(IU)
3457 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3458 *                    (RWorkspace: 0)
3459 *
3460                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
3461      $                            WORK( ITAUP ), WORK( IWORK ),
3462      $                            LWORK-IWORK+1, IERR )
3463 *
3464 *                    Generate left bidiagonalizing vectors in U
3465 *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3466 *                    (RWorkspace: 0)
3467 *
3468                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
3469      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3470                      IRWORK = IE + M
3471 *
3472 *                    Perform bidiagonal QR iteration, computing left
3473 *                    singular vectors of L in U and computing right
3474 *                    singular vectors of L in WORK(IU)
3475 *                    (CWorkspace: need M*M)
3476 *                    (RWorkspace: need BDSPAC)
3477 *
3478                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
3479      $                            WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
3480      $                            RWORK( IRWORK ), INFO )
3481 *
3482 *                    Multiply right singular vectors of L in WORK(IU) by
3483 *                    Q in VT, storing result in A
3484 *                    (CWorkspace: need M*M)
3485 *                    (RWorkspace: 0)
3486 *
3487                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
3488      $                           LDWRKU, VT, LDVT, CZERO, A, LDA )
3489 *
3490 *                    Copy right singular vectors of A from A to VT
3491 *
3492                      CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
3493 *
3494                   ELSE
3495 *
3496 *                    Insufficient workspace for a fast algorithm
3497 *
3498                      ITAU = 1
3499                      IWORK = ITAU + M
3500 *
3501 *                    Compute A=L*Q, copying result to VT
3502 *                    (CWorkspace: need 2*M, prefer M+M*NB)
3503 *                    (RWorkspace: 0)
3504 *
3505                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
3506      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3507                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
3508 *
3509 *                    Generate Q in VT
3510 *                    (CWorkspace: need M+N, prefer M+N*NB)
3511 *                    (RWorkspace: 0)
3512 *
3513                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3514      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3515 *
3516 *                    Copy L to U, zeroing out above it
3517 *
3518                      CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
3519                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
3520      $                            U( 1, 2 ), LDU )
3521                      IE = 1
3522                      ITAUQ = ITAU
3523                      ITAUP = ITAUQ + M
3524                      IWORK = ITAUP + M
3525 *
3526 *                    Bidiagonalize L in U
3527 *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3528 *                    (RWorkspace: need M)
3529 *
3530                      CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
3531      $                            WORK( ITAUQ ), WORK( ITAUP ),
3532      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3533 *
3534 *                    Multiply right bidiagonalizing vectors in U by Q
3535 *                    in VT
3536 *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3537 *                    (RWorkspace: 0)
3538 *
3539                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
3540      $                            WORK( ITAUP ), VT, LDVT,
3541      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3542 *
3543 *                    Generate left bidiagonalizing vectors in U
3544 *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
3545 *                    (RWorkspace: 0)
3546 *
3547                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
3548      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3549                      IRWORK = IE + M
3550 *
3551 *                    Perform bidiagonal QR iteration, computing left
3552 *                    singular vectors of A in U and computing right
3553 *                    singular vectors of A in VT
3554 *                    (CWorkspace: 0)
3555 *                    (RWorkspace: need BDSPAC)
3556 *
3557                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
3558      $                            LDVT, U, LDU, CDUM, 1,
3559      $                            RWORK( IRWORK ), INFO )
3560 *
3561                   END IF
3562 *
3563                END IF
3564 *
3565             END IF
3566 *
3567          ELSE
3568 *
3569 *           N .LT. MNTHR
3570 *
3571 *           Path 10t(N greater than M, but not much larger)
3572 *           Reduce to bidiagonal form without LQ decomposition
3573 *
3574             IE = 1
3575             ITAUQ = 1
3576             ITAUP = ITAUQ + M
3577             IWORK = ITAUP + M
3578 *
3579 *           Bidiagonalize A
3580 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
3581 *           (RWorkspace: M)
3582 *
3583             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
3584      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
3585      $                   IERR )
3586             IF( WNTUAS ) THEN
3587 *
3588 *              If left singular vectors desired in U, copy result to U
3589 *              and generate left bidiagonalizing vectors in U
3590 *              (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3591 *              (RWorkspace: 0)
3592 *
3593                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
3594                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
3595      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3596             END IF
3597             IF( WNTVAS ) THEN
3598 *
3599 *              If right singular vectors desired in VT, copy result to
3600 *              VT and generate right bidiagonalizing vectors in VT
3601 *              (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
3602 *              (RWorkspace: 0)
3603 *
3604                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
3605                IF( WNTVA )
3606      $            NRVT = N
3607                IF( WNTVS )
3608      $            NRVT = M
3609                CALL ZUNGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ),
3610      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3611             END IF
3612             IF( WNTUO ) THEN
3613 *
3614 *              If left singular vectors desired in A, generate left
3615 *              bidiagonalizing vectors in A
3616 *              (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3617 *              (RWorkspace: 0)
3618 *
3619                CALL ZUNGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
3620      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3621             END IF
3622             IF( WNTVO ) THEN
3623 *
3624 *              If right singular vectors desired in A, generate right
3625 *              bidiagonalizing vectors in A
3626 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
3627 *              (RWorkspace: 0)
3628 *
3629                CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
3630      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3631             END IF
3632             IRWORK = IE + M
3633             IF( WNTUAS .OR. WNTUO )
3634      $         NRU = M
3635             IF( WNTUN )
3636      $         NRU = 0
3637             IF( WNTVAS .OR. WNTVO )
3638      $         NCVT = N
3639             IF( WNTVN )
3640      $         NCVT = 0
3641             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
3642 *
3643 *              Perform bidiagonal QR iteration, if desired, computing
3644 *              left singular vectors in U and computing right singular
3645 *              vectors in VT
3646 *              (CWorkspace: 0)
3647 *              (RWorkspace: need BDSPAC)
3648 *
3649                CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
3650      $                      LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
3651      $                      INFO )
3652             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
3653 *
3654 *              Perform bidiagonal QR iteration, if desired, computing
3655 *              left singular vectors in U and computing right singular
3656 *              vectors in A
3657 *              (CWorkspace: 0)
3658 *              (RWorkspace: need BDSPAC)
3659 *
3660                CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), A,
3661      $                      LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
3662      $                      INFO )
3663             ELSE
3664 *
3665 *              Perform bidiagonal QR iteration, if desired, computing
3666 *              left singular vectors in A and computing right singular
3667 *              vectors in VT
3668 *              (CWorkspace: 0)
3669 *              (RWorkspace: need BDSPAC)
3670 *
3671                CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
3672      $                      LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
3673      $                      INFO )
3674             END IF
3675 *
3676          END IF
3677 *
3678       END IF
3679 *
3680 *     Undo scaling if necessary
3681 *
3682       IF( ISCL.EQ.1 ) THEN
3683          IF( ANRM.GT.BIGNUM )
3684      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
3685      $                   IERR )
3686          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
3687      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
3688      $                   RWORK( IE ), MINMN, IERR )
3689          IF( ANRM.LT.SMLNUM )
3690      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
3691      $                   IERR )
3692          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
3693      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
3694      $                   RWORK( IE ), MINMN, IERR )
3695       END IF
3696 *
3697 *     Return optimal workspace in WORK(1)
3698 *
3699       WORK( 1 ) = MAXWRK
3700 *
3701       RETURN
3702 *
3703 *     End of ZGESVD
3704 *
3705       END