1*> \brief \b ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLAQR4 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr4.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr4.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr4.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 22* IHIZ, Z, LDZ, WORK, LWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N 26* LOGICAL WANTT, WANTZ 27* .. 28* .. Array Arguments .. 29* COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> ZLAQR4 implements one level of recursion for ZLAQR0. 39*> It is a complete implementation of the small bulge multi-shift 40*> QR algorithm. It may be called by ZLAQR0 and, for large enough 41*> deflation window size, it may be called by ZLAQR3. This 42*> subroutine is identical to ZLAQR0 except that it calls ZLAQR2 43*> instead of ZLAQR3. 44*> 45*> ZLAQR4 computes the eigenvalues of a Hessenberg matrix H 46*> and, optionally, the matrices T and Z from the Schur decomposition 47*> H = Z T Z**H, where T is an upper triangular matrix (the 48*> Schur form), and Z is the unitary matrix of Schur vectors. 49*> 50*> Optionally Z may be postmultiplied into an input unitary 51*> matrix Q so that this routine can give the Schur factorization 52*> of a matrix A which has been reduced to the Hessenberg form H 53*> by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H. 54*> \endverbatim 55* 56* Arguments: 57* ========== 58* 59*> \param[in] WANTT 60*> \verbatim 61*> WANTT is LOGICAL 62*> = .TRUE. : the full Schur form T is required; 63*> = .FALSE.: only eigenvalues are required. 64*> \endverbatim 65*> 66*> \param[in] WANTZ 67*> \verbatim 68*> WANTZ is LOGICAL 69*> = .TRUE. : the matrix of Schur vectors Z is required; 70*> = .FALSE.: Schur vectors are not required. 71*> \endverbatim 72*> 73*> \param[in] N 74*> \verbatim 75*> N is INTEGER 76*> The order of the matrix H. N >= 0. 77*> \endverbatim 78*> 79*> \param[in] ILO 80*> \verbatim 81*> ILO is INTEGER 82*> \endverbatim 83*> 84*> \param[in] IHI 85*> \verbatim 86*> IHI is INTEGER 87*> It is assumed that H is already upper triangular in rows 88*> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1, 89*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a 90*> previous call to ZGEBAL, and then passed to ZGEHRD when the 91*> matrix output by ZGEBAL is reduced to Hessenberg form. 92*> Otherwise, ILO and IHI should be set to 1 and N, 93*> respectively. If N > 0, then 1 <= ILO <= IHI <= N. 94*> If N = 0, then ILO = 1 and IHI = 0. 95*> \endverbatim 96*> 97*> \param[in,out] H 98*> \verbatim 99*> H is COMPLEX*16 array, dimension (LDH,N) 100*> On entry, the upper Hessenberg matrix H. 101*> On exit, if INFO = 0 and WANTT is .TRUE., then H 102*> contains the upper triangular matrix T from the Schur 103*> decomposition (the Schur form). If INFO = 0 and WANT is 104*> .FALSE., then the contents of H are unspecified on exit. 105*> (The output value of H when INFO > 0 is given under the 106*> description of INFO below.) 107*> 108*> This subroutine may explicitly set H(i,j) = 0 for i > j and 109*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. 110*> \endverbatim 111*> 112*> \param[in] LDH 113*> \verbatim 114*> LDH is INTEGER 115*> The leading dimension of the array H. LDH >= max(1,N). 116*> \endverbatim 117*> 118*> \param[out] W 119*> \verbatim 120*> W is COMPLEX*16 array, dimension (N) 121*> The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored 122*> in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are 123*> stored in the same order as on the diagonal of the Schur 124*> form returned in H, with W(i) = H(i,i). 125*> \endverbatim 126*> 127*> \param[in] ILOZ 128*> \verbatim 129*> ILOZ is INTEGER 130*> \endverbatim 131*> 132*> \param[in] IHIZ 133*> \verbatim 134*> IHIZ is INTEGER 135*> Specify the rows of Z to which transformations must be 136*> applied if WANTZ is .TRUE.. 137*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. 138*> \endverbatim 139*> 140*> \param[in,out] Z 141*> \verbatim 142*> Z is COMPLEX*16 array, dimension (LDZ,IHI) 143*> If WANTZ is .FALSE., then Z is not referenced. 144*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is 145*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the 146*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). 147*> (The output value of Z when INFO > 0 is given under 148*> the description of INFO below.) 149*> \endverbatim 150*> 151*> \param[in] LDZ 152*> \verbatim 153*> LDZ is INTEGER 154*> The leading dimension of the array Z. if WANTZ is .TRUE. 155*> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1. 156*> \endverbatim 157*> 158*> \param[out] WORK 159*> \verbatim 160*> WORK is COMPLEX*16 array, dimension LWORK 161*> On exit, if LWORK = -1, WORK(1) returns an estimate of 162*> the optimal value for LWORK. 163*> \endverbatim 164*> 165*> \param[in] LWORK 166*> \verbatim 167*> LWORK is INTEGER 168*> The dimension of the array WORK. LWORK >= max(1,N) 169*> is sufficient, but LWORK typically as large as 6*N may 170*> be required for optimal performance. A workspace query 171*> to determine the optimal workspace size is recommended. 172*> 173*> If LWORK = -1, then ZLAQR4 does a workspace query. 174*> In this case, ZLAQR4 checks the input parameters and 175*> estimates the optimal workspace size for the given 176*> values of N, ILO and IHI. The estimate is returned 177*> in WORK(1). No error message related to LWORK is 178*> issued by XERBLA. Neither H nor Z are accessed. 179*> \endverbatim 180*> 181*> \param[out] INFO 182*> \verbatim 183*> INFO is INTEGER 184*> = 0: successful exit 185*> > 0: if INFO = i, ZLAQR4 failed to compute all of 186*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR 187*> and WI contain those eigenvalues which have been 188*> successfully computed. (Failures are rare.) 189*> 190*> If INFO > 0 and WANT is .FALSE., then on exit, 191*> the remaining unconverged eigenvalues are the eigen- 192*> values of the upper Hessenberg matrix rows and 193*> columns ILO through INFO of the final, output 194*> value of H. 195*> 196*> If INFO > 0 and WANTT is .TRUE., then on exit 197*> 198*> (*) (initial value of H)*U = U*(final value of H) 199*> 200*> where U is a unitary matrix. The final 201*> value of H is upper Hessenberg and triangular in 202*> rows and columns INFO+1 through IHI. 203*> 204*> If INFO > 0 and WANTZ is .TRUE., then on exit 205*> 206*> (final value of Z(ILO:IHI,ILOZ:IHIZ) 207*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U 208*> 209*> where U is the unitary matrix in (*) (regard- 210*> less of the value of WANTT.) 211*> 212*> If INFO > 0 and WANTZ is .FALSE., then Z is not 213*> accessed. 214*> \endverbatim 215* 216* Authors: 217* ======== 218* 219*> \author Univ. of Tennessee 220*> \author Univ. of California Berkeley 221*> \author Univ. of Colorado Denver 222*> \author NAG Ltd. 223* 224*> \date December 2016 225* 226*> \ingroup complex16OTHERauxiliary 227* 228*> \par Contributors: 229* ================== 230*> 231*> Karen Braman and Ralph Byers, Department of Mathematics, 232*> University of Kansas, USA 233* 234*> \par References: 235* ================ 236*> 237*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 238*> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 239*> Performance, SIAM Journal of Matrix Analysis, volume 23, pages 240*> 929--947, 2002. 241*> \n 242*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 243*> Algorithm Part II: Aggressive Early Deflation, SIAM Journal 244*> of Matrix Analysis, volume 23, pages 948--973, 2002. 245*> 246* ===================================================================== 247 SUBROUTINE ZLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 248 $ IHIZ, Z, LDZ, WORK, LWORK, INFO ) 249* 250* -- LAPACK auxiliary routine (version 3.7.0) -- 251* -- LAPACK is a software package provided by Univ. of Tennessee, -- 252* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 253* December 2016 254* 255* .. Scalar Arguments .. 256 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N 257 LOGICAL WANTT, WANTZ 258* .. 259* .. Array Arguments .. 260 COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) 261* .. 262* 263* ================================================================ 264* 265* .. Parameters .. 266* 267* ==== Matrices of order NTINY or smaller must be processed by 268* . ZLAHQR because of insufficient subdiagonal scratch space. 269* . (This is a hard limit.) ==== 270 INTEGER NTINY 271 PARAMETER ( NTINY = 11 ) 272* 273* ==== Exceptional deflation windows: try to cure rare 274* . slow convergence by varying the size of the 275* . deflation window after KEXNW iterations. ==== 276 INTEGER KEXNW 277 PARAMETER ( KEXNW = 5 ) 278* 279* ==== Exceptional shifts: try to cure rare slow convergence 280* . with ad-hoc exceptional shifts every KEXSH iterations. 281* . ==== 282 INTEGER KEXSH 283 PARAMETER ( KEXSH = 6 ) 284* 285* ==== The constant WILK1 is used to form the exceptional 286* . shifts. ==== 287 DOUBLE PRECISION WILK1 288 PARAMETER ( WILK1 = 0.75d0 ) 289 COMPLEX*16 ZERO, ONE 290 PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ), 291 $ ONE = ( 1.0d0, 0.0d0 ) ) 292 DOUBLE PRECISION TWO 293 PARAMETER ( TWO = 2.0d0 ) 294* .. 295* .. Local Scalars .. 296 COMPLEX*16 AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2 297 DOUBLE PRECISION S 298 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS, 299 $ KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS, 300 $ LWKOPT, NDEC, NDFL, NH, NHO, NIBBLE, NMIN, NS, 301 $ NSMAX, NSR, NVE, NW, NWMAX, NWR, NWUPBD 302 LOGICAL SORTED 303 CHARACTER JBCMPZ*2 304* .. 305* .. External Functions .. 306 INTEGER ILAENV 307 EXTERNAL ILAENV 308* .. 309* .. Local Arrays .. 310 COMPLEX*16 ZDUM( 1, 1 ) 311* .. 312* .. External Subroutines .. 313 EXTERNAL ZLACPY, ZLAHQR, ZLAQR2, ZLAQR5 314* .. 315* .. Intrinsic Functions .. 316 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, MAX, MIN, MOD, 317 $ SQRT 318* .. 319* .. Statement Functions .. 320 DOUBLE PRECISION CABS1 321* .. 322* .. Statement Function definitions .. 323 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 324* .. 325* .. Executable Statements .. 326 INFO = 0 327* 328* ==== Quick return for N = 0: nothing to do. ==== 329* 330 IF( N.EQ.0 ) THEN 331 WORK( 1 ) = ONE 332 RETURN 333 END IF 334* 335 IF( N.LE.NTINY ) THEN 336* 337* ==== Tiny matrices must use ZLAHQR. ==== 338* 339 LWKOPT = 1 340 IF( LWORK.NE.-1 ) 341 $ CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 342 $ IHIZ, Z, LDZ, INFO ) 343 ELSE 344* 345* ==== Use small bulge multi-shift QR with aggressive early 346* . deflation on larger-than-tiny matrices. ==== 347* 348* ==== Hope for the best. ==== 349* 350 INFO = 0 351* 352* ==== Set up job flags for ILAENV. ==== 353* 354 IF( WANTT ) THEN 355 JBCMPZ( 1: 1 ) = 'S' 356 ELSE 357 JBCMPZ( 1: 1 ) = 'E' 358 END IF 359 IF( WANTZ ) THEN 360 JBCMPZ( 2: 2 ) = 'V' 361 ELSE 362 JBCMPZ( 2: 2 ) = 'N' 363 END IF 364* 365* ==== NWR = recommended deflation window size. At this 366* . point, N .GT. NTINY = 11, so there is enough 367* . subdiagonal workspace for NWR.GE.2 as required. 368* . (In fact, there is enough subdiagonal space for 369* . NWR.GE.3.) ==== 370* 371 NWR = ILAENV( 13, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 372 NWR = MAX( 2, NWR ) 373 NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR ) 374* 375* ==== NSR = recommended number of simultaneous shifts. 376* . At this point N .GT. NTINY = 11, so there is at 377* . enough subdiagonal workspace for NSR to be even 378* . and greater than or equal to two as required. ==== 379* 380 NSR = ILAENV( 15, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 381 NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO ) 382 NSR = MAX( 2, NSR-MOD( NSR, 2 ) ) 383* 384* ==== Estimate optimal workspace ==== 385* 386* ==== Workspace query call to ZLAQR2 ==== 387* 388 CALL ZLAQR2( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ, 389 $ IHIZ, Z, LDZ, LS, LD, W, H, LDH, N, H, LDH, N, H, 390 $ LDH, WORK, -1 ) 391* 392* ==== Optimal workspace = MAX(ZLAQR5, ZLAQR2) ==== 393* 394 LWKOPT = MAX( 3*NSR / 2, INT( WORK( 1 ) ) ) 395* 396* ==== Quick return in case of workspace query. ==== 397* 398 IF( LWORK.EQ.-1 ) THEN 399 WORK( 1 ) = DCMPLX( LWKOPT, 0 ) 400 RETURN 401 END IF 402* 403* ==== ZLAHQR/ZLAQR0 crossover point ==== 404* 405 NMIN = ILAENV( 12, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 406 NMIN = MAX( NTINY, NMIN ) 407* 408* ==== Nibble crossover point ==== 409* 410 NIBBLE = ILAENV( 14, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 411 NIBBLE = MAX( 0, NIBBLE ) 412* 413* ==== Accumulate reflections during ttswp? Use block 414* . 2-by-2 structure during matrix-matrix multiply? ==== 415* 416 KACC22 = ILAENV( 16, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 417 KACC22 = MAX( 0, KACC22 ) 418 KACC22 = MIN( 2, KACC22 ) 419* 420* ==== NWMAX = the largest possible deflation window for 421* . which there is sufficient workspace. ==== 422* 423 NWMAX = MIN( ( N-1 ) / 3, LWORK / 2 ) 424 NW = NWMAX 425* 426* ==== NSMAX = the Largest number of simultaneous shifts 427* . for which there is sufficient workspace. ==== 428* 429 NSMAX = MIN( ( N+6 ) / 9, 2*LWORK / 3 ) 430 NSMAX = NSMAX - MOD( NSMAX, 2 ) 431* 432* ==== NDFL: an iteration count restarted at deflation. ==== 433* 434 NDFL = 1 435* 436* ==== ITMAX = iteration limit ==== 437* 438 ITMAX = MAX( 30, 2*KEXSH )*MAX( 10, ( IHI-ILO+1 ) ) 439* 440* ==== Last row and column in the active block ==== 441* 442 KBOT = IHI 443* 444* ==== Main Loop ==== 445* 446 DO 70 IT = 1, ITMAX 447* 448* ==== Done when KBOT falls below ILO ==== 449* 450 IF( KBOT.LT.ILO ) 451 $ GO TO 80 452* 453* ==== Locate active block ==== 454* 455 DO 10 K = KBOT, ILO + 1, -1 456 IF( H( K, K-1 ).EQ.ZERO ) 457 $ GO TO 20 458 10 CONTINUE 459 K = ILO 460 20 CONTINUE 461 KTOP = K 462* 463* ==== Select deflation window size: 464* . Typical Case: 465* . If possible and advisable, nibble the entire 466* . active block. If not, use size MIN(NWR,NWMAX) 467* . or MIN(NWR+1,NWMAX) depending upon which has 468* . the smaller corresponding subdiagonal entry 469* . (a heuristic). 470* . 471* . Exceptional Case: 472* . If there have been no deflations in KEXNW or 473* . more iterations, then vary the deflation window 474* . size. At first, because, larger windows are, 475* . in general, more powerful than smaller ones, 476* . rapidly increase the window to the maximum possible. 477* . Then, gradually reduce the window size. ==== 478* 479 NH = KBOT - KTOP + 1 480 NWUPBD = MIN( NH, NWMAX ) 481 IF( NDFL.LT.KEXNW ) THEN 482 NW = MIN( NWUPBD, NWR ) 483 ELSE 484 NW = MIN( NWUPBD, 2*NW ) 485 END IF 486 IF( NW.LT.NWMAX ) THEN 487 IF( NW.GE.NH-1 ) THEN 488 NW = NH 489 ELSE 490 KWTOP = KBOT - NW + 1 491 IF( CABS1( H( KWTOP, KWTOP-1 ) ).GT. 492 $ CABS1( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1 493 END IF 494 END IF 495 IF( NDFL.LT.KEXNW ) THEN 496 NDEC = -1 497 ELSE IF( NDEC.GE.0 .OR. NW.GE.NWUPBD ) THEN 498 NDEC = NDEC + 1 499 IF( NW-NDEC.LT.2 ) 500 $ NDEC = 0 501 NW = NW - NDEC 502 END IF 503* 504* ==== Aggressive early deflation: 505* . split workspace under the subdiagonal into 506* . - an nw-by-nw work array V in the lower 507* . left-hand-corner, 508* . - an NW-by-at-least-NW-but-more-is-better 509* . (NW-by-NHO) horizontal work array along 510* . the bottom edge, 511* . - an at-least-NW-but-more-is-better (NHV-by-NW) 512* . vertical work array along the left-hand-edge. 513* . ==== 514* 515 KV = N - NW + 1 516 KT = NW + 1 517 NHO = ( N-NW-1 ) - KT + 1 518 KWV = NW + 2 519 NVE = ( N-NW ) - KWV + 1 520* 521* ==== Aggressive early deflation ==== 522* 523 CALL ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 524 $ IHIZ, Z, LDZ, LS, LD, W, H( KV, 1 ), LDH, NHO, 525 $ H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH, WORK, 526 $ LWORK ) 527* 528* ==== Adjust KBOT accounting for new deflations. ==== 529* 530 KBOT = KBOT - LD 531* 532* ==== KS points to the shifts. ==== 533* 534 KS = KBOT - LS + 1 535* 536* ==== Skip an expensive QR sweep if there is a (partly 537* . heuristic) reason to expect that many eigenvalues 538* . will deflate without it. Here, the QR sweep is 539* . skipped if many eigenvalues have just been deflated 540* . or if the remaining active block is small. 541* 542 IF( ( LD.EQ.0 ) .OR. ( ( 100*LD.LE.NW*NIBBLE ) .AND. ( KBOT- 543 $ KTOP+1.GT.MIN( NMIN, NWMAX ) ) ) ) THEN 544* 545* ==== NS = nominal number of simultaneous shifts. 546* . This may be lowered (slightly) if ZLAQR2 547* . did not provide that many shifts. ==== 548* 549 NS = MIN( NSMAX, NSR, MAX( 2, KBOT-KTOP ) ) 550 NS = NS - MOD( NS, 2 ) 551* 552* ==== If there have been no deflations 553* . in a multiple of KEXSH iterations, 554* . then try exceptional shifts. 555* . Otherwise use shifts provided by 556* . ZLAQR2 above or from the eigenvalues 557* . of a trailing principal submatrix. ==== 558* 559 IF( MOD( NDFL, KEXSH ).EQ.0 ) THEN 560 KS = KBOT - NS + 1 561 DO 30 I = KBOT, KS + 1, -2 562 W( I ) = H( I, I ) + WILK1*CABS1( H( I, I-1 ) ) 563 W( I-1 ) = W( I ) 564 30 CONTINUE 565 ELSE 566* 567* ==== Got NS/2 or fewer shifts? Use ZLAHQR 568* . on a trailing principal submatrix to 569* . get more. (Since NS.LE.NSMAX.LE.(N+6)/9, 570* . there is enough space below the subdiagonal 571* . to fit an NS-by-NS scratch array.) ==== 572* 573 IF( KBOT-KS+1.LE.NS / 2 ) THEN 574 KS = KBOT - NS + 1 575 KT = N - NS + 1 576 CALL ZLACPY( 'A', NS, NS, H( KS, KS ), LDH, 577 $ H( KT, 1 ), LDH ) 578 CALL ZLAHQR( .false., .false., NS, 1, NS, 579 $ H( KT, 1 ), LDH, W( KS ), 1, 1, ZDUM, 580 $ 1, INF ) 581 KS = KS + INF 582* 583* ==== In case of a rare QR failure use 584* . eigenvalues of the trailing 2-by-2 585* . principal submatrix. Scale to avoid 586* . overflows, underflows and subnormals. 587* . (The scale factor S can not be zero, 588* . because H(KBOT,KBOT-1) is nonzero.) ==== 589* 590 IF( KS.GE.KBOT ) THEN 591 S = CABS1( H( KBOT-1, KBOT-1 ) ) + 592 $ CABS1( H( KBOT, KBOT-1 ) ) + 593 $ CABS1( H( KBOT-1, KBOT ) ) + 594 $ CABS1( H( KBOT, KBOT ) ) 595 AA = H( KBOT-1, KBOT-1 ) / S 596 CC = H( KBOT, KBOT-1 ) / S 597 BB = H( KBOT-1, KBOT ) / S 598 DD = H( KBOT, KBOT ) / S 599 TR2 = ( AA+DD ) / TWO 600 DET = ( AA-TR2 )*( DD-TR2 ) - BB*CC 601 RTDISC = SQRT( -DET ) 602 W( KBOT-1 ) = ( TR2+RTDISC )*S 603 W( KBOT ) = ( TR2-RTDISC )*S 604* 605 KS = KBOT - 1 606 END IF 607 END IF 608* 609 IF( KBOT-KS+1.GT.NS ) THEN 610* 611* ==== Sort the shifts (Helps a little) ==== 612* 613 SORTED = .false. 614 DO 50 K = KBOT, KS + 1, -1 615 IF( SORTED ) 616 $ GO TO 60 617 SORTED = .true. 618 DO 40 I = KS, K - 1 619 IF( CABS1( W( I ) ).LT.CABS1( W( I+1 ) ) ) 620 $ THEN 621 SORTED = .false. 622 SWAP = W( I ) 623 W( I ) = W( I+1 ) 624 W( I+1 ) = SWAP 625 END IF 626 40 CONTINUE 627 50 CONTINUE 628 60 CONTINUE 629 END IF 630 END IF 631* 632* ==== If there are only two shifts, then use 633* . only one. ==== 634* 635 IF( KBOT-KS+1.EQ.2 ) THEN 636 IF( CABS1( W( KBOT )-H( KBOT, KBOT ) ).LT. 637 $ CABS1( W( KBOT-1 )-H( KBOT, KBOT ) ) ) THEN 638 W( KBOT-1 ) = W( KBOT ) 639 ELSE 640 W( KBOT ) = W( KBOT-1 ) 641 END IF 642 END IF 643* 644* ==== Use up to NS of the the smallest magnitude 645* . shifts. If there aren't NS shifts available, 646* . then use them all, possibly dropping one to 647* . make the number of shifts even. ==== 648* 649 NS = MIN( NS, KBOT-KS+1 ) 650 NS = NS - MOD( NS, 2 ) 651 KS = KBOT - NS + 1 652* 653* ==== Small-bulge multi-shift QR sweep: 654* . split workspace under the subdiagonal into 655* . - a KDU-by-KDU work array U in the lower 656* . left-hand-corner, 657* . - a KDU-by-at-least-KDU-but-more-is-better 658* . (KDU-by-NHo) horizontal work array WH along 659* . the bottom edge, 660* . - and an at-least-KDU-but-more-is-better-by-KDU 661* . (NVE-by-KDU) vertical work WV arrow along 662* . the left-hand-edge. ==== 663* 664 KDU = 3*NS - 3 665 KU = N - KDU + 1 666 KWH = KDU + 1 667 NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1 668 KWV = KDU + 4 669 NVE = N - KDU - KWV + 1 670* 671* ==== Small-bulge multi-shift QR sweep ==== 672* 673 CALL ZLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NS, 674 $ W( KS ), H, LDH, ILOZ, IHIZ, Z, LDZ, WORK, 675 $ 3, H( KU, 1 ), LDH, NVE, H( KWV, 1 ), LDH, 676 $ NHO, H( KU, KWH ), LDH ) 677 END IF 678* 679* ==== Note progress (or the lack of it). ==== 680* 681 IF( LD.GT.0 ) THEN 682 NDFL = 1 683 ELSE 684 NDFL = NDFL + 1 685 END IF 686* 687* ==== End of main loop ==== 688 70 CONTINUE 689* 690* ==== Iteration limit exceeded. Set INFO to show where 691* . the problem occurred and exit. ==== 692* 693 INFO = KBOT 694 80 CONTINUE 695 END IF 696* 697* ==== Return the optimal value of LWORK. ==== 698* 699 WORK( 1 ) = DCMPLX( LWKOPT, 0 ) 700* 701* ==== End of ZLAQR4 ==== 702* 703 END 704