SRC/zgesdd.f

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