lapack-netlib/SRC/sgesvdx.f

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