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