1*> \brief \b DBDSQR 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DBDSQR + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, 22* LDU, C, LDC, WORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER UPLO 26* INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), 30* $ VT( LDVT, * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DBDSQR computes the singular values and, optionally, the right and/or 40*> left singular vectors from the singular value decomposition (SVD) of 41*> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit 42*> zero-shift QR algorithm. The SVD of B has the form 43*> 44*> B = Q * S * P**T 45*> 46*> where S is the diagonal matrix of singular values, Q is an orthogonal 47*> matrix of left singular vectors, and P is an orthogonal matrix of 48*> right singular vectors. If left singular vectors are requested, this 49*> subroutine actually returns U*Q instead of Q, and, if right singular 50*> vectors are requested, this subroutine returns P**T*VT instead of 51*> P**T, for given real input matrices U and VT. When U and VT are the 52*> orthogonal matrices that reduce a general matrix A to bidiagonal 53*> form: A = U*B*VT, as computed by DGEBRD, then 54*> 55*> A = (U*Q) * S * (P**T*VT) 56*> 57*> is the SVD of A. Optionally, the subroutine may also compute Q**T*C 58*> for a given real input matrix C. 59*> 60*> See "Computing Small Singular Values of Bidiagonal Matrices With 61*> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, 62*> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, 63*> no. 5, pp. 873-912, Sept 1990) and 64*> "Accurate singular values and differential qd algorithms," by 65*> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics 66*> Department, University of California at Berkeley, July 1992 67*> for a detailed description of the algorithm. 68*> \endverbatim 69* 70* Arguments: 71* ========== 72* 73*> \param[in] UPLO 74*> \verbatim 75*> UPLO is CHARACTER*1 76*> = 'U': B is upper bidiagonal; 77*> = 'L': B is lower bidiagonal. 78*> \endverbatim 79*> 80*> \param[in] N 81*> \verbatim 82*> N is INTEGER 83*> The order of the matrix B. N >= 0. 84*> \endverbatim 85*> 86*> \param[in] NCVT 87*> \verbatim 88*> NCVT is INTEGER 89*> The number of columns of the matrix VT. NCVT >= 0. 90*> \endverbatim 91*> 92*> \param[in] NRU 93*> \verbatim 94*> NRU is INTEGER 95*> The number of rows of the matrix U. NRU >= 0. 96*> \endverbatim 97*> 98*> \param[in] NCC 99*> \verbatim 100*> NCC is INTEGER 101*> The number of columns of the matrix C. NCC >= 0. 102*> \endverbatim 103*> 104*> \param[in,out] D 105*> \verbatim 106*> D is DOUBLE PRECISION array, dimension (N) 107*> On entry, the n diagonal elements of the bidiagonal matrix B. 108*> On exit, if INFO=0, the singular values of B in decreasing 109*> order. 110*> \endverbatim 111*> 112*> \param[in,out] E 113*> \verbatim 114*> E is DOUBLE PRECISION array, dimension (N-1) 115*> On entry, the N-1 offdiagonal elements of the bidiagonal 116*> matrix B. 117*> On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E 118*> will contain the diagonal and superdiagonal elements of a 119*> bidiagonal matrix orthogonally equivalent to the one given 120*> as input. 121*> \endverbatim 122*> 123*> \param[in,out] VT 124*> \verbatim 125*> VT is DOUBLE PRECISION array, dimension (LDVT, NCVT) 126*> On entry, an N-by-NCVT matrix VT. 127*> On exit, VT is overwritten by P**T * VT. 128*> Not referenced if NCVT = 0. 129*> \endverbatim 130*> 131*> \param[in] LDVT 132*> \verbatim 133*> LDVT is INTEGER 134*> The leading dimension of the array VT. 135*> LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. 136*> \endverbatim 137*> 138*> \param[in,out] U 139*> \verbatim 140*> U is DOUBLE PRECISION array, dimension (LDU, N) 141*> On entry, an NRU-by-N matrix U. 142*> On exit, U is overwritten by U * Q. 143*> Not referenced if NRU = 0. 144*> \endverbatim 145*> 146*> \param[in] LDU 147*> \verbatim 148*> LDU is INTEGER 149*> The leading dimension of the array U. LDU >= max(1,NRU). 150*> \endverbatim 151*> 152*> \param[in,out] C 153*> \verbatim 154*> C is DOUBLE PRECISION array, dimension (LDC, NCC) 155*> On entry, an N-by-NCC matrix C. 156*> On exit, C is overwritten by Q**T * C. 157*> Not referenced if NCC = 0. 158*> \endverbatim 159*> 160*> \param[in] LDC 161*> \verbatim 162*> LDC is INTEGER 163*> The leading dimension of the array C. 164*> LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. 165*> \endverbatim 166*> 167*> \param[out] WORK 168*> \verbatim 169*> WORK is DOUBLE PRECISION array, dimension (4*(N-1)) 170*> \endverbatim 171*> 172*> \param[out] INFO 173*> \verbatim 174*> INFO is INTEGER 175*> = 0: successful exit 176*> < 0: If INFO = -i, the i-th argument had an illegal value 177*> > 0: 178*> if NCVT = NRU = NCC = 0, 179*> = 1, a split was marked by a positive value in E 180*> = 2, current block of Z not diagonalized after 30*N 181*> iterations (in inner while loop) 182*> = 3, termination criterion of outer while loop not met 183*> (program created more than N unreduced blocks) 184*> else NCVT = NRU = NCC = 0, 185*> the algorithm did not converge; D and E contain the 186*> elements of a bidiagonal matrix which is orthogonally 187*> similar to the input matrix B; if INFO = i, i 188*> elements of E have not converged to zero. 189*> \endverbatim 190* 191*> \par Internal Parameters: 192* ========================= 193*> 194*> \verbatim 195*> TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8))) 196*> TOLMUL controls the convergence criterion of the QR loop. 197*> If it is positive, TOLMUL*EPS is the desired relative 198*> precision in the computed singular values. 199*> If it is negative, abs(TOLMUL*EPS*sigma_max) is the 200*> desired absolute accuracy in the computed singular 201*> values (corresponds to relative accuracy 202*> abs(TOLMUL*EPS) in the largest singular value. 203*> abs(TOLMUL) should be between 1 and 1/EPS, and preferably 204*> between 10 (for fast convergence) and .1/EPS 205*> (for there to be some accuracy in the results). 206*> Default is to lose at either one eighth or 2 of the 207*> available decimal digits in each computed singular value 208*> (whichever is smaller). 209*> 210*> MAXITR INTEGER, default = 6 211*> MAXITR controls the maximum number of passes of the 212*> algorithm through its inner loop. The algorithms stops 213*> (and so fails to converge) if the number of passes 214*> through the inner loop exceeds MAXITR*N**2. 215*> 216*> \endverbatim 217* 218*> \par Note: 219* =========== 220*> 221*> \verbatim 222*> Bug report from Cezary Dendek. 223*> On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is 224*> removed since it can overflow pretty easily (for N larger or equal 225*> than 18,919). We instead use MAXITDIVN = MAXITR*N. 226*> \endverbatim 227* 228* Authors: 229* ======== 230* 231*> \author Univ. of Tennessee 232*> \author Univ. of California Berkeley 233*> \author Univ. of Colorado Denver 234*> \author NAG Ltd. 235* 236*> \ingroup auxOTHERcomputational 237* 238* ===================================================================== 239 SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, 240 $ LDU, C, LDC, WORK, INFO ) 241* 242* -- LAPACK computational routine -- 243* -- LAPACK is a software package provided by Univ. of Tennessee, -- 244* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 245* 246* .. Scalar Arguments .. 247 CHARACTER UPLO 248 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU 249* .. 250* .. Array Arguments .. 251 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), 252 $ VT( LDVT, * ), WORK( * ) 253* .. 254* 255* ===================================================================== 256* 257* .. Parameters .. 258 DOUBLE PRECISION ZERO 259 PARAMETER ( ZERO = 0.0D0 ) 260 DOUBLE PRECISION ONE 261 PARAMETER ( ONE = 1.0D0 ) 262 DOUBLE PRECISION NEGONE 263 PARAMETER ( NEGONE = -1.0D0 ) 264 DOUBLE PRECISION HNDRTH 265 PARAMETER ( HNDRTH = 0.01D0 ) 266 DOUBLE PRECISION TEN 267 PARAMETER ( TEN = 10.0D0 ) 268 DOUBLE PRECISION HNDRD 269 PARAMETER ( HNDRD = 100.0D0 ) 270 DOUBLE PRECISION MEIGTH 271 PARAMETER ( MEIGTH = -0.125D0 ) 272 INTEGER MAXITR 273 PARAMETER ( MAXITR = 6 ) 274* .. 275* .. Local Scalars .. 276 LOGICAL LOWER, ROTATE 277 INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M, 278 $ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM 279 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, 280 $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, 281 $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA, 282 $ SN, THRESH, TOL, TOLMUL, UNFL 283* .. 284* .. External Functions .. 285 LOGICAL LSAME 286 DOUBLE PRECISION DLAMCH 287 EXTERNAL LSAME, DLAMCH 288* .. 289* .. External Subroutines .. 290 EXTERNAL DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT, 291 $ DSCAL, DSWAP, XERBLA 292* .. 293* .. Intrinsic Functions .. 294 INTRINSIC ABS, DBLE, MAX, MIN, SIGN, SQRT 295* .. 296* .. Executable Statements .. 297* 298* Test the input parameters. 299* 300 INFO = 0 301 LOWER = LSAME( UPLO, 'L' ) 302 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN 303 INFO = -1 304 ELSE IF( N.LT.0 ) THEN 305 INFO = -2 306 ELSE IF( NCVT.LT.0 ) THEN 307 INFO = -3 308 ELSE IF( NRU.LT.0 ) THEN 309 INFO = -4 310 ELSE IF( NCC.LT.0 ) THEN 311 INFO = -5 312 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. 313 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN 314 INFO = -9 315 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN 316 INFO = -11 317 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. 318 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN 319 INFO = -13 320 END IF 321 IF( INFO.NE.0 ) THEN 322 CALL XERBLA( 'DBDSQR', -INFO ) 323 RETURN 324 END IF 325 IF( N.EQ.0 ) 326 $ RETURN 327 IF( N.EQ.1 ) 328 $ GO TO 160 329* 330* ROTATE is true if any singular vectors desired, false otherwise 331* 332 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) 333* 334* If no singular vectors desired, use qd algorithm 335* 336 IF( .NOT.ROTATE ) THEN 337 CALL DLASQ1( N, D, E, WORK, INFO ) 338* 339* If INFO equals 2, dqds didn't finish, try to finish 340* 341 IF( INFO .NE. 2 ) RETURN 342 INFO = 0 343 END IF 344* 345 NM1 = N - 1 346 NM12 = NM1 + NM1 347 NM13 = NM12 + NM1 348 IDIR = 0 349* 350* Get machine constants 351* 352 EPS = DLAMCH( 'Epsilon' ) 353 UNFL = DLAMCH( 'Safe minimum' ) 354* 355* If matrix lower bidiagonal, rotate to be upper bidiagonal 356* by applying Givens rotations on the left 357* 358 IF( LOWER ) THEN 359 DO 10 I = 1, N - 1 360 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 361 D( I ) = R 362 E( I ) = SN*D( I+1 ) 363 D( I+1 ) = CS*D( I+1 ) 364 WORK( I ) = CS 365 WORK( NM1+I ) = SN 366 10 CONTINUE 367* 368* Update singular vectors if desired 369* 370 IF( NRU.GT.0 ) 371 $ CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U, 372 $ LDU ) 373 IF( NCC.GT.0 ) 374 $ CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C, 375 $ LDC ) 376 END IF 377* 378* Compute singular values to relative accuracy TOL 379* (By setting TOL to be negative, algorithm will compute 380* singular values to absolute accuracy ABS(TOL)*norm(input matrix)) 381* 382 TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) ) 383 TOL = TOLMUL*EPS 384* 385* Compute approximate maximum, minimum singular values 386* 387 SMAX = ZERO 388 DO 20 I = 1, N 389 SMAX = MAX( SMAX, ABS( D( I ) ) ) 390 20 CONTINUE 391 DO 30 I = 1, N - 1 392 SMAX = MAX( SMAX, ABS( E( I ) ) ) 393 30 CONTINUE 394 SMINL = ZERO 395 IF( TOL.GE.ZERO ) THEN 396* 397* Relative accuracy desired 398* 399 SMINOA = ABS( D( 1 ) ) 400 IF( SMINOA.EQ.ZERO ) 401 $ GO TO 50 402 MU = SMINOA 403 DO 40 I = 2, N 404 MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) ) 405 SMINOA = MIN( SMINOA, MU ) 406 IF( SMINOA.EQ.ZERO ) 407 $ GO TO 50 408 40 CONTINUE 409 50 CONTINUE 410 SMINOA = SMINOA / SQRT( DBLE( N ) ) 411 THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) ) 412 ELSE 413* 414* Absolute accuracy desired 415* 416 THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) ) 417 END IF 418* 419* Prepare for main iteration loop for the singular values 420* (MAXIT is the maximum number of passes through the inner 421* loop permitted before nonconvergence signalled.) 422* 423 MAXITDIVN = MAXITR*N 424 ITERDIVN = 0 425 ITER = -1 426 OLDLL = -1 427 OLDM = -1 428* 429* M points to last element of unconverged part of matrix 430* 431 M = N 432* 433* Begin main iteration loop 434* 435 60 CONTINUE 436* 437* Check for convergence or exceeding iteration count 438* 439 IF( M.LE.1 ) 440 $ GO TO 160 441* 442 IF( ITER.GE.N ) THEN 443 ITER = ITER - N 444 ITERDIVN = ITERDIVN + 1 445 IF( ITERDIVN.GE.MAXITDIVN ) 446 $ GO TO 200 447 END IF 448* 449* Find diagonal block of matrix to work on 450* 451 IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH ) 452 $ D( M ) = ZERO 453 SMAX = ABS( D( M ) ) 454 SMIN = SMAX 455 DO 70 LLL = 1, M - 1 456 LL = M - LLL 457 ABSS = ABS( D( LL ) ) 458 ABSE = ABS( E( LL ) ) 459 IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH ) 460 $ D( LL ) = ZERO 461 IF( ABSE.LE.THRESH ) 462 $ GO TO 80 463 SMIN = MIN( SMIN, ABSS ) 464 SMAX = MAX( SMAX, ABSS, ABSE ) 465 70 CONTINUE 466 LL = 0 467 GO TO 90 468 80 CONTINUE 469 E( LL ) = ZERO 470* 471* Matrix splits since E(LL) = 0 472* 473 IF( LL.EQ.M-1 ) THEN 474* 475* Convergence of bottom singular value, return to top of loop 476* 477 M = M - 1 478 GO TO 60 479 END IF 480 90 CONTINUE 481 LL = LL + 1 482* 483* E(LL) through E(M-1) are nonzero, E(LL-1) is zero 484* 485 IF( LL.EQ.M-1 ) THEN 486* 487* 2 by 2 block, handle separately 488* 489 CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR, 490 $ COSR, SINL, COSL ) 491 D( M-1 ) = SIGMX 492 E( M-1 ) = ZERO 493 D( M ) = SIGMN 494* 495* Compute singular vectors, if desired 496* 497 IF( NCVT.GT.0 ) 498 $ CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR, 499 $ SINR ) 500 IF( NRU.GT.0 ) 501 $ CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL ) 502 IF( NCC.GT.0 ) 503 $ CALL DROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL, 504 $ SINL ) 505 M = M - 2 506 GO TO 60 507 END IF 508* 509* If working on new submatrix, choose shift direction 510* (from larger end diagonal element towards smaller) 511* 512 IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN 513 IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN 514* 515* Chase bulge from top (big end) to bottom (small end) 516* 517 IDIR = 1 518 ELSE 519* 520* Chase bulge from bottom (big end) to top (small end) 521* 522 IDIR = 2 523 END IF 524 END IF 525* 526* Apply convergence tests 527* 528 IF( IDIR.EQ.1 ) THEN 529* 530* Run convergence test in forward direction 531* First apply standard test to bottom of matrix 532* 533 IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR. 534 $ ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN 535 E( M-1 ) = ZERO 536 GO TO 60 537 END IF 538* 539 IF( TOL.GE.ZERO ) THEN 540* 541* If relative accuracy desired, 542* apply convergence criterion forward 543* 544 MU = ABS( D( LL ) ) 545 SMINL = MU 546 DO 100 LLL = LL, M - 1 547 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN 548 E( LLL ) = ZERO 549 GO TO 60 550 END IF 551 MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) 552 SMINL = MIN( SMINL, MU ) 553 100 CONTINUE 554 END IF 555* 556 ELSE 557* 558* Run convergence test in backward direction 559* First apply standard test to top of matrix 560* 561 IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR. 562 $ ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN 563 E( LL ) = ZERO 564 GO TO 60 565 END IF 566* 567 IF( TOL.GE.ZERO ) THEN 568* 569* If relative accuracy desired, 570* apply convergence criterion backward 571* 572 MU = ABS( D( M ) ) 573 SMINL = MU 574 DO 110 LLL = M - 1, LL, -1 575 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN 576 E( LLL ) = ZERO 577 GO TO 60 578 END IF 579 MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) 580 SMINL = MIN( SMINL, MU ) 581 110 CONTINUE 582 END IF 583 END IF 584 OLDLL = LL 585 OLDM = M 586* 587* Compute shift. First, test if shifting would ruin relative 588* accuracy, and if so set the shift to zero. 589* 590 IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE. 591 $ MAX( EPS, HNDRTH*TOL ) ) THEN 592* 593* Use a zero shift to avoid loss of relative accuracy 594* 595 SHIFT = ZERO 596 ELSE 597* 598* Compute the shift from 2-by-2 block at end of matrix 599* 600 IF( IDIR.EQ.1 ) THEN 601 SLL = ABS( D( LL ) ) 602 CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R ) 603 ELSE 604 SLL = ABS( D( M ) ) 605 CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R ) 606 END IF 607* 608* Test if shift negligible, and if so set to zero 609* 610 IF( SLL.GT.ZERO ) THEN 611 IF( ( SHIFT / SLL )**2.LT.EPS ) 612 $ SHIFT = ZERO 613 END IF 614 END IF 615* 616* Increment iteration count 617* 618 ITER = ITER + M - LL 619* 620* If SHIFT = 0, do simplified QR iteration 621* 622 IF( SHIFT.EQ.ZERO ) THEN 623 IF( IDIR.EQ.1 ) THEN 624* 625* Chase bulge from top to bottom 626* Save cosines and sines for later singular vector updates 627* 628 CS = ONE 629 OLDCS = ONE 630 DO 120 I = LL, M - 1 631 CALL DLARTG( D( I )*CS, E( I ), CS, SN, R ) 632 IF( I.GT.LL ) 633 $ E( I-1 ) = OLDSN*R 634 CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) ) 635 WORK( I-LL+1 ) = CS 636 WORK( I-LL+1+NM1 ) = SN 637 WORK( I-LL+1+NM12 ) = OLDCS 638 WORK( I-LL+1+NM13 ) = OLDSN 639 120 CONTINUE 640 H = D( M )*CS 641 D( M ) = H*OLDCS 642 E( M-1 ) = H*OLDSN 643* 644* Update singular vectors 645* 646 IF( NCVT.GT.0 ) 647 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), 648 $ WORK( N ), VT( LL, 1 ), LDVT ) 649 IF( NRU.GT.0 ) 650 $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), 651 $ WORK( NM13+1 ), U( 1, LL ), LDU ) 652 IF( NCC.GT.0 ) 653 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), 654 $ WORK( NM13+1 ), C( LL, 1 ), LDC ) 655* 656* Test convergence 657* 658 IF( ABS( E( M-1 ) ).LE.THRESH ) 659 $ E( M-1 ) = ZERO 660* 661 ELSE 662* 663* Chase bulge from bottom to top 664* Save cosines and sines for later singular vector updates 665* 666 CS = ONE 667 OLDCS = ONE 668 DO 130 I = M, LL + 1, -1 669 CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R ) 670 IF( I.LT.M ) 671 $ E( I ) = OLDSN*R 672 CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) ) 673 WORK( I-LL ) = CS 674 WORK( I-LL+NM1 ) = -SN 675 WORK( I-LL+NM12 ) = OLDCS 676 WORK( I-LL+NM13 ) = -OLDSN 677 130 CONTINUE 678 H = D( LL )*CS 679 D( LL ) = H*OLDCS 680 E( LL ) = H*OLDSN 681* 682* Update singular vectors 683* 684 IF( NCVT.GT.0 ) 685 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), 686 $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) 687 IF( NRU.GT.0 ) 688 $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), 689 $ WORK( N ), U( 1, LL ), LDU ) 690 IF( NCC.GT.0 ) 691 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ), 692 $ WORK( N ), C( LL, 1 ), LDC ) 693* 694* Test convergence 695* 696 IF( ABS( E( LL ) ).LE.THRESH ) 697 $ E( LL ) = ZERO 698 END IF 699 ELSE 700* 701* Use nonzero shift 702* 703 IF( IDIR.EQ.1 ) THEN 704* 705* Chase bulge from top to bottom 706* Save cosines and sines for later singular vector updates 707* 708 F = ( ABS( D( LL ) )-SHIFT )* 709 $ ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) ) 710 G = E( LL ) 711 DO 140 I = LL, M - 1 712 CALL DLARTG( F, G, COSR, SINR, R ) 713 IF( I.GT.LL ) 714 $ E( I-1 ) = R 715 F = COSR*D( I ) + SINR*E( I ) 716 E( I ) = COSR*E( I ) - SINR*D( I ) 717 G = SINR*D( I+1 ) 718 D( I+1 ) = COSR*D( I+1 ) 719 CALL DLARTG( F, G, COSL, SINL, R ) 720 D( I ) = R 721 F = COSL*E( I ) + SINL*D( I+1 ) 722 D( I+1 ) = COSL*D( I+1 ) - SINL*E( I ) 723 IF( I.LT.M-1 ) THEN 724 G = SINL*E( I+1 ) 725 E( I+1 ) = COSL*E( I+1 ) 726 END IF 727 WORK( I-LL+1 ) = COSR 728 WORK( I-LL+1+NM1 ) = SINR 729 WORK( I-LL+1+NM12 ) = COSL 730 WORK( I-LL+1+NM13 ) = SINL 731 140 CONTINUE 732 E( M-1 ) = F 733* 734* Update singular vectors 735* 736 IF( NCVT.GT.0 ) 737 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), 738 $ WORK( N ), VT( LL, 1 ), LDVT ) 739 IF( NRU.GT.0 ) 740 $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), 741 $ WORK( NM13+1 ), U( 1, LL ), LDU ) 742 IF( NCC.GT.0 ) 743 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), 744 $ WORK( NM13+1 ), C( LL, 1 ), LDC ) 745* 746* Test convergence 747* 748 IF( ABS( E( M-1 ) ).LE.THRESH ) 749 $ E( M-1 ) = ZERO 750* 751 ELSE 752* 753* Chase bulge from bottom to top 754* Save cosines and sines for later singular vector updates 755* 756 F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT / 757 $ D( M ) ) 758 G = E( M-1 ) 759 DO 150 I = M, LL + 1, -1 760 CALL DLARTG( F, G, COSR, SINR, R ) 761 IF( I.LT.M ) 762 $ E( I ) = R 763 F = COSR*D( I ) + SINR*E( I-1 ) 764 E( I-1 ) = COSR*E( I-1 ) - SINR*D( I ) 765 G = SINR*D( I-1 ) 766 D( I-1 ) = COSR*D( I-1 ) 767 CALL DLARTG( F, G, COSL, SINL, R ) 768 D( I ) = R 769 F = COSL*E( I-1 ) + SINL*D( I-1 ) 770 D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 ) 771 IF( I.GT.LL+1 ) THEN 772 G = SINL*E( I-2 ) 773 E( I-2 ) = COSL*E( I-2 ) 774 END IF 775 WORK( I-LL ) = COSR 776 WORK( I-LL+NM1 ) = -SINR 777 WORK( I-LL+NM12 ) = COSL 778 WORK( I-LL+NM13 ) = -SINL 779 150 CONTINUE 780 E( LL ) = F 781* 782* Test convergence 783* 784 IF( ABS( E( LL ) ).LE.THRESH ) 785 $ E( LL ) = ZERO 786* 787* Update singular vectors if desired 788* 789 IF( NCVT.GT.0 ) 790 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), 791 $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) 792 IF( NRU.GT.0 ) 793 $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), 794 $ WORK( N ), U( 1, LL ), LDU ) 795 IF( NCC.GT.0 ) 796 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ), 797 $ WORK( N ), C( LL, 1 ), LDC ) 798 END IF 799 END IF 800* 801* QR iteration finished, go back and check convergence 802* 803 GO TO 60 804* 805* All singular values converged, so make them positive 806* 807 160 CONTINUE 808 DO 170 I = 1, N 809 IF( D( I ).LT.ZERO ) THEN 810 D( I ) = -D( I ) 811* 812* Change sign of singular vectors, if desired 813* 814 IF( NCVT.GT.0 ) 815 $ CALL DSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT ) 816 END IF 817 170 CONTINUE 818* 819* Sort the singular values into decreasing order (insertion sort on 820* singular values, but only one transposition per singular vector) 821* 822 DO 190 I = 1, N - 1 823* 824* Scan for smallest D(I) 825* 826 ISUB = 1 827 SMIN = D( 1 ) 828 DO 180 J = 2, N + 1 - I 829 IF( D( J ).LE.SMIN ) THEN 830 ISUB = J 831 SMIN = D( J ) 832 END IF 833 180 CONTINUE 834 IF( ISUB.NE.N+1-I ) THEN 835* 836* Swap singular values and vectors 837* 838 D( ISUB ) = D( N+1-I ) 839 D( N+1-I ) = SMIN 840 IF( NCVT.GT.0 ) 841 $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ), 842 $ LDVT ) 843 IF( NRU.GT.0 ) 844 $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 ) 845 IF( NCC.GT.0 ) 846 $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC ) 847 END IF 848 190 CONTINUE 849 GO TO 220 850* 851* Maximum number of iterations exceeded, failure to converge 852* 853 200 CONTINUE 854 INFO = 0 855 DO 210 I = 1, N - 1 856 IF( E( I ).NE.ZERO ) 857 $ INFO = INFO + 1 858 210 CONTINUE 859 220 CONTINUE 860 RETURN 861* 862* End of DBDSQR 863* 864 END 865