SRC/cgesdd.f

   1 *> \brief \b CGESDD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGESDD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesdd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesdd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesdd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
  22 *                          WORK, LWORK, RWORK, IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ
  26 *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       REAL               RWORK( * ), S( * )
  31 *       COMPLEX            A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
  32 *      $                   WORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> CGESDD computes the singular value decomposition (SVD) of a complex
  42 *> M-by-N matrix A, optionally computing the left and/or right singular
  43 *> vectors, by using divide-and-conquer method. The SVD is written
  44 *>
  45 *>      A = U * SIGMA * conjugate-transpose(V)
  46 *>
  47 *> where SIGMA is an M-by-N matrix which is zero except for its
  48 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
  49 *> V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
  50 *> are the singular values of A; they are real and non-negative, and
  51 *> are returned in descending order.  The first min(m,n) columns of
  52 *> U and V are the left and right singular vectors of A.
  53 *>
  54 *> Note that the routine returns VT = V**H, not V.
  55 *>
  56 *> The divide and conquer algorithm makes very mild assumptions about
  57 *> floating point arithmetic. It will work on machines with a guard
  58 *> digit in add/subtract, or on those binary machines without guard
  59 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
  60 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
  61 *> without guard digits, but we know of none.
  62 *> \endverbatim
  63 *
  64 *  Arguments:
  65 *  ==========
  66 *
  67 *> \param[in] JOBZ
  68 *> \verbatim
  69 *>          JOBZ is CHARACTER*1
  70 *>          Specifies options for computing all or part of the matrix U:
  71 *>          = 'A':  all M columns of U and all N rows of V**H are
  72 *>                  returned in the arrays U and VT;
  73 *>          = 'S':  the first min(M,N) columns of U and the first
  74 *>                  min(M,N) rows of V**H are returned in the arrays U
  75 *>                  and VT;
  76 *>          = 'O':  If M >= N, the first N columns of U are overwritten
  77 *>                  in the array A and all rows of V**H are returned in
  78 *>                  the array VT;
  79 *>                  otherwise, all columns of U are returned in the
  80 *>                  array U and the first M rows of V**H are overwritten
  81 *>                  in the array A;
  82 *>          = 'N':  no columns of U or rows of V**H are computed.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] M
  86 *> \verbatim
  87 *>          M is INTEGER
  88 *>          The number of rows of the input matrix A.  M >= 0.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] N
  92 *> \verbatim
  93 *>          N is INTEGER
  94 *>          The number of columns of the input matrix A.  N >= 0.
  95 *> \endverbatim
  96 *>
  97 *> \param[in,out] A
  98 *> \verbatim
  99 *>          A is COMPLEX array, dimension (LDA,N)
 100 *>          On entry, the M-by-N matrix A.
 101 *>          On exit,
 102 *>          if JOBZ = 'O',  A is overwritten with the first N columns
 103 *>                          of U (the left singular vectors, stored
 104 *>                          columnwise) if M >= N;
 105 *>                          A is overwritten with the first M rows
 106 *>                          of V**H (the right singular vectors, stored
 107 *>                          rowwise) otherwise.
 108 *>          if JOBZ .ne. 'O', the contents of A are destroyed.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] LDA
 112 *> \verbatim
 113 *>          LDA is INTEGER
 114 *>          The leading dimension of the array A.  LDA >= max(1,M).
 115 *> \endverbatim
 116 *>
 117 *> \param[out] S
 118 *> \verbatim
 119 *>          S is REAL array, dimension (min(M,N))
 120 *>          The singular values of A, sorted so that S(i) >= S(i+1).
 121 *> \endverbatim
 122 *>
 123 *> \param[out] U
 124 *> \verbatim
 125 *>          U is COMPLEX array, dimension (LDU,UCOL)
 126 *>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
 127 *>          UCOL = min(M,N) if JOBZ = 'S'.
 128 *>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
 129 *>          unitary matrix U;
 130 *>          if JOBZ = 'S', U contains the first min(M,N) columns of U
 131 *>          (the left singular vectors, stored columnwise);
 132 *>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] LDU
 136 *> \verbatim
 137 *>          LDU is INTEGER
 138 *>          The leading dimension of the array U.  LDU >= 1;
 139 *>          if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
 140 *> \endverbatim
 141 *>
 142 *> \param[out] VT
 143 *> \verbatim
 144 *>          VT is COMPLEX array, dimension (LDVT,N)
 145 *>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
 146 *>          N-by-N unitary matrix V**H;
 147 *>          if JOBZ = 'S', VT contains the first min(M,N) rows of
 148 *>          V**H (the right singular vectors, stored rowwise);
 149 *>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
 150 *> \endverbatim
 151 *>
 152 *> \param[in] LDVT
 153 *> \verbatim
 154 *>          LDVT is INTEGER
 155 *>          The leading dimension of the array VT.  LDVT >= 1;
 156 *>          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
 157 *>          if JOBZ = 'S', LDVT >= min(M,N).
 158 *> \endverbatim
 159 *>
 160 *> \param[out] WORK
 161 *> \verbatim
 162 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 163 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 164 *> \endverbatim
 165 *>
 166 *> \param[in] LWORK
 167 *> \verbatim
 168 *>          LWORK is INTEGER
 169 *>          The dimension of the array WORK. LWORK >= 1.
 170 *>          If LWORK = -1, a workspace query is assumed.  The optimal
 171 *>          size for the WORK array is calculated and stored in WORK(1),
 172 *>          and no other work except argument checking is performed.
 173 *>
 174 *>          Let mx = max(M,N) and mn = min(M,N).
 175 *>          If JOBZ = 'N', LWORK >= 2*mn + mx.
 176 *>          If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx.
 177 *>          If JOBZ = 'S', LWORK >=   mn*mn + 3*mn.
 178 *>          If JOBZ = 'A', LWORK >=   mn*mn + 2*mn + mx.
 179 *>          These are not tight minimums in all cases; see comments inside code.
 180 *>          For good performance, LWORK should generally be larger;
 181 *>          a query is recommended.
 182 *> \endverbatim
 183 *>
 184 *> \param[out] RWORK
 185 *> \verbatim
 186 *>          RWORK is REAL array, dimension (MAX(1,LRWORK))
 187 *>          Let mx = max(M,N) and mn = min(M,N).
 188 *>          If JOBZ = 'N',    LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn);
 189 *>          else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn;
 190 *>          else              LRWORK >= max( 5*mn*mn + 5*mn,
 191 *>                                           2*mx*mn + 2*mn*mn + mn ).
 192 *> \endverbatim
 193 *>
 194 *> \param[out] IWORK
 195 *> \verbatim
 196 *>          IWORK is INTEGER array, dimension (8*min(M,N))
 197 *> \endverbatim
 198 *>
 199 *> \param[out] INFO
 200 *> \verbatim
 201 *>          INFO is INTEGER
 202 *>          = 0:  successful exit.
 203 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 204 *>          > 0:  The updating process of SBDSDC did not converge.
 205 *> \endverbatim
 206 *
 207 *  Authors:
 208 *  ========
 209 *
 210 *> \author Univ. of Tennessee
 211 *> \author Univ. of California Berkeley
 212 *> \author Univ. of Colorado Denver
 213 *> \author NAG Ltd.
 214 *
 215 *> \date June 2016
 216 *
 217 *> \ingroup complexGEsing
 218 *
 219 *> \par Contributors:
 220 *  ==================
 221 *>
 222 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 223 *>     California at Berkeley, USA
 224 *>
 225 *  =====================================================================
 226       SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 227      $                   WORK, LWORK, RWORK, IWORK, INFO )
 228       implicit none
 229 *
 230 *  -- LAPACK driver routine (version 3.6.1) --
 231 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 232 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 233 *     June 2016
 234 *
 235 *     .. Scalar Arguments ..
 236       CHARACTER          JOBZ
 237       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
 238 *     ..
 239 *     .. Array Arguments ..
 240       INTEGER            IWORK( * )
 241       REAL               RWORK( * ), S( * )
 242       COMPLEX            A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
 243      $                   WORK( * )
 244 *     ..
 245 *
 246 *  =====================================================================
 247 *
 248 *     .. Parameters ..
 249       COMPLEX            CZERO, CONE
 250       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 251      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
 252       REAL               ZERO, ONE
 253       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 254 *     ..
 255 *     .. Local Scalars ..
 256       LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
 257       INTEGER            BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
 258      $                   ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
 259      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
 260      $                   MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
 261       INTEGER            LWORK_CGEBRD_MN, LWORK_CGEBRD_MM,
 262      $                   LWORK_CGEBRD_NN, LWORK_CGELQF_MN,
 263      $                   LWORK_CGEQRF_MN,
 264      $                   LWORK_CUNGBR_P_MN, LWORK_CUNGBR_P_NN,
 265      $                   LWORK_CUNGBR_Q_MN, LWORK_CUNGBR_Q_MM,
 266      $                   LWORK_CUNGLQ_MN, LWORK_CUNGLQ_NN,
 267      $                   LWORK_CUNGQR_MM, LWORK_CUNGQR_MN,
 268      $                   LWORK_CUNMBR_PRC_MM, LWORK_CUNMBR_QLN_MM,
 269      $                   LWORK_CUNMBR_PRC_MN, LWORK_CUNMBR_QLN_MN,
 270      $                   LWORK_CUNMBR_PRC_NN, LWORK_CUNMBR_QLN_NN
 271       REAL   ANRM, BIGNUM, EPS, SMLNUM
 272 *     ..
 273 *     .. Local Arrays ..
 274       INTEGER            IDUM( 1 )
 275       REAL               DUM( 1 )
 276       COMPLEX            CDUM( 1 )
 277 *     ..
 278 *     .. External Subroutines ..
 279       EXTERNAL           CGEBRD, CGELQF, CGEMM, CGEQRF, CLACP2, CLACPY,
 280      $                   CLACRM, CLARCM, CLASCL, CLASET, CUNGBR, CUNGLQ,
 281      $                   CUNGQR, CUNMBR, SBDSDC, SLASCL, XERBLA
 282 *     ..
 283 *     .. External Functions ..
 284       LOGICAL            LSAME
 285       REAL               SLAMCH, CLANGE
 286       EXTERNAL           LSAME, SLAMCH, CLANGE
 287 *     ..
 288 *     .. Intrinsic Functions ..
 289       INTRINSIC          INT, MAX, MIN, SQRT
 290 *     ..
 291 *     .. Executable Statements ..
 292 *
 293 *     Test the input arguments
 294 *
 295       INFO   = 0
 296       MINMN  = MIN( M, N )
 297       MNTHR1 = INT( MINMN*17.0E0 / 9.0E0 )
 298       MNTHR2 = INT( MINMN*5.0E0 / 3.0E0 )
 299       WNTQA  = LSAME( JOBZ, 'A' )
 300       WNTQS  = LSAME( JOBZ, 'S' )
 301       WNTQAS = WNTQA .OR. WNTQS
 302       WNTQO  = LSAME( JOBZ, 'O' )
 303       WNTQN  = LSAME( JOBZ, 'N' )
 304       LQUERY = ( LWORK.EQ.-1 )
 305       MINWRK = 1
 306       MAXWRK = 1
 307 *
 308       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
 309          INFO = -1
 310       ELSE IF( M.LT.0 ) THEN
 311          INFO = -2
 312       ELSE IF( N.LT.0 ) THEN
 313          INFO = -3
 314       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 315          INFO = -5
 316       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
 317      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
 318          INFO = -8
 319       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
 320      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
 321      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
 322          INFO = -10
 323       END IF
 324 *
 325 *     Compute workspace
 326 *       Note: Comments in the code beginning "Workspace:" describe the
 327 *       minimal amount of workspace allocated at that point in the code,
 328 *       as well as the preferred amount for good performance.
 329 *       CWorkspace refers to complex workspace, and RWorkspace to
 330 *       real workspace. NB refers to the optimal block size for the
 331 *       immediately following subroutine, as returned by ILAENV.)
 332 *
 333       IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN
 334          IF( M.GE.N ) THEN
 335 *
 336 *           There is no complex work space needed for bidiagonal SVD
 337 *           The real work space needed for bidiagonal SVD (sbdsdc) is
 338 *           BDSPAC = 3*N*N + 4*N for singular values and vectors;
 339 *           BDSPAC = 4*N         for singular values only;
 340 *           not including e, RU, and RVT matrices.
 341 *
 342 *           Compute space preferred for each routine
 343             CALL CGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
 344      $                   CDUM(1), CDUM(1), -1, IERR )
 345             LWORK_CGEBRD_MN = INT( CDUM(1) )
 346 *
 347             CALL CGEBRD( N, N, CDUM(1), N, DUM(1), DUM(1), CDUM(1),
 348      $                   CDUM(1), CDUM(1), -1, IERR )
 349             LWORK_CGEBRD_NN = INT( CDUM(1) )
 350 *
 351             CALL CGEQRF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR )
 352             LWORK_CGEQRF_MN = INT( CDUM(1) )
 353 *
 354             CALL CUNGBR( 'P', N, N, N, CDUM(1), N, CDUM(1), CDUM(1),
 355      $                   -1, IERR )
 356             LWORK_CUNGBR_P_NN = INT( CDUM(1) )
 357 *
 358             CALL CUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
 359      $                   -1, IERR )
 360             LWORK_CUNGBR_Q_MM = INT( CDUM(1) )
 361 *
 362             CALL CUNGBR( 'Q', M, N, N, CDUM(1), M, CDUM(1), CDUM(1),
 363      $                   -1, IERR )
 364             LWORK_CUNGBR_Q_MN = INT( CDUM(1) )
 365 *
 366             CALL CUNGQR( M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
 367      $                   -1, IERR )
 368             LWORK_CUNGQR_MM = INT( CDUM(1) )
 369 *
 370             CALL CUNGQR( M, N, N, CDUM(1), M, CDUM(1), CDUM(1),
 371      $                   -1, IERR )
 372             LWORK_CUNGQR_MN = INT( CDUM(1) )
 373 *
 374             CALL CUNMBR( 'P', 'R', 'C', N, N, N, CDUM(1), N, CDUM(1),
 375      $                   CDUM(1), N, CDUM(1), -1, IERR )
 376             LWORK_CUNMBR_PRC_NN = INT( CDUM(1) )
 377 *
 378             CALL CUNMBR( 'Q', 'L', 'N', M, M, N, CDUM(1), M, CDUM(1),
 379      $                   CDUM(1), M, CDUM(1), -1, IERR )
 380             LWORK_CUNMBR_QLN_MM = INT( CDUM(1) )
 381 *
 382             CALL CUNMBR( 'Q', 'L', 'N', M, N, N, CDUM(1), M, CDUM(1),
 383      $                   CDUM(1), M, CDUM(1), -1, IERR )
 384             LWORK_CUNMBR_QLN_MN = INT( CDUM(1) )
 385 *
 386             CALL CUNMBR( 'Q', 'L', 'N', N, N, N, CDUM(1), N, CDUM(1),
 387      $                   CDUM(1), N, CDUM(1), -1, IERR )
 388             LWORK_CUNMBR_QLN_NN = INT( CDUM(1) )
 389 *
 390             IF( M.GE.MNTHR1 ) THEN
 391                IF( WNTQN ) THEN
 392 *
 393 *                 Path 1 (M >> N, JOBZ='N')
 394 *
 395                   MAXWRK = N + LWORK_CGEQRF_MN
 396                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CGEBRD_NN )
 397                   MINWRK = 3*N
 398                ELSE IF( WNTQO ) THEN
 399 *
 400 *                 Path 2 (M >> N, JOBZ='O')
 401 *
 402                   WRKBL = N + LWORK_CGEQRF_MN
 403                   WRKBL = MAX( WRKBL,   N + LWORK_CUNGQR_MN )
 404                   WRKBL = MAX( WRKBL, 2*N + LWORK_CGEBRD_NN )
 405                   WRKBL = MAX( WRKBL, 2*N + LWORK_CUNMBR_QLN_NN )
 406                   WRKBL = MAX( WRKBL, 2*N + LWORK_CUNMBR_PRC_NN )
 407                   MAXWRK = M*N + N*N + WRKBL
 408                   MINWRK = 2*N*N + 3*N
 409                ELSE IF( WNTQS ) THEN
 410 *
 411 *                 Path 3 (M >> N, JOBZ='S')
 412 *
 413                   WRKBL = N + LWORK_CGEQRF_MN
 414                   WRKBL = MAX( WRKBL,   N + LWORK_CUNGQR_MN )
 415                   WRKBL = MAX( WRKBL, 2*N + LWORK_CGEBRD_NN )
 416                   WRKBL = MAX( WRKBL, 2*N + LWORK_CUNMBR_QLN_NN )
 417                   WRKBL = MAX( WRKBL, 2*N + LWORK_CUNMBR_PRC_NN )
 418                   MAXWRK = N*N + WRKBL
 419                   MINWRK = N*N + 3*N
 420                ELSE IF( WNTQA ) THEN
 421 *
 422 *                 Path 4 (M >> N, JOBZ='A')
 423 *
 424                   WRKBL = N + LWORK_CGEQRF_MN
 425                   WRKBL = MAX( WRKBL,   N + LWORK_CUNGQR_MM )
 426                   WRKBL = MAX( WRKBL, 2*N + LWORK_CGEBRD_NN )
 427                   WRKBL = MAX( WRKBL, 2*N + LWORK_CUNMBR_QLN_NN )
 428                   WRKBL = MAX( WRKBL, 2*N + LWORK_CUNMBR_PRC_NN )
 429                   MAXWRK = N*N + WRKBL
 430                   MINWRK = N*N + MAX( 3*N, N + M )
 431                END IF
 432             ELSE IF( M.GE.MNTHR2 ) THEN
 433 *
 434 *              Path 5 (M >> N, but not as much as MNTHR1)
 435 *
 436                MAXWRK = 2*N + LWORK_CGEBRD_MN
 437                MINWRK = 2*N + M
 438                IF( WNTQO ) THEN
 439 *                 Path 5o (M >> N, JOBZ='O')
 440                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR_P_NN )
 441                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR_Q_MN )
 442                   MAXWRK = MAXWRK + M*N
 443                   MINWRK = MINWRK + N*N
 444                ELSE IF( WNTQS ) THEN
 445 *                 Path 5s (M >> N, JOBZ='S')
 446                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR_P_NN )
 447                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR_Q_MN )
 448                ELSE IF( WNTQA ) THEN
 449 *                 Path 5a (M >> N, JOBZ='A')
 450                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR_P_NN )
 451                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR_Q_MM )
 452                END IF
 453             ELSE
 454 *
 455 *              Path 6 (M >= N, but not much larger)
 456 *
 457                MAXWRK = 2*N + LWORK_CGEBRD_MN
 458                MINWRK = 2*N + M
 459                IF( WNTQO ) THEN
 460 *                 Path 6o (M >= N, JOBZ='O')
 461                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR_PRC_NN )
 462                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR_QLN_MN )
 463                   MAXWRK = MAXWRK + M*N
 464                   MINWRK = MINWRK + N*N
 465                ELSE IF( WNTQS ) THEN
 466 *                 Path 6s (M >= N, JOBZ='S')
 467                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR_QLN_MN )
 468                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR_PRC_NN )
 469                ELSE IF( WNTQA ) THEN
 470 *                 Path 6a (M >= N, JOBZ='A')
 471                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR_QLN_MM )
 472                   MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR_PRC_NN )
 473                END IF
 474             END IF
 475          ELSE
 476 *
 477 *           There is no complex work space needed for bidiagonal SVD
 478 *           The real work space needed for bidiagonal SVD (sbdsdc) is
 479 *           BDSPAC = 3*M*M + 4*M for singular values and vectors;
 480 *           BDSPAC = 4*M         for singular values only;
 481 *           not including e, RU, and RVT matrices.
 482 *
 483 *           Compute space preferred for each routine
 484             CALL CGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
 485      $                   CDUM(1), CDUM(1), -1, IERR )
 486             LWORK_CGEBRD_MN = INT( CDUM(1) )
 487 *
 488             CALL CGEBRD( M, M, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
 489      $                   CDUM(1), CDUM(1), -1, IERR )
 490             LWORK_CGEBRD_MM = INT( CDUM(1) )
 491 *
 492             CALL CGELQF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR )
 493             LWORK_CGELQF_MN = INT( CDUM(1) )
 494 *
 495             CALL CUNGBR( 'P', M, N, M, CDUM(1), M, CDUM(1), CDUM(1),
 496      $                   -1, IERR )
 497             LWORK_CUNGBR_P_MN = INT( CDUM(1) )
 498 *
 499             CALL CUNGBR( 'P', N, N, M, CDUM(1), N, CDUM(1), CDUM(1),
 500      $                   -1, IERR )
 501             LWORK_CUNGBR_P_NN = INT( CDUM(1) )
 502 *
 503             CALL CUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
 504      $                   -1, IERR )
 505             LWORK_CUNGBR_Q_MM = INT( CDUM(1) )
 506 *
 507             CALL CUNGLQ( M, N, M, CDUM(1), M, CDUM(1), CDUM(1),
 508      $                   -1, IERR )
 509             LWORK_CUNGLQ_MN = INT( CDUM(1) )
 510 *
 511             CALL CUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1),
 512      $                   -1, IERR )
 513             LWORK_CUNGLQ_NN = INT( CDUM(1) )
 514 *
 515             CALL CUNMBR( 'P', 'R', 'C', M, M, M, CDUM(1), M, CDUM(1),
 516      $                   CDUM(1), M, CDUM(1), -1, IERR )
 517             LWORK_CUNMBR_PRC_MM = INT( CDUM(1) )
 518 *
 519             CALL CUNMBR( 'P', 'R', 'C', M, N, M, CDUM(1), M, CDUM(1),
 520      $                   CDUM(1), M, CDUM(1), -1, IERR )
 521             LWORK_CUNMBR_PRC_MN = INT( CDUM(1) )
 522 *
 523             CALL CUNMBR( 'P', 'R', 'C', N, N, M, CDUM(1), N, CDUM(1),
 524      $                   CDUM(1), N, CDUM(1), -1, IERR )
 525             LWORK_CUNMBR_PRC_NN = INT( CDUM(1) )
 526 *
 527             CALL CUNMBR( 'Q', 'L', 'N', M, M, M, CDUM(1), M, CDUM(1),
 528      $                   CDUM(1), M, CDUM(1), -1, IERR )
 529             LWORK_CUNMBR_QLN_MM = INT( CDUM(1) )
 530 *
 531             IF( N.GE.MNTHR1 ) THEN
 532                IF( WNTQN ) THEN
 533 *
 534 *                 Path 1t (N >> M, JOBZ='N')
 535 *
 536                   MAXWRK = M + LWORK_CGELQF_MN
 537                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CGEBRD_MM )
 538                   MINWRK = 3*M
 539                ELSE IF( WNTQO ) THEN
 540 *
 541 *                 Path 2t (N >> M, JOBZ='O')
 542 *
 543                   WRKBL = M + LWORK_CGELQF_MN
 544                   WRKBL = MAX( WRKBL,   M + LWORK_CUNGLQ_MN )
 545                   WRKBL = MAX( WRKBL, 2*M + LWORK_CGEBRD_MM )
 546                   WRKBL = MAX( WRKBL, 2*M + LWORK_CUNMBR_QLN_MM )
 547                   WRKBL = MAX( WRKBL, 2*M + LWORK_CUNMBR_PRC_MM )
 548                   MAXWRK = M*N + M*M + WRKBL
 549                   MINWRK = 2*M*M + 3*M
 550                ELSE IF( WNTQS ) THEN
 551 *
 552 *                 Path 3t (N >> M, JOBZ='S')
 553 *
 554                   WRKBL = M + LWORK_CGELQF_MN
 555                   WRKBL = MAX( WRKBL,   M + LWORK_CUNGLQ_MN )
 556                   WRKBL = MAX( WRKBL, 2*M + LWORK_CGEBRD_MM )
 557                   WRKBL = MAX( WRKBL, 2*M + LWORK_CUNMBR_QLN_MM )
 558                   WRKBL = MAX( WRKBL, 2*M + LWORK_CUNMBR_PRC_MM )
 559                   MAXWRK = M*M + WRKBL
 560                   MINWRK = M*M + 3*M
 561                ELSE IF( WNTQA ) THEN
 562 *
 563 *                 Path 4t (N >> M, JOBZ='A')
 564 *
 565                   WRKBL = M + LWORK_CGELQF_MN
 566                   WRKBL = MAX( WRKBL,   M + LWORK_CUNGLQ_NN )
 567                   WRKBL = MAX( WRKBL, 2*M + LWORK_CGEBRD_MM )
 568                   WRKBL = MAX( WRKBL, 2*M + LWORK_CUNMBR_QLN_MM )
 569                   WRKBL = MAX( WRKBL, 2*M + LWORK_CUNMBR_PRC_MM )
 570                   MAXWRK = M*M + WRKBL
 571                   MINWRK = M*M + MAX( 3*M, M + N )
 572                END IF
 573             ELSE IF( N.GE.MNTHR2 ) THEN
 574 *
 575 *              Path 5t (N >> M, but not as much as MNTHR1)
 576 *
 577                MAXWRK = 2*M + LWORK_CGEBRD_MN
 578                MINWRK = 2*M + N
 579                IF( WNTQO ) THEN
 580 *                 Path 5to (N >> M, JOBZ='O')
 581                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR_Q_MM )
 582                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR_P_MN )
 583                   MAXWRK = MAXWRK + M*N
 584                   MINWRK = MINWRK + M*M
 585                ELSE IF( WNTQS ) THEN
 586 *                 Path 5ts (N >> M, JOBZ='S')
 587                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR_Q_MM )
 588                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR_P_MN )
 589                ELSE IF( WNTQA ) THEN
 590 *                 Path 5ta (N >> M, JOBZ='A')
 591                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR_Q_MM )
 592                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR_P_NN )
 593                END IF
 594             ELSE
 595 *
 596 *              Path 6t (N > M, but not much larger)
 597 *
 598                MAXWRK = 2*M + LWORK_CGEBRD_MN
 599                MINWRK = 2*M + N
 600                IF( WNTQO ) THEN
 601 *                 Path 6to (N > M, JOBZ='O')
 602                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR_QLN_MM )
 603                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR_PRC_MN )
 604                   MAXWRK = MAXWRK + M*N
 605                   MINWRK = MINWRK + M*M
 606                ELSE IF( WNTQS ) THEN
 607 *                 Path 6ts (N > M, JOBZ='S')
 608                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR_QLN_MM )
 609                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR_PRC_MN )
 610                ELSE IF( WNTQA ) THEN
 611 *                 Path 6ta (N > M, JOBZ='A')
 612                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR_QLN_MM )
 613                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR_PRC_NN )
 614                END IF
 615             END IF
 616          END IF
 617          MAXWRK = MAX( MAXWRK, MINWRK )
 618       END IF
 619       IF( INFO.EQ.0 ) THEN
 620          WORK( 1 ) = MAXWRK
 621          IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
 622             INFO = -12
 623          END IF
 624       END IF
 625 *
 626       IF( INFO.NE.0 ) THEN
 627          CALL XERBLA( 'CGESDD', -INFO )
 628          RETURN
 629       ELSE IF( LQUERY ) THEN
 630          RETURN
 631       END IF
 632 *
 633 *     Quick return if possible
 634 *
 635       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 636          RETURN
 637       END IF
 638 *
 639 *     Get machine constants
 640 *
 641       EPS = SLAMCH( 'P' )
 642       SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
 643       BIGNUM = ONE / SMLNUM
 644 *
 645 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 646 *
 647       ANRM = CLANGE( 'M', M, N, A, LDA, DUM )
 648       ISCL = 0
 649       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 650          ISCL = 1
 651          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
 652       ELSE IF( ANRM.GT.BIGNUM ) THEN
 653          ISCL = 1
 654          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
 655       END IF
 656 *
 657       IF( M.GE.N ) THEN
 658 *
 659 *        A has at least as many rows as columns. If A has sufficiently
 660 *        more rows than columns, first reduce using the QR
 661 *        decomposition (if sufficient workspace available)
 662 *
 663          IF( M.GE.MNTHR1 ) THEN
 664 *
 665             IF( WNTQN ) THEN
 666 *
 667 *              Path 1 (M >> N, JOBZ='N')
 668 *              No singular vectors to be computed
 669 *
 670                ITAU = 1
 671                NWORK = ITAU + N
 672 *
 673 *              Compute A=Q*R
 674 *              CWorkspace: need   N [tau] + N    [work]
 675 *              CWorkspace: prefer N [tau] + N*NB [work]
 676 *              RWorkspace: need   0
 677 *
 678                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 679      $                      LWORK-NWORK+1, IERR )
 680 *
 681 *              Zero out below R
 682 *
 683                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
 684      $                      LDA )
 685                IE = 1
 686                ITAUQ = 1
 687                ITAUP = ITAUQ + N
 688                NWORK = ITAUP + N
 689 *
 690 *              Bidiagonalize R in A
 691 *              CWorkspace: need   2*N [tauq, taup] + N      [work]
 692 *              CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work]
 693 *              RWorkspace: need   N [e]
 694 *
 695                CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
 696      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 697      $                      IERR )
 698                NRWORK = IE + N
 699 *
 700 *              Perform bidiagonal SVD, compute singular values only
 701 *              CWorkspace: need   0
 702 *              RWorkspace: need   N [e] + BDSPAC
 703 *
 704                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
 705      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
 706 *
 707             ELSE IF( WNTQO ) THEN
 708 *
 709 *              Path 2 (M >> N, JOBZ='O')
 710 *              N left singular vectors to be overwritten on A and
 711 *              N right singular vectors to be computed in VT
 712 *
 713                IU = 1
 714 *
 715 *              WORK(IU) is N by N
 716 *
 717                LDWRKU = N
 718                IR = IU + LDWRKU*N
 719                IF( LWORK .GE. M*N + N*N + 3*N ) THEN
 720 *
 721 *                 WORK(IR) is M by N
 722 *
 723                   LDWRKR = M
 724                ELSE
 725                   LDWRKR = ( LWORK - N*N - 3*N ) / N
 726                END IF
 727                ITAU = IR + LDWRKR*N
 728                NWORK = ITAU + N
 729 *
 730 *              Compute A=Q*R
 731 *              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work]
 732 *              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
 733 *              RWorkspace: need   0
 734 *
 735                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 736      $                      LWORK-NWORK+1, IERR )
 737 *
 738 *              Copy R to WORK( IR ), zeroing out below it
 739 *
 740                CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
 741                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
 742      $                      LDWRKR )
 743 *
 744 *              Generate Q in A
 745 *              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work]
 746 *              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
 747 *              RWorkspace: need   0
 748 *
 749                CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
 750      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 751                IE = 1
 752                ITAUQ = ITAU
 753                ITAUP = ITAUQ + N
 754                NWORK = ITAUP + N
 755 *
 756 *              Bidiagonalize R in WORK(IR)
 757 *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N      [work]
 758 *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
 759 *              RWorkspace: need   N [e]
 760 *
 761                CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
 762      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
 763      $                      LWORK-NWORK+1, IERR )
 764 *
 765 *              Perform bidiagonal SVD, computing left singular vectors
 766 *              of R in WORK(IRU) and computing right singular vectors
 767 *              of R in WORK(IRVT)
 768 *              CWorkspace: need   0
 769 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 770 *
 771                IRU = IE + N
 772                IRVT = IRU + N*N
 773                NRWORK = IRVT + N*N
 774                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
 775      $                      N, RWORK( IRVT ), N, DUM, IDUM,
 776      $                      RWORK( NRWORK ), IWORK, INFO )
 777 *
 778 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 779 *              Overwrite WORK(IU) by the left singular vectors of R
 780 *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work]
 781 *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
 782 *              RWorkspace: need   0
 783 *
 784                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
 785      $                      LDWRKU )
 786                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
 787      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
 788      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 789 *
 790 *              Copy real matrix RWORK(IRVT) to complex matrix VT
 791 *              Overwrite VT by the right singular vectors of R
 792 *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work]
 793 *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
 794 *              RWorkspace: need   0
 795 *
 796                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
 797                CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
 798      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 799      $                      LWORK-NWORK+1, IERR )
 800 *
 801 *              Multiply Q in A by left singular vectors of R in
 802 *              WORK(IU), storing result in WORK(IR) and copying to A
 803 *              CWorkspace: need   N*N [U] + N*N [R]
 804 *              CWorkspace: prefer N*N [U] + M*N [R]
 805 *              RWorkspace: need   0
 806 *
 807                DO 10 I = 1, M, LDWRKR
 808                   CHUNK = MIN( M-I+1, LDWRKR )
 809                   CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
 810      $                        LDA, WORK( IU ), LDWRKU, CZERO,
 811      $                        WORK( IR ), LDWRKR )
 812                   CALL CLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
 813      $                         A( I, 1 ), LDA )
 814    10          CONTINUE
 815 *
 816             ELSE IF( WNTQS ) THEN
 817 *
 818 *              Path 3 (M >> N, JOBZ='S')
 819 *              N left singular vectors to be computed in U and
 820 *              N right singular vectors to be computed in VT
 821 *
 822                IR = 1
 823 *
 824 *              WORK(IR) is N by N
 825 *
 826                LDWRKR = N
 827                ITAU = IR + LDWRKR*N
 828                NWORK = ITAU + N
 829 *
 830 *              Compute A=Q*R
 831 *              CWorkspace: need   N*N [R] + N [tau] + N    [work]
 832 *              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
 833 *              RWorkspace: need   0
 834 *
 835                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 836      $                      LWORK-NWORK+1, IERR )
 837 *
 838 *              Copy R to WORK(IR), zeroing out below it
 839 *
 840                CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
 841                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
 842      $                      LDWRKR )
 843 *
 844 *              Generate Q in A
 845 *              CWorkspace: need   N*N [R] + N [tau] + N    [work]
 846 *              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
 847 *              RWorkspace: need   0
 848 *
 849                CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
 850      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 851                IE = 1
 852                ITAUQ = ITAU
 853                ITAUP = ITAUQ + N
 854                NWORK = ITAUP + N
 855 *
 856 *              Bidiagonalize R in WORK(IR)
 857 *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N      [work]
 858 *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
 859 *              RWorkspace: need   N [e]
 860 *
 861                CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
 862      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
 863      $                      LWORK-NWORK+1, IERR )
 864 *
 865 *              Perform bidiagonal SVD, computing left singular vectors
 866 *              of bidiagonal matrix in RWORK(IRU) and computing right
 867 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 868 *              CWorkspace: need   0
 869 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 870 *
 871                IRU = IE + N
 872                IRVT = IRU + N*N
 873                NRWORK = IRVT + N*N
 874                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
 875      $                      N, RWORK( IRVT ), N, DUM, IDUM,
 876      $                      RWORK( NRWORK ), IWORK, INFO )
 877 *
 878 *              Copy real matrix RWORK(IRU) to complex matrix U
 879 *              Overwrite U by left singular vectors of R
 880 *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work]
 881 *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
 882 *              RWorkspace: need   0
 883 *
 884                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
 885                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
 886      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
 887      $                      LWORK-NWORK+1, IERR )
 888 *
 889 *              Copy real matrix RWORK(IRVT) to complex matrix VT
 890 *              Overwrite VT by right singular vectors of R
 891 *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work]
 892 *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
 893 *              RWorkspace: need   0
 894 *
 895                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
 896                CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
 897      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 898      $                      LWORK-NWORK+1, IERR )
 899 *
 900 *              Multiply Q in A by left singular vectors of R in
 901 *              WORK(IR), storing result in U
 902 *              CWorkspace: need   N*N [R]
 903 *              RWorkspace: need   0
 904 *
 905                CALL CLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
 906                CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
 907      $                     LDWRKR, CZERO, U, LDU )
 908 *
 909             ELSE IF( WNTQA ) THEN
 910 *
 911 *              Path 4 (M >> N, JOBZ='A')
 912 *              M left singular vectors to be computed in U and
 913 *              N right singular vectors to be computed in VT
 914 *
 915                IU = 1
 916 *
 917 *              WORK(IU) is N by N
 918 *
 919                LDWRKU = N
 920                ITAU = IU + LDWRKU*N
 921                NWORK = ITAU + N
 922 *
 923 *              Compute A=Q*R, copying result to U
 924 *              CWorkspace: need   N*N [U] + N [tau] + N    [work]
 925 *              CWorkspace: prefer N*N [U] + N [tau] + N*NB [work]
 926 *              RWorkspace: need   0
 927 *
 928                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 929      $                      LWORK-NWORK+1, IERR )
 930                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
 931 *
 932 *              Generate Q in U
 933 *              CWorkspace: need   N*N [U] + N [tau] + M    [work]
 934 *              CWorkspace: prefer N*N [U] + N [tau] + M*NB [work]
 935 *              RWorkspace: need   0
 936 *
 937                CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
 938      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 939 *
 940 *              Produce R in A, zeroing out below it
 941 *
 942                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
 943      $                      LDA )
 944                IE = 1
 945                ITAUQ = ITAU
 946                ITAUP = ITAUQ + N
 947                NWORK = ITAUP + N
 948 *
 949 *              Bidiagonalize R in A
 950 *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N      [work]
 951 *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work]
 952 *              RWorkspace: need   N [e]
 953 *
 954                CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
 955      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 956      $                      IERR )
 957                IRU = IE + N
 958                IRVT = IRU + N*N
 959                NRWORK = IRVT + N*N
 960 *
 961 *              Perform bidiagonal SVD, computing left singular vectors
 962 *              of bidiagonal matrix in RWORK(IRU) and computing right
 963 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 964 *              CWorkspace: need   0
 965 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 966 *
 967                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
 968      $                      N, RWORK( IRVT ), N, DUM, IDUM,
 969      $                      RWORK( NRWORK ), IWORK, INFO )
 970 *
 971 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 972 *              Overwrite WORK(IU) by left singular vectors of R
 973 *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work]
 974 *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
 975 *              RWorkspace: need   0
 976 *
 977                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
 978      $                      LDWRKU )
 979                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
 980      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
 981      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 982 *
 983 *              Copy real matrix RWORK(IRVT) to complex matrix VT
 984 *              Overwrite VT by right singular vectors of R
 985 *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work]
 986 *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
 987 *              RWorkspace: need   0
 988 *
 989                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
 990                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
 991      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 992      $                      LWORK-NWORK+1, IERR )
 993 *
 994 *              Multiply Q in U by left singular vectors of R in
 995 *              WORK(IU), storing result in A
 996 *              CWorkspace: need   N*N [U]
 997 *              RWorkspace: need   0
 998 *
 999                CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
