zgesdd.f (flens/lapack/interface/ref

   1       SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,

   2      $                   LWORK, RWORK, IWORK, INFO )

   3 *

   4 *  -- LAPACK driver routine (version 3.2.2) --

   5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --

   6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

   7 *     June 2010

   8 *     8-15-00:  Improve consistency of WS calculations (eca)

   9 *

  10 *     .. Scalar Arguments ..

  11       CHARACTER          JOBZ

  12       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N

  13 *     ..

  14 *     .. Array Arguments ..

  15       INTEGER            IWORK( * )

  16       DOUBLE PRECISION   RWORK( * ), S( * )

  17       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),

  18      $                   WORK( * )

  19 *     ..

  20 *

  21 *  Purpose

  22 *  =======

  23 *

  24 *  ZGESDD computes the singular value decomposition (SVD) of a complex

  25 *  M-by-N matrix A, optionally computing the left and/or right singular

  26 *  vectors, by using divide-and-conquer method. The SVD is written

  27 *

  28 *       A = U * SIGMA * conjugate-transpose(V)

  29 *

  30 *  where SIGMA is an M-by-N matrix which is zero except for its

  31 *  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and

  32 *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA

  33 *  are the singular values of A; they are real and non-negative, and

  34 *  are returned in descending order.  The first min(m,n) columns of

  35 *  U and V are the left and right singular vectors of A.

  36 *

  37 *  Note that the routine returns VT = V**H, not V.

  38 *

  39 *  The divide and conquer algorithm makes very mild assumptions about

  40 *  floating point arithmetic. It will work on machines with a guard

  41 *  digit in add/subtract, or on those binary machines without guard

  42 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or

  43 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines

  44 *  without guard digits, but we know of none.

  45 *

  46 *  Arguments

  47 *  =========

  48 *

  49 *  JOBZ    (input) CHARACTER*1

  50 *          Specifies options for computing all or part of the matrix U:

  51 *          = 'A':  all M columns of U and all N rows of V**H are

  52 *                  returned in the arrays U and VT;

  53 *          = 'S':  the first min(M,N) columns of U and the first

  54 *                  min(M,N) rows of V**H are returned in the arrays U

  55 *                  and VT;

  56 *          = 'O':  If M >= N, the first N columns of U are overwritten

  57 *                  in the array A and all rows of V**H are returned in

  58 *                  the array VT;

  59 *                  otherwise, all columns of U are returned in the

  60 *                  array U and the first M rows of V**H are overwritten

  61 *                  in the array A;

  62 *          = 'N':  no columns of U or rows of V**H are computed.

  63 *

  64 *  M       (input) INTEGER

  65 *          The number of rows of the input matrix A.  M >= 0.

  66 *

  67 *  N       (input) INTEGER

  68 *          The number of columns of the input matrix A.  N >= 0.

  69 *

  70 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)

  71 *          On entry, the M-by-N matrix A.

  72 *          On exit,

  73 *          if JOBZ = 'O',  A is overwritten with the first N columns

  74 *                          of U (the left singular vectors, stored

  75 *                          columnwise) if M >= N;

  76 *                          A is overwritten with the first M rows

  77 *                          of V**H (the right singular vectors, stored

  78 *                          rowwise) otherwise.

  79 *          if JOBZ .ne. 'O', the contents of A are destroyed.

  80 *

  81 *  LDA     (input) INTEGER

  82 *          The leading dimension of the array A.  LDA >= max(1,M).

  83 *

  84 *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))

  85 *          The singular values of A, sorted so that S(i) >= S(i+1).

  86 *

  87 *  U       (output) COMPLEX*16 array, dimension (LDU,UCOL)

  88 *          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;

  89 *          UCOL = min(M,N) if JOBZ = 'S'.

  90 *          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M

  91 *          unitary matrix U;

  92 *          if JOBZ = 'S', U contains the first min(M,N) columns of U

  93 *          (the left singular vectors, stored columnwise);

  94 *          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.

  95 *

  96 *  LDU     (input) INTEGER

  97 *          The leading dimension of the array U.  LDU >= 1; if

  98 *          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.

  99 *

 100 *  VT      (output) COMPLEX*16 array, dimension (LDVT,N)

 101 *          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the

 102 *          N-by-N unitary matrix V**H;

 103 *          if JOBZ = 'S', VT contains the first min(M,N) rows of

 104 *          V**H (the right singular vectors, stored rowwise);

 105 *          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.

 106 *

 107 *  LDVT    (input) INTEGER

 108 *          The leading dimension of the array VT.  LDVT >= 1; if

 109 *          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;

 110 *          if JOBZ = 'S', LDVT >= min(M,N).

 111 *

 112 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))

 113 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

 114 *

 115 *  LWORK   (input) INTEGER

 116 *          The dimension of the array WORK. LWORK >= 1.

 117 *          if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N).

 118 *          if JOBZ = 'O',

 119 *                LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).

 120 *          if JOBZ = 'S' or 'A',

 121 *                LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N).

 122 *          For good performance, LWORK should generally be larger.

 123 *

 124 *          If LWORK = -1, a workspace query is assumed.  The optimal

 125 *          size for the WORK array is calculated and stored in WORK(1),

 126 *          and no other work except argument checking is performed.

 127 *

 128 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))

 129 *          If JOBZ = 'N', LRWORK >= 5*min(M,N).

 130 *          Otherwise,

 131 *          LRWORK >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1)

 132 *

 133 *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))

 134 *

 135 *  INFO    (output) INTEGER

 136 *          = 0:  successful exit.

 137 *          < 0:  if INFO = -i, the i-th argument had an illegal value.

 138 *          > 0:  The updating process of DBDSDC did not converge.

 139 *

 140 *  Further Details

 141 *  ===============

 142 *

 143 *  Based on contributions by

 144 *     Ming Gu and Huan Ren, Computer Science Division, University of

 145 *     California at Berkeley, USA

 146 *

 147 *  =====================================================================

 148 *

 149 *     .. Parameters ..

 150       INTEGER            LQUERV

 151       PARAMETER          ( LQUERV = -1 )

 152       COMPLEX*16         CZERO, CONE

 153       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),

 154      $                   CONE = ( 1.0D+0, 0.0D+0 ) )

 155       DOUBLE PRECISION   ZERO, ONE

 156       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )

 157 *     ..

 158 *     .. Local Scalars ..

 159       LOGICAL            WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS

 160       INTEGER            BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,

 161      $                   ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,

 162      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,

 163      $                   MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL

 164       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM

 165 *     ..

 166 *     .. Local Arrays ..

 167       INTEGER            IDUM( 1 )

 168       DOUBLE PRECISION   DUM( 1 )

 169 *     ..

 170 *     .. External Subroutines ..

 171       EXTERNAL           DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM,

 172      $                   ZGEQRF, ZLACP2, ZLACPY, ZLACRM, ZLARCM, ZLASCL,

 173      $                   ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR

 174 *     ..

 175 *     .. External Functions ..

 176       LOGICAL            LSAME

 177       INTEGER            ILAENV

 178       DOUBLE PRECISION   DLAMCH, ZLANGE

 179       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE

 180 *     ..

 181 *     .. Intrinsic Functions ..

 182       INTRINSIC          INT, MAX, MIN, SQRT

 183 *     ..

 184 *     .. Executable Statements ..

 185 *

 186 *     Test the input arguments

 187 *

 188       INFO = 0

 189       MINMN = MIN( M, N )

 190       MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )

 191       MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )

 192       WNTQA = LSAME( JOBZ, 'A' )

 193       WNTQS = LSAME( JOBZ, 'S' )

 194       WNTQAS = WNTQA .OR. WNTQS

 195       WNTQO = LSAME( JOBZ, 'O' )

 196       WNTQN = LSAME( JOBZ, 'N' )

 197       MINWRK = 1

 198       MAXWRK = 1

 199 *

 200       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN

 201          INFO = -1

 202       ELSE IF( M.LT.0 ) THEN

 203          INFO = -2

 204       ELSE IF( N.LT.0 ) THEN

 205          INFO = -3

 206       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN

 207          INFO = -5

 208       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.

 209      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN

 210          INFO = -8

 211       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.

 212      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.

 213      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN

 214          INFO = -10

 215       END IF

 216 *

 217 *     Compute workspace

 218 *      (Note: Comments in the code beginning "Workspace:" describe the

 219 *       minimal amount of workspace needed at that point in the code,

 220 *       as well as the preferred amount for good performance.

 221 *       CWorkspace refers to complex workspace, and RWorkspace to

 222 *       real workspace. NB refers to the optimal block size for the

 223 *       immediately following subroutine, as returned by ILAENV.)

 224 *

 225       IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN

 226          IF( M.GE.N ) THEN

 227 *

 228 *           There is no complex work space needed for bidiagonal SVD

 229 *           The real work space needed for bidiagonal SVD is BDSPAC

 230 *           for computing singular values and singular vectors; BDSPAN

 231 *           for computing singular values only.

 232 *           BDSPAC = 5*N*N + 7*N

 233 *           BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8))

 234 *

 235             IF( M.GE.MNTHR1 ) THEN

 236                IF( WNTQN ) THEN

 237 *

 238 *                 Path 1 (M much larger than N, JOBZ='N')

 239 *

 240                   MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1,

 241      $                     -1 )

 242                   MAXWRK = MAX( MAXWRK, 2*N+2*N*

 243      $                     ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )

 244                   MINWRK = 3*N

 245                ELSE IF( WNTQO ) THEN

 246 *

 247 *                 Path 2 (M much larger than N, JOBZ='O')

 248 *

 249                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )

 250                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,

 251      $                    N, N, -1 ) )

 252                   WRKBL = MAX( WRKBL, 2*N+2*N*

 253      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )

 254                   WRKBL = MAX( WRKBL, 2*N+N*

 255      $                    ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )

 256                   WRKBL = MAX( WRKBL, 2*N+N*

 257      $                    ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )

 258                   MAXWRK = M*N + N*N + WRKBL

 259                   MINWRK = 2*N*N + 3*N

 260                ELSE IF( WNTQS ) THEN

 261 *

 262 *                 Path 3 (M much larger than N, JOBZ='S')

 263 *

 264                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )

 265                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,

 266      $                    N, N, -1 ) )

 267                   WRKBL = MAX( WRKBL, 2*N+2*N*

 268      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )

 269                   WRKBL = MAX( WRKBL, 2*N+N*

 270      $                    ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )

 271                   WRKBL = MAX( WRKBL, 2*N+N*

 272      $                    ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )

 273                   MAXWRK = N*N + WRKBL

 274                   MINWRK = N*N + 3*N

 275                ELSE IF( WNTQA ) THEN

 276 *

 277 *                 Path 4 (M much larger than N, JOBZ='A')

 278 *

 279                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )

 280                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,

 281      $                    M, N, -1 ) )

 282                   WRKBL = MAX( WRKBL, 2*N+2*N*

 283      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )

 284                   WRKBL = MAX( WRKBL, 2*N+N*

 285      $                    ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )

 286                   WRKBL = MAX( WRKBL, 2*N+N*

 287      $                    ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )

 288                   MAXWRK = N*N + WRKBL

 289                   MINWRK = N*N + 2*N + M

 290                END IF

 291             ELSE IF( M.GE.MNTHR2 ) THEN

 292 *

 293 *              Path 5 (M much larger than N, but not as much as MNTHR1)

 294 *

 295                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,

 296      $                  -1, -1 )

 297                MINWRK = 2*N + M

 298                IF( WNTQO ) THEN

 299                   MAXWRK = MAX( MAXWRK, 2*N+N*

 300      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )

 301                   MAXWRK = MAX( MAXWRK, 2*N+N*

 302      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )

 303                   MAXWRK = MAXWRK + M*N

 304                   MINWRK = MINWRK + N*N

 305                ELSE IF( WNTQS ) THEN

 306                   MAXWRK = MAX( MAXWRK, 2*N+N*

 307      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )

 308                   MAXWRK = MAX( MAXWRK, 2*N+N*

 309      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )

 310                ELSE IF( WNTQA ) THEN

 311                   MAXWRK = MAX( MAXWRK, 2*N+N*

 312      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )

 313                   MAXWRK = MAX( MAXWRK, 2*N+M*

 314      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )

 315                END IF

 316             ELSE

 317 *

 318 *              Path 6 (M at least N, but not much larger)

 319 *

 320                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,

 321      $                  -1, -1 )

 322                MINWRK = 2*N + M

 323                IF( WNTQO ) THEN

 324                   MAXWRK = MAX( MAXWRK, 2*N+N*

 325      $                     ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )

 326                   MAXWRK = MAX( MAXWRK, 2*N+N*

 327      $                     ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) )

 328                   MAXWRK = MAXWRK + M*N

 329                   MINWRK = MINWRK + N*N

 330                ELSE IF( WNTQS ) THEN

 331                   MAXWRK = MAX( MAXWRK, 2*N+N*

 332      $                     ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )

 333                   MAXWRK = MAX( MAXWRK, 2*N+N*

 334      $                     ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) )

 335                ELSE IF( WNTQA ) THEN

 336                   MAXWRK = MAX( MAXWRK, 2*N+N*

 337      $                     ILAENV( 1, 'ZUNGBR', 'PRC', N, N, N, -1 ) )

 338                   MAXWRK = MAX( MAXWRK, 2*N+M*

 339      $                     ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )

 340                END IF

 341             END IF

 342          ELSE

 343 *

 344 *           There is no complex work space needed for bidiagonal SVD

 345 *           The real work space needed for bidiagonal SVD is BDSPAC

 346 *           for computing singular values and singular vectors; BDSPAN

 347 *           for computing singular values only.

 348 *           BDSPAC = 5*M*M + 7*M

 349 *           BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8))

 350 *

 351             IF( N.GE.MNTHR1 ) THEN

 352                IF( WNTQN ) THEN

 353 *

 354 *                 Path 1t (N much larger than M, JOBZ='N')

 355 *

 356                   MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,

 357      $                     -1 )

 358                   MAXWRK = MAX( MAXWRK, 2*M+2*M*

 359      $                     ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )

 360                   MINWRK = 3*M

 361                ELSE IF( WNTQO ) THEN

 362 *

 363 *                 Path 2t (N much larger than M, JOBZ='O')

 364 *

 365                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )

 366                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,

 367      $                    N, M, -1 ) )

 368                   WRKBL = MAX( WRKBL, 2*M+2*M*

 369      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )

 370                   WRKBL = MAX( WRKBL, 2*M+M*

 371      $                    ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )

 372                   WRKBL = MAX( WRKBL, 2*M+M*

 373      $                    ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )

 374                   MAXWRK = M*N + M*M + WRKBL

 375                   MINWRK = 2*M*M + 3*M

 376                ELSE IF( WNTQS ) THEN

 377 *

 378 *                 Path 3t (N much larger than M, JOBZ='S')

 379 *

 380                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )

 381                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,

 382      $                    N, M, -1 ) )

 383                   WRKBL = MAX( WRKBL, 2*M+2*M*

 384      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )

 385                   WRKBL = MAX( WRKBL, 2*M+M*

 386      $                    ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )

 387                   WRKBL = MAX( WRKBL, 2*M+M*

 388      $                    ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )

 389                   MAXWRK = M*M + WRKBL

 390                   MINWRK = M*M + 3*M

 391                ELSE IF( WNTQA ) THEN

 392 *

 393 *                 Path 4t (N much larger than M, JOBZ='A')

 394 *

 395                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )

 396                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,

 397      $                    N, M, -1 ) )

 398                   WRKBL = MAX( WRKBL, 2*M+2*M*

 399      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )

 400                   WRKBL = MAX( WRKBL, 2*M+M*

 401      $                    ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )

 402                   WRKBL = MAX( WRKBL, 2*M+M*

 403      $                    ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )

 404                   MAXWRK = M*M + WRKBL

 405                   MINWRK = M*M + 2*M + N

 406                END IF

 407             ELSE IF( N.GE.MNTHR2 ) THEN

 408 *

 409 *              Path 5t (N much larger than M, but not as much as MNTHR1)

 410 *

 411                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,

 412      $                  -1, -1 )

 413                MINWRK = 2*M + N

 414                IF( WNTQO ) THEN

 415                   MAXWRK = MAX( MAXWRK, 2*M+M*

 416      $                     ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )

 417                   MAXWRK = MAX( MAXWRK, 2*M+M*

 418      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )

 419                   MAXWRK = MAXWRK + M*N

 420                   MINWRK = MINWRK + M*M

 421                ELSE IF( WNTQS ) THEN

 422                   MAXWRK = MAX( MAXWRK, 2*M+M*

 423      $                     ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )

 424                   MAXWRK = MAX( MAXWRK, 2*M+M*

 425      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )

 426                ELSE IF( WNTQA ) THEN

 427                   MAXWRK = MAX( MAXWRK, 2*M+N*

 428      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) )

 429                   MAXWRK = MAX( MAXWRK, 2*M+M*

 430      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )

 431                END IF

 432             ELSE

 433 *

 434 *              Path 6t (N greater than M, but not much larger)

 435 *

 436                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,

 437      $                  -1, -1 )

 438                MINWRK = 2*M + N

 439                IF( WNTQO ) THEN

 440                   MAXWRK = MAX( MAXWRK, 2*M+M*

 441      $                     ILAENV( 1, 'ZUNMBR', 'PRC', M, N, M, -1 ) )

 442                   MAXWRK = MAX( MAXWRK, 2*M+M*

 443      $                     ILAENV( 1, 'ZUNMBR', 'QLN', M, M, N, -1 ) )

 444                   MAXWRK = MAXWRK + M*N

 445                   MINWRK = MINWRK + M*M

 446                ELSE IF( WNTQS ) THEN

 447                   MAXWRK = MAX( MAXWRK, 2*M+M*

 448      $                     ILAENV( 1, 'ZUNGBR', 'PRC', M, N, M, -1 ) )

 449                   MAXWRK = MAX( MAXWRK, 2*M+M*

 450      $                     ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )

 451                ELSE IF( WNTQA ) THEN

 452                   MAXWRK = MAX( MAXWRK, 2*M+N*

 453      $                     ILAENV( 1, 'ZUNGBR', 'PRC', N, N, M, -1 ) )

 454                   MAXWRK = MAX( MAXWRK, 2*M+M*

 455      $                     ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )

 456                END IF

 457             END IF

 458          END IF

 459          MAXWRK = MAX( MAXWRK, MINWRK )

 460       END IF

 461       IF( INFO.EQ.0 ) THEN

 462          WORK( 1 ) = MAXWRK

 463          IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )

 464      $      INFO = -13

 465       END IF

 466 *

 467 *     Quick returns

 468 *

 469       IF( INFO.NE.0 ) THEN

 470          CALL XERBLA( 'ZGESDD', -INFO )

 471          RETURN

 472       END IF

 473       IF( LWORK.EQ.LQUERV )

 474      $   RETURN

 475       IF( M.EQ.0 .OR. N.EQ.0 ) THEN

 476          RETURN

 477       END IF

 478 *

 479 *     Get machine constants

 480 *

 481       EPS = DLAMCH( 'P' )

 482       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS

 483       BIGNUM = ONE / SMLNUM

 484 *

 485 *     Scale A if max element outside range [SMLNUM,BIGNUM]

 486 *

 487       ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )

 488       ISCL = 0

 489       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN

 490          ISCL = 1

 491          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )

 492       ELSE IF( ANRM.GT.BIGNUM ) THEN

 493          ISCL = 1

 494          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )

 495       END IF

 496 *

 497       IF( M.GE.N ) THEN

 498 *

 499 *        A has at least as many rows as columns. If A has sufficiently

 500 *        more rows than columns, first reduce using the QR

 501 *        decomposition (if sufficient workspace available)

 502 *

 503          IF( M.GE.MNTHR1 ) THEN

 504 *

 505             IF( WNTQN ) THEN

 506 *

 507 *              Path 1 (M much larger than N, JOBZ='N')

 508 *              No singular vectors to be computed

 509 *

 510                ITAU = 1

 511                NWORK = ITAU + N

 512 *

 513 *              Compute A=Q*R

 514 *              (CWorkspace: need 2*N, prefer N+N*NB)

 515 *              (RWorkspace: need 0)

 516 *

 517                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

 518      $                      LWORK-NWORK+1, IERR )

 519 *

 520 *              Zero out below R

 521 *

 522                CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),

 523      $                      LDA )

 524                IE = 1

 525                ITAUQ = 1

 526                ITAUP = ITAUQ + N

 527                NWORK = ITAUP + N

 528 *

 529 *              Bidiagonalize R in A

 530 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)

 531 *              (RWorkspace: need N)

 532 *

 533                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

 534      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

 535      $                      IERR )

 536                NRWORK = IE + N

 537 *

 538 *              Perform bidiagonal SVD, compute singular values only

 539 *              (CWorkspace: 0)

 540 *              (RWorkspace: need BDSPAN)

 541 *

 542                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,

 543      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )

 544 *

 545             ELSE IF( WNTQO ) THEN

 546 *

 547 *              Path 2 (M much larger than N, JOBZ='O')

 548 *              N left singular vectors to be overwritten on A and

 549 *              N right singular vectors to be computed in VT

 550 *

 551                IU = 1

 552 *

 553 *              WORK(IU) is N by N

 554 *

 555                LDWRKU = N

 556                IR = IU + LDWRKU*N

 557                IF( LWORK.GE.M*N+N*N+3*N ) THEN

 558 *

 559 *                 WORK(IR) is M by N

 560 *

 561                   LDWRKR = M

 562                ELSE

 563                   LDWRKR = ( LWORK-N*N-3*N ) / N

 564                END IF

 565                ITAU = IR + LDWRKR*N

 566                NWORK = ITAU + N

 567 *

 568 *              Compute A=Q*R

 569 *              (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)

 570 *              (RWorkspace: 0)

 571 *

 572                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

 573      $                      LWORK-NWORK+1, IERR )

 574 *

 575 *              Copy R to WORK( IR ), zeroing out below it

 576 *

 577                CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )

 578                CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),

 579      $                      LDWRKR )

 580 *

 581 *              Generate Q in A

 582 *              (CWorkspace: need 2*N, prefer N+N*NB)

 583 *              (RWorkspace: 0)

 584 *

 585                CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),

 586      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 587                IE = 1

 588                ITAUQ = ITAU

 589                ITAUP = ITAUQ + N

 590                NWORK = ITAUP + N

 591 *

 592 *              Bidiagonalize R in WORK(IR)

 593 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)

 594 *              (RWorkspace: need N)

 595 *

 596                CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),

 597      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),

 598      $                      LWORK-NWORK+1, IERR )

 599 *

 600 *              Perform bidiagonal SVD, computing left singular vectors

 601 *              of R in WORK(IRU) and computing right singular vectors

 602 *              of R in WORK(IRVT)

 603 *              (CWorkspace: need 0)

 604 *              (RWorkspace: need BDSPAC)

 605 *

 606                IRU = IE + N

 607                IRVT = IRU + N*N

 608                NRWORK = IRVT + N*N

 609                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

 610      $                      N, RWORK( IRVT ), N, DUM, IDUM,

 611      $                      RWORK( NRWORK ), IWORK, INFO )

 612 *

 613 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)

 614 *              Overwrite WORK(IU) by the left singular vectors of R

 615 *              (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)

 616 *              (RWorkspace: 0)

 617 *

 618                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),

 619      $                      LDWRKU )

 620                CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,

 621      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,

 622      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 623 *

 624 *              Copy real matrix RWORK(IRVT) to complex matrix VT

 625 *              Overwrite VT by the right singular vectors of R

 626 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)

 627 *              (RWorkspace: 0)

 628 *

 629                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )

 630                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,

 631      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

 632      $                      LWORK-NWORK+1, IERR )

 633 *

 634 *              Multiply Q in A by left singular vectors of R in

 635 *              WORK(IU), storing result in WORK(IR) and copying to A

 636 *              (CWorkspace: need 2*N*N, prefer N*N+M*N)

 637 *              (RWorkspace: 0)

 638 *

 639                DO 10 I = 1, M, LDWRKR

 640                   CHUNK = MIN( M-I+1, LDWRKR )

 641                   CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),

 642      $                        LDA, WORK( IU ), LDWRKU, CZERO,

 643      $                        WORK( IR ), LDWRKR )

 644                   CALL ZLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,

 645      $                         A( I, 1 ), LDA )

 646    10          CONTINUE

 647 *

 648             ELSE IF( WNTQS ) THEN

 649 *

 650 *              Path 3 (M much larger than N, JOBZ='S')

 651 *              N left singular vectors to be computed in U and

 652 *              N right singular vectors to be computed in VT

 653 *

 654                IR = 1

 655 *

 656 *              WORK(IR) is N by N

 657 *

 658                LDWRKR = N

 659                ITAU = IR + LDWRKR*N

 660                NWORK = ITAU + N

 661 *

 662 *              Compute A=Q*R

 663 *              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)

 664 *              (RWorkspace: 0)

 665 *

 666                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

 667      $                      LWORK-NWORK+1, IERR )

 668 *

 669 *              Copy R to WORK(IR), zeroing out below it

 670 *

 671                CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )

 672                CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),

 673      $                      LDWRKR )

 674 *

 675 *              Generate Q in A

 676 *              (CWorkspace: need 2*N, prefer N+N*NB)

 677 *              (RWorkspace: 0)

 678 *

 679                CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),

 680      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 681                IE = 1

 682                ITAUQ = ITAU

 683                ITAUP = ITAUQ + N

 684                NWORK = ITAUP + N

 685 *

 686 *              Bidiagonalize R in WORK(IR)

 687 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)

 688 *              (RWorkspace: need N)

 689 *

 690                CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),

 691      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),

 692      $                      LWORK-NWORK+1, IERR )

 693 *

 694 *              Perform bidiagonal SVD, computing left singular vectors

 695 *              of bidiagonal matrix in RWORK(IRU) and computing right

 696 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

 697 *              (CWorkspace: need 0)

 698 *              (RWorkspace: need BDSPAC)

 699 *

 700                IRU = IE + N

 701                IRVT = IRU + N*N

 702                NRWORK = IRVT + N*N

 703                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

 704      $                      N, RWORK( IRVT ), N, DUM, IDUM,

 705      $                      RWORK( NRWORK ), IWORK, INFO )

 706 *

 707 *              Copy real matrix RWORK(IRU) to complex matrix U

 708 *              Overwrite U by left singular vectors of R

 709 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)

 710 *              (RWorkspace: 0)

 711 *

 712                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )

 713                CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,

 714      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

 715      $                      LWORK-NWORK+1, IERR )

 716 *

 717 *              Copy real matrix RWORK(IRVT) to complex matrix VT

 718 *              Overwrite VT by right singular vectors of R

 719 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)

 720 *              (RWorkspace: 0)

 721 *

 722                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )

 723                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,

 724      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

 725      $                      LWORK-NWORK+1, IERR )

 726 *

 727 *              Multiply Q in A by left singular vectors of R in

 728 *              WORK(IR), storing result in U

 729 *              (CWorkspace: need N*N)

 730 *              (RWorkspace: 0)

 731 *

 732                CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )

 733                CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),

 734      $                     LDWRKR, CZERO, U, LDU )

 735 *

 736             ELSE IF( WNTQA ) THEN

 737 *

 738 *              Path 4 (M much larger than N, JOBZ='A')

 739 *              M left singular vectors to be computed in U and

 740 *              N right singular vectors to be computed in VT

 741 *

 742                IU = 1

 743 *

 744 *              WORK(IU) is N by N

 745 *

 746                LDWRKU = N

 747                ITAU = IU + LDWRKU*N

 748                NWORK = ITAU + N

 749 *

 750 *              Compute A=Q*R, copying result to U

 751 *              (CWorkspace: need 2*N, prefer N+N*NB)

 752 *              (RWorkspace: 0)

 753 *

 754                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

 755      $                      LWORK-NWORK+1, IERR )

 756                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )

 757 *

 758 *              Generate Q in U

 759 *              (CWorkspace: need N+M, prefer N+M*NB)

 760 *              (RWorkspace: 0)

 761 *

 762                CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),

 763      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 764 *

 765 *              Produce R in A, zeroing out below it

 766 *

 767                CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),

 768      $                      LDA )

 769                IE = 1

 770                ITAUQ = ITAU

 771                ITAUP = ITAUQ + N

 772                NWORK = ITAUP + N

 773 *

 774 *              Bidiagonalize R in A

 775 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)

 776 *              (RWorkspace: need N)

 777 *

 778                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

 779      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

 780      $                      IERR )

 781                IRU = IE + N

 782                IRVT = IRU + N*N

 783                NRWORK = IRVT + N*N

 784 *

 785 *              Perform bidiagonal SVD, computing left singular vectors

 786 *              of bidiagonal matrix in RWORK(IRU) and computing right

 787 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

 788 *              (CWorkspace: need 0)

 789 *              (RWorkspace: need BDSPAC)

 790 *

 791                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

 792      $                      N, RWORK( IRVT ), N, DUM, IDUM,

 793      $                      RWORK( NRWORK ), IWORK, INFO )

 794 *

 795 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)

 796 *              Overwrite WORK(IU) by left singular vectors of R

 797 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)

 798 *              (RWorkspace: 0)

 799 *

 800                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),

 801      $                      LDWRKU )

 802                CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,

 803      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,

 804      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 805 *

 806 *              Copy real matrix RWORK(IRVT) to complex matrix VT

 807 *              Overwrite VT by right singular vectors of R

 808 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)

 809 *              (RWorkspace: 0)

 810 *

 811                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )

 812                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,

 813      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

 814      $                      LWORK-NWORK+1, IERR )

 815 *

 816 *              Multiply Q in U by left singular vectors of R in

 817 *              WORK(IU), storing result in A

 818 *              (CWorkspace: need N*N)

 819 *              (RWorkspace: 0)

 820 *

 821                CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),

 822      $                     LDWRKU, CZERO, A, LDA )

 823 *

 824 *              Copy left singular vectors of A from A to U

 825 *

 826                CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )

 827 *

 828             END IF

 829 *

 830          ELSE IF( M.GE.MNTHR2 ) THEN

 831 *

 832 *           MNTHR2 <= M < MNTHR1

 833 *

 834 *           Path 5 (M much larger than N, but not as much as MNTHR1)

 835 *           Reduce to bidiagonal form without QR decomposition, use

 836 *           ZUNGBR and matrix multiplication to compute singular vectors

 837 *

 838             IE = 1

 839             NRWORK = IE + N

 840             ITAUQ = 1

 841             ITAUP = ITAUQ + N

 842             NWORK = ITAUP + N

 843 *

 844 *           Bidiagonalize A

 845 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)

 846 *           (RWorkspace: need N)

 847 *

 848             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

 849      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

 850      $                   IERR )

 851             IF( WNTQN ) THEN

 852 *

 853 *              Compute singular values only

 854 *              (Cworkspace: 0)

 855 *              (Rworkspace: need BDSPAN)

 856 *

 857                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,

 858      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )

 859             ELSE IF( WNTQO ) THEN

 860                IU = NWORK

 861                IRU = NRWORK

 862                IRVT = IRU + N*N

 863                NRWORK = IRVT + N*N

 864 *

 865 *              Copy A to VT, generate P**H

 866 *              (Cworkspace: need 2*N, prefer N+N*NB)

 867 *              (Rworkspace: 0)

 868 *

 869                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )

 870                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),

 871      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 872 *

 873 *              Generate Q in A

 874 *              (CWorkspace: need 2*N, prefer N+N*NB)

 875 *              (RWorkspace: 0)

 876 *

 877                CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),

 878      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 879 *

 880                IF( LWORK.GE.M*N+3*N ) THEN

 881 *

 882 *                 WORK( IU ) is M by N

 883 *

 884                   LDWRKU = M

 885                ELSE

 886 *

 887 *                 WORK(IU) is LDWRKU by N

 888 *

 889                   LDWRKU = ( LWORK-3*N ) / N

 890                END IF

 891                NWORK = IU + LDWRKU*N

 892 *

 893 *              Perform bidiagonal SVD, computing left singular vectors

 894 *              of bidiagonal matrix in RWORK(IRU) and computing right

 895 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

 896 *              (CWorkspace: need 0)

 897 *              (RWorkspace: need BDSPAC)

 898 *

 899                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

 900      $                      N, RWORK( IRVT ), N, DUM, IDUM,

 901      $                      RWORK( NRWORK ), IWORK, INFO )

 902 *

 903 *              Multiply real matrix RWORK(IRVT) by P**H in VT,

 904 *              storing the result in WORK(IU), copying to VT

 905 *              (Cworkspace: need 0)

 906 *              (Rworkspace: need 3*N*N)

 907 *

 908                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,

 909      $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )

 910                CALL ZLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )

 911 *

 912 *              Multiply Q in A by real matrix RWORK(IRU), storing the

 913 *              result in WORK(IU), copying to A

 914 *              (CWorkspace: need N*N, prefer M*N)

 915 *              (Rworkspace: need 3*N*N, prefer N*N+2*M*N)

 916 *

 917                NRWORK = IRVT

 918                DO 20 I = 1, M, LDWRKU

 919                   CHUNK = MIN( M-I+1, LDWRKU )

 920                   CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),

 921      $                         N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )

 922                   CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,

 923      $                         A( I, 1 ), LDA )

 924    20          CONTINUE

 925 *

 926             ELSE IF( WNTQS ) THEN

 927 *

 928 *              Copy A to VT, generate P**H

 929 *              (Cworkspace: need 2*N, prefer N+N*NB)

 930 *              (Rworkspace: 0)

 931 *

 932                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )

 933                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),

 934      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 935 *

 936 *              Copy A to U, generate Q

 937 *              (Cworkspace: need 2*N, prefer N+N*NB)

 938 *              (Rworkspace: 0)

 939 *

 940                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )

 941                CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),

 942      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 943 *

 944 *              Perform bidiagonal SVD, computing left singular vectors

 945 *              of bidiagonal matrix in RWORK(IRU) and computing right

 946 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

 947 *              (CWorkspace: need 0)

 948 *              (RWorkspace: need BDSPAC)

 949 *

 950                IRU = NRWORK

 951                IRVT = IRU + N*N

 952                NRWORK = IRVT + N*N

 953                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

 954      $                      N, RWORK( IRVT ), N, DUM, IDUM,

 955      $                      RWORK( NRWORK ), IWORK, INFO )

 956 *

 957 *              Multiply real matrix RWORK(IRVT) by P**H in VT,

 958 *              storing the result in A, copying to VT

 959 *              (Cworkspace: need 0)

 960 *              (Rworkspace: need 3*N*N)

 961 *

 962                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,

 963      $                      RWORK( NRWORK ) )

 964                CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )

 965 *

 966 *              Multiply Q in U by real matrix RWORK(IRU), storing the

 967 *              result in A, copying to U

 968 *              (CWorkspace: need 0)

 969 *              (Rworkspace: need N*N+2*M*N)

 970 *

 971                NRWORK = IRVT

 972                CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,

 973      $                      RWORK( NRWORK ) )

 974                CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )

 975             ELSE

 976 *

 977 *              Copy A to VT, generate P**H

 978 *              (Cworkspace: need 2*N, prefer N+N*NB)

 979 *              (Rworkspace: 0)

 980 *

 981                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )

 982                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),

 983      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 984 *

 985 *              Copy A to U, generate Q

 986 *              (Cworkspace: need 2*N, prefer N+N*NB)

 987 *              (Rworkspace: 0)

 988 *

 989                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )

 990                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),

 991      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

 992 *

 993 *              Perform bidiagonal SVD, computing left singular vectors

 994 *              of bidiagonal matrix in RWORK(IRU) and computing right

 995 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

 996 *              (CWorkspace: need 0)

 997 *              (RWorkspace: need BDSPAC)

 998 *

 999                IRU = NRWORK

