SRC/zgesvdx.f

   1 *> \brief <b> ZGESVDX computes the singular value decomposition (SVD) for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGESVDX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *     SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
  22 *    $                    IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
  23 *    $                    LWORK, RWORK, IWORK, INFO )
  24 *
  25 *
  26 *     .. Scalar Arguments ..
  27 *      CHARACTER          JOBU, JOBVT, RANGE
  28 *      INTEGER            IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
  29 *      DOUBLE PRECISION   VL, VU
  30 *     ..
  31 *     .. Array Arguments ..
  32 *      INTEGER            IWORK( * )
  33 *      DOUBLE PRECISION   S( * ), RWORK( * )
  34 *      COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
  35 *     $                   WORK( * )
  36 *     ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *>  ZGESVDX computes the singular value decomposition (SVD) of a complex
  45 *>  M-by-N matrix A, optionally computing the left and/or right singular
  46 *>  vectors. The SVD is written
  47 *>
  48 *>      A = U * SIGMA * transpose(V)
  49 *>
  50 *>  where SIGMA is an M-by-N matrix which is zero except for its
  51 *>  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
  52 *>  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
  53 *>  are the singular values of A; they are real and non-negative, and
  54 *>  are returned in descending order.  The first min(m,n) columns of
  55 *>  U and V are the left and right singular vectors of A.
  56 *>
  57 *>  ZGESVDX uses an eigenvalue problem for obtaining the SVD, which
  58 *>  allows for the computation of a subset of singular values and
  59 *>  vectors. See DBDSVDX for details.
  60 *>
  61 *>  Note that the routine returns V**T, not V.
  62 *> \endverbatim
  63 *
  64 *  Arguments:
  65 *  ==========
  66 *
  67 *> \param[in] JOBU
  68 *> \verbatim
  69 *>          JOBU is CHARACTER*1
  70 *>          Specifies options for computing all or part of the matrix U:
  71 *>          = 'V':  the first min(m,n) columns of U (the left singular
  72 *>                  vectors) or as specified by RANGE are returned in
  73 *>                  the array U;
  74 *>          = 'N':  no columns of U (no left singular vectors) are
  75 *>                  computed.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] JOBVT
  79 *> \verbatim
  80 *>          JOBVT is CHARACTER*1
  81 *>           Specifies options for computing all or part of the matrix
  82 *>           V**T:
  83 *>           = 'V':  the first min(m,n) rows of V**T (the right singular
  84 *>                   vectors) or as specified by RANGE are returned in
  85 *>                   the array VT;
  86 *>           = 'N':  no rows of V**T (no right singular vectors) are
  87 *>                   computed.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] RANGE
  91 *> \verbatim
  92 *>          RANGE is CHARACTER*1
  93 *>          = 'A': all singular values will be found.
  94 *>          = 'V': all singular values in the half-open interval (VL,VU]
  95 *>                 will be found.
  96 *>          = 'I': the IL-th through IU-th singular values will be found.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] M
 100 *> \verbatim
 101 *>          M is INTEGER
 102 *>          The number of rows of the input matrix A.  M >= 0.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] N
 106 *> \verbatim
 107 *>          N is INTEGER
 108 *>          The number of columns of the input matrix A.  N >= 0.
 109 *> \endverbatim
 110 *>
 111 *> \param[in,out] A
 112 *> \verbatim
 113 *>          A is COMPLEX*16 array, dimension (LDA,N)
 114 *>          On entry, the M-by-N matrix A.
 115 *>          On exit, the contents of A are destroyed.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] LDA
 119 *> \verbatim
 120 *>          LDA is INTEGER
 121 *>          The leading dimension of the array A.  LDA >= max(1,M).
 122 *> \endverbatim
 123 *>
 124 *> \param[in] VL
 125 *> \verbatim
 126 *>          VL is DOUBLE PRECISION
 127 *>          If RANGE='V', the lower bound of the interval to
 128 *>          be searched for singular values. VU > VL.
 129 *>          Not referenced if RANGE = 'A' or 'I'.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] VU
 133 *> \verbatim
 134 *>          VU is DOUBLE PRECISION
 135 *>          If RANGE='V', the upper bound of the interval to
 136 *>          be searched for singular values. VU > VL.
 137 *>          Not referenced if RANGE = 'A' or 'I'.
 138 *> \endverbatim
 139 *>
 140 *> \param[in] IL
 141 *> \verbatim
 142 *>          IL is INTEGER
 143 *>          If RANGE='I', the index of the
 144 *>          smallest singular value to be returned.
 145 *>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
 146 *>          Not referenced if RANGE = 'A' or 'V'.
 147 *> \endverbatim
 148 *>
 149 *> \param[in] IU
 150 *> \verbatim
 151 *>          IU is INTEGER
 152 *>          If RANGE='I', the index of the
 153 *>          largest singular value to be returned.
 154 *>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
 155 *>          Not referenced if RANGE = 'A' or 'V'.
 156 *> \endverbatim
 157 *>
 158 *> \param[out] NS
 159 *> \verbatim
 160 *>          NS is INTEGER
 161 *>          The total number of singular values found,
 162 *>          0 <= NS <= min(M,N).
 163 *>          If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
 164 *> \endverbatim
 165 *>
 166 *> \param[out] S
 167 *> \verbatim
 168 *>          S is DOUBLE PRECISION array, dimension (min(M,N))
 169 *>          The singular values of A, sorted so that S(i) >= S(i+1).
 170 *> \endverbatim
 171 *>
 172 *> \param[out] U
 173 *> \verbatim
 174 *>          U is COMPLEX*16 array, dimension (LDU,UCOL)
 175 *>          If JOBU = 'V', U contains columns of U (the left singular
 176 *>          vectors, stored columnwise) as specified by RANGE; if
 177 *>          JOBU = 'N', U is not referenced.
 178 *>          Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
 179 *>          the exact value of NS is not known in advance and an upper
 180 *>          bound must be used.
 181 *> \endverbatim
 182 *>
 183 *> \param[in] LDU
 184 *> \verbatim
 185 *>          LDU is INTEGER
 186 *>          The leading dimension of the array U.  LDU >= 1; if
 187 *>          JOBU = 'V', LDU >= M.
 188 *> \endverbatim
 189 *>
 190 *> \param[out] VT
 191 *> \verbatim
 192 *>          VT is COMPLEX*16 array, dimension (LDVT,N)
 193 *>          If JOBVT = 'V', VT contains the rows of V**T (the right singular
 194 *>          vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
 195 *>          VT is not referenced.
 196 *>          Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
 197 *>          the exact value of NS is not known in advance and an upper
 198 *>          bound must be used.
 199 *> \endverbatim
 200 *>
 201 *> \param[in] LDVT
 202 *> \verbatim
 203 *>          LDVT is INTEGER
 204 *>          The leading dimension of the array VT.  LDVT >= 1; if
 205 *>          JOBVT = 'V', LDVT >= NS (see above).
 206 *> \endverbatim
 207 *>
 208 *> \param[out] WORK
 209 *> \verbatim
 210 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 211 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
 212 *> \endverbatim
 213 *>
 214 *> \param[in] LWORK
 215 *> \verbatim
 216 *>          LWORK is INTEGER
 217 *>          The dimension of the array WORK.
 218 *>          LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
 219 *>          comments inside the code):
 220 *>             - PATH 1  (M much larger than N)
 221 *>             - PATH 1t (N much larger than M)
 222 *>          LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
 223 *>          For good performance, LWORK should generally be larger.
 224 *>
 225 *>          If LWORK = -1, then a workspace query is assumed; the routine
 226 *>          only calculates the optimal size of the WORK array, returns
 227 *>          this value as the first entry of the WORK array, and no error
 228 *>          message related to LWORK is issued by XERBLA.
 229 *> \endverbatim
 230 *>
 231 *> \param[out] RWORK
 232 *> \verbatim
 233 *>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
 234 *>          LRWORK >= MIN(M,N)*(MIN(M,N)*2+15*MIN(M,N)).
 235 *> \endverbatim
 236 *>
 237 *> \param[out] IWORK
 238 *> \verbatim
 239 *>          IWORK is INTEGER array, dimension (12*MIN(M,N))
 240 *>          If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
 241 *>          then IWORK contains the indices of the eigenvectors that failed
 242 *>          to converge in DBDSVDX/DSTEVX.
 243 *> \endverbatim
 244 *>
 245 *> \param[out] INFO
 246 *> \verbatim
 247 *>     INFO is INTEGER
 248 *>           = 0:  successful exit
 249 *>           < 0:  if INFO = -i, the i-th argument had an illegal value
 250 *>           > 0:  if INFO = i, then i eigenvectors failed to converge
 251 *>                 in DBDSVDX/DSTEVX.
 252 *>                 if INFO = N*2 + 1, an internal error occurred in
 253 *>                 DBDSVDX
 254 *> \endverbatim
 255 *
 256 *  Authors:
 257 *  ========
 258 *
 259 *> \author Univ. of Tennessee
 260 *> \author Univ. of California Berkeley
 261 *> \author Univ. of Colorado Denver
 262 *> \author NAG Ltd.
 263 *
 264 *> \date June 2016
 265 *
 266 *> \ingroup complex16GEsing
 267 *
 268 *  =====================================================================
 269       SUBROUTINE ZGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
 270      $                    IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
 271      $                    LWORK, RWORK, IWORK, INFO )
 272 *
 273 *  -- LAPACK driver routine (version 3.6.1) --
 274 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 275 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 276 *     June 2016
 277 *
 278 *     .. Scalar Arguments ..
 279       CHARACTER          JOBU, JOBVT, RANGE
 280       INTEGER            IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
 281       DOUBLE PRECISION   VL, VU
 282 *     ..
 283 *     .. Array Arguments ..
 284       INTEGER            IWORK( * )
 285       DOUBLE PRECISION   S( * ), RWORK( * )
 286       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
 287      $                   WORK( * )
 288 *     ..
 289 *
 290 *  =====================================================================
 291 *
 292 *     .. Parameters ..
 293       COMPLEX*16         CZERO, CONE
 294       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
 295      $                   CONE = ( 1.0D0, 0.0D0 ) )
 296       DOUBLE PRECISION   ZERO, ONE
 297       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 298 *     ..
 299 *     .. Local Scalars ..
 300       CHARACTER          JOBZ, RNGTGK
 301       LOGICAL            ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
 302       INTEGER            I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
 303      $                   ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ,
 304      $                   IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR
 305       DOUBLE PRECISION   ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
 306 *     ..
 307 *     .. Local Arrays ..
 308       DOUBLE PRECISION   DUM( 1 )
 309 *     ..
 310 *     .. External Subroutines ..
 311       EXTERNAL           ZGEBRD, ZGELQF, ZGEQRF, ZLASCL, ZLASET,
 312      $                   DLASCL, XERBLA
 313 *     ..
 314 *     .. External Functions ..
 315       LOGICAL            LSAME
 316       INTEGER            ILAENV
 317       DOUBLE PRECISION   DLAMCH, ZLANGE
 318       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
 319 *     ..
 320 *     .. Intrinsic Functions ..
 321       INTRINSIC          MAX, MIN, SQRT
 322 *     ..
 323 *     .. Executable Statements ..
 324 *
 325 *     Test the input arguments.
 326 *
 327       NS = 0
 328       INFO = 0
 329       ABSTOL = 2*DLAMCH('S')
 330       LQUERY = ( LWORK.EQ.-1 )
 331       MINMN = MIN( M, N )
 332
 333       WANTU = LSAME( JOBU, 'V' )
 334       WANTVT = LSAME( JOBVT, 'V' )
 335       IF( WANTU .OR. WANTVT ) THEN
 336          JOBZ = 'V'
 337       ELSE
 338          JOBZ = 'N'
 339       END IF
 340       ALLS = LSAME( RANGE, 'A' )
 341       VALS = LSAME( RANGE, 'V' )
 342       INDS = LSAME( RANGE, 'I' )
 343 *
 344       INFO = 0
 345       IF( .NOT.LSAME( JOBU, 'V' ) .AND.
 346      $    .NOT.LSAME( JOBU, 'N' ) ) THEN
 347          INFO = -1
 348       ELSE IF( .NOT.LSAME( JOBVT, 'V' ) .AND.
 349      $         .NOT.LSAME( JOBVT, 'N' ) ) THEN
 350          INFO = -2
 351       ELSE IF( .NOT.( ALLS .OR. VALS .OR. INDS ) ) THEN
 352          INFO = -3
 353       ELSE IF( M.LT.0 ) THEN
 354          INFO = -4
 355       ELSE IF( N.LT.0 ) THEN
 356          INFO = -5
 357       ELSE IF( M.GT.LDA ) THEN
 358          INFO = -7
 359       ELSE IF( MINMN.GT.0 ) THEN
 360          IF( VALS ) THEN
 361             IF( VL.LT.ZERO ) THEN
 362                INFO = -8
 363             ELSE IF( VU.LE.VL ) THEN
 364                INFO = -9
 365             END IF
 366          ELSE IF( INDS ) THEN
 367             IF( IL.LT.1 .OR. IL.GT.MAX( 1, MINMN ) ) THEN
 368                INFO = -10
 369             ELSE IF( IU.LT.MIN( MINMN, IL ) .OR. IU.GT.MINMN ) THEN
 370                INFO = -11
 371             END IF
 372          END IF
 373          IF( INFO.EQ.0 ) THEN
 374             IF( WANTU .AND. LDU.LT.M ) THEN
 375                INFO = -15
 376             ELSE IF( WANTVT ) THEN
 377                IF( INDS ) THEN
 378                    IF( LDVT.LT.IU-IL+1 ) THEN
 379                        INFO = -17
 380                    END IF
 381                ELSE IF( LDVT.LT.MINMN ) THEN
 382                    INFO = -17
 383                END IF
 384             END IF
 385          END IF
 386       END IF
 387 *
 388 *     Compute workspace
 389 *     (Note: Comments in the code beginning "Workspace:" describe the
 390 *     minimal amount of workspace needed at that point in the code,
 391 *     as well as the preferred amount for good performance.
 392 *     NB refers to the optimal block size for the immediately
 393 *     following subroutine, as returned by ILAENV.)
 394 *
 395       IF( INFO.EQ.0 ) THEN
 396          MINWRK = 1
 397          MAXWRK = 1
 398          IF( MINMN.GT.0 ) THEN
 399             IF( M.GE.N ) THEN
 400                MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
 401                IF( M.GE.MNTHR ) THEN
 402 *
 403 *                 Path 1 (M much larger than N)
 404 *
 405                   MINWRK = N*(N+5)
 406                   MAXWRK = N + N*ILAENV(1,'ZGEQRF',' ',M,N,-1,-1)
 407                   MAXWRK = MAX(MAXWRK,
 408      $                     N*N+2*N+2*N*ILAENV(1,'ZGEBRD',' ',N,N,-1,-1))
 409                   IF (WANTU .OR. WANTVT) THEN
 410                      MAXWRK = MAX(MAXWRK,
 411      $                       N*N+2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1))
 412                   END IF
 413                ELSE
 414 *
 415 *                 Path 2 (M at least N, but not much larger)
 416 *
 417                   MINWRK = 3*N + M
 418                   MAXWRK = 2*N + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1)
 419                   IF (WANTU .OR. WANTVT) THEN
 420                      MAXWRK = MAX(MAXWRK,
 421      $                        2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1))
 422                   END IF
 423                END IF
 424             ELSE
 425                MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
 426                IF( N.GE.MNTHR ) THEN
 427 *
 428 *                 Path 1t (N much larger than M)
 429 *
 430                   MINWRK = M*(M+5)
 431                   MAXWRK = M + M*ILAENV(1,'ZGELQF',' ',M,N,-1,-1)
 432                   MAXWRK = MAX(MAXWRK,
 433      $                     M*M+2*M+2*M*ILAENV(1,'ZGEBRD',' ',M,M,-1,-1))
 434                   IF (WANTU .OR. WANTVT) THEN
 435                      MAXWRK = MAX(MAXWRK,
 436      $                       M*M+2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1))
 437                   END IF
 438                ELSE
 439 *
 440 *                 Path 2t (N greater than M, but not much larger)
 441 *
 442 *
 443                   MINWRK = 3*M + N
 444                   MAXWRK = 2*M + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1)
 445                   IF (WANTU .OR. WANTVT) THEN
 446                      MAXWRK = MAX(MAXWRK,
 447      $                        2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1))
 448                   END IF
 449                END IF
 450             END IF
 451          END IF
 452          MAXWRK = MAX( MAXWRK, MINWRK )
 453          WORK( 1 ) = DCMPLX( DBLE( MAXWRK ), ZERO )
 454 *
 455          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 456             INFO = -19
 457          END IF
 458       END IF
 459 *
 460       IF( INFO.NE.0 ) THEN
 461          CALL XERBLA( 'ZGESVDX', -INFO )
 462          RETURN
 463       ELSE IF( LQUERY ) THEN
 464          RETURN
 465       END IF
 466 *
 467 *     Quick return if possible
 468 *
 469       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 470          RETURN
 471       END IF
 472 *
 473 *     Set singular values indices accord to RANGE='A'.
 474 *
 475       IF( ALLS ) THEN
 476          RNGTGK = 'I'
 477          ILTGK = 1
 478          IUTGK = MIN( M, N )
 479       ELSE IF( INDS ) THEN
 480          RNGTGK = 'I'
 481          ILTGK = IL
 482          IUTGK = IU
 483       ELSE
 484          RNGTGK = 'V'
 485          ILTGK = 0
 486          IUTGK = 0
 487       END IF
 488 *
 489 *     Get machine constants
 490 *
 491       EPS = DLAMCH( 'P' )
 492       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
 493       BIGNUM = ONE / SMLNUM
 494 *
 495 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 496 *
 497       ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
 498       ISCL = 0
 499       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 500          ISCL = 1
 501          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
 502       ELSE IF( ANRM.GT.BIGNUM ) THEN
 503          ISCL = 1
 504          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
 505       END IF
 506 *
 507       IF( M.GE.N ) THEN
 508 *
 509 *        A has at least as many rows as columns. If A has sufficiently
 510 *        more rows than columns, first reduce A using the QR
 511 *        decomposition.
 512 *
 513          IF( M.GE.MNTHR ) THEN
 514 *
 515 *           Path 1 (M much larger than N):
 516 *           A = Q * R = Q * ( QB * B * PB**T )
 517 *                     = Q * ( QB * ( UB * S * VB**T ) * PB**T )
 518 *           U = Q * QB * UB; V**T = VB**T * PB**T
 519 *
 520 *           Compute A=Q*R
 521 *           (Workspace: need 2*N, prefer N+N*NB)
 522 *
 523             ITAU = 1
 524             ITEMP = ITAU + N
 525             CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
 526      $                   LWORK-ITEMP+1, INFO )
 527 *
 528 *           Copy R into WORK and bidiagonalize it:
 529 *           (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB)
 530 *
 531             IQRF = ITEMP
 532             ITAUQ = ITEMP + N*N
 533             ITAUP = ITAUQ + N
 534             ITEMP = ITAUP + N
 535             ID = 1
 536             IE = ID + N
 537             ITGKZ = IE + N
 538             CALL ZLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )
 539             CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
 540      $                   WORK( IQRF+1 ), N )
 541             CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ),
 542      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
 543      $                   WORK( ITEMP ), LWORK-ITEMP+1, INFO )
 544             ITEMPR = ITGKZ + N*(N*2+1)
 545 *
 546 *           Solve eigenvalue problem TGK*Z=Z*S.
 547 *           (Workspace: need 2*N*N+14*N)
 548 *
 549             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
 550      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
 551      $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
 552      $                    IWORK, INFO)
 553 *
 554 *           If needed, compute left singular vectors.
 555 *
 556             IF( WANTU ) THEN
 557                K = ITGKZ
 558                DO I = 1, NS
 559                   DO J = 1, N
 560                      U( J, I ) = DCMPLX( RWORK( K ), ZERO )
 561                      K = K + 1
 562                   END DO
 563                   K = K + N
 564                END DO
 565                CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
 566 *
 567 *              Call ZUNMBR to compute QB*UB.
 568 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 569 *
 570                CALL ZUNMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,
 571      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
 572      $                      LWORK-ITEMP+1, INFO )
 573 *
 574 *              Call ZUNMQR to compute Q*(QB*UB).
 575 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 576 *
 577                CALL ZUNMQR( 'L', 'N', M, NS, N, A, LDA,
 578      $                      WORK( ITAU ), U, LDU, WORK( ITEMP ),
 579      $                      LWORK-ITEMP+1, INFO )
 580             END IF
 581 *
 582 *           If needed, compute right singular vectors.
 583 *
 584             IF( WANTVT) THEN
 585                K = ITGKZ + N
 586                DO I = 1, NS
 587                   DO J = 1, N
 588                      VT( I, J ) = DCMPLX( RWORK( K ), ZERO )
 589                      K = K + 1
 590                   END DO
 591                   K = K + N
 592                END DO
 593 *
 594 *              Call ZUNMBR to compute VB**T * PB**T
 595 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 596 *
 597                CALL ZUNMBR( 'P', 'R', 'C', NS, N, N, WORK( IQRF ), N,
 598      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
 599      $                      LWORK-ITEMP+1, INFO )
 600             END IF
 601          ELSE
 602 *
 603 *           Path 2 (M at least N, but not much larger)
 604 *           Reduce A to bidiagonal form without QR decomposition
 605 *           A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
 606 *           U = QB * UB; V**T = VB**T * PB**T
 607 *
 608 *           Bidiagonalize A
 609 *           (Workspace: need 2*N+M, prefer 2*N+(M+N)*NB)
 610 *
 611             ITAUQ = 1
 612             ITAUP = ITAUQ + N
 613             ITEMP = ITAUP + N
 614             ID = 1
 615             IE = ID + N
 616             ITGKZ = IE + N
 617             CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),
 618      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
 619      $                   LWORK-ITEMP+1, INFO )
 620             ITEMPR = ITGKZ + N*(N*2+1)
 621 *
 622 *           Solve eigenvalue problem TGK*Z=Z*S.
 623 *           (Workspace: need 2*N*N+14*N)
 624 *
 625             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
 626      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
 627      $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
 628      $                    IWORK, INFO)
 629 *
 630 *           If needed, compute left singular vectors.
 631 *
 632             IF( WANTU ) THEN
 633                K = ITGKZ
 634                DO I = 1, NS
 635                   DO J = 1, N
 636                      U( J, I ) = DCMPLX( RWORK( K ), ZERO )
 637                      K = K + 1
 638                   END DO
 639                   K = K + N
 640                END DO
 641                CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
 642 *
 643 *              Call ZUNMBR to compute QB*UB.
 644 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 645 *
 646                CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
 647      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
 648      $                      LWORK-ITEMP+1, IERR )
 649             END IF
 650 *
 651 *           If needed, compute right singular vectors.
 652 *
 653             IF( WANTVT) THEN
 654                K = ITGKZ + N
 655                DO I = 1, NS
 656                   DO J = 1, N
 657                      VT( I, J ) = DCMPLX( RWORK( K ), ZERO )
 658                      K = K + 1
 659                   END DO
 660                   K = K + N
 661                END DO
 662 *
 663 *              Call ZUNMBR to compute VB**T * PB**T
 664 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 665 *
 666                CALL ZUNMBR( 'P', 'R', 'C', NS, N, N, A, LDA,
 667      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
 668      $                      LWORK-ITEMP+1, IERR )
 669             END IF
 670          END IF
 671       ELSE
 672 *
 673 *        A has more columns than rows. If A has sufficiently more
 674 *        columns than rows, first reduce A using the LQ decomposition.
 675 *
 676          IF( N.GE.MNTHR ) THEN
 677 *
 678 *           Path 1t (N much larger than M):
 679 *           A = L * Q = ( QB * B * PB**T ) * Q
 680 *                     = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
 681 *           U = QB * UB ; V**T = VB**T * PB**T * Q
 682 *
 683 *           Compute A=L*Q
 684 *           (Workspace: need 2*M, prefer M+M*NB)
 685 *
 686             ITAU = 1
 687             ITEMP = ITAU + M
 688             CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
 689      $                   LWORK-ITEMP+1, INFO )
 690
 691 *           Copy L into WORK and bidiagonalize it:
 692 *           (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB)
 693 *
 694             ILQF = ITEMP
 695             ITAUQ = ILQF + M*M
 696             ITAUP = ITAUQ + M
 697             ITEMP = ITAUP + M
 698             ID = 1
 699             IE = ID + M
 700             ITGKZ = IE + M
 701             CALL ZLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )
 702             CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
 703      $                   WORK( ILQF+M ), M )
 704             CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ),
 705      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
 706      $                   WORK( ITEMP ), LWORK-ITEMP+1, INFO )
 707             ITEMPR = ITGKZ + M*(M*2+1)
 708 *
 709 *           Solve eigenvalue problem TGK*Z=Z*S.
 710 *           (Workspace: need 2*M*M+14*M)
 711 *
 712             CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ),
 713      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
 714      $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),
 715      $                    IWORK, INFO)
 716 *
 717 *           If needed, compute left singular vectors.
 718 *
 719             IF( WANTU ) THEN
 720                K = ITGKZ
 721                DO I = 1, NS
 722                   DO J = 1, M
 723                      U( J, I ) = DCMPLX( RWORK( K ), ZERO )
 724                      K = K + 1
 725                   END DO
 726                   K = K + M
 727                END DO
 728 *
 729 *              Call ZUNMBR to compute QB*UB.
 730 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 731 *
 732                CALL ZUNMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,
 733      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
 734      $                      LWORK-ITEMP+1, INFO )
 735             END IF
 736 *
 737 *           If needed, compute right singular vectors.
 738 *
 739             IF( WANTVT) THEN
 740                K = ITGKZ + M
 741                DO I = 1, NS
 742                   DO J = 1, M
 743                      VT( I, J ) = DCMPLX( RWORK( K ), ZERO )
 744                      K = K + 1
 745                   END DO
 746                   K = K + M
 747                END DO
 748                CALL ZLASET( 'A', NS, N-M, CZERO, CZERO,
 749      $                      VT( 1,M+1 ), LDVT )
 750 *
 751 *              Call ZUNMBR to compute (VB**T)*(PB**T)
 752 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 753 *
 754                CALL ZUNMBR( 'P', 'R', 'C', NS, M, M, WORK( ILQF ), M,
 755      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
 756      $                      LWORK-ITEMP+1, INFO )
 757 *
 758 *              Call ZUNMLQ to compute ((VB**T)*(PB**T))*Q.
 759 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 760 *
 761                CALL ZUNMLQ( 'R', 'N', NS, N, M, A, LDA,
 762      $                      WORK( ITAU ), VT, LDVT, WORK( ITEMP ),
 763      $                      LWORK-ITEMP+1, INFO )
 764             END IF
 765          ELSE
 766 *
 767 *           Path 2t (N greater than M, but not much larger)
 768 *           Reduce to bidiagonal form without LQ decomposition
 769 *           A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
 770 *           U = QB * UB; V**T = VB**T * PB**T
 771 *
 772 *           Bidiagonalize A
 773 *           (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB)
 774 *
 775             ITAUQ = 1
 776             ITAUP = ITAUQ + M
 777             ITEMP = ITAUP + M
 778             ID = 1
 779             IE = ID + M
 780             ITGKZ = IE + M
 781             CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),
 782      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
 783      $                   LWORK-ITEMP+1, INFO )
 784             ITEMPR = ITGKZ + M*(M*2+1)
 785 *
 786 *           Solve eigenvalue problem TGK*Z=Z*S.
 787 *           (Workspace: need 2*M*M+14*M)
 788 *
 789             CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ),
 790      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
 791      $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),
 792      $                    IWORK, INFO)
 793 *
 794 *           If needed, compute left singular vectors.
 795 *
 796             IF( WANTU ) THEN
 797                K = ITGKZ
 798                DO I = 1, NS
 799                   DO J = 1, M
 800                      U( J, I ) = DCMPLX( RWORK( K ), ZERO )
 801                      K = K + 1
 802                   END DO
 803                   K = K + M
 804                END DO
 805 *
 806 *              Call ZUNMBR to compute QB*UB.
 807 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 808 *
 809                CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
 810      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
 811      $                      LWORK-ITEMP+1, INFO )
 812             END IF
 813 *
 814 *           If needed, compute right singular vectors.
 815 *
 816             IF( WANTVT) THEN
 817                K = ITGKZ + M
 818                DO I = 1, NS
 819                   DO J = 1, M
 820                      VT( I, J ) = DCMPLX( RWORK( K ), ZERO )
 821                      K = K + 1
 822                   END DO
 823                   K = K + M
 824                END DO
 825                CALL ZLASET( 'A', NS, N-M, CZERO, CZERO,
 826      $                      VT( 1,M+1 ), LDVT )
 827 *
 828 *              Call ZUNMBR to compute VB**T * PB**T
 829 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 830 *
 831                CALL ZUNMBR( 'P', 'R', 'C', NS, N, M, A, LDA,
 832      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
 833      $                      LWORK-ITEMP+1, INFO )
 834             END IF
 835          END IF
 836       END IF
 837 *
 838 *     Undo scaling if necessary
 839 *
 840       IF( ISCL.EQ.1 ) THEN
 841          IF( ANRM.GT.BIGNUM )
 842      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1,
 843      $                   S, MINMN, INFO )
 844          IF( ANRM.LT.SMLNUM )
 845      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1,
 846      $                   S, MINMN, INFO )
 847       END IF
 848 *
 849 *     Return optimal workspace in WORK(1)
 850 *
 851       WORK( 1 ) = DCMPLX( DBLE( MAXWRK ), ZERO )
 852 *
 853       RETURN
 854 *
 855 *     End of ZGESVDX
 856 *
 857       END