SRC/cgesvd.f

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