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