lapack-netlib/SRC/cgesvdx.f

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