1*> \brief \b DLAQR5 performs a single small-bulge multi-shift QR sweep. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLAQR5 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr5.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr5.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr5.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, 22* SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, 23* LDU, NV, WV, LDWV, NH, WH, LDWH ) 24* 25* .. Scalar Arguments .. 26* INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, 27* $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV 28* LOGICAL WANTT, WANTZ 29* .. 30* .. Array Arguments .. 31* DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), U( LDU, * ), 32* $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ), 33* $ Z( LDZ, * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> DLAQR5, called by DLAQR0, performs a 43*> single small-bulge multi-shift QR sweep. 44*> \endverbatim 45* 46* Arguments: 47* ========== 48* 49*> \param[in] WANTT 50*> \verbatim 51*> WANTT is LOGICAL 52*> WANTT = .true. if the quasi-triangular Schur factor 53*> is being computed. WANTT is set to .false. otherwise. 54*> \endverbatim 55*> 56*> \param[in] WANTZ 57*> \verbatim 58*> WANTZ is LOGICAL 59*> WANTZ = .true. if the orthogonal Schur factor is being 60*> computed. WANTZ is set to .false. otherwise. 61*> \endverbatim 62*> 63*> \param[in] KACC22 64*> \verbatim 65*> KACC22 is INTEGER with value 0, 1, or 2. 66*> Specifies the computation mode of far-from-diagonal 67*> orthogonal updates. 68*> = 0: DLAQR5 does not accumulate reflections and does not 69*> use matrix-matrix multiply to update far-from-diagonal 70*> matrix entries. 71*> = 1: DLAQR5 accumulates reflections and uses matrix-matrix 72*> multiply to update the far-from-diagonal matrix entries. 73*> = 2: Same as KACC22 = 1. This option used to enable exploiting 74*> the 2-by-2 structure during matrix multiplications, but 75*> this is no longer supported. 76*> \endverbatim 77*> 78*> \param[in] N 79*> \verbatim 80*> N is INTEGER 81*> N is the order of the Hessenberg matrix H upon which this 82*> subroutine operates. 83*> \endverbatim 84*> 85*> \param[in] KTOP 86*> \verbatim 87*> KTOP is INTEGER 88*> \endverbatim 89*> 90*> \param[in] KBOT 91*> \verbatim 92*> KBOT is INTEGER 93*> These are the first and last rows and columns of an 94*> isolated diagonal block upon which the QR sweep is to be 95*> applied. It is assumed without a check that 96*> either KTOP = 1 or H(KTOP,KTOP-1) = 0 97*> and 98*> either KBOT = N or H(KBOT+1,KBOT) = 0. 99*> \endverbatim 100*> 101*> \param[in] NSHFTS 102*> \verbatim 103*> NSHFTS is INTEGER 104*> NSHFTS gives the number of simultaneous shifts. NSHFTS 105*> must be positive and even. 106*> \endverbatim 107*> 108*> \param[in,out] SR 109*> \verbatim 110*> SR is DOUBLE PRECISION array, dimension (NSHFTS) 111*> \endverbatim 112*> 113*> \param[in,out] SI 114*> \verbatim 115*> SI is DOUBLE PRECISION array, dimension (NSHFTS) 116*> SR contains the real parts and SI contains the imaginary 117*> parts of the NSHFTS shifts of origin that define the 118*> multi-shift QR sweep. On output SR and SI may be 119*> reordered. 120*> \endverbatim 121*> 122*> \param[in,out] H 123*> \verbatim 124*> H is DOUBLE PRECISION array, dimension (LDH,N) 125*> On input H contains a Hessenberg matrix. On output a 126*> multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied 127*> to the isolated diagonal block in rows and columns KTOP 128*> through KBOT. 129*> \endverbatim 130*> 131*> \param[in] LDH 132*> \verbatim 133*> LDH is INTEGER 134*> LDH is the leading dimension of H just as declared in the 135*> calling procedure. LDH >= MAX(1,N). 136*> \endverbatim 137*> 138*> \param[in] ILOZ 139*> \verbatim 140*> ILOZ is INTEGER 141*> \endverbatim 142*> 143*> \param[in] IHIZ 144*> \verbatim 145*> IHIZ is INTEGER 146*> Specify the rows of Z to which transformations must be 147*> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N 148*> \endverbatim 149*> 150*> \param[in,out] Z 151*> \verbatim 152*> Z is DOUBLE PRECISION array, dimension (LDZ,IHIZ) 153*> If WANTZ = .TRUE., then the QR Sweep orthogonal 154*> similarity transformation is accumulated into 155*> Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. 156*> If WANTZ = .FALSE., then Z is unreferenced. 157*> \endverbatim 158*> 159*> \param[in] LDZ 160*> \verbatim 161*> LDZ is INTEGER 162*> LDA is the leading dimension of Z just as declared in 163*> the calling procedure. LDZ >= N. 164*> \endverbatim 165*> 166*> \param[out] V 167*> \verbatim 168*> V is DOUBLE PRECISION array, dimension (LDV,NSHFTS/2) 169*> \endverbatim 170*> 171*> \param[in] LDV 172*> \verbatim 173*> LDV is INTEGER 174*> LDV is the leading dimension of V as declared in the 175*> calling procedure. LDV >= 3. 176*> \endverbatim 177*> 178*> \param[out] U 179*> \verbatim 180*> U is DOUBLE PRECISION array, dimension (LDU,2*NSHFTS) 181*> \endverbatim 182*> 183*> \param[in] LDU 184*> \verbatim 185*> LDU is INTEGER 186*> LDU is the leading dimension of U just as declared in the 187*> in the calling subroutine. LDU >= 2*NSHFTS. 188*> \endverbatim 189*> 190*> \param[in] NV 191*> \verbatim 192*> NV is INTEGER 193*> NV is the number of rows in WV agailable for workspace. 194*> NV >= 1. 195*> \endverbatim 196*> 197*> \param[out] WV 198*> \verbatim 199*> WV is DOUBLE PRECISION array, dimension (LDWV,2*NSHFTS) 200*> \endverbatim 201*> 202*> \param[in] LDWV 203*> \verbatim 204*> LDWV is INTEGER 205*> LDWV is the leading dimension of WV as declared in the 206*> in the calling subroutine. LDWV >= NV. 207*> \endverbatim 208* 209*> \param[in] NH 210*> \verbatim 211*> NH is INTEGER 212*> NH is the number of columns in array WH available for 213*> workspace. NH >= 1. 214*> \endverbatim 215*> 216*> \param[out] WH 217*> \verbatim 218*> WH is DOUBLE PRECISION array, dimension (LDWH,NH) 219*> \endverbatim 220*> 221*> \param[in] LDWH 222*> \verbatim 223*> LDWH is INTEGER 224*> Leading dimension of WH just as declared in the 225*> calling procedure. LDWH >= 2*NSHFTS. 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 doubleOTHERauxiliary 237* 238*> \par Contributors: 239* ================== 240*> 241*> Karen Braman and Ralph Byers, Department of Mathematics, 242*> University of Kansas, USA 243*> 244*> Lars Karlsson, Daniel Kressner, and Bruno Lang 245*> 246*> Thijs Steel, Department of Computer science, 247*> KU Leuven, Belgium 248* 249*> \par References: 250* ================ 251*> 252*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 253*> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 254*> Performance, SIAM Journal of Matrix Analysis, volume 23, pages 255*> 929--947, 2002. 256*> 257*> Lars Karlsson, Daniel Kressner, and Bruno Lang, Optimally packed 258*> chains of bulges in multishift QR algorithms. 259*> ACM Trans. Math. Softw. 40, 2, Article 12 (February 2014). 260*> 261* ===================================================================== 262 SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, 263 $ SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, 264 $ LDU, NV, WV, LDWV, NH, WH, LDWH ) 265 IMPLICIT NONE 266* 267* -- LAPACK auxiliary routine -- 268* -- LAPACK is a software package provided by Univ. of Tennessee, -- 269* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 270* 271* .. Scalar Arguments .. 272 INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, 273 $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV 274 LOGICAL WANTT, WANTZ 275* .. 276* .. Array Arguments .. 277 DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), U( LDU, * ), 278 $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ), 279 $ Z( LDZ, * ) 280* .. 281* 282* ================================================================ 283* .. Parameters .. 284 DOUBLE PRECISION ZERO, ONE 285 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 ) 286* .. 287* .. Local Scalars .. 288 DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM, 289 $ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, TST1, TST2, 290 $ ULP 291 INTEGER I, I2, I4, INCOL, J, JBOT, JCOL, JLEN, 292 $ JROW, JTOP, K, K1, KDU, KMS, KRCOL, 293 $ M, M22, MBOT, MTOP, NBMPS, NDCOL, 294 $ NS, NU 295 LOGICAL ACCUM, BMP22 296* .. 297* .. External Functions .. 298 DOUBLE PRECISION DLAMCH 299 EXTERNAL DLAMCH 300* .. 301* .. Intrinsic Functions .. 302* 303 INTRINSIC ABS, DBLE, MAX, MIN, MOD 304* .. 305* .. Local Arrays .. 306 DOUBLE PRECISION VT( 3 ) 307* .. 308* .. External Subroutines .. 309 EXTERNAL DGEMM, DLABAD, DLACPY, DLAQR1, DLARFG, DLASET, 310 $ DTRMM 311* .. 312* .. Executable Statements .. 313* 314* ==== If there are no shifts, then there is nothing to do. ==== 315* 316 IF( NSHFTS.LT.2 ) 317 $ RETURN 318* 319* ==== If the active block is empty or 1-by-1, then there 320* . is nothing to do. ==== 321* 322 IF( KTOP.GE.KBOT ) 323 $ RETURN 324* 325* ==== Shuffle shifts into pairs of real shifts and pairs 326* . of complex conjugate shifts assuming complex 327* . conjugate shifts are already adjacent to one 328* . another. ==== 329* 330 DO 10 I = 1, NSHFTS - 2, 2 331 IF( SI( I ).NE.-SI( I+1 ) ) THEN 332* 333 SWAP = SR( I ) 334 SR( I ) = SR( I+1 ) 335 SR( I+1 ) = SR( I+2 ) 336 SR( I+2 ) = SWAP 337* 338 SWAP = SI( I ) 339 SI( I ) = SI( I+1 ) 340 SI( I+1 ) = SI( I+2 ) 341 SI( I+2 ) = SWAP 342 END IF 343 10 CONTINUE 344* 345* ==== NSHFTS is supposed to be even, but if it is odd, 346* . then simply reduce it by one. The shuffle above 347* . ensures that the dropped shift is real and that 348* . the remaining shifts are paired. ==== 349* 350 NS = NSHFTS - MOD( NSHFTS, 2 ) 351* 352* ==== Machine constants for deflation ==== 353* 354 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) 355 SAFMAX = ONE / SAFMIN 356 CALL DLABAD( SAFMIN, SAFMAX ) 357 ULP = DLAMCH( 'PRECISION' ) 358 SMLNUM = SAFMIN*( DBLE( N ) / ULP ) 359* 360* ==== Use accumulated reflections to update far-from-diagonal 361* . entries ? ==== 362* 363 ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 ) 364* 365* ==== clear trash ==== 366* 367 IF( KTOP+2.LE.KBOT ) 368 $ H( KTOP+2, KTOP ) = ZERO 369* 370* ==== NBMPS = number of 2-shift bulges in the chain ==== 371* 372 NBMPS = NS / 2 373* 374* ==== KDU = width of slab ==== 375* 376 KDU = 4*NBMPS 377* 378* ==== Create and chase chains of NBMPS bulges ==== 379* 380 DO 180 INCOL = KTOP - 2*NBMPS + 1, KBOT - 2, 2*NBMPS 381* 382* JTOP = Index from which updates from the right start. 383* 384 IF( ACCUM ) THEN 385 JTOP = MAX( KTOP, INCOL ) 386 ELSE IF( WANTT ) THEN 387 JTOP = 1 388 ELSE 389 JTOP = KTOP 390 END IF 391* 392 NDCOL = INCOL + KDU 393 IF( ACCUM ) 394 $ CALL DLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU ) 395* 396* ==== Near-the-diagonal bulge chase. The following loop 397* . performs the near-the-diagonal part of a small bulge 398* . multi-shift QR sweep. Each 4*NBMPS column diagonal 399* . chunk extends from column INCOL to column NDCOL 400* . (including both column INCOL and column NDCOL). The 401* . following loop chases a 2*NBMPS+1 column long chain of 402* . NBMPS bulges 2*NBMPS columns to the right. (INCOL 403* . may be less than KTOP and and NDCOL may be greater than 404* . KBOT indicating phantom columns from which to chase 405* . bulges before they are actually introduced or to which 406* . to chase bulges beyond column KBOT.) ==== 407* 408 DO 145 KRCOL = INCOL, MIN( INCOL+2*NBMPS-1, KBOT-2 ) 409* 410* ==== Bulges number MTOP to MBOT are active double implicit 411* . shift bulges. There may or may not also be small 412* . 2-by-2 bulge, if there is room. The inactive bulges 413* . (if any) must wait until the active bulges have moved 414* . down the diagonal to make room. The phantom matrix 415* . paradigm described above helps keep track. ==== 416* 417 MTOP = MAX( 1, ( KTOP-KRCOL ) / 2+1 ) 418 MBOT = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 2 ) 419 M22 = MBOT + 1 420 BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+2*( M22-1 ) ).EQ. 421 $ ( KBOT-2 ) 422* 423* ==== Generate reflections to chase the chain right 424* . one column. (The minimum value of K is KTOP-1.) ==== 425* 426 IF ( BMP22 ) THEN 427* 428* ==== Special case: 2-by-2 reflection at bottom treated 429* . separately ==== 430* 431 K = KRCOL + 2*( M22-1 ) 432 IF( K.EQ.KTOP-1 ) THEN 433 CALL DLAQR1( 2, H( K+1, K+1 ), LDH, SR( 2*M22-1 ), 434 $ SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ), 435 $ V( 1, M22 ) ) 436 BETA = V( 1, M22 ) 437 CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) ) 438 ELSE 439 BETA = H( K+1, K ) 440 V( 2, M22 ) = H( K+2, K ) 441 CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) ) 442 H( K+1, K ) = BETA 443 H( K+2, K ) = ZERO 444 END IF 445 446* 447* ==== Perform update from right within 448* . computational window. ==== 449* 450 DO 30 J = JTOP, MIN( KBOT, K+3 ) 451 REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )* 452 $ H( J, K+2 ) ) 453 H( J, K+1 ) = H( J, K+1 ) - REFSUM 454 H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M22 ) 455 30 CONTINUE 456* 457* ==== Perform update from left within 458* . computational window. ==== 459* 460 IF( ACCUM ) THEN 461 JBOT = MIN( NDCOL, KBOT ) 462 ELSE IF( WANTT ) THEN 463 JBOT = N 464 ELSE 465 JBOT = KBOT 466 END IF 467 DO 40 J = K+1, JBOT 468 REFSUM = V( 1, M22 )*( H( K+1, J )+V( 2, M22 )* 469 $ H( K+2, J ) ) 470 H( K+1, J ) = H( K+1, J ) - REFSUM 471 H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 ) 472 40 CONTINUE 473* 474* ==== The following convergence test requires that 475* . the tradition small-compared-to-nearby-diagonals 476* . criterion and the Ahues & Tisseur (LAWN 122, 1997) 477* . criteria both be satisfied. The latter improves 478* . accuracy in some examples. Falling back on an 479* . alternate convergence criterion when TST1 or TST2 480* . is zero (as done here) is traditional but probably 481* . unnecessary. ==== 482* 483 IF( K.GE.KTOP ) THEN 484 IF( H( K+1, K ).NE.ZERO ) THEN 485 TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) 486 IF( TST1.EQ.ZERO ) THEN 487 IF( K.GE.KTOP+1 ) 488 $ TST1 = TST1 + ABS( H( K, K-1 ) ) 489 IF( K.GE.KTOP+2 ) 490 $ TST1 = TST1 + ABS( H( K, K-2 ) ) 491 IF( K.GE.KTOP+3 ) 492 $ TST1 = TST1 + ABS( H( K, K-3 ) ) 493 IF( K.LE.KBOT-2 ) 494 $ TST1 = TST1 + ABS( H( K+2, K+1 ) ) 495 IF( K.LE.KBOT-3 ) 496 $ TST1 = TST1 + ABS( H( K+3, K+1 ) ) 497 IF( K.LE.KBOT-4 ) 498 $ TST1 = TST1 + ABS( H( K+4, K+1 ) ) 499 END IF 500 IF( ABS( H( K+1, K ) ) 501 $ .LE.MAX( SMLNUM, ULP*TST1 ) ) THEN 502 H12 = MAX( ABS( H( K+1, K ) ), 503 $ ABS( H( K, K+1 ) ) ) 504 H21 = MIN( ABS( H( K+1, K ) ), 505 $ ABS( H( K, K+1 ) ) ) 506 H11 = MAX( ABS( H( K+1, K+1 ) ), 507 $ ABS( H( K, K )-H( K+1, K+1 ) ) ) 508 H22 = MIN( ABS( H( K+1, K+1 ) ), 509 $ ABS( H( K, K )-H( K+1, K+1 ) ) ) 510 SCL = H11 + H12 511 TST2 = H22*( H11 / SCL ) 512* 513 IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. 514 $ MAX( SMLNUM, ULP*TST2 ) ) THEN 515 H( K+1, K ) = ZERO 516 END IF 517 END IF 518 END IF 519 END IF 520* 521* ==== Accumulate orthogonal transformations. ==== 522* 523 IF( ACCUM ) THEN 524 KMS = K - INCOL 525 DO 50 J = MAX( 1, KTOP-INCOL ), KDU 526 REFSUM = V( 1, M22 )*( U( J, KMS+1 )+ 527 $ V( 2, M22 )*U( J, KMS+2 ) ) 528 U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM 529 U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M22 ) 530 50 CONTINUE 531 ELSE IF( WANTZ ) THEN 532 DO 60 J = ILOZ, IHIZ 533 REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )* 534 $ Z( J, K+2 ) ) 535 Z( J, K+1 ) = Z( J, K+1 ) - REFSUM 536 Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M22 ) 537 60 CONTINUE 538 END IF 539 END IF 540* 541* ==== Normal case: Chain of 3-by-3 reflections ==== 542* 543 DO 80 M = MBOT, MTOP, -1 544 K = KRCOL + 2*( M-1 ) 545 IF( K.EQ.KTOP-1 ) THEN 546 CALL DLAQR1( 3, H( KTOP, KTOP ), LDH, SR( 2*M-1 ), 547 $ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ), 548 $ V( 1, M ) ) 549 ALPHA = V( 1, M ) 550 CALL DLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) ) 551 ELSE 552* 553* ==== Perform delayed transformation of row below 554* . Mth bulge. Exploit fact that first two elements 555* . of row are actually zero. ==== 556* 557 REFSUM = V( 1, M )*V( 3, M )*H( K+3, K+2 ) 558 H( K+3, K ) = -REFSUM 559 H( K+3, K+1 ) = -REFSUM*V( 2, M ) 560 H( K+3, K+2 ) = H( K+3, K+2 ) - REFSUM*V( 3, M ) 561* 562* ==== Calculate reflection to move 563* . Mth bulge one step. ==== 564* 565 BETA = H( K+1, K ) 566 V( 2, M ) = H( K+2, K ) 567 V( 3, M ) = H( K+3, K ) 568 CALL DLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) ) 569* 570* ==== A Bulge may collapse because of vigilant 571* . deflation or destructive underflow. In the 572* . underflow case, try the two-small-subdiagonals 573* . trick to try to reinflate the bulge. ==== 574* 575 IF( H( K+3, K ).NE.ZERO .OR. H( K+3, K+1 ).NE. 576 $ ZERO .OR. H( K+3, K+2 ).EQ.ZERO ) THEN 577* 578* ==== Typical case: not collapsed (yet). ==== 579* 580 H( K+1, K ) = BETA 581 H( K+2, K ) = ZERO 582 H( K+3, K ) = ZERO 583 ELSE 584* 585* ==== Atypical case: collapsed. Attempt to 586* . reintroduce ignoring H(K+1,K) and H(K+2,K). 587* . If the fill resulting from the new 588* . reflector is too large, then abandon it. 589* . Otherwise, use the new one. ==== 590* 591 CALL DLAQR1( 3, H( K+1, K+1 ), LDH, SR( 2*M-1 ), 592 $ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ), 593 $ VT ) 594 ALPHA = VT( 1 ) 595 CALL DLARFG( 3, ALPHA, VT( 2 ), 1, VT( 1 ) ) 596 REFSUM = VT( 1 )*( H( K+1, K )+VT( 2 )* 597 $ H( K+2, K ) ) 598* 599 IF( ABS( H( K+2, K )-REFSUM*VT( 2 ) )+ 600 $ ABS( REFSUM*VT( 3 ) ).GT.ULP* 601 $ ( ABS( H( K, K ) )+ABS( H( K+1, 602 $ K+1 ) )+ABS( H( K+2, K+2 ) ) ) ) THEN 603* 604* ==== Starting a new bulge here would 605* . create non-negligible fill. Use 606* . the old one with trepidation. ==== 607* 608 H( K+1, K ) = BETA 609 H( K+2, K ) = ZERO 610 H( K+3, K ) = ZERO 611 ELSE 612* 613* ==== Starting a new bulge here would 614* . create only negligible fill. 615* . Replace the old reflector with 616* . the new one. ==== 617* 618 H( K+1, K ) = H( K+1, K ) - REFSUM 619 H( K+2, K ) = ZERO 620 H( K+3, K ) = ZERO 621 V( 1, M ) = VT( 1 ) 622 V( 2, M ) = VT( 2 ) 623 V( 3, M ) = VT( 3 ) 624 END IF 625 END IF 626 END IF 627* 628* ==== Apply reflection from the right and 629* . the first column of update from the left. 630* . These updates are required for the vigilant 631* . deflation check. We still delay most of the 632* . updates from the left for efficiency. ==== 633* 634 DO 70 J = JTOP, MIN( KBOT, K+3 ) 635 REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )* 636 $ H( J, K+2 )+V( 3, M )*H( J, K+3 ) ) 637 H( J, K+1 ) = H( J, K+1 ) - REFSUM 638 H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M ) 639 H( J, K+3 ) = H( J, K+3 ) - REFSUM*V( 3, M ) 640 70 CONTINUE 641* 642* ==== Perform update from left for subsequent 643* . column. ==== 644* 645 REFSUM = V( 1, M )*( H( K+1, K+1 )+V( 2, M )* 646 $ H( K+2, K+1 )+V( 3, M )*H( K+3, K+1 ) ) 647 H( K+1, K+1 ) = H( K+1, K+1 ) - REFSUM 648 H( K+2, K+1 ) = H( K+2, K+1 ) - REFSUM*V( 2, M ) 649 H( K+3, K+1 ) = H( K+3, K+1 ) - REFSUM*V( 3, M ) 650* 651* ==== The following convergence test requires that 652* . the tradition small-compared-to-nearby-diagonals 653* . criterion and the Ahues & Tisseur (LAWN 122, 1997) 654* . criteria both be satisfied. The latter improves 655* . accuracy in some examples. Falling back on an 656* . alternate convergence criterion when TST1 or TST2 657* . is zero (as done here) is traditional but probably 658* . unnecessary. ==== 659* 660 IF( K.LT.KTOP) 661 $ CYCLE 662 IF( H( K+1, K ).NE.ZERO ) THEN 663 TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) 664 IF( TST1.EQ.ZERO ) THEN 665 IF( K.GE.KTOP+1 ) 666 $ TST1 = TST1 + ABS( H( K, K-1 ) ) 667 IF( K.GE.KTOP+2 ) 668 $ TST1 = TST1 + ABS( H( K, K-2 ) ) 669 IF( K.GE.KTOP+3 ) 670 $ TST1 = TST1 + ABS( H( K, K-3 ) ) 671 IF( K.LE.KBOT-2 ) 672 $ TST1 = TST1 + ABS( H( K+2, K+1 ) ) 673 IF( K.LE.KBOT-3 ) 674 $ TST1 = TST1 + ABS( H( K+3, K+1 ) ) 675 IF( K.LE.KBOT-4 ) 676 $ TST1 = TST1 + ABS( H( K+4, K+1 ) ) 677 END IF 678 IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) 679 $ THEN 680 H12 = MAX( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) 681 H21 = MIN( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) 682 H11 = MAX( ABS( H( K+1, K+1 ) ), 683 $ ABS( H( K, K )-H( K+1, K+1 ) ) ) 684 H22 = MIN( ABS( H( K+1, K+1 ) ), 685 $ ABS( H( K, K )-H( K+1, K+1 ) ) ) 686 SCL = H11 + H12 687 TST2 = H22*( H11 / SCL ) 688* 689 IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. 690 $ MAX( SMLNUM, ULP*TST2 ) ) THEN 691 H( K+1, K ) = ZERO 692 END IF 693 END IF 694 END IF 695 80 CONTINUE 696* 697* ==== Multiply H by reflections from the left ==== 698* 699 IF( ACCUM ) THEN 700 JBOT = MIN( NDCOL, KBOT ) 701 ELSE IF( WANTT ) THEN 702 JBOT = N 703 ELSE 704 JBOT = KBOT 705 END IF 706* 707 DO 100 M = MBOT, MTOP, -1 708 K = KRCOL + 2*( M-1 ) 709 DO 90 J = MAX( KTOP, KRCOL + 2*M ), JBOT 710 REFSUM = V( 1, M )*( H( K+1, J )+V( 2, M )* 711 $ H( K+2, J )+V( 3, M )*H( K+3, J ) ) 712 H( K+1, J ) = H( K+1, J ) - REFSUM 713 H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M ) 714 H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M ) 715 90 CONTINUE 716 100 CONTINUE 717* 718* ==== Accumulate orthogonal transformations. ==== 719* 720 IF( ACCUM ) THEN 721* 722* ==== Accumulate U. (If needed, update Z later 723* . with an efficient matrix-matrix 724* . multiply.) ==== 725* 726 DO 120 M = MBOT, MTOP, -1 727 K = KRCOL + 2*( M-1 ) 728 KMS = K - INCOL 729 I2 = MAX( 1, KTOP-INCOL ) 730 I2 = MAX( I2, KMS-(KRCOL-INCOL)+1 ) 731 I4 = MIN( KDU, KRCOL + 2*( MBOT-1 ) - INCOL + 5 ) 732 DO 110 J = I2, I4 733 REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )* 734 $ U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) ) 735 U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM 736 U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M ) 737 U( J, KMS+3 ) = U( J, KMS+3 ) - REFSUM*V( 3, M ) 738 110 CONTINUE 739 120 CONTINUE 740 ELSE IF( WANTZ ) THEN 741* 742* ==== U is not accumulated, so update Z 743* . now by multiplying by reflections 744* . from the right. ==== 745* 746 DO 140 M = MBOT, MTOP, -1 747 K = KRCOL + 2*( M-1 ) 748 DO 130 J = ILOZ, IHIZ 749 REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )* 750 $ Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) ) 751 Z( J, K+1 ) = Z( J, K+1 ) - REFSUM 752 Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M ) 753 Z( J, K+3 ) = Z( J, K+3 ) - REFSUM*V( 3, M ) 754 130 CONTINUE 755 140 CONTINUE 756 END IF 757* 758* ==== End of near-the-diagonal bulge chase. ==== 759* 760 145 CONTINUE 761* 762* ==== Use U (if accumulated) to update far-from-diagonal 763* . entries in H. If required, use U to update Z as 764* . well. ==== 765* 766 IF( ACCUM ) THEN 767 IF( WANTT ) THEN 768 JTOP = 1 769 JBOT = N 770 ELSE 771 JTOP = KTOP 772 JBOT = KBOT 773 END IF 774 K1 = MAX( 1, KTOP-INCOL ) 775 NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1 776* 777* ==== Horizontal Multiply ==== 778* 779 DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH 780 JLEN = MIN( NH, JBOT-JCOL+1 ) 781 CALL DGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ), 782 $ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH, 783 $ LDWH ) 784 CALL DLACPY( 'ALL', NU, JLEN, WH, LDWH, 785 $ H( INCOL+K1, JCOL ), LDH ) 786 150 CONTINUE 787* 788* ==== Vertical multiply ==== 789* 790 DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV 791 JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW ) 792 CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE, 793 $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ), 794 $ LDU, ZERO, WV, LDWV ) 795 CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV, 796 $ H( JROW, INCOL+K1 ), LDH ) 797 160 CONTINUE 798* 799* ==== Z multiply (also vertical) ==== 800* 801 IF( WANTZ ) THEN 802 DO 170 JROW = ILOZ, IHIZ, NV 803 JLEN = MIN( NV, IHIZ-JROW+1 ) 804 CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE, 805 $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ), 806 $ LDU, ZERO, WV, LDWV ) 807 CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV, 808 $ Z( JROW, INCOL+K1 ), LDZ ) 809 170 CONTINUE 810 END IF 811 END IF 812 180 CONTINUE 813* 814* ==== End of DLAQR5 ==== 815* 816 END 817