1000                IRVT = IRU + N*N

1001                NRWORK = IRVT + N*N

1002                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

1003      $                      N, RWORK( IRVT ), N, DUM, IDUM,

1004      $                      RWORK( NRWORK ), IWORK, INFO )

1005 *

1006 *              Multiply real matrix RWORK(IRVT) by P**H in VT,

1007 *              storing the result in A, copying to VT

1008 *              (Cworkspace: need 0)

1009 *              (Rworkspace: need 3*N*N)

1010 *

1011                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,

1012      $                      RWORK( NRWORK ) )

1013                CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )

1014 *

1015 *              Multiply Q in U by real matrix RWORK(IRU), storing the

1016 *              result in A, copying to U

1017 *              (CWorkspace: 0)

1018 *              (Rworkspace: need 3*N*N)

1019 *

1020                NRWORK = IRVT

1021                CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,

1022      $                      RWORK( NRWORK ) )

1023                CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )

1024             END IF

1025 *

1026          ELSE

1027 *

1028 *           M .LT. MNTHR2

1029 *

1030 *           Path 6 (M at least N, but not much larger)

1031 *           Reduce to bidiagonal form without QR decomposition

1032 *           Use ZUNMBR to compute singular vectors

1033 *

1034             IE = 1

1035             NRWORK = IE + N

