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*> \date December 2016 187* 188*> \ingroup doubleOTHERauxiliary 189* 190*> \par Further Details: 191* ===================== 192*> 193*> \verbatim 194*> 195*> 02-96 Based on modifications by 196*> David Day, Sandia National Laboratory, USA 197*> 198*> 12-04 Further modifications by 199*> Ralph Byers, University of Kansas, USA 200*> This is a modified version of DLAHQR from LAPACK version 3.0. 201*> It is (1) more robust against overflow and underflow and 202*> (2) adopts the more conservative Ahues & Tisseur stopping 203*> criterion (LAWN 122, 1997). 204*> \endverbatim 205*> 206* ===================================================================== 207 SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, 208 $ ILOZ, IHIZ, Z, LDZ, INFO ) 209* 210* -- LAPACK auxiliary routine (version 3.7.0) -- 211* -- LAPACK is a software package provided by Univ. of Tennessee, -- 212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 213* December 2016 214* 215* .. Scalar Arguments .. 216 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N 217 LOGICAL WANTT, WANTZ 218* .. 219* .. Array Arguments .. 220 DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) 221* .. 222* 223* ========================================================= 224* 225* .. Parameters .. 226 DOUBLE PRECISION ZERO, ONE, TWO 227 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0 ) 228 DOUBLE PRECISION DAT1, DAT2 229 PARAMETER ( DAT1 = 3.0d0 / 4.0d0, DAT2 = -0.4375d0 ) 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* .. 238* .. Local Arrays .. 239 DOUBLE PRECISION V( 3 ) 240* .. 241* .. External Functions .. 242 DOUBLE PRECISION DLAMCH 243 EXTERNAL DLAMCH 244* .. 245* .. External Subroutines .. 246 EXTERNAL DCOPY, DLABAD, DLANV2, DLARFG, DROT 247* .. 248* .. Intrinsic Functions .. 249 INTRINSIC ABS, DBLE, MAX, MIN, SQRT 250* .. 251* .. Executable Statements .. 252* 253 INFO = 0 254* 255* Quick return if possible 256* 257 IF( N.EQ.0 ) 258 $ RETURN 259 IF( ILO.EQ.IHI ) THEN 260 WR( ILO ) = H( ILO, ILO ) 261 WI( ILO ) = ZERO 262 RETURN 263 END IF 264* 265* ==== clear out the trash ==== 266 DO 10 J = ILO, IHI - 3 267 H( J+2, J ) = ZERO 268 H( J+3, J ) = ZERO 269 10 CONTINUE 270 IF( ILO.LE.IHI-2 ) 271 $ H( IHI, IHI-2 ) = ZERO 272* 273 NH = IHI - ILO + 1 274 NZ = IHIZ - ILOZ + 1 275* 276* Set machine-dependent constants for the stopping criterion. 277* 278 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) 279 SAFMAX = ONE / SAFMIN 280 CALL DLABAD( SAFMIN, SAFMAX ) 281 ULP = DLAMCH( 'PRECISION' ) 282 SMLNUM = SAFMIN*( DBLE( NH ) / ULP ) 283* 284* I1 and I2 are the indices of the first row and last column of H 285* to which transformations must be applied. If eigenvalues only are 286* being computed, I1 and I2 are set inside the main loop. 287* 288 IF( WANTT ) THEN 289 I1 = 1 290 I2 = N 291 END IF 292* 293* ITMAX is the total number of QR iterations allowed. 294* 295 ITMAX = 30 * MAX( 10, NH ) 296* 297* The main loop begins here. I is the loop index and decreases from 298* IHI to ILO in steps of 1 or 2. Each iteration of the loop works 299* with the active submatrix in rows and columns L to I. 300* Eigenvalues I+1 to IHI have already converged. Either L = ILO or 301* H(L,L-1) is negligible so that the matrix splits. 302* 303 I = IHI 304 20 CONTINUE 305 L = ILO 306 IF( I.LT.ILO ) 307 $ GO TO 160 308* 309* Perform QR iterations on rows and columns ILO to I until a 310* submatrix of order 1 or 2 splits off at the bottom because a 311* subdiagonal element has become negligible. 312* 313 DO 140 ITS = 0, ITMAX 314* 315* Look for a single small subdiagonal element. 316* 317 DO 30 K = I, L + 1, -1 318 IF( ABS( H( K, K-1 ) ).LE.SMLNUM ) 319 $ GO TO 40 320 TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) ) 321 IF( TST.EQ.ZERO ) THEN 322 IF( K-2.GE.ILO ) 323 $ TST = TST + ABS( H( K-1, K-2 ) ) 324 IF( K+1.LE.IHI ) 325 $ TST = TST + ABS( H( K+1, K ) ) 326 END IF 327* ==== The following is a conservative small subdiagonal 328* . deflation criterion due to Ahues & Tisseur (LAWN 122, 329* . 1997). It has better mathematical foundation and 330* . improves accuracy in some cases. ==== 331 IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN 332 AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) 333 BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) 334 AA = MAX( ABS( H( K, K ) ), 335 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) 336 BB = MIN( ABS( H( K, K ) ), 337 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) 338 S = AA + AB 339 IF( BA*( AB / S ).LE.MAX( SMLNUM, 340 $ ULP*( BB*( AA / S ) ) ) )GO TO 40 341 END IF 342 30 CONTINUE 343 40 CONTINUE 344 L = K 345 IF( L.GT.ILO ) THEN 346* 347* H(L,L-1) is negligible 348* 349 H( L, L-1 ) = ZERO 350 END IF 351* 352* Exit from loop if a submatrix of order 1 or 2 has split off. 353* 354 IF( L.GE.I-1 ) 355 $ GO TO 150 356* 357* Now the active submatrix is in rows and columns L to I. If 358* eigenvalues only are being computed, only the active submatrix 359* need be transformed. 360* 361 IF( .NOT.WANTT ) THEN 362 I1 = L 363 I2 = I 364 END IF 365* 366 IF( ITS.EQ.10 ) THEN 367* 368* Exceptional shift. 369* 370 S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) ) 371 H11 = DAT1*S + H( L, L ) 372 H12 = DAT2*S 373 H21 = S 374 H22 = H11 375 ELSE IF( ITS.EQ.20 ) THEN 376* 377* Exceptional shift. 378* 379 S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) 380 H11 = DAT1*S + H( I, I ) 381 H12 = DAT2*S 382 H21 = S 383 H22 = H11 384 ELSE 385* 386* Prepare to use Francis' double shift 387* (i.e. 2nd degree generalized Rayleigh quotient) 388* 389 H11 = H( I-1, I-1 ) 390 H21 = H( I, I-1 ) 391 H12 = H( I-1, I ) 392 H22 = H( I, I ) 393 END IF 394 S = ABS( H11 ) + ABS( H12 ) + ABS( H21 ) + ABS( H22 ) 395 IF( S.EQ.ZERO ) THEN 396 RT1R = ZERO 397 RT1I = ZERO 398 RT2R = ZERO 399 RT2I = ZERO 400 ELSE 401 H11 = H11 / S 402 H21 = H21 / S 403 H12 = H12 / S 404 H22 = H22 / S 405 TR = ( H11+H22 ) / TWO 406 DET = ( H11-TR )*( H22-TR ) - H12*H21 407 RTDISC = SQRT( ABS( DET ) ) 408 IF( DET.GE.ZERO ) THEN 409* 410* ==== complex conjugate shifts ==== 411* 412 RT1R = TR*S 413 RT2R = RT1R 414 RT1I = RTDISC*S 415 RT2I = -RT1I 416 ELSE 417* 418* ==== real shifts (use only one of them) ==== 419* 420 RT1R = TR + RTDISC 421 RT2R = TR - RTDISC 422 IF( ABS( RT1R-H22 ).LE.ABS( RT2R-H22 ) ) THEN 423 RT1R = RT1R*S 424 RT2R = RT1R 425 ELSE 426 RT2R = RT2R*S 427 RT1R = RT2R 428 END IF 429 RT1I = ZERO 430 RT2I = ZERO 431 END IF 432 END IF 433* 434* Look for two consecutive small subdiagonal elements. 435* 436 DO 50 M = I - 2, L, -1 437* Determine the effect of starting the double-shift QR 438* iteration at row M, and see if this would make H(M,M-1) 439* negligible. (The following uses scaling to avoid 440* overflows and most underflows.) 441* 442 H21S = H( M+1, M ) 443 S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S ) 444 H21S = H( M+1, M ) / S 445 V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )* 446 $ ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S ) 447 V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R ) 448 V( 3 ) = H21S*H( M+2, M+1 ) 449 S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) ) 450 V( 1 ) = V( 1 ) / S 451 V( 2 ) = V( 2 ) / S 452 V( 3 ) = V( 3 ) / S 453 IF( M.EQ.L ) 454 $ GO TO 60 455 IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE. 456 $ ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M, 457 $ M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60 458 50 CONTINUE 459 60 CONTINUE 460* 461* Double-shift QR step 462* 463 DO 130 K = M, I - 1 464* 465* The first iteration of this loop determines a reflection G 466* from the vector V and applies it from left and right to H, 467* thus creating a nonzero bulge below the subdiagonal. 468* 469* Each subsequent iteration determines a reflection G to 470* restore the Hessenberg form in the (K-1)th column, and thus 471* chases the bulge one step toward the bottom of the active 472* submatrix. NR is the order of G. 473* 474 NR = MIN( 3, I-K+1 ) 475 IF( K.GT.M ) 476 $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 ) 477 CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) 478 IF( K.GT.M ) THEN 479 H( K, K-1 ) = V( 1 ) 480 H( K+1, K-1 ) = ZERO 481 IF( K.LT.I-1 ) 482 $ H( K+2, K-1 ) = ZERO 483 ELSE IF( M.GT.L ) THEN 484* ==== Use the following instead of 485* . H( K, K-1 ) = -H( K, K-1 ) to 486* . avoid a bug when v(2) and v(3) 487* . underflow. ==== 488 H( K, K-1 ) = H( K, K-1 )*( ONE-T1 ) 489 END IF 490 V2 = V( 2 ) 491 T2 = T1*V2 492 IF( NR.EQ.3 ) THEN 493 V3 = V( 3 ) 494 T3 = T1*V3 495* 496* Apply G from the left to transform the rows of the matrix 497* in columns K to I2. 498* 499 DO 70 J = K, I2 500 SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J ) 501 H( K, J ) = H( K, J ) - SUM*T1 502 H( K+1, J ) = H( K+1, J ) - SUM*T2 503 H( K+2, J ) = H( K+2, J ) - SUM*T3 504 70 CONTINUE 505* 506* Apply G from the right to transform the columns of the 507* matrix in rows I1 to min(K+3,I). 508* 509 DO 80 J = I1, MIN( K+3, I ) 510 SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 ) 511 H( J, K ) = H( J, K ) - SUM*T1 512 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 513 H( J, K+2 ) = H( J, K+2 ) - SUM*T3 514 80 CONTINUE 515* 516 IF( WANTZ ) THEN 517* 518* Accumulate transformations in the matrix Z 519* 520 DO 90 J = ILOZ, IHIZ 521 SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 ) 522 Z( J, K ) = Z( J, K ) - SUM*T1 523 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 524 Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3 525 90 CONTINUE 526 END IF 527 ELSE IF( NR.EQ.2 ) THEN 528* 529* Apply G from the left to transform the rows of the matrix 530* in columns K to I2. 531* 532 DO 100 J = K, I2 533 SUM = H( K, J ) + V2*H( K+1, J ) 534 H( K, J ) = H( K, J ) - SUM*T1 535 H( K+1, J ) = H( K+1, J ) - SUM*T2 536 100 CONTINUE 537* 538* Apply G from the right to transform the columns of the 539* matrix in rows I1 to min(K+3,I). 540* 541 DO 110 J = I1, I 542 SUM = H( J, K ) + V2*H( J, K+1 ) 543 H( J, K ) = H( J, K ) - SUM*T1 544 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 545 110 CONTINUE 546* 547 IF( WANTZ ) THEN 548* 549* Accumulate transformations in the matrix Z 550* 551 DO 120 J = ILOZ, IHIZ 552 SUM = Z( J, K ) + V2*Z( J, K+1 ) 553 Z( J, K ) = Z( J, K ) - SUM*T1 554 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 555 120 CONTINUE 556 END IF 557 END IF 558 130 CONTINUE 559* 560 140 CONTINUE 561* 562* Failure to converge in remaining number of iterations 563* 564 INFO = I 565 RETURN 566* 567 150 CONTINUE 568* 569 IF( L.EQ.I ) THEN 570* 571* H(I,I-1) is negligible: one eigenvalue has converged. 572* 573 WR( I ) = H( I, I ) 574 WI( I ) = ZERO 575 ELSE IF( L.EQ.I-1 ) THEN 576* 577* H(I-1,I-2) is negligible: a pair of eigenvalues have converged. 578* 579* Transform the 2-by-2 submatrix to standard Schur form, 580* and compute and store the eigenvalues. 581* 582 CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), 583 $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ), 584 $ CS, SN ) 585* 586 IF( WANTT ) THEN 587* 588* Apply the transformation to the rest of H. 589* 590 IF( I2.GT.I ) 591 $ CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, 592 $ CS, SN ) 593 CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN ) 594 END IF 595 IF( WANTZ ) THEN 596* 597* Apply the transformation to Z. 598* 599 CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN ) 600 END IF 601 END IF 602* 603* return to start of the main loop with new value of I. 604* 605 I = L - 1 606 GO TO 20 607* 608 160 CONTINUE 609 RETURN 610* 611* End of DLAHQR 612* 613 END 614