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