1*> \brief <b> DGELSD computes the minimum-norm solution to a linear least squares problem 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 DGELSD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 22* WORK, LWORK, IWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 26* DOUBLE PRECISION RCOND 27* .. 28* .. Array Arguments .. 29* INTEGER IWORK( * ) 30* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DGELSD computes the minimum-norm solution to a real linear least 40*> squares problem: 41*> minimize 2-norm(| b - A*x |) 42*> using the singular value decomposition (SVD) of A. A is an M-by-N 43*> matrix which may be rank-deficient. 44*> 45*> Several right hand side vectors b and solution vectors x can be 46*> handled in a single call; they are stored as the columns of the 47*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution 48*> matrix X. 49*> 50*> The problem is solved in three steps: 51*> (1) Reduce the coefficient matrix A to bidiagonal form with 52*> Householder transformations, reducing the original problem 53*> into a "bidiagonal least squares problem" (BLS) 54*> (2) Solve the BLS using a divide and conquer approach. 55*> (3) Apply back all the Householder transformations to solve 56*> the original least squares problem. 57*> 58*> The effective rank of A is determined by treating as zero those 59*> singular values which are less than RCOND times the largest singular 60*> value. 61*> 62*> The divide and conquer algorithm makes very mild assumptions about 63*> floating point arithmetic. It will work on machines with a guard 64*> digit in add/subtract, or on those binary machines without guard 65*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 66*> Cray-2. It could conceivably fail on hexadecimal or decimal machines 67*> without guard digits, but we know of none. 68*> \endverbatim 69* 70* Arguments: 71* ========== 72* 73*> \param[in] M 74*> \verbatim 75*> M is INTEGER 76*> The number of rows of A. M >= 0. 77*> \endverbatim 78*> 79*> \param[in] N 80*> \verbatim 81*> N is INTEGER 82*> The number of columns of A. N >= 0. 83*> \endverbatim 84*> 85*> \param[in] NRHS 86*> \verbatim 87*> NRHS is INTEGER 88*> The number of right hand sides, i.e., the number of columns 89*> of the matrices B and X. NRHS >= 0. 90*> \endverbatim 91*> 92*> \param[in,out] A 93*> \verbatim 94*> A is DOUBLE PRECISION array, dimension (LDA,N) 95*> On entry, the M-by-N matrix A. 96*> On exit, A has been destroyed. 97*> \endverbatim 98*> 99*> \param[in] LDA 100*> \verbatim 101*> LDA is INTEGER 102*> The leading dimension of the array A. LDA >= max(1,M). 103*> \endverbatim 104*> 105*> \param[in,out] B 106*> \verbatim 107*> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 108*> On entry, the M-by-NRHS right hand side matrix B. 109*> On exit, B is overwritten by the N-by-NRHS solution 110*> matrix X. If m >= n and RANK = n, the residual 111*> sum-of-squares for the solution in the i-th column is given 112*> by the sum of squares of elements n+1:m in that column. 113*> \endverbatim 114*> 115*> \param[in] LDB 116*> \verbatim 117*> LDB is INTEGER 118*> The leading dimension of the array B. LDB >= max(1,max(M,N)). 119*> \endverbatim 120*> 121*> \param[out] S 122*> \verbatim 123*> S is DOUBLE PRECISION array, dimension (min(M,N)) 124*> The singular values of A in decreasing order. 125*> The condition number of A in the 2-norm = S(1)/S(min(m,n)). 126*> \endverbatim 127*> 128*> \param[in] RCOND 129*> \verbatim 130*> RCOND is DOUBLE PRECISION 131*> RCOND is used to determine the effective rank of A. 132*> Singular values S(i) <= RCOND*S(1) are treated as zero. 133*> If RCOND < 0, machine precision is used instead. 134*> \endverbatim 135*> 136*> \param[out] RANK 137*> \verbatim 138*> RANK is INTEGER 139*> The effective rank of A, i.e., the number of singular values 140*> which are greater than RCOND*S(1). 141*> \endverbatim 142*> 143*> \param[out] WORK 144*> \verbatim 145*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 146*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 147*> \endverbatim 148*> 149*> \param[in] LWORK 150*> \verbatim 151*> LWORK is INTEGER 152*> The dimension of the array WORK. LWORK must be at least 1. 153*> The exact minimum amount of workspace needed depends on M, 154*> N and NRHS. As long as LWORK is at least 155*> 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2, 156*> if M is greater than or equal to N or 157*> 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2, 158*> if M is less than N, the code will execute correctly. 159*> SMLSIZ is returned by ILAENV and is equal to the maximum 160*> size of the subproblems at the bottom of the computation 161*> tree (usually about 25), and 162*> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 163*> For good performance, LWORK should generally be larger. 164*> 165*> If LWORK = -1, then a workspace query is assumed; the routine 166*> only calculates the optimal size of the WORK array, returns 167*> this value as the first entry of the WORK array, and no error 168*> message related to LWORK is issued by XERBLA. 169*> \endverbatim 170*> 171*> \param[out] IWORK 172*> \verbatim 173*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 174*> LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN), 175*> where MINMN = MIN( M,N ). 176*> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. 177*> \endverbatim 178*> 179*> \param[out] INFO 180*> \verbatim 181*> INFO is INTEGER 182*> = 0: successful exit 183*> < 0: if INFO = -i, the i-th argument had an illegal value. 184*> > 0: the algorithm for computing the SVD failed to converge; 185*> if INFO = i, i off-diagonal elements of an intermediate 186*> bidiagonal form did not converge to zero. 187*> \endverbatim 188* 189* Authors: 190* ======== 191* 192*> \author Univ. of Tennessee 193*> \author Univ. of California Berkeley 194*> \author Univ. of Colorado Denver 195*> \author NAG Ltd. 196* 197*> \ingroup doubleGEsolve 198* 199*> \par Contributors: 200* ================== 201*> 202*> Ming Gu and Ren-Cang Li, Computer Science Division, University of 203*> California at Berkeley, USA \n 204*> Osni Marques, LBNL/NERSC, USA \n 205* 206* ===================================================================== 207 SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 208 $ WORK, LWORK, IWORK, INFO ) 209* 210* -- LAPACK driver routine -- 211* -- LAPACK is a software package provided by Univ. of Tennessee, -- 212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 213* 214* .. Scalar Arguments .. 215 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 216 DOUBLE PRECISION RCOND 217* .. 218* .. Array Arguments .. 219 INTEGER IWORK( * ) 220 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) 221* .. 222* 223* ===================================================================== 224* 225* .. Parameters .. 226 DOUBLE PRECISION ZERO, ONE, TWO 227 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 228* .. 229* .. Local Scalars .. 230 LOGICAL LQUERY 231 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ, 232 $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK, 233 $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD 234 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM 235* .. 236* .. External Subroutines .. 237 EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD, 238 $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA 239* .. 240* .. External Functions .. 241 INTEGER ILAENV 242 DOUBLE PRECISION DLAMCH, DLANGE 243 EXTERNAL ILAENV, DLAMCH, DLANGE 244* .. 245* .. Intrinsic Functions .. 246 INTRINSIC DBLE, INT, LOG, MAX, MIN 247* .. 248* .. Executable Statements .. 249* 250* Test the input arguments. 251* 252 INFO = 0 253 MINMN = MIN( M, N ) 254 MAXMN = MAX( M, N ) 255 MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 ) 256 LQUERY = ( LWORK.EQ.-1 ) 257 IF( M.LT.0 ) THEN 258 INFO = -1 259 ELSE IF( N.LT.0 ) THEN 260 INFO = -2 261 ELSE IF( NRHS.LT.0 ) THEN 262 INFO = -3 263 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 264 INFO = -5 265 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 266 INFO = -7 267 END IF 268* 269 SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 ) 270* 271* Compute workspace. 272* (Note: Comments in the code beginning "Workspace:" describe the 273* minimal amount of workspace needed at that point in the code, 274* as well as the preferred amount for good performance. 275* NB refers to the optimal block size for the immediately 276* following subroutine, as returned by ILAENV.) 277* 278 MINWRK = 1 279 LIWORK = 1 280 MINMN = MAX( 1, MINMN ) 281 NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) / 282 $ LOG( TWO ) ) + 1, 0 ) 283* 284 IF( INFO.EQ.0 ) THEN 285 MAXWRK = 0 286 LIWORK = 3*MINMN*NLVL + 11*MINMN 287 MM = M 288 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 289* 290* Path 1a - overdetermined, with many more rows than columns. 291* 292 MM = N 293 MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N, 294 $ -1, -1 ) ) 295 MAXWRK = MAX( MAXWRK, N+NRHS* 296 $ ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, -1 ) ) 297 END IF 298 IF( M.GE.N ) THEN 299* 300* Path 1 - overdetermined or exactly determined. 301* 302 MAXWRK = MAX( MAXWRK, 3*N+( MM+N )* 303 $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) ) 304 MAXWRK = MAX( MAXWRK, 3*N+NRHS* 305 $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) ) 306 MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* 307 $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) ) 308 WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2 309 MAXWRK = MAX( MAXWRK, 3*N+WLALSD ) 310 MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD ) 311 END IF 312 IF( N.GT.M ) THEN 313 WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2 314 IF( N.GE.MNTHR ) THEN 315* 316* Path 2a - underdetermined, with many more columns 317* than rows. 318* 319 MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 320 MAXWRK = MAX( MAXWRK, M*M+4*M+2*M* 321 $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) 322 MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS* 323 $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) ) 324 MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )* 325 $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) ) 326 IF( NRHS.GT.1 ) THEN 327 MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS ) 328 ELSE 329 MAXWRK = MAX( MAXWRK, M*M+2*M ) 330 END IF 331 MAXWRK = MAX( MAXWRK, M+NRHS* 332 $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) ) 333 MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD ) 334! XXX: Ensure the Path 2a case below is triggered. The workspace 335! calculation should use queries for all routines eventually. 336 MAXWRK = MAX( MAXWRK, 337 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) 338 ELSE 339* 340* Path 2 - remaining underdetermined cases. 341* 342 MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N, 343 $ -1, -1 ) 344 MAXWRK = MAX( MAXWRK, 3*M+NRHS* 345 $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) ) 346 MAXWRK = MAX( MAXWRK, 3*M+M* 347 $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) ) 348 MAXWRK = MAX( MAXWRK, 3*M+WLALSD ) 349 END IF 350 MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD ) 351 END IF 352 MINWRK = MIN( MINWRK, MAXWRK ) 353 WORK( 1 ) = MAXWRK 354 IWORK( 1 ) = LIWORK 355 356 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 357 INFO = -12 358 END IF 359 END IF 360* 361 IF( INFO.NE.0 ) THEN 362 CALL XERBLA( 'DGELSD', -INFO ) 363 RETURN 364 ELSE IF( LQUERY ) THEN 365 GO TO 10 366 END IF 367* 368* Quick return if possible. 369* 370 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 371 RANK = 0 372 RETURN 373 END IF 374* 375* Get machine parameters. 376* 377 EPS = DLAMCH( 'P' ) 378 SFMIN = DLAMCH( 'S' ) 379 SMLNUM = SFMIN / EPS 380 BIGNUM = ONE / SMLNUM 381 CALL DLABAD( SMLNUM, BIGNUM ) 382* 383* Scale A if max entry outside range [SMLNUM,BIGNUM]. 384* 385 ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) 386 IASCL = 0 387 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 388* 389* Scale matrix norm up to SMLNUM. 390* 391 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 392 IASCL = 1 393 ELSE IF( ANRM.GT.BIGNUM ) THEN 394* 395* Scale matrix norm down to BIGNUM. 396* 397 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 398 IASCL = 2 399 ELSE IF( ANRM.EQ.ZERO ) THEN 400* 401* Matrix all zero. Return zero solution. 402* 403 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 404 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 405 RANK = 0 406 GO TO 10 407 END IF 408* 409* Scale B if max entry outside range [SMLNUM,BIGNUM]. 410* 411 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) 412 IBSCL = 0 413 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 414* 415* Scale matrix norm up to SMLNUM. 416* 417 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 418 IBSCL = 1 419 ELSE IF( BNRM.GT.BIGNUM ) THEN 420* 421* Scale matrix norm down to BIGNUM. 422* 423 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 424 IBSCL = 2 425 END IF 426* 427* If M < N make sure certain entries of B are zero. 428* 429 IF( M.LT.N ) 430 $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 431* 432* Overdetermined case. 433* 434 IF( M.GE.N ) THEN 435* 436* Path 1 - overdetermined or exactly determined. 437* 438 MM = M 439 IF( M.GE.MNTHR ) THEN 440* 441* Path 1a - overdetermined, with many more rows than columns. 442* 443 MM = N 444 ITAU = 1 445 NWORK = ITAU + N 446* 447* Compute A=Q*R. 448* (Workspace: need 2*N, prefer N+N*NB) 449* 450 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 451 $ LWORK-NWORK+1, INFO ) 452* 453* Multiply B by transpose(Q). 454* (Workspace: need N+NRHS, prefer N+NRHS*NB) 455* 456 CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, 457 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 458* 459* Zero out below R. 460* 461 IF( N.GT.1 ) THEN 462 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 463 END IF 464 END IF 465* 466 IE = 1 467 ITAUQ = IE + N 468 ITAUP = ITAUQ + N 469 NWORK = ITAUP + N 470* 471* Bidiagonalize R in A. 472* (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) 473* 474 CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 475 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+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 DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 482 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 483* 484* Solve the bidiagonal least squares problem. 485* 486 CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB, 487 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 488 IF( INFO.NE.0 ) THEN 489 GO TO 10 490 END IF 491* 492* Multiply B by right bidiagonalizing vectors of R. 493* 494 CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), 495 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 496* 497 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 498 $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN 499* 500* Path 2a - underdetermined, with many more columns than rows 501* and sufficient workspace for an efficient algorithm. 502* 503 LDWORK = M 504 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 505 $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA 506 ITAU = 1 507 NWORK = M + 1 508* 509* Compute A=L*Q. 510* (Workspace: need 2*M, prefer M+M*NB) 511* 512 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 513 $ LWORK-NWORK+1, INFO ) 514 IL = NWORK 515* 516* Copy L to WORK(IL), zeroing out above its diagonal. 517* 518 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 519 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), 520 $ LDWORK ) 521 IE = IL + LDWORK*M 522 ITAUQ = IE + M 523 ITAUP = ITAUQ + M 524 NWORK = ITAUP + M 525* 526* Bidiagonalize L in WORK(IL). 527* (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) 528* 529 CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), 530 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 531 $ LWORK-NWORK+1, INFO ) 532* 533* Multiply B by transpose of left bidiagonalizing vectors of L. 534* (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 535* 536 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, 537 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), 538 $ LWORK-NWORK+1, INFO ) 539* 540* Solve the bidiagonal least squares problem. 541* 542 CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 543 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 544 IF( INFO.NE.0 ) THEN 545 GO TO 10 546 END IF 547* 548* Multiply B by right bidiagonalizing vectors of L. 549* 550 CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, 551 $ WORK( ITAUP ), B, LDB, WORK( NWORK ), 552 $ LWORK-NWORK+1, INFO ) 553* 554* Zero out below first M rows of B. 555* 556 CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) 557 NWORK = ITAU + M 558* 559* Multiply transpose(Q) by B. 560* (Workspace: need M+NRHS, prefer M+NRHS*NB) 561* 562 CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, 563 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 564* 565 ELSE 566* 567* Path 2 - remaining underdetermined cases. 568* 569 IE = 1 570 ITAUQ = IE + M 571 ITAUP = ITAUQ + M 572 NWORK = ITAUP + M 573* 574* Bidiagonalize A. 575* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) 576* 577 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 578 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 579 $ INFO ) 580* 581* Multiply B by transpose of left bidiagonalizing vectors. 582* (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) 583* 584 CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), 585 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 586* 587* Solve the bidiagonal least squares problem. 588* 589 CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, 590 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) 591 IF( INFO.NE.0 ) THEN 592 GO TO 10 593 END IF 594* 595* Multiply B by right bidiagonalizing vectors of A. 596* 597 CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), 598 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 599* 600 END IF 601* 602* Undo scaling. 603* 604 IF( IASCL.EQ.1 ) THEN 605 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 606 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 607 $ INFO ) 608 ELSE IF( IASCL.EQ.2 ) THEN 609 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 610 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 611 $ INFO ) 612 END IF 613 IF( IBSCL.EQ.1 ) THEN 614 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 615 ELSE IF( IBSCL.EQ.2 ) THEN 616 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 617 END IF 618* 619 10 CONTINUE 620 WORK( 1 ) = MAXWRK 621 IWORK( 1 ) = LIWORK 622 RETURN 623* 624* End of DGELSD 625* 626 END 627