1036             ITAUQ = 1

1037             ITAUP = ITAUQ + N

1038             NWORK = ITAUP + N

1039 *

1040 *           Bidiagonalize A

1041 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)

1042 *           (RWorkspace: need N)

1043 *

1044             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

1045      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

1046      $                   IERR )

1047             IF( WNTQN ) THEN

1048 *

1049 *              Compute singular values only

1050 *              (Cworkspace: 0)

1051 *              (Rworkspace: need BDSPAN)

1052 *

1053                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,

1054      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )

1055             ELSE IF( WNTQO ) THEN

1056                IU = NWORK

1057                IRU = NRWORK

1058                IRVT = IRU + N*N

1059                NRWORK = IRVT + N*N

1060                IF( LWORK.GE.M*N+3*N ) THEN

1061 *

1062 *                 WORK( IU ) is M by N

1063 *

1064                   LDWRKU = M

1065                ELSE

1066 *

1067 *                 WORK( IU ) is LDWRKU by N

1068 *

1069                   LDWRKU = ( LWORK-3*N ) / N

1070                END IF

1071                NWORK = IU + LDWRKU*N

1072 *

1073 *              Perform bidiagonal SVD, computing left singular vectors

1074 *              of bidiagonal matrix in RWORK(IRU) and computing right

