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