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