1075 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1076 *              (CWorkspace: need 0)

1077 *              (RWorkspace: need BDSPAC)

1078 *

1079                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

1080      $                      N, RWORK( IRVT ), N, DUM, IDUM,

1081      $                      RWORK( NRWORK ), IWORK, INFO )

1082 *

1083 *              Copy real matrix RWORK(IRVT) to complex matrix VT

1084 *              Overwrite VT by right singular vectors of A

1085 *              (Cworkspace: need 2*N, prefer N+N*NB)

1086 *              (Rworkspace: need 0)

1087 *

1088                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )

1089                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,

1090      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

1091      $                      LWORK-NWORK+1, IERR )

1092 *

1093                IF( LWORK.GE.M*N+3*N ) THEN

1094 *

1095 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)

1096 *              Overwrite WORK(IU) by left singular vectors of A, copying

1097 *              to A

1098 *              (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)

1099 *              (Rworkspace: need 0)

1100 *

1101                   CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),

1102      $                         LDWRKU )

1103                   CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),

1104      $                         LDWRKU )

1105                   CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,

1106      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,

1107      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )

1108                   CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )

1109                ELSE

1110 *

1111 *                 Generate Q in A

1112 *                 (Cworkspace: need 2*N, prefer N+N*NB)

1113 *                 (Rworkspace: need 0)

