SRC/dlahqr.f

   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 *>          .GT. 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 .GT. 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 thorugh INFO of the final, output
 164 *>                  value of H.
 165 *>
 166 *>                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
 167 *>          (*)       (initial value of H)*U  = U*(final value of H)
 168 *>                  where U is an orthognal 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 .GT. 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 November 2015
 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.6.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 *     November 2015
 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