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