1114 *

1115                   CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),

1116      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )

1117 *

1118 *                 Multiply Q in A by real matrix RWORK(IRU), storing the

1119 *                 result in WORK(IU), copying to A

1120 *                 (CWorkspace: need N*N, prefer M*N)

1121 *                 (Rworkspace: need 3*N*N, prefer N*N+2*M*N)

1122 *

1123                   NRWORK = IRVT

1124                   DO 30 I = 1, M, LDWRKU

1125                      CHUNK = MIN( M-I+1, LDWRKU )

1126                      CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA,

1127      $                            RWORK( IRU ), N, WORK( IU ), LDWRKU,

1128      $                            RWORK( NRWORK ) )

1129                      CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,

1130      $                            A( I, 1 ), LDA )

1131    30             CONTINUE

1132                END IF

1133 *

1134             ELSE IF( WNTQS ) THEN

1135 *

1136 *              Perform bidiagonal SVD, computing left singular vectors

1137 *              of bidiagonal matrix in RWORK(IRU) and computing right

1138 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1139 *              (CWorkspace: need 0)

1140 *              (RWorkspace: need BDSPAC)

1141 *

1142                IRU = NRWORK

1143                IRVT = IRU + N*N

1144                NRWORK = IRVT + N*N

1145                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

1146      $                      N, RWORK( IRVT ), N, DUM, IDUM,

1147      $                      RWORK( NRWORK ), IWORK, INFO )

1148 *

1149 *              Copy real matrix RWORK(IRU) to complex matrix U

1150 *              Overwrite U by left singular vectors of A

1151 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)

1152 *              (RWorkspace: 0)

1153 *

1154                CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )

1155                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )

1156                CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,

1157      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1158      $                      LWORK-NWORK+1, IERR )

1159 *

1160 *              Copy real matrix RWORK(IRVT) to complex matrix VT

1161 *              Overwrite VT by right singular vectors of A

1162 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)

1163 *              (RWorkspace: 0)

1164 *

1165                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )

1166                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,

1167      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

1168      $                      LWORK-NWORK+1, IERR )

1169             ELSE

1170 *

1171 *              Perform bidiagonal SVD, computing left singular vectors

1172 *              of bidiagonal matrix in RWORK(IRU) and computing right

1173 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1174 *              (CWorkspace: need 0)

1175 *              (RWorkspace: need BDSPAC)

1176 *

1177                IRU = NRWORK

1178                IRVT = IRU + N*N

1179                NRWORK = IRVT + N*N

1180                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),

1181      $                      N, RWORK( IRVT ), N, DUM, IDUM,

1182      $                      RWORK( NRWORK ), IWORK, INFO )

1183 *

1184 *              Set the right corner of U to identity matrix

1185 *

1186                CALL ZLASET( 'F', M, M, CZERO, CZERO, U, LDU )

