1*> \brief \b DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLAHQR + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahqr.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahqr.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahqr.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, 22* ILOZ, IHIZ, Z, LDZ, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N 26* LOGICAL WANTT, WANTZ 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DLAHQR is an auxiliary routine called by DHSEQR to update the 39*> eigenvalues and Schur decomposition already computed by DHSEQR, by 40*> dealing with the Hessenberg submatrix in rows and columns ILO to 41*> IHI. 42*> \endverbatim 43* 44* Arguments: 45* ========== 46* 47*> \param[in] WANTT 48*> \verbatim 49*> WANTT is LOGICAL 50*> = .TRUE. : the full Schur form T is required; 51*> = .FALSE.: only eigenvalues are required. 52*> \endverbatim 53*> 54*> \param[in] WANTZ 55*> \verbatim 56*> WANTZ is LOGICAL 57*> = .TRUE. : the matrix of Schur vectors Z is required; 58*> = .FALSE.: Schur vectors are not required. 59*> \endverbatim 60*> 61*> \param[in] N 62*> \verbatim 63*> N is INTEGER 64*> The order of the matrix H. N >= 0. 65*> \endverbatim 66*> 67*> \param[in] ILO 68*> \verbatim 69*> ILO is INTEGER 70*> \endverbatim 71*> 72*> \param[in] IHI 73*> \verbatim 74*> IHI is INTEGER 75*> It is assumed that H is already upper quasi-triangular in 76*> rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless 77*> ILO = 1). DLAHQR works primarily with the Hessenberg 78*> submatrix in rows and columns ILO to IHI, but applies 79*> transformations to all of H if WANTT is .TRUE.. 80*> 1 <= ILO <= max(1,IHI); IHI <= N. 81*> \endverbatim 82*> 83*> \param[in,out] H 84*> \verbatim 85*> H is DOUBLE PRECISION array, dimension (LDH,N) 86*> On entry, the upper Hessenberg matrix H. 87*> On exit, if INFO is zero and if WANTT is .TRUE., H is upper 88*> quasi-triangular in rows and columns ILO:IHI, with any 89*> 2-by-2 diagonal blocks in standard form. If INFO is zero 90*> and WANTT is .FALSE., the contents of H are unspecified on 91*> exit. The output state of H if INFO is nonzero is given 92*> below under the description of INFO. 93*> \endverbatim 94*> 95*> \param[in] LDH 96*> \verbatim 97*> LDH is INTEGER 98*> The leading dimension of the array H. LDH >= max(1,N). 99*> \endverbatim 100*> 101*> \param[out] WR 102*> \verbatim 103*> WR is DOUBLE PRECISION array, dimension (N) 104*> \endverbatim 105*> 106*> \param[out] WI 107*> \verbatim 108*> WI is DOUBLE PRECISION array, dimension (N) 109*> The real and imaginary parts, respectively, of the computed 110*> eigenvalues ILO to IHI are stored in the corresponding 111*> elements of WR and WI. If two eigenvalues are computed as a 112*> complex conjugate pair, they are stored in consecutive 113*> elements of WR and WI, say the i-th and (i+1)th, with 114*> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the 115*> eigenvalues are stored in the same order as on the diagonal 116*> of the Schur form returned in H, with WR(i) = H(i,i), and, if 117*> H(i:i+1,i:i+1) is a 2-by-2 diagonal block, 118*> WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(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 DOUBLE PRECISION array, dimension (LDZ,N) 137*> If WANTZ is .TRUE., on entry Z must contain the current 138*> matrix Z of transformations accumulated by DHSEQR, and on 139*> exit Z has been updated; transformations are applied only to 140*> the submatrix Z(ILOZ:IHIZ,ILO:IHI). 141*> If WANTZ is .FALSE., Z is not referenced. 142*> \endverbatim 143*> 144*> \param[in] LDZ 145*> \verbatim 146*> LDZ is INTEGER 147*> The leading dimension of the array Z. LDZ >= max(1,N). 148*> \endverbatim 149*> 150*> \param[out] INFO 151*> \verbatim 152*> INFO is INTEGER 153*> = 0: successful exit 154*> > 0: If INFO = i, DLAHQR failed to compute all the 155*> eigenvalues ILO to IHI in a total of 30 iterations 156*> per eigenvalue; elements i+1:ihi of WR and WI 157*> contain those eigenvalues which have been 158*> successfully computed. 159*> 160*> If INFO > 0 and WANTT is .FALSE., then on exit, 161*> the remaining unconverged eigenvalues are the 162*> eigenvalues of the upper Hessenberg matrix rows 163*> and columns ILO through INFO of the final, output 164*> value of H. 165*> 166*> If INFO > 0 and WANTT is .TRUE., then on exit 167*> (*) (initial value of H)*U = U*(final value of H) 168*> where U is an orthogonal matrix. The final 169*> value of H is upper Hessenberg and triangular in 170*> rows and columns INFO+1 through IHI. 171*> 172*> If INFO > 0 and WANTZ is .TRUE., then on exit 173*> (final value of Z) = (initial value of Z)*U 174*> where U is the orthogonal matrix in (*) 175*> (regardless of the value of WANTT.) 176*> \endverbatim 177* 178* Authors: 179* ======== 180* 181*> \author Univ. of Tennessee 182*> \author Univ. of California Berkeley 183*> \author Univ. of Colorado Denver 184*> \author NAG Ltd. 185* 186*> \ingroup doubleOTHERauxiliary 187* 188*> \par Further Details: 189* ===================== 190*> 191*> \verbatim 192*> 193*> 02-96 Based on modifications by 194*> David Day, Sandia National Laboratory, USA 195*> 196*> 12-04 Further modifications by 197*> Ralph Byers, University of Kansas, USA 198*> This is a modified version of DLAHQR from LAPACK version 3.0. 199*> It is (1) more robust against overflow and underflow and 200*> (2) adopts the more conservative Ahues & Tisseur stopping 201*> criterion (LAWN 122, 1997). 202*> \endverbatim 203*> 204* ===================================================================== 205 SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, 206 $ ILOZ, IHIZ, Z, LDZ, INFO ) 207 IMPLICIT NONE 208* 209* -- LAPACK auxiliary routine -- 210* -- LAPACK is a software package provided by Univ. of Tennessee, -- 211* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 212* 213* .. Scalar Arguments .. 214 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N 215 LOGICAL WANTT, WANTZ 216* .. 217* .. Array Arguments .. 218 DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) 219* .. 220* 221* ========================================================= 222* 223* .. Parameters .. 224 DOUBLE PRECISION ZERO, ONE, TWO 225 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0 ) 226 DOUBLE PRECISION DAT1, DAT2 227 PARAMETER ( DAT1 = 3.0d0 / 4.0d0, DAT2 = -0.4375d0 ) 228 INTEGER KEXSH 229 PARAMETER ( KEXSH = 10 ) 230* .. 231* .. Local Scalars .. 232 DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S, 233 $ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX, 234 $ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST, 235 $ ULP, V2, V3 236 INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ, 237 $ KDEFL 238* .. 239* .. Local Arrays .. 240 DOUBLE PRECISION V( 3 ) 241* .. 242* .. External Functions .. 243 DOUBLE PRECISION DLAMCH 244 EXTERNAL DLAMCH 245* .. 246* .. External Subroutines .. 247 EXTERNAL DCOPY, DLABAD, DLANV2, DLARFG, DROT 248* .. 249* .. Intrinsic Functions .. 250 INTRINSIC ABS, DBLE, MAX, MIN, SQRT 251* .. 252* .. Executable Statements .. 253* 254 INFO = 0 255* 256* Quick return if possible 257* 258 IF( N.EQ.0 ) 259 $ RETURN 260 IF( ILO.EQ.IHI ) THEN 261 WR( ILO ) = H( ILO, ILO ) 262 WI( ILO ) = ZERO 263 RETURN 264 END IF 265* 266* ==== clear out the trash ==== 267 DO 10 J = ILO, IHI - 3 268 H( J+2, J ) = ZERO 269 H( J+3, J ) = ZERO 270 10 CONTINUE 271 IF( ILO.LE.IHI-2 ) 272 $ H( IHI, IHI-2 ) = ZERO 273* 274 NH = IHI - ILO + 1 275 NZ = IHIZ - ILOZ + 1 276* 277* Set machine-dependent constants for the stopping criterion. 278* 279 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) 280 SAFMAX = ONE / SAFMIN 281 CALL DLABAD( SAFMIN, SAFMAX ) 282 ULP = DLAMCH( 'PRECISION' ) 283 SMLNUM = SAFMIN*( DBLE( NH ) / ULP ) 284* 285* I1 and I2 are the indices of the first row and last column of H 286* to which transformations must be applied. If eigenvalues only are 287* being computed, I1 and I2 are set inside the main loop. 288* 289 IF( WANTT ) THEN 290 I1 = 1 291 I2 = N 292 END IF 293* 294* ITMAX is the total number of QR iterations allowed. 295* 296 ITMAX = 30 * MAX( 10, NH ) 297* 298* KDEFL counts the number of iterations since a deflation 299* 300 KDEFL = 0 301* 302* The main loop begins here. I is the loop index and decreases from 303* IHI to ILO in steps of 1 or 2. Each iteration of the loop works 304* with the active submatrix in rows and columns L to I. 305* Eigenvalues I+1 to IHI have already converged. Either L = ILO or 306* H(L,L-1) is negligible so that the matrix splits. 307* 308 I = IHI 309 20 CONTINUE 310 L = ILO 311 IF( I.LT.ILO ) 312 $ GO TO 160 313* 314* Perform QR iterations on rows and columns ILO to I until a 315* submatrix of order 1 or 2 splits off at the bottom because a 316* subdiagonal element has become negligible. 317* 318 DO 140 ITS = 0, ITMAX 319* 320* Look for a single small subdiagonal element. 321* 322 DO 30 K = I, L + 1, -1 323 IF( ABS( H( K, K-1 ) ).LE.SMLNUM ) 324 $ GO TO 40 325 TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) ) 326 IF( TST.EQ.ZERO ) THEN 327 IF( K-2.GE.ILO ) 328 $ TST = TST + ABS( H( K-1, K-2 ) ) 329 IF( K+1.LE.IHI ) 330 $ TST = TST + ABS( H( K+1, K ) ) 331 END IF 332* ==== The following is a conservative small subdiagonal 333* . deflation criterion due to Ahues & Tisseur (LAWN 122, 334* . 1997). It has better mathematical foundation and 335* . improves accuracy in some cases. ==== 336 IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN 337 AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) 338 BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) 339 AA = MAX( ABS( H( K, K ) ), 340 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) 341 BB = MIN( ABS( H( K, K ) ), 342 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) 343 S = AA + AB 344 IF( BA*( AB / S ).LE.MAX( SMLNUM, 345 $ ULP*( BB*( AA / S ) ) ) )GO TO 40 346 END IF 347 30 CONTINUE 348 40 CONTINUE 349 L = K 350 IF( L.GT.ILO ) THEN 351* 352* H(L,L-1) is negligible 353* 354 H( L, L-1 ) = ZERO 355 END IF 356* 357* Exit from loop if a submatrix of order 1 or 2 has split off. 358* 359 IF( L.GE.I-1 ) 360 $ GO TO 150 361 KDEFL = KDEFL + 1 362* 363* Now the active submatrix is in rows and columns L to I. If 364* eigenvalues only are being computed, only the active submatrix 365* need be transformed. 366* 367 IF( .NOT.WANTT ) THEN 368 I1 = L 369 I2 = I 370 END IF 371* 372 IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN 373* 374* Exceptional shift. 375* 376 S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) 377 H11 = DAT1*S + H( I, I ) 378 H12 = DAT2*S 379 H21 = S 380 H22 = H11 381 ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN 382* 383* Exceptional shift. 384* 385 S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) ) 386 H11 = DAT1*S + H( L, L ) 387 H12 = DAT2*S 388 H21 = S 389 H22 = H11 390 ELSE 391* 392* Prepare to use Francis' double shift 393* (i.e. 2nd degree generalized Rayleigh quotient) 394* 395 H11 = H( I-1, I-1 ) 396 H21 = H( I, I-1 ) 397 H12 = H( I-1, I ) 398 H22 = H( I, I ) 399 END IF 400 S = ABS( H11 ) + ABS( H12 ) + ABS( H21 ) + ABS( H22 ) 401 IF( S.EQ.ZERO ) THEN 402 RT1R = ZERO 403 RT1I = ZERO 404 RT2R = ZERO 405 RT2I = ZERO 406 ELSE 407 H11 = H11 / S 408 H21 = H21 / S 409 H12 = H12 / S 410 H22 = H22 / S 411 TR = ( H11+H22 ) / TWO 412 DET = ( H11-TR )*( H22-TR ) - H12*H21 413 RTDISC = SQRT( ABS( DET ) ) 414 IF( DET.GE.ZERO ) THEN 415* 416* ==== complex conjugate shifts ==== 417* 418 RT1R = TR*S 419 RT2R = RT1R 420 RT1I = RTDISC*S 421 RT2I = -RT1I 422 ELSE 423* 424* ==== real shifts (use only one of them) ==== 425* 426 RT1R = TR + RTDISC 427 RT2R = TR - RTDISC 428 IF( ABS( RT1R-H22 ).LE.ABS( RT2R-H22 ) ) THEN 429 RT1R = RT1R*S 430 RT2R = RT1R 431 ELSE 432 RT2R = RT2R*S 433 RT1R = RT2R 434 END IF 435 RT1I = ZERO 436 RT2I = ZERO 437 END IF 438 END IF 439* 440* Look for two consecutive small subdiagonal elements. 441* 442 DO 50 M = I - 2, L, -1 443* Determine the effect of starting the double-shift QR 444* iteration at row M, and see if this would make H(M,M-1) 445* negligible. (The following uses scaling to avoid 446* overflows and most underflows.) 447* 448 H21S = H( M+1, M ) 449 S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S ) 450 H21S = H( M+1, M ) / S 451 V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )* 452 $ ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S ) 453 V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R ) 454 V( 3 ) = H21S*H( M+2, M+1 ) 455 S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) ) 456 V( 1 ) = V( 1 ) / S 457 V( 2 ) = V( 2 ) / S 458 V( 3 ) = V( 3 ) / S 459 IF( M.EQ.L ) 460 $ GO TO 60 461 IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE. 462 $ ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M, 463 $ M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60 464 50 CONTINUE 465 60 CONTINUE 466* 467* Double-shift QR step 468* 469 DO 130 K = M, I - 1 470* 471* The first iteration of this loop determines a reflection G 472* from the vector V and applies it from left and right to H, 473* thus creating a nonzero bulge below the subdiagonal. 474* 475* Each subsequent iteration determines a reflection G to 476* restore the Hessenberg form in the (K-1)th column, and thus 477* chases the bulge one step toward the bottom of the active 478* submatrix. NR is the order of G. 479* 480 NR = MIN( 3, I-K+1 ) 481 IF( K.GT.M ) 482 $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 ) 483 CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) 484 IF( K.GT.M ) THEN 485 H( K, K-1 ) = V( 1 ) 486 H( K+1, K-1 ) = ZERO 487 IF( K.LT.I-1 ) 488 $ H( K+2, K-1 ) = ZERO 489 ELSE IF( M.GT.L ) THEN 490* ==== Use the following instead of 491* . H( K, K-1 ) = -H( K, K-1 ) to 492* . avoid a bug when v(2) and v(3) 493* . underflow. ==== 494 H( K, K-1 ) = H( K, K-1 )*( ONE-T1 ) 495 END IF 496 V2 = V( 2 ) 497 T2 = T1*V2 498 IF( NR.EQ.3 ) THEN 499 V3 = V( 3 ) 500 T3 = T1*V3 501* 502* Apply G from the left to transform the rows of the matrix 503* in columns K to I2. 504* 505 DO 70 J = K, I2 506 SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J ) 507 H( K, J ) = H( K, J ) - SUM*T1 508 H( K+1, J ) = H( K+1, J ) - SUM*T2 509 H( K+2, J ) = H( K+2, J ) - SUM*T3 510 70 CONTINUE 511* 512* Apply G from the right to transform the columns of the 513* matrix in rows I1 to min(K+3,I). 514* 515 DO 80 J = I1, MIN( K+3, I ) 516 SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 ) 517 H( J, K ) = H( J, K ) - SUM*T1 518 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 519 H( J, K+2 ) = H( J, K+2 ) - SUM*T3 520 80 CONTINUE 521* 522 IF( WANTZ ) THEN 523* 524* Accumulate transformations in the matrix Z 525* 526 DO 90 J = ILOZ, IHIZ 527 SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 ) 528 Z( J, K ) = Z( J, K ) - SUM*T1 529 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 530 Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3 531 90 CONTINUE 532 END IF 533 ELSE IF( NR.EQ.2 ) THEN 534* 535* Apply G from the left to transform the rows of the matrix 536* in columns K to I2. 537* 538 DO 100 J = K, I2 539 SUM = H( K, J ) + V2*H( K+1, J ) 540 H( K, J ) = H( K, J ) - SUM*T1 541 H( K+1, J ) = H( K+1, J ) - SUM*T2 542 100 CONTINUE 543* 544* Apply G from the right to transform the columns of the 545* matrix in rows I1 to min(K+3,I). 546* 547 DO 110 J = I1, I 548 SUM = H( J, K ) + V2*H( J, K+1 ) 549 H( J, K ) = H( J, K ) - SUM*T1 550 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 551 110 CONTINUE 552* 553 IF( WANTZ ) THEN 554* 555* Accumulate transformations in the matrix Z 556* 557 DO 120 J = ILOZ, IHIZ 558 SUM = Z( J, K ) + V2*Z( J, K+1 ) 559 Z( J, K ) = Z( J, K ) - SUM*T1 560 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 561 120 CONTINUE 562 END IF 563 END IF 564 130 CONTINUE 565* 566 140 CONTINUE 567* 568* Failure to converge in remaining number of iterations 569* 570 INFO = I 571 RETURN 572* 573 150 CONTINUE 574* 575 IF( L.EQ.I ) THEN 576* 577* H(I,I-1) is negligible: one eigenvalue has converged. 578* 579 WR( I ) = H( I, I ) 580 WI( I ) = ZERO 581 ELSE IF( L.EQ.I-1 ) THEN 582* 583* H(I-1,I-2) is negligible: a pair of eigenvalues have converged. 584* 585* Transform the 2-by-2 submatrix to standard Schur form, 586* and compute and store the eigenvalues. 587* 588 CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), 589 $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ), 590 $ CS, SN ) 591* 592 IF( WANTT ) THEN 593* 594* Apply the transformation to the rest of H. 595* 596 IF( I2.GT.I ) 597 $ CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, 598 $ CS, SN ) 599 CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN ) 600 END IF 601 IF( WANTZ ) THEN 602* 603* Apply the transformation to Z. 604* 605 CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN ) 606 END IF 607 END IF 608* reset deflation counter 609 KDEFL = 0 610* 611* return to start of the main loop with new value of I. 612* 613 I = L - 1 614 GO TO 20 615* 616 160 CONTINUE 617 RETURN 618* 619* End of DLAHQR 620* 621 END 622