1*> \brief <b> SGELSS solves overdetermined or underdetermined systems 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 SGELSS + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelss.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelss.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelss.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 22* WORK, LWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 26* REAL RCOND 27* .. 28* .. Array Arguments .. 29* REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> SGELSS computes the minimum norm solution to a real linear least 39*> squares problem: 40*> 41*> Minimize 2-norm(| b - A*x |). 42*> 43*> using the singular value decomposition (SVD) of A. A is an M-by-N 44*> matrix which may be rank-deficient. 45*> 46*> Several right hand side vectors b and solution vectors x can be 47*> handled in a single call; they are stored as the columns of the 48*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix 49*> X. 50*> 51*> The effective rank of A is determined by treating as zero those 52*> singular values which are less than RCOND times the largest singular 53*> value. 54*> \endverbatim 55* 56* Arguments: 57* ========== 58* 59*> \param[in] M 60*> \verbatim 61*> M is INTEGER 62*> The number of rows of the matrix A. M >= 0. 63*> \endverbatim 64*> 65*> \param[in] N 66*> \verbatim 67*> N is INTEGER 68*> The number of columns of the matrix A. N >= 0. 69*> \endverbatim 70*> 71*> \param[in] NRHS 72*> \verbatim 73*> NRHS is INTEGER 74*> The number of right hand sides, i.e., the number of columns 75*> of the matrices B and X. NRHS >= 0. 76*> \endverbatim 77*> 78*> \param[in,out] A 79*> \verbatim 80*> A is REAL array, dimension (LDA,N) 81*> On entry, the M-by-N matrix A. 82*> On exit, the first min(m,n) rows of A are overwritten with 83*> its right singular vectors, stored rowwise. 84*> \endverbatim 85*> 86*> \param[in] LDA 87*> \verbatim 88*> LDA is INTEGER 89*> The leading dimension of the array A. LDA >= max(1,M). 90*> \endverbatim 91*> 92*> \param[in,out] B 93*> \verbatim 94*> B is REAL array, dimension (LDB,NRHS) 95*> On entry, the M-by-NRHS right hand side matrix B. 96*> On exit, B is overwritten by the N-by-NRHS solution 97*> matrix X. If m >= n and RANK = n, the residual 98*> sum-of-squares for the solution in the i-th column is given 99*> by the sum of squares of elements n+1:m in that column. 100*> \endverbatim 101*> 102*> \param[in] LDB 103*> \verbatim 104*> LDB is INTEGER 105*> The leading dimension of the array B. LDB >= max(1,max(M,N)). 106*> \endverbatim 107*> 108*> \param[out] S 109*> \verbatim 110*> S is REAL array, dimension (min(M,N)) 111*> The singular values of A in decreasing order. 112*> The condition number of A in the 2-norm = S(1)/S(min(m,n)). 113*> \endverbatim 114*> 115*> \param[in] RCOND 116*> \verbatim 117*> RCOND is REAL 118*> RCOND is used to determine the effective rank of A. 119*> Singular values S(i) <= RCOND*S(1) are treated as zero. 120*> If RCOND < 0, machine precision is used instead. 121*> \endverbatim 122*> 123*> \param[out] RANK 124*> \verbatim 125*> RANK is INTEGER 126*> The effective rank of A, i.e., the number of singular values 127*> which are greater than RCOND*S(1). 128*> \endverbatim 129*> 130*> \param[out] WORK 131*> \verbatim 132*> WORK is REAL array, dimension (MAX(1,LWORK)) 133*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 134*> \endverbatim 135*> 136*> \param[in] LWORK 137*> \verbatim 138*> LWORK is INTEGER 139*> The dimension of the array WORK. LWORK >= 1, and also: 140*> LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS ) 141*> For good performance, LWORK should generally be larger. 142*> 143*> If LWORK = -1, then a workspace query is assumed; the routine 144*> only calculates the optimal size of the WORK array, returns 145*> this value as the first entry of the WORK array, and no error 146*> message related to LWORK is issued by XERBLA. 147*> \endverbatim 148*> 149*> \param[out] INFO 150*> \verbatim 151*> INFO is INTEGER 152*> = 0: successful exit 153*> < 0: if INFO = -i, the i-th argument had an illegal value. 154*> > 0: the algorithm for computing the SVD failed to converge; 155*> if INFO = i, i off-diagonal elements of an intermediate 156*> bidiagonal form did not converge to zero. 157*> \endverbatim 158* 159* Authors: 160* ======== 161* 162*> \author Univ. of Tennessee 163*> \author Univ. of California Berkeley 164*> \author Univ. of Colorado Denver 165*> \author NAG Ltd. 166* 167*> \date December 2016 168* 169*> \ingroup realGEsolve 170* 171* ===================================================================== 172 SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 173 $ WORK, LWORK, INFO ) 174* 175* -- LAPACK driver routine (version 3.7.0) -- 176* -- LAPACK is a software package provided by Univ. of Tennessee, -- 177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 178* December 2016 179* 180* .. Scalar Arguments .. 181 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 182 REAL RCOND 183* .. 184* .. Array Arguments .. 185 REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 186* .. 187* 188* ===================================================================== 189* 190* .. Parameters .. 191 REAL ZERO, ONE 192 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 193* .. 194* .. Local Scalars .. 195 LOGICAL LQUERY 196 INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL, 197 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, 198 $ MAXWRK, MINMN, MINWRK, MM, MNTHR 199 INTEGER LWORK_SGEQRF, LWORK_SORMQR, LWORK_SGEBRD, 200 $ LWORK_SORMBR, LWORK_SORGBR, LWORK_SORMLQ 201 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR 202* .. 203* .. Local Arrays .. 204 REAL DUM( 1 ) 205* .. 206* .. External Subroutines .. 207 EXTERNAL SBDSQR, SCOPY, SGEBRD, SGELQF, SGEMM, SGEMV, 208 $ SGEQRF, SLABAD, SLACPY, SLASCL, SLASET, SORGBR, 209 $ SORMBR, SORMLQ, SORMQR, SRSCL, XERBLA 210* .. 211* .. External Functions .. 212 INTEGER ILAENV 213 REAL SLAMCH, SLANGE 214 EXTERNAL ILAENV, SLAMCH, SLANGE 215* .. 216* .. Intrinsic Functions .. 217 INTRINSIC MAX, MIN 218* .. 219* .. Executable Statements .. 220* 221* Test the input arguments 222* 223 INFO = 0 224 MINMN = MIN( M, N ) 225 MAXMN = MAX( M, N ) 226 LQUERY = ( LWORK.EQ.-1 ) 227 IF( M.LT.0 ) THEN 228 INFO = -1 229 ELSE IF( N.LT.0 ) THEN 230 INFO = -2 231 ELSE IF( NRHS.LT.0 ) THEN 232 INFO = -3 233 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 234 INFO = -5 235 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 236 INFO = -7 237 END IF 238* 239* Compute workspace 240* (Note: Comments in the code beginning "Workspace:" describe the 241* minimal amount of workspace needed at that point in the code, 242* as well as the preferred amount for good performance. 243* NB refers to the optimal block size for the immediately 244* following subroutine, as returned by ILAENV.) 245* 246 IF( INFO.EQ.0 ) THEN 247 MINWRK = 1 248 MAXWRK = 1 249 IF( MINMN.GT.0 ) THEN 250 MM = M 251 MNTHR = ILAENV( 6, 'SGELSS', ' ', M, N, NRHS, -1 ) 252 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 253* 254* Path 1a - overdetermined, with many more rows than 255* columns 256* 257* Compute space needed for SGEQRF 258 CALL SGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) 259 LWORK_SGEQRF=DUM(1) 260* Compute space needed for SORMQR 261 CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B, 262 $ LDB, DUM(1), -1, INFO ) 263 LWORK_SORMQR=DUM(1) 264 MM = N 265 MAXWRK = MAX( MAXWRK, N + LWORK_SGEQRF ) 266 MAXWRK = MAX( MAXWRK, N + LWORK_SORMQR ) 267 END IF 268 IF( M.GE.N ) THEN 269* 270* Path 1 - overdetermined or exactly determined 271* 272* Compute workspace needed for SBDSQR 273* 274 BDSPAC = MAX( 1, 5*N ) 275* Compute space needed for SGEBRD 276 CALL SGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), 277 $ DUM(1), DUM(1), -1, INFO ) 278 LWORK_SGEBRD=DUM(1) 279* Compute space needed for SORMBR 280 CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1), 281 $ B, LDB, DUM(1), -1, INFO ) 282 LWORK_SORMBR=DUM(1) 283* Compute space needed for SORGBR 284 CALL SORGBR( 'P', N, N, N, A, LDA, DUM(1), 285 $ DUM(1), -1, INFO ) 286 LWORK_SORGBR=DUM(1) 287* Compute total workspace needed 288 MAXWRK = MAX( MAXWRK, 3*N + LWORK_SGEBRD ) 289 MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORMBR ) 290 MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORGBR ) 291 MAXWRK = MAX( MAXWRK, BDSPAC ) 292 MAXWRK = MAX( MAXWRK, N*NRHS ) 293 MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC ) 294 MAXWRK = MAX( MINWRK, MAXWRK ) 295 END IF 296 IF( N.GT.M ) THEN 297* 298* Compute workspace needed for SBDSQR 299* 300 BDSPAC = MAX( 1, 5*M ) 301 MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC ) 302 IF( N.GE.MNTHR ) THEN 303* 304* Path 2a - underdetermined, with many more columns 305* than rows 306* 307* Compute space needed for SGEBRD 308 CALL SGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), 309 $ DUM(1), DUM(1), -1, INFO ) 310 LWORK_SGEBRD=DUM(1) 311* Compute space needed for SORMBR 312 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, 313 $ DUM(1), B, LDB, DUM(1), -1, INFO ) 314 LWORK_SORMBR=DUM(1) 315* Compute space needed for SORGBR 316 CALL SORGBR( 'P', M, M, M, A, LDA, DUM(1), 317 $ DUM(1), -1, INFO ) 318 LWORK_SORGBR=DUM(1) 319* Compute space needed for SORMLQ 320 CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1), 321 $ B, LDB, DUM(1), -1, INFO ) 322 LWORK_SORMLQ=DUM(1) 323* Compute total workspace needed 324 MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, 325 $ -1 ) 326 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SGEBRD ) 327 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORMBR ) 328 MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORGBR ) 329 MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC ) 330 IF( NRHS.GT.1 ) THEN 331 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 332 ELSE 333 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 334 END IF 335 MAXWRK = MAX( MAXWRK, M + LWORK_SORMLQ ) 336 ELSE 337* 338* Path 2 - underdetermined 339* 340* Compute space needed for SGEBRD 341 CALL SGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), 342 $ DUM(1), DUM(1), -1, INFO ) 343 LWORK_SGEBRD=DUM(1) 344* Compute space needed for SORMBR 345 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA, 346 $ DUM(1), B, LDB, DUM(1), -1, INFO ) 347 LWORK_SORMBR=DUM(1) 348* Compute space needed for SORGBR 349 CALL SORGBR( 'P', M, N, M, A, LDA, DUM(1), 350 $ DUM(1), -1, INFO ) 351 LWORK_SORGBR=DUM(1) 352 MAXWRK = 3*M + LWORK_SGEBRD 353 MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORMBR ) 354 MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORGBR ) 355 MAXWRK = MAX( MAXWRK, BDSPAC ) 356 MAXWRK = MAX( MAXWRK, N*NRHS ) 357 END IF 358 END IF 359 MAXWRK = MAX( MINWRK, MAXWRK ) 360 END IF 361 WORK( 1 ) = MAXWRK 362* 363 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) 364 $ INFO = -12 365 END IF 366* 367 IF( INFO.NE.0 ) THEN 368 CALL XERBLA( 'SGELSS', -INFO ) 369 RETURN 370 ELSE IF( LQUERY ) THEN 371 RETURN 372 END IF 373* 374* Quick return if possible 375* 376 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 377 RANK = 0 378 RETURN 379 END IF 380* 381* Get machine parameters 382* 383 EPS = SLAMCH( 'P' ) 384 SFMIN = SLAMCH( 'S' ) 385 SMLNUM = SFMIN / EPS 386 BIGNUM = ONE / SMLNUM 387 CALL SLABAD( SMLNUM, BIGNUM ) 388* 389* Scale A if max element outside range [SMLNUM,BIGNUM] 390* 391 ANRM = SLANGE( 'M', M, N, A, LDA, WORK ) 392 IASCL = 0 393 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 394* 395* Scale matrix norm up to SMLNUM 396* 397 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 398 IASCL = 1 399 ELSE IF( ANRM.GT.BIGNUM ) THEN 400* 401* Scale matrix norm down to BIGNUM 402* 403 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 404 IASCL = 2 405 ELSE IF( ANRM.EQ.ZERO ) THEN 406* 407* Matrix all zero. Return zero solution. 408* 409 CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 410 CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN ) 411 RANK = 0 412 GO TO 70 413 END IF 414* 415* Scale B if max element outside range [SMLNUM,BIGNUM] 416* 417 BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK ) 418 IBSCL = 0 419 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 420* 421* Scale matrix norm up to SMLNUM 422* 423 CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 424 IBSCL = 1 425 ELSE IF( BNRM.GT.BIGNUM ) THEN 426* 427* Scale matrix norm down to BIGNUM 428* 429 CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 430 IBSCL = 2 431 END IF 432* 433* Overdetermined case 434* 435 IF( M.GE.N ) THEN 436* 437* Path 1 - overdetermined or exactly determined 438* 439 MM = M 440 IF( M.GE.MNTHR ) THEN 441* 442* Path 1a - overdetermined, with many more rows than columns 443* 444 MM = N 445 ITAU = 1 446 IWORK = ITAU + N 447* 448* Compute A=Q*R 449* (Workspace: need 2*N, prefer N+N*NB) 450* 451 CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 452 $ LWORK-IWORK+1, INFO ) 453* 454* Multiply B by transpose(Q) 455* (Workspace: need N+NRHS, prefer N+NRHS*NB) 456* 457 CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 458 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 459* 460* Zero out below R 461* 462 IF( N.GT.1 ) 463 $ CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 464 END IF 465* 466 IE = 1 467 ITAUQ = IE + N 468 ITAUP = ITAUQ + N 469 IWORK = ITAUP + N 470* 471* Bidiagonalize R in A 472* (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 473* 474 CALL SGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 475 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 476 $ INFO ) 477* 478* Multiply B by transpose of left bidiagonalizing vectors of R 479* (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) 480* 481 CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 482 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 483* 484* Generate right bidiagonalizing vectors of R in A 485* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) 486* 487 CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 488 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 489 IWORK = IE + N 490* 491* Perform bidiagonal QR iteration 492* multiply B by transpose of left singular vectors 493* compute right singular vectors in A 494* (Workspace: need BDSPAC) 495* 496 CALL SBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, 497 $ 1, B, LDB, WORK( IWORK ), INFO ) 498 IF( INFO.NE.0 ) 499 $ GO TO 70 500* 501* Multiply B by reciprocals of singular values 502* 503 THR = MAX( RCOND*S( 1 ), SFMIN ) 504 IF( RCOND.LT.ZERO ) 505 $ THR = MAX( EPS*S( 1 ), SFMIN ) 506 RANK = 0 507 DO 10 I = 1, N 508 IF( S( I ).GT.THR ) THEN 509 CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 510 RANK = RANK + 1 511 ELSE 512 CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 513 END IF 514 10 CONTINUE 515* 516* Multiply B by right singular vectors 517* (Workspace: need N, prefer N*NRHS) 518* 519 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 520 CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO, 521 $ WORK, LDB ) 522 CALL SLACPY( 'G', N, NRHS, WORK, LDB, B, LDB ) 523 ELSE IF( NRHS.GT.1 ) THEN 524 CHUNK = LWORK / N 525 DO 20 I = 1, NRHS, CHUNK 526 BL = MIN( NRHS-I+1, CHUNK ) 527 CALL SGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ), 528 $ LDB, ZERO, WORK, N ) 529 CALL SLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB ) 530 20 CONTINUE 531 ELSE 532 CALL SGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 533 CALL SCOPY( N, WORK, 1, B, 1 ) 534 END IF 535* 536 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 537 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN 538* 539* Path 2a - underdetermined, with many more columns than rows 540* and sufficient workspace for an efficient algorithm 541* 542 LDWORK = M 543 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 544 $ M*LDA+M+M*NRHS ) )LDWORK = LDA 545 ITAU = 1 546 IWORK = M + 1 547* 548* Compute A=L*Q 549* (Workspace: need 2*M, prefer M+M*NB) 550* 551 CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 552 $ LWORK-IWORK+1, INFO ) 553 IL = IWORK 554* 555* Copy L to WORK(IL), zeroing out above it 556* 557 CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 558 CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 559 $ LDWORK ) 560 IE = IL + LDWORK*M 561 ITAUQ = IE + M 562 ITAUP = ITAUQ + M 563 IWORK = ITAUP + M 564* 565* Bidiagonalize L in WORK(IL) 566* (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 567* 568 CALL SGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 569 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ), 570 $ LWORK-IWORK+1, INFO ) 571* 572* Multiply B by transpose of left bidiagonalizing vectors of L 573* (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 574* 575 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 576 $ WORK( ITAUQ ), B, LDB, WORK( IWORK ), 577 $ LWORK-IWORK+1, INFO ) 578* 579* Generate right bidiagonalizing vectors of R in WORK(IL) 580* (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB) 581* 582 CALL SORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ), 583 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 584 IWORK = IE + M 585* 586* Perform bidiagonal QR iteration, 587* computing right singular vectors of L in WORK(IL) and 588* multiplying B by transpose of left singular vectors 589* (Workspace: need M*M+M+BDSPAC) 590* 591 CALL SBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ), 592 $ LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO ) 593 IF( INFO.NE.0 ) 594 $ GO TO 70 595* 596* Multiply B by reciprocals of singular values 597* 598 THR = MAX( RCOND*S( 1 ), SFMIN ) 599 IF( RCOND.LT.ZERO ) 600 $ THR = MAX( EPS*S( 1 ), SFMIN ) 601 RANK = 0 602 DO 30 I = 1, M 603 IF( S( I ).GT.THR ) THEN 604 CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 605 RANK = RANK + 1 606 ELSE 607 CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 608 END IF 609 30 CONTINUE 610 IWORK = IE 611* 612* Multiply B by right singular vectors of L in WORK(IL) 613* (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS) 614* 615 IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN 616 CALL SGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK, 617 $ B, LDB, ZERO, WORK( IWORK ), LDB ) 618 CALL SLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB ) 619 ELSE IF( NRHS.GT.1 ) THEN 620 CHUNK = ( LWORK-IWORK+1 ) / M 621 DO 40 I = 1, NRHS, CHUNK 622 BL = MIN( NRHS-I+1, CHUNK ) 623 CALL SGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK, 624 $ B( 1, I ), LDB, ZERO, WORK( IWORK ), M ) 625 CALL SLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ), 626 $ LDB ) 627 40 CONTINUE 628 ELSE 629 CALL SGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ), 630 $ 1, ZERO, WORK( IWORK ), 1 ) 631 CALL SCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 ) 632 END IF 633* 634* Zero out below first M rows of B 635* 636 CALL SLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 637 IWORK = ITAU + M 638* 639* Multiply transpose(Q) by B 640* (Workspace: need M+NRHS, prefer M+NRHS*NB) 641* 642 CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 643 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 644* 645 ELSE 646* 647* Path 2 - remaining underdetermined cases 648* 649 IE = 1 650 ITAUQ = IE + M 651 ITAUP = ITAUQ + M 652 IWORK = ITAUP + M 653* 654* Bidiagonalize A 655* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 656* 657 CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 658 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 659 $ INFO ) 660* 661* Multiply B by transpose of left bidiagonalizing vectors 662* (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 663* 664 CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 665 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 666* 667* Generate right bidiagonalizing vectors in A 668* (Workspace: need 4*M, prefer 3*M+M*NB) 669* 670 CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 671 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 672 IWORK = IE + M 673* 674* Perform bidiagonal QR iteration, 675* computing right singular vectors of A in A and 676* multiplying B by transpose of left singular vectors 677* (Workspace: need BDSPAC) 678* 679 CALL SBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, 680 $ 1, B, LDB, WORK( IWORK ), INFO ) 681 IF( INFO.NE.0 ) 682 $ GO TO 70 683* 684* Multiply B by reciprocals of singular values 685* 686 THR = MAX( RCOND*S( 1 ), SFMIN ) 687 IF( RCOND.LT.ZERO ) 688 $ THR = MAX( EPS*S( 1 ), SFMIN ) 689 RANK = 0 690 DO 50 I = 1, M 691 IF( S( I ).GT.THR ) THEN 692 CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 693 RANK = RANK + 1 694 ELSE 695 CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 696 END IF 697 50 CONTINUE 698* 699* Multiply B by right singular vectors of A 700* (Workspace: need N, prefer N*NRHS) 701* 702 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 703 CALL SGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO, 704 $ WORK, LDB ) 705 CALL SLACPY( 'F', N, NRHS, WORK, LDB, B, LDB ) 706 ELSE IF( NRHS.GT.1 ) THEN 707 CHUNK = LWORK / N 708 DO 60 I = 1, NRHS, CHUNK 709 BL = MIN( NRHS-I+1, CHUNK ) 710 CALL SGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ), 711 $ LDB, ZERO, WORK, N ) 712 CALL SLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB ) 713 60 CONTINUE 714 ELSE 715 CALL SGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 ) 716 CALL SCOPY( N, WORK, 1, B, 1 ) 717 END IF 718 END IF 719* 720* Undo scaling 721* 722 IF( IASCL.EQ.1 ) THEN 723 CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 724 CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 725 $ INFO ) 726 ELSE IF( IASCL.EQ.2 ) THEN 727 CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 728 CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 729 $ INFO ) 730 END IF 731 IF( IBSCL.EQ.1 ) THEN 732 CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 733 ELSE IF( IBSCL.EQ.2 ) THEN 734 CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 735 END IF 736* 737 70 CONTINUE 738 WORK( 1 ) = MAXWRK 739 RETURN 740* 741* End of SGELSS 742* 743 END 744