1187                IF( M.GT.N ) THEN

1188                   CALL ZLASET( 'F', M-N, M-N, CZERO, CONE,

1189      $                         U( N+1, N+1 ), LDU )

1190                END IF

1191 *

1192 *              Copy real matrix RWORK(IRU) to complex matrix U

1193 *              Overwrite U by left singular vectors of A

1194 *              (CWorkspace: need 2*N+M, prefer 2*N+M*NB)

1195 *              (RWorkspace: 0)

1196 *

1197                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )

1198                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,

1199      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1200      $                      LWORK-NWORK+1, IERR )

1201 *

1202 *              Copy real matrix RWORK(IRVT) to complex matrix VT

1203 *              Overwrite VT by right singular vectors of A

1204 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)

1205 *              (RWorkspace: 0)

1206 *

1207                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )

1208                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,

1209      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

1210      $                      LWORK-NWORK+1, IERR )

1211             END IF

1212 *

1213          END IF

1214 *

1215       ELSE

1216 *

1217 *        A has more columns than rows. If A has sufficiently more

1218 *        columns than rows, first reduce using the LQ decomposition (if

1219 *        sufficient workspace available)

1220 *

1221          IF( N.GE.MNTHR1 ) THEN

1222 *

1223             IF( WNTQN ) THEN

1224 *

1225 *              Path 1t (N much larger than M, JOBZ='N')

1226 *              No singular vectors to be computed

1227 *

1228                ITAU = 1

1229                NWORK = ITAU + M

1230 *

1231 *              Compute A=L*Q

1232 *              (CWorkspace: need 2*M, prefer M+M*NB)

1233 *              (RWorkspace: 0)

1234 *

1235                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

1236      $                      LWORK-NWORK+1, IERR )

1237 *

1238 *              Zero out above L

1239 *

1240                CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),

1241      $                      LDA )

1242                IE = 1

1243                ITAUQ = 1

1244                ITAUP = ITAUQ + M

1245                NWORK = ITAUP + M

1246 *

1247 *              Bidiagonalize L in A

1248 *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)

1249 *              (RWorkspace: need M)

1250 *

1251                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

1252      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

1253      $                      IERR )

1254                NRWORK = IE + M

1255 *

1256 *              Perform bidiagonal SVD, compute singular values only

1257 *              (CWorkspace: 0)

1258 *              (RWorkspace: need BDSPAN)

1259 *

1260                CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,

1261      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )

1262 *

1263             ELSE IF( WNTQO ) THEN

1264 *

1265 *              Path 2t (N much larger than M, JOBZ='O')

1266 *              M right singular vectors to be overwritten on A and

1267 *              M left singular vectors to be computed in U

1268 *

1269                IVT = 1

1270                LDWKVT = M

1271 *

1272 *              WORK(IVT) is M by M

1273 *

1274                IL = IVT + LDWKVT*M

1275                IF( LWORK.GE.M*N+M*M+3*M ) THEN

1276 *

1277 *                 WORK(IL) M by N

1278 *

1279                   LDWRKL = M

1280                   CHUNK = N

1281                ELSE

1282 *

1283 *                 WORK(IL) is M by CHUNK

1284 *

1285                   LDWRKL = M

1286                   CHUNK = ( LWORK-M*M-3*M ) / M

1287                END IF

1288                ITAU = IL + LDWRKL*CHUNK

1289                NWORK = ITAU + M

1290 *

1291 *              Compute A=L*Q

1292 *              (CWorkspace: need 2*M, prefer M+M*NB)

1293 *              (RWorkspace: 0)

1294 *

1295                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

1296      $                      LWORK-NWORK+1, IERR )

1297 *

1298 *              Copy L to WORK(IL), zeroing about above it

1299 *

1300                CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )

1301                CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,

1302      $                      WORK( IL+LDWRKL ), LDWRKL )

1303 *

1304 *              Generate Q in A

1305 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)

1306 *              (RWorkspace: 0)

1307 *

1308                CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),

1309      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1310                IE = 1

1311                ITAUQ = ITAU

1312                ITAUP = ITAUQ + M

1313                NWORK = ITAUP + M

1314 *

1315 *              Bidiagonalize L in WORK(IL)

1316 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)

1317 *              (RWorkspace: need M)

1318 *

1319                CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),

1320      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),

1321      $                      LWORK-NWORK+1, IERR )

1322 *

1323 *              Perform bidiagonal SVD, computing left singular vectors

1324 *              of bidiagonal matrix in RWORK(IRU) and computing right

1325 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1326 *              (CWorkspace: need 0)

1327 *              (RWorkspace: need BDSPAC)

1328 *

1329                IRU = IE + M

1330                IRVT = IRU + M*M

1331                NRWORK = IRVT + M*M

1332                CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1333      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1334      $                      RWORK( NRWORK ), IWORK, INFO )

1335 *

1336 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)

1337 *              Overwrite WORK(IU) by the left singular vectors of L

1338 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)

1339 *              (RWorkspace: 0)

1340 *

1341                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )

1342                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,

1343      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1344      $                      LWORK-NWORK+1, IERR )

1345 *

1346 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)

1347 *              Overwrite WORK(IVT) by the right singular vectors of L

1348 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)

1349 *              (RWorkspace: 0)

1350 *

1351                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),

1352      $                      LDWKVT )

1353                CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,

1354      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,

1355      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1356 *

1357 *              Multiply right singular vectors of L in WORK(IL) by Q

1358 *              in A, storing result in WORK(IL) and copying to A

1359 *              (CWorkspace: need 2*M*M, prefer M*M+M*N))

1360 *              (RWorkspace: 0)

1361 *

1362                DO 40 I = 1, N, CHUNK

1363                   BLK = MIN( N-I+1, CHUNK )

1364                   CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,

1365      $                        A( 1, I ), LDA, CZERO, WORK( IL ),

1366      $                        LDWRKL )

1367                   CALL ZLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,

1368      $                         A( 1, I ), LDA )

1369    40          CONTINUE

1370 *

1371             ELSE IF( WNTQS ) THEN

1372 *

1373 *             Path 3t (N much larger than M, JOBZ='S')

1374 *             M right singular vectors to be computed in VT and

1375 *             M left singular vectors to be computed in U

1376 *

1377                IL = 1

1378 *

1379 *              WORK(IL) is M by M

1380 *

1381                LDWRKL = M

1382                ITAU = IL + LDWRKL*M

1383                NWORK = ITAU + M

1384 *

1385 *              Compute A=L*Q

1386 *              (CWorkspace: need 2*M, prefer M+M*NB)

1387 *              (RWorkspace: 0)

1388 *

1389                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

1390      $                      LWORK-NWORK+1, IERR )

1391 *

1392 *              Copy L to WORK(IL), zeroing out above it

1393 *

1394                CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )

1395                CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,

1396      $                      WORK( IL+LDWRKL ), LDWRKL )

1397 *

1398 *              Generate Q in A

1399 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)

1400 *              (RWorkspace: 0)

1401 *

1402                CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),

1403      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1404                IE = 1

1405                ITAUQ = ITAU

1406                ITAUP = ITAUQ + M

1407                NWORK = ITAUP + M

1408 *

1409 *              Bidiagonalize L in WORK(IL)

1410 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)

1411 *              (RWorkspace: need M)

1412 *

1413                CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),

1414      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),

1415      $                      LWORK-NWORK+1, IERR )

1416 *

1417 *              Perform bidiagonal SVD, computing left singular vectors

1418 *              of bidiagonal matrix in RWORK(IRU) and computing right

1419 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1420 *              (CWorkspace: need 0)

1421 *              (RWorkspace: need BDSPAC)

1422 *

1423                IRU = IE + M

1424                IRVT = IRU + M*M

1425                NRWORK = IRVT + M*M

1426                CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1427      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1428      $                      RWORK( NRWORK ), IWORK, INFO )

1429 *

1430 *              Copy real matrix RWORK(IRU) to complex matrix U

1431 *              Overwrite U by left singular vectors of L

1432 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)

1433 *              (RWorkspace: 0)

1434 *

1435                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )

1436                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,

1437      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1438      $                      LWORK-NWORK+1, IERR )

1439 *

1440 *              Copy real matrix RWORK(IRVT) to complex matrix VT

1441 *              Overwrite VT by left singular vectors of L

1442 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)

1443 *              (RWorkspace: 0)

1444 *

1445                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )

1446                CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,

1447      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

1448      $                      LWORK-NWORK+1, IERR )

1449 *

1450 *              Copy VT to WORK(IL), multiply right singular vectors of L