1000      $                     LDWRKU, CZERO, A, LDA )
1001 *
1002 *              Copy left singular vectors of A from A to U
1003 *
1004                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1005 *
1006             END IF
1007 *
1008          ELSE IF( M.GE.MNTHR2 ) THEN
1009 *
1010 *           MNTHR2 <= M < MNTHR1
1011 *
1012 *           Path 5 (M >> N, but not as much as MNTHR1)
1013 *           Reduce to bidiagonal form without QR decomposition, use
1014 *           CUNGBR and matrix multiplication to compute singular vectors
1015 *
1016             IE = 1
1017             NRWORK = IE + N
1018             ITAUQ = 1
1019             ITAUP = ITAUQ + N
1020             NWORK = ITAUP + N
1021 *
1022 *           Bidiagonalize A
1023 *           CWorkspace: need   2*N [tauq, taup] + M        [work]
1024 *           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1025 *           RWorkspace: need   N [e]
1026 *
1027             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1028      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1029      $                   IERR )
1030             IF( WNTQN ) THEN
1031 *
1032 *              Path 5n (M >> N, JOBZ='N')
1033 *              Compute singular values only
1034 *              CWorkspace: need   0
1035 *              RWorkspace: need   N [e] + BDSPAC
1036 *
1037                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1,DUM,1,
1038      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1039             ELSE IF( WNTQO ) THEN
1040                IU = NWORK
1041                IRU = NRWORK
1042                IRVT = IRU + N*N
1043                NRWORK = IRVT + N*N
1044 *
1045 *              Path 5o (M >> N, JOBZ='O')
1046 *              Copy A to VT, generate P**H
1047 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1048 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1049 *              RWorkspace: need   0
1050 *
1051                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
1052                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1053      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1054 *
1055 *              Generate Q in A
1056 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1057 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1058 *              RWorkspace: need   0
1059 *
1060                CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
1061      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1062 *
1063                IF( LWORK .GE. M*N + 3*N ) THEN
1064 *
1065 *                 WORK( IU ) is M by N
1066 *
1067                   LDWRKU = M
1068                ELSE
1069 *
1070 *                 WORK(IU) is LDWRKU by N
1071 *
1072                   LDWRKU = ( LWORK - 3*N ) / N
1073                END IF
1074                NWORK = IU + LDWRKU*N
1075 *
1076 *              Perform bidiagonal SVD, computing left singular vectors
1077 *              of bidiagonal matrix in RWORK(IRU) and computing right
1078 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1079 *              CWorkspace: need   0
1080 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1081 *
1082                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1083      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1084      $                      RWORK( NRWORK ), IWORK, INFO )
1085 *
1086 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1087 *              storing the result in WORK(IU), copying to VT
1088 *              CWorkspace: need   2*N [tauq, taup] + N*N [U]
1089 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1090 *
1091                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
1092      $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )
1093                CALL CLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )
1094 *
1095 *              Multiply Q in A by real matrix RWORK(IRU), storing the
1096 *              result in WORK(IU), copying to A
1097 *              CWorkspace: need   2*N [tauq, taup] + N*N [U]
1098 *              CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1099 *              RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork]
1100 *              RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1101 *
1102                NRWORK = IRVT
1103                DO 20 I = 1, M, LDWRKU
1104                   CHUNK = MIN( M-I+1, LDWRKU )
1105                   CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),
1106      $                         N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )
1107                   CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
1108      $                         A( I, 1 ), LDA )
1109    20          CONTINUE
1110 *
1111             ELSE IF( WNTQS ) THEN
1112 *
1113 *              Path 5s (M >> N, JOBZ='S')
1114 *              Copy A to VT, generate P**H
1115 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1116 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1117 *              RWorkspace: need   0
1118 *
1119                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
1120                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1121      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1122 *
1123 *              Copy A to U, generate Q
1124 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1125 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1126 *              RWorkspace: need   0
1127 *
1128                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1129                CALL CUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
1130      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1131 *
1132 *              Perform bidiagonal SVD, computing left singular vectors
1133 *              of bidiagonal matrix in RWORK(IRU) and computing right
1134 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1135 *              CWorkspace: need   0
1136 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1137 *
1138                IRU = NRWORK
1139                IRVT = IRU + N*N
1140                NRWORK = IRVT + N*N
1141                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1142      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1143      $                      RWORK( NRWORK ), IWORK, INFO )
1144 *
1145 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1146 *              storing the result in A, copying to VT
1147 *              CWorkspace: need   0
1148 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1149 *
1150                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
1151      $                      RWORK( NRWORK ) )
1152                CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
1153 *
1154 *              Multiply Q in U by real matrix RWORK(IRU), storing the
1155 *              result in A, copying to U
1156 *              CWorkspace: need   0
1157 *              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1158 *
1159                NRWORK = IRVT
1160                CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
1161      $                      RWORK( NRWORK ) )
1162                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1163             ELSE
1164 *
1165 *              Path 5a (M >> N, JOBZ='A')
1166 *              Copy A to VT, generate P**H
1167 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1168 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1169 *              RWorkspace: need   0
1170 *
1171                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
1172                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1173      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1174 *
1175 *              Copy A to U, generate Q
1176 *              CWorkspace: need   2*N [tauq, taup] + M    [work]
1177 *              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1178 *              RWorkspace: need   0
1179 *
1180                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1181                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1182      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1183 *
1184 *              Perform bidiagonal SVD, computing left singular vectors
1185 *              of bidiagonal matrix in RWORK(IRU) and computing right
1186 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1187 *              CWorkspace: need   0
1188 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1189 *
1190                IRU = NRWORK
1191                IRVT = IRU + N*N
1192                NRWORK = IRVT + N*N
1193                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1194      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1195      $                      RWORK( NRWORK ), IWORK, INFO )
1196 *
1197 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1198 *              storing the result in A, copying to VT
1199 *              CWorkspace: need   0
1200 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1201 *
1202                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
1203      $                      RWORK( NRWORK ) )
1204                CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
1205 *
1206 *              Multiply Q in U by real matrix RWORK(IRU), storing the
1207 *              result in A, copying to U
1208 *              CWorkspace: need   0
1209 *              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1210 *
1211                NRWORK = IRVT
1212                CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
1213      $                      RWORK( NRWORK ) )
1214                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1215             END IF
1216 *
1217          ELSE
1218 *
1219 *           M .LT. MNTHR2
1220 *
1221 *           Path 6 (M >= N, but not much larger)
1222 *           Reduce to bidiagonal form without QR decomposition
1223 *           Use CUNMBR to compute singular vectors
1224 *
1225             IE = 1
1226             NRWORK = IE + N
1227             ITAUQ = 1
1228             ITAUP = ITAUQ + N
1229             NWORK = ITAUP + N
1230 *
1231 *           Bidiagonalize A
1232 *           CWorkspace: need   2*N [tauq, taup] + M        [work]
1233 *           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1234 *           RWorkspace: need   N [e]
1235 *
1236             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1237      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1238      $                   IERR )
1239             IF( WNTQN ) THEN
1240 *
1241 *              Path 6n (M >= N, JOBZ='N')
1242 *              Compute singular values only
1243 *              CWorkspace: need   0
1244 *              RWorkspace: need   N [e] + BDSPAC
1245 *
1246                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
1247      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1248             ELSE IF( WNTQO ) THEN
1249                IU = NWORK
1250                IRU = NRWORK
1251                IRVT = IRU + N*N
1252                NRWORK = IRVT + N*N
1253                IF( LWORK .GE. M*N + 3*N ) THEN
1254 *
1255 *                 WORK( IU ) is M by N
1256 *
1257                   LDWRKU = M
1258                ELSE
1259 *
1260 *                 WORK( IU ) is LDWRKU by N
1261 *
1262                   LDWRKU = ( LWORK - 3*N ) / N
1263                END IF
1264                NWORK = IU + LDWRKU*N
1265 *
1266 *              Path 6o (M >= N, JOBZ='O')
1267 *              Perform bidiagonal SVD, computing left singular vectors
1268 *              of bidiagonal matrix in RWORK(IRU) and computing right
1269 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1270 *              CWorkspace: need   0
1271 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1272 *
1273                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1274      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1275      $                      RWORK( NRWORK ), IWORK, INFO )
1276 *
1277 *              Copy real matrix RWORK(IRVT) to complex matrix VT
1278 *              Overwrite VT by right singular vectors of A
1279 *              CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work]
1280 *              CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1281 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
1282 *
1283                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1284                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1285      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1286      $                      LWORK-NWORK+1, IERR )
1287 *
1288                IF( LWORK .GE. M*N + 3*N ) THEN
1289 *
1290 *                 Path 6o-fast
1291 *                 Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1292 *                 Overwrite WORK(IU) by left singular vectors of A, copying
1293 *                 to A
1294 *                 CWorkspace: need   2*N [tauq, taup] + M*N [U] + N    [work]
1295 *                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work]
1296 *                 RWorkspace: need   N [e] + N*N [RU]
1297 *
1298                   CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
1299      $                         LDWRKU )
1300                   CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
1301      $                         LDWRKU )
1302                   CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1303      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
1304      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1305                   CALL CLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
1306                ELSE
1307 *
1308 *                 Path 6o-slow
1309 *                 Generate Q in A
1310 *                 CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work]
1311 *                 CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1312 *                 RWorkspace: need   0
1313 *
1314                   CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
1315      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1316 *
1317 *                 Multiply Q in A by real matrix RWORK(IRU), storing the
1318 *                 result in WORK(IU), copying to A
1319 *                 CWorkspace: need   2*N [tauq, taup] + N*N [U]
1320 *                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1321 *                 RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork]
1322 *                 RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1323 *
1324                   NRWORK = IRVT
1325                   DO 30 I = 1, M, LDWRKU
1326                      CHUNK = MIN( M-I+1, LDWRKU )
1327                      CALL CLACRM( CHUNK, N, A( I, 1 ), LDA,
1328      $                            RWORK( IRU ), N, WORK( IU ), LDWRKU,
1329      $                            RWORK( NRWORK ) )
1330                      CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
1331      $                            A( I, 1 ), LDA )
1332    30             CONTINUE
1333                END IF
1334 *
1335             ELSE IF( WNTQS ) THEN
1336 *
1337 *              Path 6s (M >= N, JOBZ='S')
1338 *              Perform bidiagonal SVD, computing left singular vectors
1339 *              of bidiagonal matrix in RWORK(IRU) and computing right
1340 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1341 *              CWorkspace: need   0
1342 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1343 *
1344                IRU = NRWORK
1345                IRVT = IRU + N*N
1346                NRWORK = IRVT + N*N
1347                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1348      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1349      $                      RWORK( NRWORK ), IWORK, INFO )
1350 *
1351 *              Copy real matrix RWORK(IRU) to complex matrix U
1352 *              Overwrite U by left singular vectors of A
1353 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1354 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1355 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
1356 *
1357                CALL CLASET( 'F', M, N, CZERO, CZERO, U, LDU )
1358                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1359                CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1360      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1361      $                      LWORK-NWORK+1, IERR )
1362 *
1363 *              Copy real matrix RWORK(IRVT) to complex matrix VT
1364 *              Overwrite VT by right singular vectors of A
1365 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1366 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1367 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
1368 *
1369                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1370                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1371      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1372      $                      LWORK-NWORK+1, IERR )
1373             ELSE
1374 *
1375 *              Path 6a (M >= N, JOBZ='A')
1376 *              Perform bidiagonal SVD, computing left singular vectors
1377 *              of bidiagonal matrix in RWORK(IRU) and computing right
1378 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1379 *              CWorkspace: need   0
1380 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1381 *
1382                IRU = NRWORK
1383                IRVT = IRU + N*N
1384                NRWORK = IRVT + N*N
1385                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1386      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1387      $                      RWORK( NRWORK ), IWORK, INFO )
1388 *
1389 *              Set the right corner of U to identity matrix
1390 *
1391                CALL CLASET( 'F', M, M, CZERO, CZERO, U, LDU )
1392                IF( M.GT.N ) THEN
1393                   CALL CLASET( 'F', M-N, M-N, CZERO, CONE,
1394      $                         U( N+1, N+1 ), LDU )
1395                END IF
1396 *
1397 *              Copy real matrix RWORK(IRU) to complex matrix U
1398 *              Overwrite U by left singular vectors of A
1399 *              CWorkspace: need   2*N [tauq, taup] + M    [work]
1400 *              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1401 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
1402 *
1403                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1404                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1405      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1406      $                      LWORK-NWORK+1, IERR )
1407 *
1408 *              Copy real matrix RWORK(IRVT) to complex matrix VT
1409 *              Overwrite VT by right singular vectors of A
1410 *              CWorkspace: need   2*N [tauq, taup] + N    [work]
1411 *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1412 *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
1413 *
1414                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1415                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1416      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1417      $                      LWORK-NWORK+1, IERR )
1418             END IF
1419 *
1420          END IF
1421 *
1422       ELSE
1423 *
1424 *        A has more columns than rows. If A has sufficiently more
1425 *        columns than rows, first reduce using the LQ decomposition (if
1426 *        sufficient workspace available)
1427 *
1428          IF( N.GE.MNTHR1 ) THEN
1429 *
1430             IF( WNTQN ) THEN
1431 *
1432 *              Path 1t (N >> M, JOBZ='N')
1433 *              No singular vectors to be computed
1434 *
1435                ITAU = 1
1436                NWORK = ITAU + M
1437 *
1438 *              Compute A=L*Q
1439 *              CWorkspace: need   M [tau] + M    [work]
1440 *              CWorkspace: prefer M [tau] + M*NB [work]
1441 *              RWorkspace: need   0
1442 *
1443                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1444      $                      LWORK-NWORK+1, IERR )
1445 *
1446 *              Zero out above L
1447 *
1448                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1449      $                      LDA )
1450                IE = 1
1451                ITAUQ = 1
1452                ITAUP = ITAUQ + M
1453                NWORK = ITAUP + M
1454 *
1455 *              Bidiagonalize L in A
1456 *              CWorkspace: need   2*M [tauq, taup] + M      [work]
1457 *              CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work]
1458 *              RWorkspace: need   M [e]
1459 *
1460                CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1461      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1462      $                      IERR )
1463                NRWORK = IE + M
1464 *
1465 *              Perform bidiagonal SVD, compute singular values only
1466 *              CWorkspace: need   0
1467 *              RWorkspace: need   M [e] + BDSPAC
1468 *
1469                CALL SBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
1470      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1471 *
1472             ELSE IF( WNTQO ) THEN
1473 *
1474 *              Path 2t (N >> M, JOBZ='O')
1475 *              M right singular vectors to be overwritten on A and
1476 *              M left singular vectors to be computed in U
1477 *
1478                IVT = 1
1479                LDWKVT = M
1480 *
1481 *              WORK(IVT) is M by M
1482 *
1483                IL = IVT + LDWKVT*M
1484                IF( LWORK .GE. M*N + M*M + 3*M ) THEN
1485 *
1486 *                 WORK(IL) M by N
1487 *
1488                   LDWRKL = M
1489                   CHUNK = N
1490                ELSE
1491 *
1492 *                 WORK(IL) is M by CHUNK
1493 *
1494                   LDWRKL = M
1495                   CHUNK = ( LWORK - M*M - 3*M ) / M
1496                END IF
1497                ITAU = IL + LDWRKL*CHUNK
1498                NWORK = ITAU + M
1499 *
1500 *              Compute A=L*Q
1501 *              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
1502 *              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1503 *              RWorkspace: need   0
1504 *
1505                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1506      $                      LWORK-NWORK+1, IERR )
1507 *
1508 *              Copy L to WORK(IL), zeroing about above it
1509 *
1510                CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1511                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
1512      $                      WORK( IL+LDWRKL ), LDWRKL )
1513 *
1514 *              Generate Q in A
1515 *              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
1516 *              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1517 *              RWorkspace: need   0
1518 *
1519                CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1520      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1521                IE = 1
1522                ITAUQ = ITAU
1523                ITAUP = ITAUQ + M
1524                NWORK = ITAUP + M
1525 *
1526 *              Bidiagonalize L in WORK(IL)
1527 *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M      [work]
1528 *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1529 *              RWorkspace: need   M [e]
1530 *
1531                CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1532      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1533      $                      LWORK-NWORK+1, IERR )
1534 *
1535 *              Perform bidiagonal SVD, computing left singular vectors
1536 *              of bidiagonal matrix in RWORK(IRU) and computing right
1537 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1538 *              CWorkspace: need   0
1539 *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1540 *
1541                IRU = IE + M
1542                IRVT = IRU + M*M
1543                NRWORK = IRVT + M*M
1544                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1545      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1546      $                      RWORK( NRWORK ), IWORK, INFO )
1547 *
1548 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1549 *              Overwrite WORK(IU) by the left singular vectors of L
1550 *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work]
1551 *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1552 *              RWorkspace: need   0
1553 *
1554                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1555                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1556      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1557      $                      LWORK-NWORK+1, IERR )
1558 *
1559 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1560 *              Overwrite WORK(IVT) by the right singular vectors of L
1561 *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work]
1562 *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1563 *              RWorkspace: need   0
1564 *
1565                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1566      $                      LDWKVT )
1567                CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1568      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1569      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1570 *
1571 *              Multiply right singular vectors of L in WORK(IL) by Q
1572 *              in A, storing result in WORK(IL) and copying to A
1573 *              CWorkspace: need   M*M [VT] + M*M [L]
1574 *              CWorkspace: prefer M*M [VT] + M*N [L]
1575 *              RWorkspace: need   0
1576 *
1577                DO 40 I = 1, N, CHUNK
1578                   BLK = MIN( N-I+1, CHUNK )
1579                   CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,
1580      $                        A( 1, I ), LDA, CZERO, WORK( IL ),
1581      $                        LDWRKL )
1582                   CALL CLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1583      $                         A( 1, I ), LDA )
1584    40          CONTINUE
1585 *
1586             ELSE IF( WNTQS ) THEN
1587 *
1588 *              Path 3t (N >> M, JOBZ='S')
1589 *              M right singular vectors to be computed in VT and
1590 *              M left singular vectors to be computed in U
1591 *
1592                IL = 1
1593 *
1594 *              WORK(IL) is M by M
1595 *
1596                LDWRKL = M
1597                ITAU = IL + LDWRKL*M
1598                NWORK = ITAU + M
1599 *
1600 *              Compute A=L*Q
1601 *              CWorkspace: need   M*M [L] + M [tau] + M    [work]
1602 *              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1603 *              RWorkspace: need   0
1604 *
1605                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1606      $                      LWORK-NWORK+1, IERR )
1607 *
1608 *              Copy L to WORK(IL), zeroing out above it
1609 *
1610                CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1611                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
1612      $                      WORK( IL+LDWRKL ), LDWRKL )
1613 *
1614 *              Generate Q in A
1615 *              CWorkspace: need   M*M [L] + M [tau] + M    [work]
1616 *              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1617 *              RWorkspace: need   0
1618 *
1619                CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1620      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1621                IE = 1
1622                ITAUQ = ITAU
1623                ITAUP = ITAUQ + M
1624                NWORK = ITAUP + M
1625 *
1626 *              Bidiagonalize L in WORK(IL)
1627 *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M      [work]
1628 *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1629 *              RWorkspace: need   M [e]
1630 *
1631                CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1632      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1633      $                      LWORK-NWORK+1, IERR )
1634 *
1635 *              Perform bidiagonal SVD, computing left singular vectors
1636 *              of bidiagonal matrix in RWORK(IRU) and computing right
1637 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1638 *              CWorkspace: need   0
1639 *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1640 *
1641                IRU = IE + M
1642                IRVT = IRU + M*M
1643                NRWORK = IRVT + M*M
1644                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1645      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1646      $                      RWORK( NRWORK ), IWORK, INFO )
1647 *
1648 *              Copy real matrix RWORK(IRU) to complex matrix U
1649 *              Overwrite U by left singular vectors of L
1650 *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work]
1651 *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1652 *              RWorkspace: need   0
1653 *
1654                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1655                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1656      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1657      $                      LWORK-NWORK+1, IERR )
1658 *
1659 *              Copy real matrix RWORK(IRVT) to complex matrix VT
1660 *              Overwrite VT by left singular vectors of L
1661 *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work]
1662 *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1663 *              RWorkspace: need   0
1664 *
1665                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1666                CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1667      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1668      $                      LWORK-NWORK+1, IERR )
1669 *
1670 *              Copy VT to WORK(IL), multiply right singular vectors of L
1671 *              in WORK(IL) by Q in A, storing result in VT
1672 *              CWorkspace: need   M*M [L]
1673 *              RWorkspace: need   0
1674 *
1675                CALL CLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1676                CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
1677      $                     A, LDA, CZERO, VT, LDVT )
1678 *
1679             ELSE IF( WNTQA ) THEN
1680 *
1681 *              Path 4t (N >> M, JOBZ='A')
1682 *              N right singular vectors to be computed in VT and
1683 *              M left singular vectors to be computed in U
1684 *
1685                IVT = 1
1686 *
1687 *              WORK(IVT) is M by M
1688 *
1689                LDWKVT = M
1690                ITAU = IVT + LDWKVT*M
1691                NWORK = ITAU + M
1692 *
1693 *              Compute A=L*Q, copying result to VT
1694 *              CWorkspace: need   M*M [VT] + M [tau] + M    [work]
1695 *              CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work]
1696 *              RWorkspace: need   0
1697 *
1698                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1699      $                      LWORK-NWORK+1, IERR )
1700                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1701 *
1702 *              Generate Q in VT
1703 *              CWorkspace: need   M*M [VT] + M [tau] + N    [work]
1704 *              CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work]
1705 *              RWorkspace: need   0
1706 *
1707                CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1708      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1709 *
1710 *              Produce L in A, zeroing out above it
1711 *
1712                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1713      $                      LDA )
1714                IE = 1
1715                ITAUQ = ITAU
1716                ITAUP = ITAUQ + M
1717                NWORK = ITAUP + M
1718 *
1719 *              Bidiagonalize L in A
1720 *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M      [work]
1721 *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work]
1722 *              RWorkspace: need   M [e]
1723 *
1724                CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1725      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1726      $                      IERR )
1727 *
1728 *              Perform bidiagonal SVD, computing left singular vectors
1729 *              of bidiagonal matrix in RWORK(IRU) and computing right
1730 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1731 *              CWorkspace: need   0
1732 *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1733 *
1734                IRU = IE + M
1735                IRVT = IRU + M*M
1736                NRWORK = IRVT + M*M
1737                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1738      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1739      $                      RWORK( NRWORK ), IWORK, INFO )
1740 *
1741 *              Copy real matrix RWORK(IRU) to complex matrix U
1742 *              Overwrite U by left singular vectors of L
1743 *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work]
1744 *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1745 *              RWorkspace: need   0
1746 *
1747                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1748                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1749      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1750      $                      LWORK-NWORK+1, IERR )
1751 *
1752 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1753 *              Overwrite WORK(IVT) by right singular vectors of L
1754 *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work]
1755 *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1756 *              RWorkspace: need   0
1757 *
1758                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1759      $                      LDWKVT )
1760                CALL CUNMBR( 'P', 'R', 'C', M, M, M, A, LDA,
1761      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1762      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1763 *
1764 *              Multiply right singular vectors of L in WORK(IVT) by
1765 *              Q in VT, storing result in A
1766 *              CWorkspace: need   M*M [VT]
1767 *              RWorkspace: need   0
1768 *
1769                CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,
1770      $                     VT, LDVT, CZERO, A, LDA )
1771 *
1772 *              Copy right singular vectors of A from A to VT
1773 *
1774                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1775 *
1776             END IF
1777 *
1778          ELSE IF( N.GE.MNTHR2 ) THEN
1779 *
1780 *           MNTHR2 <= N < MNTHR1
1781 *
1782 *           Path 5t (N >> M, but not as much as MNTHR1)
1783 *           Reduce to bidiagonal form without QR decomposition, use
1784 *           CUNGBR and matrix multiplication to compute singular vectors
1785 *
1786             IE = 1
1787             NRWORK = IE + M
1788             ITAUQ = 1
1789             ITAUP = ITAUQ + M
1790             NWORK = ITAUP + M
1791 *
1792 *           Bidiagonalize A
1793 *           CWorkspace: need   2*M [tauq, taup] + N        [work]
1794 *           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
1795 *           RWorkspace: need   M [e]
1796 *
1797             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1798      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1799      $                   IERR )
1800 *
1801             IF( WNTQN ) THEN
1802 *
1803 *              Path 5tn (N >> M, JOBZ='N')
1804 *              Compute singular values only
1805 *              CWorkspace: need   0
1806 *              RWorkspace: need   M [e] + BDSPAC
1807 *
1808                CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
1809      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1810             ELSE IF( WNTQO ) THEN
1811                IRVT = NRWORK
1812                IRU = IRVT + M*M
1813                NRWORK = IRU + M*M
1814                IVT = NWORK
1815 *
1816 *              Path 5to (N >> M, JOBZ='O')
1817 *              Copy A to U, generate Q
1818 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
1819 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1820 *              RWorkspace: need   0
1821 *
1822                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1823                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1824      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1825 *
1826 *              Generate P**H in A
1827 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
1828 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1829 *              RWorkspace: need   0
1830 *
1831                CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1832      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1833 *
1834                LDWKVT = M
1835                IF( LWORK .GE. M*N + 3*M ) THEN
1836 *
1837 *                 WORK( IVT ) is M by N
1838 *
1839                   NWORK = IVT + LDWKVT*N
1840                   CHUNK = N
1841                ELSE
1842 *
1843 *                 WORK( IVT ) is M by CHUNK
1844 *
1845                   CHUNK = ( LWORK - 3*M ) / M
1846                   NWORK = IVT + LDWKVT*CHUNK
1847                END IF
1848 *
1849 *              Perform bidiagonal SVD, computing left singular vectors
1850 *              of bidiagonal matrix in RWORK(IRU) and computing right
1851 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1852 *              CWorkspace: need   0
1853 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1854 *
1855                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1856      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1857      $                      RWORK( NRWORK ), IWORK, INFO )
1858 *
1859 *              Multiply Q in U by real matrix RWORK(IRVT)
1860 *              storing the result in WORK(IVT), copying to U
1861 *              CWorkspace: need   2*M [tauq, taup] + M*M [VT]
1862 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1863 *
1864                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
1865      $                      LDWKVT, RWORK( NRWORK ) )
1866                CALL CLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU )
1867 *
1868 *              Multiply RWORK(IRVT) by P**H in A, storing the
1869 *              result in WORK(IVT), copying to A
1870 *              CWorkspace: need   2*M [tauq, taup] + M*M [VT]
1871 *              CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
1872 *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork]
1873 *              RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1874 *
1875                NRWORK = IRU
1876                DO 50 I = 1, N, CHUNK
1877                   BLK = MIN( N-I+1, CHUNK )
1878                   CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA,
1879      $                         WORK( IVT ), LDWKVT, RWORK( NRWORK ) )
1880                   CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
1881      $                         A( 1, I ), LDA )
1882    50          CONTINUE
1883             ELSE IF( WNTQS ) THEN
1884 *
1885 *              Path 5ts (N >> M, JOBZ='S')
1886 *              Copy A to U, generate Q
1887 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
1888 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1889 *              RWorkspace: need   0
1890 *
1891                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1892                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1893      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1894 *
1895 *              Copy A to VT, generate P**H
1896 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
1897 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1898 *              RWorkspace: need   0
1899 *
1900                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1901                CALL CUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
1902      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1903 *
1904 *              Perform bidiagonal SVD, computing left singular vectors
1905 *              of bidiagonal matrix in RWORK(IRU) and computing right
1906 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1907 *              CWorkspace: need   0
1908 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1909 *
1910                IRVT = NRWORK
1911                IRU = IRVT + M*M
1912                NRWORK = IRU + M*M
1913                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1914      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1915      $                      RWORK( NRWORK ), IWORK, INFO )
1916 *
1917 *              Multiply Q in U by real matrix RWORK(IRU), storing the
1918 *              result in A, copying to U
1919 *              CWorkspace: need   0
1920 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1921 *
1922                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1923      $                      RWORK( NRWORK ) )
1924                CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
1925 *
1926 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1927 *              storing the result in A, copying to VT
1928 *              CWorkspace: need   0
1929 *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1930 *
1931                NRWORK = IRU
1932                CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1933      $                      RWORK( NRWORK ) )
1934                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1935             ELSE
1936 *
1937 *              Path 5ta (N >> M, JOBZ='A')
1938 *              Copy A to U, generate Q
1939 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
1940 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1941 *              RWorkspace: need   0
1942 *
1943                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1944                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1945      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1946 *
1947 *              Copy A to VT, generate P**H
1948 *              CWorkspace: need   2*M [tauq, taup] + N    [work]
1949 *              CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
1950 *              RWorkspace: need   0
1951 *
1952                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1953                CALL CUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
1954      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1955 *
1956 *              Perform bidiagonal SVD, computing left singular vectors
1957 *              of bidiagonal matrix in RWORK(IRU) and computing right
1958 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1959 *              CWorkspace: need   0
1960 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1961 *
1962                IRVT = NRWORK
1963                IRU = IRVT + M*M
1964                NRWORK = IRU + M*M
1965                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1966      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1967      $                      RWORK( NRWORK ), IWORK, INFO )
1968 *
1969 *              Multiply Q in U by real matrix RWORK(IRU), storing the
1970 *              result in A, copying to U
1971 *              CWorkspace: need   0
1972 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1973 *
1974                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1975      $                      RWORK( NRWORK ) )
1976                CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
1977 *
1978 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1979 *              storing the result in A, copying to VT
1980 *              CWorkspace: need   0
1981 *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1982 *
1983                NRWORK = IRU
1984                CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1985      $                      RWORK( NRWORK ) )
1986                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1987             END IF
1988 *
1989          ELSE
1990 *
1991 *           N .LT. MNTHR2
1992 *
1993 *           Path 6t (N > M, but not much larger)
1994 *           Reduce to bidiagonal form without LQ decomposition
1995 *           Use CUNMBR to compute singular vectors
1996 *
1997             IE = 1
1998             NRWORK = IE + M
1999             ITAUQ = 1
2000             ITAUP = ITAUQ + M
2001             NWORK = ITAUP + M
2002 *
2003 *           Bidiagonalize A
2004 *           CWorkspace: need   2*M [tauq, taup] + N        [work]
2005 *           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
2006 *           RWorkspace: need   M [e]
2007 *
2008             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
2009      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
2010      $                   IERR )
2011             IF( WNTQN ) THEN
2012 *
2013 *              Path 6tn (N > M, JOBZ='N')
2014 *              Compute singular values only
2015 *              CWorkspace: need   0
2016 *              RWorkspace: need   M [e] + BDSPAC
2017 *
2018                CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
2019      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
2020             ELSE IF( WNTQO ) THEN
2021 *              Path 6to (N > M, JOBZ='O')
2022                LDWKVT = M
2023                IVT = NWORK
2024                IF( LWORK .GE. M*N + 3*M ) THEN
2025 *
2026 *                 WORK( IVT ) is M by N
2027 *
2028                   CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ),
2029      $                         LDWKVT )
2030                   NWORK = IVT + LDWKVT*N
2031                ELSE
2032 *
2033 *                 WORK( IVT ) is M by CHUNK
2034 *
2035                   CHUNK = ( LWORK - 3*M ) / M
2036                   NWORK = IVT + LDWKVT*CHUNK
2037                END IF
2038 *
2039 *              Perform bidiagonal SVD, computing left singular vectors
2040 *              of bidiagonal matrix in RWORK(IRU) and computing right
2041 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
2042 *              CWorkspace: need   0
2043 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2044 *
2045                IRVT = NRWORK
2046                IRU = IRVT + M*M
2047                NRWORK = IRU + M*M
2048                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
2049      $                      M, RWORK( IRVT ), M, DUM, IDUM,
2050      $                      RWORK( NRWORK ), IWORK, INFO )
2051 *
2052 *              Copy real matrix RWORK(IRU) to complex matrix U
2053 *              Overwrite U by left singular vectors of A
2054 *              CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work]
2055 *              CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2056 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
2057 *
2058                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
2059                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
2060      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
2061      $                      LWORK-NWORK+1, IERR )
2062 *
2063                IF( LWORK .GE. M*N + 3*M ) THEN
2064 *
2065 *                 Path 6to-fast
2066 *                 Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
2067 *                 Overwrite WORK(IVT) by right singular vectors of A,
2068 *                 copying to A
2069 *                 CWorkspace: need   2*M [tauq, taup] + M*N [VT] + M    [work]
2070 *                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work]
2071 *                 RWorkspace: need   M [e] + M*M [RVT]
2072 *
2073                   CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
2074      $                         LDWKVT )
2075                   CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
2076      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
2077      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
2078                   CALL CLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
2079                ELSE
2080 *
2081 *                 Path 6to-slow
2082 *                 Generate P**H in A
2083 *                 CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work]
2084 *                 CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2085 *                 RWorkspace: need   0
2086 *
2087                   CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
2088      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
2089 *
2090 *                 Multiply Q in A by real matrix RWORK(IRU), storing the
2091 *                 result in WORK(IU), copying to A
2092 *                 CWorkspace: need   2*M [tauq, taup] + M*M [VT]
2093 *                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
2094 *                 RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork]
2095 *                 RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
2096 *
2097                   NRWORK = IRU
2098                   DO 60 I = 1, N, CHUNK
2099                      BLK = MIN( N-I+1, CHUNK )
2100                      CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ),
2101      $                            LDA, WORK( IVT ), LDWKVT,
2102      $                            RWORK( NRWORK ) )
2103                      CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
2104      $                            A( 1, I ), LDA )
2105    60             CONTINUE
2106                END IF
2107             ELSE IF( WNTQS ) THEN
2108 *
2109 *              Path 6ts (N > M, JOBZ='S')
2110 *              Perform bidiagonal SVD, computing left singular vectors
2111 *              of bidiagonal matrix in RWORK(IRU) and computing right
2112 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
2113 *              CWorkspace: need   0
2114 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2115 *
2116                IRVT = NRWORK
2117                IRU = IRVT + M*M
2118                NRWORK = IRU + M*M
2119                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
2120      $                      M, RWORK( IRVT ), M, DUM, IDUM,
2121      $                      RWORK( NRWORK ), IWORK, INFO )
2122 *
2123 *              Copy real matrix RWORK(IRU) to complex matrix U
2124 *              Overwrite U by left singular vectors of A
2125 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
2126 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2127 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
2128 *
2129                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
2130                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
2131      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
2132      $                      LWORK-NWORK+1, IERR )
2133 *
2134 *              Copy real matrix RWORK(IRVT) to complex matrix VT
2135 *              Overwrite VT by right singular vectors of A
2136 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
2137 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2138 *              RWorkspace: need   M [e] + M*M [RVT]
2139 *
2140                CALL CLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
2141                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
2142                CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
2143      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
2144      $                      LWORK-NWORK+1, IERR )
2145             ELSE
2146 *
2147 *              Path 6ta (N > M, JOBZ='A')
2148 *              Perform bidiagonal SVD, computing left singular vectors
2149 *              of bidiagonal matrix in RWORK(IRU) and computing right
2150 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
2151 *              CWorkspace: need   0
2152 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2153 *
2154                IRVT = NRWORK
2155                IRU = IRVT + M*M
2156                NRWORK = IRU + M*M
2157 *
2158                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
2159      $                      M, RWORK( IRVT ), M, DUM, IDUM,
2160      $                      RWORK( NRWORK ), IWORK, INFO )
2161 *
2162 *              Copy real matrix RWORK(IRU) to complex matrix U
2163 *              Overwrite U by left singular vectors of A
2164 *              CWorkspace: need   2*M [tauq, taup] + M    [work]
2165 *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2166 *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
2167 *
2168                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
2169                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
2170      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
2171      $                      LWORK-NWORK+1, IERR )
2172 *
2173 *              Set all of VT to identity matrix
2174 *
2175                CALL CLASET( 'F', N, N, CZERO, CONE, VT, LDVT )
2176 *
2177 *              Copy real matrix RWORK(IRVT) to complex matrix VT
2178 *              Overwrite VT by right singular vectors of A
2179 *              CWorkspace: need   2*M [tauq, taup] + N    [work]
2180 *              CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
2181 *              RWorkspace: need   M [e] + M*M [RVT]
2182 *
2183                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
2184                CALL CUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,
2185      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
2186      $                      LWORK-NWORK+1, IERR )
2187             END IF
2188 *
2189          END IF
2190 *
2191       END IF
2192 *
2193 *     Undo scaling if necessary
2194 *
2195       IF( ISCL.EQ.1 ) THEN
2196          IF( ANRM.GT.BIGNUM )
2197      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
2198      $                   IERR )
2199          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
2200      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
2201      $                   RWORK( IE ), MINMN, IERR )
2202          IF( ANRM.LT.SMLNUM )
2203      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
2204      $                   IERR )
2205          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
2206      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
2207      $                   RWORK( IE ), MINMN, IERR )
2208       END IF
2209 *
2210 *     Return optimal workspace in WORK(1)
2211 *
2212       WORK( 1 ) = MAXWRK
2213 *
2214       RETURN
2215 *
2216 *     End of CGESDD
2217 *
2218       END