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*> \ingroup realGEsing 258* 259* ===================================================================== 260 SUBROUTINE SGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, 261 $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK, 262 $ LWORK, IWORK, INFO ) 263* 264* -- LAPACK driver routine -- 265* -- LAPACK is a software package provided by Univ. of Tennessee, -- 266* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 267* 268* .. Scalar Arguments .. 269 CHARACTER JOBU, JOBVT, RANGE 270 INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS 271 REAL VL, VU 272* .. 273* .. Array Arguments .. 274 INTEGER IWORK( * ) 275 REAL A( LDA, * ), S( * ), U( LDU, * ), 276 $ VT( LDVT, * ), WORK( * ) 277* .. 278* 279* ===================================================================== 280* 281* .. Parameters .. 282 REAL ZERO, ONE 283 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 284* .. 285* .. Local Scalars .. 286 CHARACTER JOBZ, RNGTGK 287 LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT 288 INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL, 289 $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK, 290 $ J, MAXWRK, MINMN, MINWRK, MNTHR 291 REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM 292* .. 293* .. Local Arrays .. 294 REAL DUM( 1 ) 295* .. 296* .. External Subroutines .. 297 EXTERNAL SBDSVDX, SGEBRD, SGELQF, SGEQRF, SLACPY, 298 $ SLASCL, SLASET, SORMBR, SORMLQ, SORMQR, 299 $ SCOPY, XERBLA 300* .. 301* .. External Functions .. 302 LOGICAL LSAME 303 INTEGER ILAENV 304 REAL SLAMCH, SLANGE 305 EXTERNAL LSAME, ILAENV, SLAMCH, SLANGE 306* .. 307* .. Intrinsic Functions .. 308 INTRINSIC MAX, MIN, SQRT 309* .. 310* .. Executable Statements .. 311* 312* Test the input arguments. 313* 314 NS = 0 315 INFO = 0 316 ABSTOL = 2*SLAMCH('S') 317 LQUERY = ( LWORK.EQ.-1 ) 318 MINMN = MIN( M, N ) 319 320 WANTU = LSAME( JOBU, 'V' ) 321 WANTVT = LSAME( JOBVT, 'V' ) 322 IF( WANTU .OR. WANTVT ) THEN 323 JOBZ = 'V' 324 ELSE 325 JOBZ = 'N' 326 END IF 327 ALLS = LSAME( RANGE, 'A' ) 328 VALS = LSAME( RANGE, 'V' ) 329 INDS = LSAME( RANGE, 'I' ) 330* 331 INFO = 0 332 IF( .NOT.LSAME( JOBU, 'V' ) .AND. 333 $ .NOT.LSAME( JOBU, 'N' ) ) THEN 334 INFO = -1 335 ELSE IF( .NOT.LSAME( JOBVT, 'V' ) .AND. 336 $ .NOT.LSAME( JOBVT, 'N' ) ) THEN 337 INFO = -2 338 ELSE IF( .NOT.( ALLS .OR. VALS .OR. INDS ) ) THEN 339 INFO = -3 340 ELSE IF( M.LT.0 ) THEN 341 INFO = -4 342 ELSE IF( N.LT.0 ) THEN 343 INFO = -5 344 ELSE IF( M.GT.LDA ) THEN 345 INFO = -7 346 ELSE IF( MINMN.GT.0 ) THEN 347 IF( VALS ) THEN 348 IF( VL.LT.ZERO ) THEN 349 INFO = -8 350 ELSE IF( VU.LE.VL ) THEN 351 INFO = -9 352 END IF 353 ELSE IF( INDS ) THEN 354 IF( IL.LT.1 .OR. IL.GT.MAX( 1, MINMN ) ) THEN 355 INFO = -10 356 ELSE IF( IU.LT.MIN( MINMN, IL ) .OR. IU.GT.MINMN ) THEN 357 INFO = -11 358 END IF 359 END IF 360 IF( INFO.EQ.0 ) THEN 361 IF( WANTU .AND. LDU.LT.M ) THEN 362 INFO = -15 363 ELSE IF( WANTVT ) THEN 364 IF( INDS ) THEN 365 IF( LDVT.LT.IU-IL+1 ) THEN 366 INFO = -17 367 END IF 368 ELSE IF( LDVT.LT.MINMN ) THEN 369 INFO = -17 370 END IF 371 END IF 372 END IF 373 END IF 374* 375* Compute workspace 376* (Note: Comments in the code beginning "Workspace:" describe the 377* minimal amount of workspace needed at that point in the code, 378* as well as the preferred amount for good performance. 379* NB refers to the optimal block size for the immediately 380* following subroutine, as returned by ILAENV.) 381* 382 IF( INFO.EQ.0 ) THEN 383 MINWRK = 1 384 MAXWRK = 1 385 IF( MINMN.GT.0 ) THEN 386 IF( M.GE.N ) THEN 387 MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 ) 388 IF( M.GE.MNTHR ) THEN 389* 390* Path 1 (M much larger than N) 391* 392 MAXWRK = N + 393 $ N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) 394 MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N* 395 $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) 396 IF (WANTU) THEN 397 MAXWRK = MAX(MAXWRK,N*(N*3+6)+N* 398 $ ILAENV( 1, 'SORMQR', ' ', N, N, -1, -1 ) ) 399 END IF 400 IF (WANTVT) THEN 401 MAXWRK = MAX(MAXWRK,N*(N*3+6)+N* 402 $ ILAENV( 1, 'SORMLQ', ' ', N, N, -1, -1 ) ) 403 END IF 404 MINWRK = N*(N*3+20) 405 ELSE 406* 407* Path 2 (M at least N, but not much larger) 408* 409 MAXWRK = 4*N + ( M+N )* 410 $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) 411 IF (WANTU) THEN 412 MAXWRK = MAX(MAXWRK,N*(N*2+5)+N* 413 $ ILAENV( 1, 'SORMQR', ' ', N, N, -1, -1 ) ) 414 END IF 415 IF (WANTVT) THEN 416 MAXWRK = MAX(MAXWRK,N*(N*2+5)+N* 417 $ ILAENV( 1, 'SORMLQ', ' ', N, N, -1, -1 ) ) 418 END IF 419 MINWRK = MAX(N*(N*2+19),4*N+M) 420 END IF 421 ELSE 422 MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 ) 423 IF( N.GE.MNTHR ) THEN 424* 425* Path 1t (N much larger than M) 426* 427 MAXWRK = M + 428 $ M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) 429 MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M* 430 $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) 431 IF (WANTU) THEN 432 MAXWRK = MAX(MAXWRK,M*(M*3+6)+M* 433 $ ILAENV( 1, 'SORMQR', ' ', M, M, -1, -1 ) ) 434 END IF 435 IF (WANTVT) THEN 436 MAXWRK = MAX(MAXWRK,M*(M*3+6)+M* 437 $ ILAENV( 1, 'SORMLQ', ' ', M, M, -1, -1 ) ) 438 END IF 439 MINWRK = M*(M*3+20) 440 ELSE 441* 442* Path 2t (N at least M, but not much larger) 443* 444 MAXWRK = 4*M + ( M+N )* 445 $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) 446 IF (WANTU) THEN 447 MAXWRK = MAX(MAXWRK,M*(M*2+5)+M* 448 $ ILAENV( 1, 'SORMQR', ' ', M, M, -1, -1 ) ) 449 END IF 450 IF (WANTVT) THEN 451 MAXWRK = MAX(MAXWRK,M*(M*2+5)+M* 452 $ ILAENV( 1, 'SORMLQ', ' ', M, M, -1, -1 ) ) 453 END IF 454 MINWRK = MAX(M*(M*2+19),4*M+N) 455 END IF 456 END IF 457 END IF 458 MAXWRK = MAX( MAXWRK, MINWRK ) 459 WORK( 1 ) = REAL( MAXWRK ) 460* 461 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 462 INFO = -19 463 END IF 464 END IF 465* 466 IF( INFO.NE.0 ) THEN 467 CALL XERBLA( 'SGESVDX', -INFO ) 468 RETURN 469 ELSE IF( LQUERY ) THEN 470 RETURN 471 END IF 472* 473* Quick return if possible 474* 475 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 476 RETURN 477 END IF 478* 479* Set singular values indices accord to RANGE. 480* 481 IF( ALLS ) THEN 482 RNGTGK = 'I' 483 ILTGK = 1 484 IUTGK = MIN( M, N ) 485 ELSE IF( INDS ) THEN 486 RNGTGK = 'I' 487 ILTGK = IL 488 IUTGK = IU 489 ELSE 490 RNGTGK = 'V' 491 ILTGK = 0 492 IUTGK = 0 493 END IF 494* 495* Get machine constants 496* 497 EPS = SLAMCH( 'P' ) 498 SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS 499 BIGNUM = ONE / SMLNUM 500* 501* Scale A if max element outside range [SMLNUM,BIGNUM] 502* 503 ANRM = SLANGE( 'M', M, N, A, LDA, DUM ) 504 ISCL = 0 505 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 506 ISCL = 1 507 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 508 ELSE IF( ANRM.GT.BIGNUM ) THEN 509 ISCL = 1 510 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 511 END IF 512* 513 IF( M.GE.N ) THEN 514* 515* A has at least as many rows as columns. If A has sufficiently 516* more rows than columns, first reduce A using the QR 517* decomposition. 518* 519 IF( M.GE.MNTHR ) THEN 520* 521* Path 1 (M much larger than N): 522* A = Q * R = Q * ( QB * B * PB**T ) 523* = Q * ( QB * ( UB * S * VB**T ) * PB**T ) 524* U = Q * QB * UB; V**T = VB**T * PB**T 525* 526* Compute A=Q*R 527* (Workspace: need 2*N, prefer N+N*NB) 528* 529 ITAU = 1 530 ITEMP = ITAU + N 531 CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ), 532 $ LWORK-ITEMP+1, INFO ) 533* 534* Copy R into WORK and bidiagonalize it: 535* (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB) 536* 537 IQRF = ITEMP 538 ID = IQRF + N*N 539 IE = ID + N 540 ITAUQ = IE + N 541 ITAUP = ITAUQ + N 542 ITEMP = ITAUP + N 543 CALL SLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N ) 544 CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IQRF+1 ), N ) 545 CALL SGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ), 546 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), 547 $ LWORK-ITEMP+1, INFO ) 548* 549* Solve eigenvalue problem TGK*Z=Z*S. 550* (Workspace: need 14*N + 2*N*(N+1)) 551* 552 ITGKZ = ITEMP 553 ITEMP = ITGKZ + N*(N*2+1) 554 CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ), 555 $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), 556 $ N*2, WORK( ITEMP ), IWORK, INFO) 557* 558* If needed, compute left singular vectors. 559* 560 IF( WANTU ) THEN 561 J = ITGKZ 562 DO I = 1, NS 563 CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) 564 J = J + N*2 565 END DO 566 CALL SLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) 567* 568* Call SORMBR to compute QB*UB. 569* (Workspace in WORK( ITEMP ): need N, prefer N*NB) 570* 571 CALL SORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N, 572 $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), 573 $ LWORK-ITEMP+1, INFO ) 574* 575* Call SORMQR to compute Q*(QB*UB). 576* (Workspace in WORK( ITEMP ): need N, prefer N*NB) 577* 578 CALL SORMQR( 'L', 'N', M, NS, N, A, LDA, 579 $ WORK( ITAU ), U, LDU, WORK( ITEMP ), 580 $ LWORK-ITEMP+1, INFO ) 581 END IF 582* 583* If needed, compute right singular vectors. 584* 585 IF( WANTVT) THEN 586 J = ITGKZ + N 587 DO I = 1, NS 588 CALL SCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT ) 589 J = J + N*2 590 END DO 591* 592* Call SORMBR to compute VB**T * PB**T 593* (Workspace in WORK( ITEMP ): need N, prefer N*NB) 594* 595 CALL SORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N, 596 $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), 597 $ LWORK-ITEMP+1, INFO ) 598 END IF 599 ELSE 600* 601* Path 2 (M at least N, but not much larger) 602* Reduce A to bidiagonal form without QR decomposition 603* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T 604* U = QB * UB; V**T = VB**T * PB**T 605* 606* Bidiagonalize A 607* (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB) 608* 609 ID = 1 610 IE = ID + N 611 ITAUQ = IE + N 612 ITAUP = ITAUQ + N 613 ITEMP = ITAUP + N 614 CALL SGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ), 615 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), 616 $ LWORK-ITEMP+1, INFO ) 617* 618* Solve eigenvalue problem TGK*Z=Z*S. 619* (Workspace: need 14*N + 2*N*(N+1)) 620* 621 ITGKZ = ITEMP 622 ITEMP = ITGKZ + N*(N*2+1) 623 CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ), 624 $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), 625 $ N*2, WORK( ITEMP ), IWORK, INFO) 626* 627* If needed, compute left singular vectors. 628* 629 IF( WANTU ) THEN 630 J = ITGKZ 631 DO I = 1, NS 632 CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) 633 J = J + N*2 634 END DO 635 CALL SLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) 636* 637* Call SORMBR to compute QB*UB. 638* (Workspace in WORK( ITEMP ): need N, prefer N*NB) 639* 640 CALL SORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, 641 $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), 642 $ LWORK-ITEMP+1, IERR ) 643 END IF 644* 645* If needed, compute right singular vectors. 646* 647 IF( WANTVT) THEN 648 J = ITGKZ + N 649 DO I = 1, NS 650 CALL SCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT ) 651 J = J + N*2 652 END DO 653* 654* Call SORMBR to compute VB**T * PB**T 655* (Workspace in WORK( ITEMP ): need N, prefer N*NB) 656* 657 CALL SORMBR( 'P', 'R', 'T', NS, N, N, A, LDA, 658 $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), 659 $ LWORK-ITEMP+1, IERR ) 660 END IF 661 END IF 662 ELSE 663* 664* A has more columns than rows. If A has sufficiently more 665* columns than rows, first reduce A using the LQ decomposition. 666* 667 IF( N.GE.MNTHR ) THEN 668* 669* Path 1t (N much larger than M): 670* A = L * Q = ( QB * B * PB**T ) * Q 671* = ( QB * ( UB * S * VB**T ) * PB**T ) * Q 672* U = QB * UB ; V**T = VB**T * PB**T * Q 673* 674* Compute A=L*Q 675* (Workspace: need 2*M, prefer M+M*NB) 676* 677 ITAU = 1 678 ITEMP = ITAU + M 679 CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ), 680 $ LWORK-ITEMP+1, INFO ) 681 682* Copy L into WORK and bidiagonalize it: 683* (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB) 684* 685 ILQF = ITEMP 686 ID = ILQF + M*M 687 IE = ID + M 688 ITAUQ = IE + M 689 ITAUP = ITAUQ + M 690 ITEMP = ITAUP + M 691 CALL SLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M ) 692 CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( ILQF+M ), M ) 693 CALL SGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ), 694 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), 695 $ LWORK-ITEMP+1, INFO ) 696* 697* Solve eigenvalue problem TGK*Z=Z*S. 698* (Workspace: need 2*M*M+14*M) 699* 700 ITGKZ = ITEMP 701 ITEMP = ITGKZ + M*(M*2+1) 702 CALL SBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ), 703 $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), 704 $ M*2, WORK( ITEMP ), IWORK, INFO) 705* 706* If needed, compute left singular vectors. 707* 708 IF( WANTU ) THEN 709 J = ITGKZ 710 DO I = 1, NS 711 CALL SCOPY( M, WORK( J ), 1, U( 1,I ), 1 ) 712 J = J + M*2 713 END DO 714* 715* Call SORMBR to compute QB*UB. 716* (Workspace in WORK( ITEMP ): need M, prefer M*NB) 717* 718 CALL SORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M, 719 $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), 720 $ LWORK-ITEMP+1, INFO ) 721 END IF 722* 723* If needed, compute right singular vectors. 724* 725 IF( WANTVT) THEN 726 J = ITGKZ + M 727 DO I = 1, NS 728 CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) 729 J = J + M*2 730 END DO 731 CALL SLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) 732* 733* Call SORMBR to compute (VB**T)*(PB**T) 734* (Workspace in WORK( ITEMP ): need M, prefer M*NB) 735* 736 CALL SORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M, 737 $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), 738 $ LWORK-ITEMP+1, INFO ) 739* 740* Call SORMLQ to compute ((VB**T)*(PB**T))*Q. 741* (Workspace in WORK( ITEMP ): need M, prefer M*NB) 742* 743 CALL SORMLQ( 'R', 'N', NS, N, M, A, LDA, 744 $ WORK( ITAU ), VT, LDVT, WORK( ITEMP ), 745 $ LWORK-ITEMP+1, INFO ) 746 END IF 747 ELSE 748* 749* Path 2t (N greater than M, but not much larger) 750* Reduce to bidiagonal form without LQ decomposition 751* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T 752* U = QB * UB; V**T = VB**T * PB**T 753* 754* Bidiagonalize A 755* (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB) 756* 757 ID = 1 758 IE = ID + M 759 ITAUQ = IE + M 760 ITAUP = ITAUQ + M 761 ITEMP = ITAUP + M 762 CALL SGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ), 763 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), 764 $ LWORK-ITEMP+1, INFO ) 765* 766* Solve eigenvalue problem TGK*Z=Z*S. 767* (Workspace: need 2*M*M+14*M) 768* 769 ITGKZ = ITEMP 770 ITEMP = ITGKZ + M*(M*2+1) 771 CALL SBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ), 772 $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), 773 $ M*2, WORK( ITEMP ), IWORK, INFO) 774* 775* If needed, compute left singular vectors. 776* 777 IF( WANTU ) THEN 778 J = ITGKZ 779 DO I = 1, NS 780 CALL SCOPY( M, WORK( J ), 1, U( 1,I ), 1 ) 781 J = J + M*2 782 END DO 783* 784* Call SORMBR to compute QB*UB. 785* (Workspace in WORK( ITEMP ): need M, prefer M*NB) 786* 787 CALL SORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, 788 $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), 789 $ LWORK-ITEMP+1, INFO ) 790 END IF 791* 792* If needed, compute right singular vectors. 793* 794 IF( WANTVT) THEN 795 J = ITGKZ + M 796 DO I = 1, NS 797 CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) 798 J = J + M*2 799 END DO 800 CALL SLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) 801* 802* Call SORMBR to compute VB**T * PB**T 803* (Workspace in WORK( ITEMP ): need M, prefer M*NB) 804* 805 CALL SORMBR( 'P', 'R', 'T', NS, N, M, A, LDA, 806 $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), 807 $ LWORK-ITEMP+1, INFO ) 808 END IF 809 END IF 810 END IF 811* 812* Undo scaling if necessary 813* 814 IF( ISCL.EQ.1 ) THEN 815 IF( ANRM.GT.BIGNUM ) 816 $ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, 817 $ S, MINMN, INFO ) 818 IF( ANRM.LT.SMLNUM ) 819 $ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, 820 $ S, MINMN, INFO ) 821 END IF 822* 823* Return optimal workspace in WORK(1) 824* 825 WORK( 1 ) = REAL( MAXWRK ) 826* 827 RETURN 828* 829* End of SGESVDX 830* 831 END 832