1451 *              in WORK(IL) by Q in A, storing result in VT

1452 *              (CWorkspace: need M*M)

1453 *              (RWorkspace: 0)

1454 *

1455                CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )

1456                CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,

1457      $                     A, LDA, CZERO, VT, LDVT )

1458 *

1459             ELSE IF( WNTQA ) THEN

1460 *

1461 *              Path 9t (N much larger than M, JOBZ='A')

1462 *              N right singular vectors to be computed in VT and

1463 *              M left singular vectors to be computed in U

1464 *

1465                IVT = 1

1466 *

1467 *              WORK(IVT) is M by M

1468 *

1469                LDWKVT = M

1470                ITAU = IVT + LDWKVT*M

1471                NWORK = ITAU + M

1472 *

1473 *              Compute A=L*Q, copying result to VT

1474 *              (CWorkspace: need 2*M, prefer M+M*NB)

1475 *              (RWorkspace: 0)

1476 *

1477                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),

1478      $                      LWORK-NWORK+1, IERR )

1479                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )

1480 *

1481 *              Generate Q in VT

1482 *              (CWorkspace: need M+N, prefer M+N*NB)

1483 *              (RWorkspace: 0)

1484 *

1485                CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),

1486      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1487 *

1488 *              Produce L in A, zeroing out above it

1489 *

1490                CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),

1491      $                      LDA )

1492                IE = 1

1493                ITAUQ = ITAU

1494                ITAUP = ITAUQ + M

1495                NWORK = ITAUP + M

1496 *

1497 *              Bidiagonalize L in A

1498 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)

1499 *              (RWorkspace: need M)

1500 *

1501                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

1502      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

1503      $                      IERR )

1504 *

1505 *              Perform bidiagonal SVD, computing left singular vectors

1506 *              of bidiagonal matrix in RWORK(IRU) and computing right

1507 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1508 *              (CWorkspace: need 0)

1509 *              (RWorkspace: need BDSPAC)

1510 *

1511                IRU = IE + M

1512                IRVT = IRU + M*M

1513                NRWORK = IRVT + M*M

1514                CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1515      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1516      $                      RWORK( NRWORK ), IWORK, INFO )

1517 *

1518 *              Copy real matrix RWORK(IRU) to complex matrix U

1519 *              Overwrite U by left singular vectors of L

1520 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)

1521 *              (RWorkspace: 0)

1522 *

1523                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )

1524                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,

1525      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1526      $                      LWORK-NWORK+1, IERR )

1527 *

1528 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)

1529 *              Overwrite WORK(IVT) by right singular vectors of L

1530 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)

1531 *              (RWorkspace: 0)

1532 *

1533                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),

1534      $                      LDWKVT )

1535                CALL ZUNMBR( 'P', 'R', 'C', M, M, M, A, LDA,

1536      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,

1537      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1538 *

1539 *              Multiply right singular vectors of L in WORK(IVT) by

1540 *              Q in VT, storing result in A

1541 *              (CWorkspace: need M*M)

1542 *              (RWorkspace: 0)

1543 *

1544                CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,

1545      $                     VT, LDVT, CZERO, A, LDA )

1546 *

1547 *              Copy right singular vectors of A from A to VT

1548 *

1549                CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )

1550 *

1551             END IF

1552 *

1553          ELSE IF( N.GE.MNTHR2 ) THEN

1554 *

1555 *           MNTHR2 <= N < MNTHR1

1556 *

1557 *           Path 5t (N much larger than M, but not as much as MNTHR1)

1558 *           Reduce to bidiagonal form without QR decomposition, use

1559 *           ZUNGBR and matrix multiplication to compute singular vectors

1560 *

1561 *

1562             IE = 1

1563             NRWORK = IE + M

1564             ITAUQ = 1

1565             ITAUP = ITAUQ + M

1566             NWORK = ITAUP + M

1567 *

1568 *           Bidiagonalize A

1569 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)

1570 *           (RWorkspace: M)

1571 *

1572             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

1573      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

1574      $                   IERR )

1575 *

1576             IF( WNTQN ) THEN

1577 *

1578 *              Compute singular values only

1579 *              (Cworkspace: 0)

1580 *              (Rworkspace: need BDSPAN)

1581 *

1582                CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,

1583      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )

1584             ELSE IF( WNTQO ) THEN

1585                IRVT = NRWORK

1586                IRU = IRVT + M*M

1587                NRWORK = IRU + M*M

1588                IVT = NWORK

1589 *

1590 *              Copy A to U, generate Q

1591 *              (Cworkspace: need 2*M, prefer M+M*NB)

1592 *              (Rworkspace: 0)

1593 *

1594                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )

1595                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),

1596      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1597 *

1598 *              Generate P**H in A

1599 *              (Cworkspace: need 2*M, prefer M+M*NB)

1600 *              (Rworkspace: 0)

1601 *

1602                CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),

1603      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1604 *

1605                LDWKVT = M

1606                IF( LWORK.GE.M*N+3*M ) THEN

1607 *

1608 *                 WORK( IVT ) is M by N

1609 *

1610                   NWORK = IVT + LDWKVT*N

1611                   CHUNK = N

1612                ELSE

1613 *

1614 *                 WORK( IVT ) is M by CHUNK

1615 *

1616                   CHUNK = ( LWORK-3*M ) / M

1617                   NWORK = IVT + LDWKVT*CHUNK

1618                END IF

1619 *

1620 *              Perform bidiagonal SVD, computing left singular vectors

1621 *              of bidiagonal matrix in RWORK(IRU) and computing right

1622 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1623 *              (CWorkspace: need 0)

1624 *              (RWorkspace: need BDSPAC)

1625 *

1626                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1627      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1628      $                      RWORK( NRWORK ), IWORK, INFO )

1629 *

1630 *              Multiply Q in U by real matrix RWORK(IRVT)

1631 *              storing the result in WORK(IVT), copying to U

1632 *              (Cworkspace: need 0)

1633 *              (Rworkspace: need 2*M*M)

1634 *

1635                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),

1636      $                      LDWKVT, RWORK( NRWORK ) )

1637                CALL ZLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU )

1638 *

1639 *              Multiply RWORK(IRVT) by P**H in A, storing the

1640 *              result in WORK(IVT), copying to A

1641 *              (CWorkspace: need M*M, prefer M*N)

1642 *              (Rworkspace: need 2*M*M, prefer 2*M*N)

1643 *

1644                NRWORK = IRU

1645                DO 50 I = 1, N, CHUNK

1646                   BLK = MIN( N-I+1, CHUNK )

1647                   CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA,

1648      $                         WORK( IVT ), LDWKVT, RWORK( NRWORK ) )

1649                   CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,

1650      $                         A( 1, I ), LDA )

1651    50          CONTINUE

1652             ELSE IF( WNTQS ) THEN

1653 *

1654 *              Copy A to U, generate Q

1655 *              (Cworkspace: need 2*M, prefer M+M*NB)

1656 *              (Rworkspace: 0)

1657 *

1658                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )

1659                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),

1660      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1661 *

1662 *              Copy A to VT, generate P**H

1663 *              (Cworkspace: need 2*M, prefer M+M*NB)

1664 *              (Rworkspace: 0)

1665 *

1666                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )

1667                CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),

1668      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1669 *

1670 *              Perform bidiagonal SVD, computing left singular vectors

1671 *              of bidiagonal matrix in RWORK(IRU) and computing right

1672 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1673 *              (CWorkspace: need 0)

1674 *              (RWorkspace: need BDSPAC)

1675 *

1676                IRVT = NRWORK

1677                IRU = IRVT + M*M

1678                NRWORK = IRU + M*M

1679                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1680      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1681      $                      RWORK( NRWORK ), IWORK, INFO )

1682 *

1683 *              Multiply Q in U by real matrix RWORK(IRU), storing the

1684 *              result in A, copying to U

1685 *              (CWorkspace: need 0)

1686 *              (Rworkspace: need 3*M*M)

1687 *

1688                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,

1689      $                      RWORK( NRWORK ) )

1690                CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )

1691 *

1692 *              Multiply real matrix RWORK(IRVT) by P**H in VT,

1693 *              storing the result in A, copying to VT

1694 *              (Cworkspace: need 0)

1695 *              (Rworkspace: need M*M+2*M*N)

1696 *

1697                NRWORK = IRU

1698                CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,

1699      $                      RWORK( NRWORK ) )

1700                CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )

1701             ELSE

1702 *

1703 *              Copy A to U, generate Q

1704 *              (Cworkspace: need 2*M, prefer M+M*NB)

1705 *              (Rworkspace: 0)

1706 *

1707                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )

1708                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),

