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