1709      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1710 *

1711 *              Copy A to VT, generate P**H

1712 *              (Cworkspace: need 2*M, prefer M+M*NB)

1713 *              (Rworkspace: 0)

1714 *

1715                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )

1716                CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),

1717      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )

1718 *

1719 *              Perform bidiagonal SVD, computing left singular vectors

1720 *              of bidiagonal matrix in RWORK(IRU) and computing right

1721 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1722 *              (CWorkspace: need 0)

1723 *              (RWorkspace: need BDSPAC)

1724 *

1725                IRVT = NRWORK

1726                IRU = IRVT + M*M

1727                NRWORK = IRU + M*M

1728                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1729      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1730      $                      RWORK( NRWORK ), IWORK, INFO )

1731 *

1732 *              Multiply Q in U by real matrix RWORK(IRU), storing the

1733 *              result in A, copying to U

1734 *              (CWorkspace: need 0)

1735 *              (Rworkspace: need 3*M*M)

1736 *

1737                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,

1738      $                      RWORK( NRWORK ) )

1739                CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )

1740 *

1741 *              Multiply real matrix RWORK(IRVT) by P**H in VT,

1742 *              storing the result in A, copying to VT

1743 *              (Cworkspace: need 0)

1744 *              (Rworkspace: need M*M+2*M*N)

1745 *

1746                CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,

1747      $                      RWORK( NRWORK ) )

1748                CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )

1749             END IF

1750 *

1751          ELSE

1752 *

1753 *           N .LT. MNTHR2

1754 *

1755 *           Path 6t (N greater than M, but not much larger)

1756 *           Reduce to bidiagonal form without LQ decomposition

1757 *           Use ZUNMBR to compute singular vectors

1758 *

1759             IE = 1

1760             NRWORK = IE + M

1761             ITAUQ = 1

1762             ITAUP = ITAUQ + M

1763             NWORK = ITAUP + M

1764 *

1765 *           Bidiagonalize A

1766 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)

1767 *           (RWorkspace: M)

1768 *

1769             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),

1770      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,

1771      $                   IERR )

1772             IF( WNTQN ) THEN

1773 *

1774 *              Compute singular values only

1775 *              (Cworkspace: 0)

1776 *              (Rworkspace: need BDSPAN)

1777 *

1778                CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,

1779      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )

1780             ELSE IF( WNTQO ) THEN

1781                LDWKVT = M

1782                IVT = NWORK

1783                IF( LWORK.GE.M*N+3*M ) THEN

1784 *

1785 *                 WORK( IVT ) is M by N

1786 *

1787                   CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ),

1788      $                         LDWKVT )

1789                   NWORK = IVT + LDWKVT*N

1790                ELSE

1791 *

1792 *                 WORK( IVT ) is M by CHUNK

1793 *

1794                   CHUNK = ( LWORK-3*M ) / M

1795                   NWORK = IVT + LDWKVT*CHUNK

1796                END IF

1797 *

1798 *              Perform bidiagonal SVD, computing left singular vectors

1799 *              of bidiagonal matrix in RWORK(IRU) and computing right

1800 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1801 *              (CWorkspace: need 0)

1802 *              (RWorkspace: need BDSPAC)

1803 *

1804                IRVT = NRWORK

1805                IRU = IRVT + M*M

1806                NRWORK = IRU + M*M

1807                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1808      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1809      $                      RWORK( NRWORK ), IWORK, INFO )

1810 *

1811 *              Copy real matrix RWORK(IRU) to complex matrix U

1812 *              Overwrite U by left singular vectors of A

1813 *              (Cworkspace: need 2*M, prefer M+M*NB)

1814 *              (Rworkspace: need 0)

1815 *

1816                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )

1817                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,

1818      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1819      $                      LWORK-NWORK+1, IERR )

1820 *

1821                IF( LWORK.GE.M*N+3*M ) THEN

1822 *

1823 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)

1824 *              Overwrite WORK(IVT) by right singular vectors of A,

1825 *              copying to A

1826 *              (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB)

1827 *              (Rworkspace: need 0)

1828 *

1829                   CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),

1830      $                         LDWKVT )

1831                   CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,

1832      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,

1833      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )

1834                   CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )

1835                ELSE

1836 *

1837 *                 Generate P**H in A

1838 *                 (Cworkspace: need 2*M, prefer M+M*NB)

1839 *                 (Rworkspace: need 0)

1840 *

1841                   CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),

1842      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )

1843 *

1844 *                 Multiply Q in A by real matrix RWORK(IRU), storing the

1845 *                 result in WORK(IU), copying to A

1846 *                 (CWorkspace: need M*M, prefer M*N)

1847 *                 (Rworkspace: need 3*M*M, prefer M*M+2*M*N)

1848 *

1849                   NRWORK = IRU

1850                   DO 60 I = 1, N, CHUNK

1851                      BLK = MIN( N-I+1, CHUNK )

1852                      CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ),

1853      $                            LDA, WORK( IVT ), LDWKVT,

1854      $                            RWORK( NRWORK ) )

1855                      CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,

1856      $                            A( 1, I ), LDA )

1857    60             CONTINUE

1858                END IF

1859             ELSE IF( WNTQS ) THEN

1860 *

1861 *              Perform bidiagonal SVD, computing left singular vectors

1862 *              of bidiagonal matrix in RWORK(IRU) and computing right

1863 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1864 *              (CWorkspace: need 0)

1865 *              (RWorkspace: need BDSPAC)

1866 *

1867                IRVT = NRWORK

1868                IRU = IRVT + M*M

1869                NRWORK = IRU + M*M

1870                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1871      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1872      $                      RWORK( NRWORK ), IWORK, INFO )

1873 *

1874 *              Copy real matrix RWORK(IRU) to complex matrix U

1875 *              Overwrite U by left singular vectors of A

1876 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)

1877 *              (RWorkspace: M*M)

1878 *

1879                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )

1880                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,

1881      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1882      $                      LWORK-NWORK+1, IERR )

1883 *

1884 *              Copy real matrix RWORK(IRVT) to complex matrix VT

1885 *              Overwrite VT by right singular vectors of A

1886 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)

1887 *              (RWorkspace: M*M)

1888 *

1889                CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )

1890                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )

1891                CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,

1892      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

1893      $                      LWORK-NWORK+1, IERR )

1894             ELSE

1895 *

1896 *              Perform bidiagonal SVD, computing left singular vectors

1897 *              of bidiagonal matrix in RWORK(IRU) and computing right

1898 *              singular vectors of bidiagonal matrix in RWORK(IRVT)

1899 *              (CWorkspace: need 0)

1900 *              (RWorkspace: need BDSPAC)

1901 *

1902                IRVT = NRWORK

1903                IRU = IRVT + M*M

1904                NRWORK = IRU + M*M

1905 *

1906                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),

1907      $                      M, RWORK( IRVT ), M, DUM, IDUM,

1908      $                      RWORK( NRWORK ), IWORK, INFO )

1909 *

1910 *              Copy real matrix RWORK(IRU) to complex matrix U

1911 *              Overwrite U by left singular vectors of A

1912 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)

1913 *              (RWorkspace: M*M)

1914 *

1915                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )

1916                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,

1917      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),

1918      $                      LWORK-NWORK+1, IERR )

1919 *

1920 *              Set all of VT to identity matrix

1921 *

1922                CALL ZLASET( 'F', N, N, CZERO, CONE, VT, LDVT )

1923 *

1924 *              Copy real matrix RWORK(IRVT) to complex matrix VT

1925 *              Overwrite VT by right singular vectors of A

1926 *              (CWorkspace: need 2*M+N, prefer 2*M+N*NB)

1927 *              (RWorkspace: M*M)

1928 *

1929                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )

1930                CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,

1931      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),

1932      $                      LWORK-NWORK+1, IERR )

1933             END IF

1934 *

1935          END IF

1936 *

1937       END IF

1938 *

1939 *     Undo scaling if necessary

1940 *

1941       IF( ISCL.EQ.1 ) THEN

1942          IF( ANRM.GT.BIGNUM )

1943      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,

1944      $                   IERR )

1945          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )

1946      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,

1947      $                   RWORK( IE ), MINMN, IERR )

1948          IF( ANRM.LT.SMLNUM )

1949      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,

1950      $                   IERR )

1951          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )

1952      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,

1953      $                   RWORK( IE ), MINMN, IERR )

1954       END IF

1955 *

1956 *     Return optimal workspace in WORK(1)

1957 *

1958       WORK( 1 ) = MAXWRK

1959 *

1960       RETURN

1961 *

1962 *     End of ZGESDD

1963 *

1964       END