*     =========================================================
*
*     Purpose
*     =======
*
*     SLAHQR is an auxiliary routine called by SHSEQR to update the
*     eigenvalues and Schur decomposition already computed by SHSEQR, by
*     dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
*
*     Arguments
*     =========
*
*     WANTT   (input) LOGICAL
*          = .TRUE. : the full Schur form T is required;
*          = .FALSE.: only eigenvalues are required.
*
*     WANTZ   (input) LOGICAL
*          = .TRUE. : the matrix of Schur vectors Z is required;
*          = .FALSE.: Schur vectors are not required.
*
*     N       (input) INTEGER
*          The order of the matrix H.  N >= 0.
*
*     ILO     (input) INTEGER
*     IHI     (input) INTEGER
*          It is assumed that H is already upper quasi-triangular in
*          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
*          ILO = 1). SLAHQR works primarily with the Hessenberg
*          submatrix in rows and columns ILO to IHI, but applies
*          transformations to all of H if WANTT is .TRUE..
*          1 <= ILO <= max(1,IHI); IHI <= N.
*
*     H       (input/output) REAL array, dimension (LDH,N)
*          On entry, the upper Hessenberg matrix H.
*          On exit, if INFO is zero and if WANTT is .TRUE., H is upper
*          quasi-triangular in rows and columns ILO:IHI, with any
*          2-by-2 diagonal blocks in standard form. If INFO is zero
*          and WANTT is .FALSE., the contents of H are unspecified on
*          exit.  The output state of H if INFO is nonzero is given
*          below under the description of INFO.
*
*     LDH     (input) INTEGER
*          The leading dimension of the array H. LDH >= max(1,N).
*
*     WR      (output) REAL array, dimension (N)
*     WI      (output) REAL array, dimension (N)
*          The real and imaginary parts, respectively, of the computed
*          eigenvalues ILO to IHI are stored in the corresponding
*          elements of WR and WI. If two eigenvalues are computed as a
*          complex conjugate pair, they are stored in consecutive
*          elements of WR and WI, say the i-th and (i+1)th, with
*          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
*          eigenvalues are stored in the same order as on the diagonal
*          of the Schur form returned in H, with WR(i) = H(i,i), and, if
*          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
*          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
*
*     ILOZ    (input) INTEGER
*     IHIZ    (input) INTEGER
*          Specify the rows of Z to which transformations must be
*          applied if WANTZ is .TRUE..
*          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*
*     Z       (input/output) REAL array, dimension (LDZ,N)
*          If WANTZ is .TRUE., on entry Z must contain the current
*          matrix Z of transformations accumulated by DHSEQR, and on
*          exit Z has been updated; transformations are applied only to
*          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
*          If WANTZ is .FALSE., Z is not referenced.
*
*     LDZ     (input) INTEGER
*          The leading dimension of the array Z. LDZ >= max(1,N).
*
*     INFO    (output) INTEGER
*           =   0: successful exit
*          .GT. 0: If INFO = i, SLAHQR failed to compute all the
*                  eigenvalues ILO to IHI in a total of 30 iterations
*                  per eigenvalue; elements i+1:ihi of WR and WI
*                  contain those eigenvalues which have been
*                  successfully computed.
*
*                  If INFO .GT. 0 and WANTT is .FALSE., then on exit,
*                  the remaining unconverged eigenvalues are the
*                  eigenvalues of the upper Hessenberg matrix rows
*                  and columns ILO thorugh INFO of the final, output
*                  value of H.
*
*                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
*          (*)       (initial value of H)*U  = U*(final value of H)
*                  where U is an orthognal matrix.    The final
*                  value of H is upper Hessenberg and triangular in
*                  rows and columns INFO+1 through IHI.
*
*     Further Details
*     ===============
*
*     02-96 Based on modifications by
*     David Day, Sandia National Laboratory, USA
*
*     12-04 Further modifications by
*     Ralph Byers, University of Kansas, USA
*     This is a modified version of SLAHQR from LAPACK version 3.0.
*     It is (1) more robust against overflow and underflow and
*     (2) adopts the more conservative Ahues & Tisseur stopping
*     criterion (LAWN 122, 1997).
*
*     =========================================================
*
      SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
     $                   ILOZ, IHIZ, Z, LDZ, INFO )
      IMPLICIT NONE
      INTEGER         IHI, ILO, IHIZ, ILOZ, INFO, LDH, LDZ, N
      INTEGER         I, ITS, J, K, L, M, NH
      LOGICAL         WANTT, WANTZ
      REAL            H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ), V( 3 )
      INTEGER         ITMAX
      PARAMETER       ( ITMAX = 30 )
      REAL            AA, BB, DET, H11, H12, H21, H21S, H22, H31S,
     $                AB, BA, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SMLNUM,
     $                CS, SN, CTEMP, STEMP, SUM, TR, TST,
     $                ULP, V1, V2, V3, BETA
      REAL            ZERO, ONE, TWO, DAT1, DAT2
      PARAMETER       ( ZERO = 0.0e0, ONE = 1.0e0, TWO = 2.0e0 )
      PARAMETER       ( DAT1 = 3.0e0 / 4.0e0, DAT2 = -0.4375e0 )
*
      INFO = 0
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
      IF( ILO.EQ.IHI ) THEN
         WR( ILO ) = H( ILO, ILO )
         WI( ILO ) = ZERO
         RETURN
      END IF
*
*     ==== clear out the trash ====
      DO 10 J = ILO, IHI - 3
         H( J+2, J ) = ZERO
         H( J+3, J ) = ZERO
   10 CONTINUE
      IF( ILO.LE.IHI-2 )
     $   H( IHI, IHI-2 ) = ZERO
*
      NH = IHI - ILO + 1
*
*     Set machine-dependent constants for the stopping criterion.
*
      ULP = TWO * EPSILON( ZERO )
      SMLNUM = TINY( ZERO )*( REAL( NH ) / ULP )
*
*     I1 and I2 are the indices of the first row and last column of H
*     to which transformations must be applied. If eigenvalues only are
*     being computed, I1 and I2 are set inside the main loop.
*
*     I1 = 1
*     I2 = N
*
*     The main loop begins here. I is the loop index and decreases from
*     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
*     with the active submatrix in rows and columns L to I.
*     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
*     H(L,L-1) is negligible so that the matrix splits.
*
      I = IHI
   20 CONTINUE
      L = ILO
      IF( I.LT.ILO )
     $   GO TO 160
*
*     Perform QR iterations on rows and columns ILO to I until a
*     submatrix of order 1 or 2 splits off at the bottom because a
*     subdiagonal element has become negligible.
*
      DO 140 ITS = 0, ITMAX
*
*        Look for a single small subdiagonal element.
*
         DO 30 K = I, L + 1, -1
            IF( ABS( H( K, K-1 ) ).LE.SMLNUM )
     $         GO TO 40
            TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
            IF( TST.EQ.ZERO ) THEN
               IF( K-2.GE.ILO )
     $            TST = TST + ABS( H( K-1, K-2 ) )
               IF( K.LT.IHI )
     $            TST = TST + ABS( H( K+1, K ) )
*               TST = SLANHS1( I-L+1, H( L, L ), LDH )
            END IF
*           ==== The following is a conservative small subdiagonal
*           .    deflation  criterion due to Ahues & Tisseur (LAWN 122,
*           .    1997). It has better mathematical foundation and
*           .    improves accuracy in some cases.  ====
            IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN
               AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
               BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
               AA = MAX( ABS( H( K, K ) ),
     $              ABS( H( K-1, K-1 )-H( K, K ) ) )
               BB = MIN( ABS( H( K, K ) ),
     $              ABS( H( K-1, K-1 )-H( K, K ) ) )
               S = AA + AB
               IF( BA*( AB / S ).LE.MAX( SMLNUM,
     $             ULP*( BB*( AA / S ) ) ) )GO TO 40
            END IF
   30    CONTINUE
   40    CONTINUE
         L = K
         IF( L.GT.ILO ) THEN
*           H(L,L-1) is negligible
            H( L, L-1 ) = ZERO
         END IF
*
*        Exit from loop if a submatrix of order 1 or 2 has split off.
*
         IF( L.GE.I-1 )
     $      GO TO 150
*
*        Now the active submatrix is in rows and columns L to I.
*
         IF( ITS.EQ.10 ) THEN
*           Exceptional shift.
            S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
            RT1R = DAT1*S + H( L, L )
            RT1I = S*SQRT( -DAT2 )
            RT2R =  RT1R
         ELSE IF( ITS.EQ.20 ) THEN
*           Exceptional shift.
            S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
            RT1R = DAT1*S + H( I, I )
            RT1I = S*SQRT( -DAT2 )
            RT2R =  RT1R
         ELSE
*           Prepare to use Francis' double shift
*           (i.e. 2nd degree generalized Rayleigh quotient)
            CALL SLANVE2( H( I-1, I-1 ), LDH,
     $                    RT1R, RT1I, RT2R, RT2I )
         END IF
*
*        Look for two consecutive small subdiagonal elements.
*
         DO 50 M = I - 2, L, -1
*           Determine the effect of starting the double-shift QR
*           iteration at row M, and see if this would make H(M,M-1)
*           negligible.  (The following uses scaling to avoid
*           overflows and most underflows.)
*
            H21S = H( M+1, M )
            H31S = H( M+2, M )
            S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + 
     $          ABS( H21S ) + ABS( H31S )
            H21S = H( M+1, M ) / S
            H31S = H( M+2, M ) / S
            V( 1 ) = H21S*H( M, M+1 ) + H31S*H( M, M+2 )
     $             + ( H( M, M )-RT1R )*( ( H( M, M )-RT2R ) / S )
     $             + ( ( RT1I*RT1I ) / S )
            V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R )
     $             + H31S*H( M+1, M+2 )
            V( 3 ) = H31S*( H( M, M )+H( M+2, M+2 )-RT1R-RT2R )
     $             + H21S*H( M+2, M+1 )
            S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) )
            V( 1 ) = V( 1 ) / S
            V( 2 ) = V( 2 ) / S
            V( 3 ) = V( 3 ) / S
            IF( M.EQ.L )
     $         GO TO 60
            IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE.
     $          ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M,
     $          M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60
   50    CONTINUE
   60    CONTINUE
*
*        Double-shift QR step
*
         DO 130 K = M, I - 2
*
*           The first iteration of this loop determines a reflection G
*           from the vector V and applies it from left and right to H,
*           thus creating a nonzero bulge below the subdiagonal.
*
*           Each subsequent iteration determines a reflection G to
*           restore the Hessenberg form in the (K-1)th column, and thus
*           chases the bulge one step toward the bottom of the active submatrix.
*
            IF( K.GT.M ) THEN
               V( 1 ) = H( K,   K-1 )
               V( 2 ) = H( K+1, K-1 )
               V( 3 ) = H( K+2, K-1 )
            END IF
            BETA = SIGN(SQRT(V(1)*V(1) + V(2)*V(2) + V(3)*V(3)), V(1))
            V1 = V( 1 ) + BETA
            V2 = V( 2 ) / V1
            V3 = V( 3 ) / V1
            V1 = V1 / BETA 
            IF( K.GT.M ) THEN
               H( K,   K-1 ) = -BETA
               H( K+1, K-1 ) =  ZERO
               H( K+2, K-1 ) =  ZERO
            ELSE IF( M.GT.L ) THEN
*               ==== Use the following instead of H( K, K-1 ) = -H( K, K-1 )
*               .    to avoid a bug when v(2) and v(3) underflow. ====
               H( K, K-1 ) = H( K, K-1 )*( ONE-V1 )
            END IF
*           Apply G from the left to transform the rows of the matrix
*           in columns K to N.
            DO 70 J = K, NOT
               SUM = V1*( H(K, J) + V2*H(K+1, J) + V3*H(K+2, J) )
               H( K,   J ) = H( K,   J ) - SUM
               H( K+1, J ) = H( K+1, J ) - SUM*V2
               H( K+2, J ) = H( K+2, J ) - SUM*V3
   70       CONTINUE
*           Apply G from the right to transform the columns of the
*           matrix in rows 1 to min(K+3,I).
            DO 80 J = 1, MIN( K+3, I )
               SUM = V1*(1(J, K) + V2*H(J, K+1) + V3*H(J, K+2)  )
               H( J, K   ) = H( J, K   ) - SUM
               H( J, K+1 ) = H( J, K+1 ) - SUM*V2
               H( J, K+2 ) = H( J, K+2 ) - SUM*V3
   80       CONTINUE
*
            IF( WANTZ ) THEN
*              Accumulate transformations in the matrix Z
               DO 90 J = ILOZ, IHIZ
                  SUM = V1*( Z(J, K) + V2*Z(J, K+1) + V3*Z(J, K+2) )
                  Z( J, K   ) = Z( J, K   ) - SUM
                  Z( J, K+1 ) = Z( J, K+1 ) - SUM*V2
                  Z( J, K+2 ) = Z( J, K+2 ) - SUM*V3
   90          CONTINUE
            END IF
  130    CONTINUE
         IF( I-1.GE.M ) THEN
            V1 = H( I-1, I-2 )
            V2 = H( I,   I-2 )
            BETA = SIGN( SQRT( V1*V1 + V2*V2 ), V1 )
            V1 = V1 + BETA
            V2 = V2 / V1
            V1 = V1 / BETA 
            H( I-1, I-2 ) = -BETA
            H( I,   I-2 ) =  ZERO
*           Apply G from the left to transform the rows of the matrix
*           in columns I-1 to I.
            DO 100 J = I-1, N
               SUM = V1*( H( I-1, J ) + V2*H( I, J ) )
               H( I-1, J ) = H( I-1, J ) - SUM
               H( I,   J ) = H( I,   J ) - SUM*V2
  100       CONTINUE
*           Apply G from the right to transform the columns of the
*           matrix in rows 1 to I.
            DO 110 J = 1, I
               SUM = V1*( H( J, I-1 ) + V2*H( J, I ) )
               H( J, I-1 ) = H( J, I-1 ) - SUM
               H( J, I   ) = H( J, I   ) - SUM*V2
  110       CONTINUE
            IF( WANTZ ) THEN
*              Accumulate transformations in the matrix Z
               DO 120 J = ILOZ, IHIZ
                  SUM = V1*( Z( J, I-1 ) + V2*Z( J, I ) )
                  Z( J, I-1 ) = Z( J, I-1 ) - SUM
                  Z( J, I   ) = Z( J, I   ) - SUM*V2
  120          CONTINUE
            END IF
         END IF
*
  140 CONTINUE
*
*     Failure to converge in remaining number of iterations
*
      INFO = I
      RETURN
*
  150 CONTINUE
*
      IF( L.EQ.I ) THEN
*        H(I,I-1) is negligible: one eigenvalue has converged.
         WR( I ) = H( I, I )
         WI( I ) = ZERO
         ND = 1
      ELSE IF( L.EQ.I-1 ) THEN
*        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
         CALL SLANVE2( H( I-1, I-1 ), LDH,
     $                 WR( I-1 ), WI( I-1 ), WR( I ), WI( I ) )
         ND = 2
      ELSE
         ND = 0
      END IF 
*
*     Transform the 2-by-2 submatrix to standard Schur form,
*
      IF( L.EQ.I-1 ) THEN
         CALL SLANVR2( H( I-1, I-1 )-H( I, I   ),
     $                 H( I-1, I   )+H( I, I-1 ), CS, SN )
         DO J = I-1, N
            CTEMP = CS*H( I-1, J ) + SN*H( I,   J )
            STEMP = CS*H( I,   J ) - SN*H( I-1, J )
            H( I-1, J ) = CTEMP
            H( I,   J ) = STEMP
         END DO
         DO J = 1, I
            CTEMP = CS*H( J, I-1 ) + SN*H( J, I   )
            STEMP = CS*H( J, I   ) - SN*H( J, I-1 )
            H( J, I-1 ) = CTEMP
            H( J, I   ) = STEMP
         END DO
*
         IF( WANTZ ) THEN
*           Apply the transformation to Z.
            DO J = ILOZ, IHIZ
               CTEMP = CS*Z( J, I-1 ) + SN*Z( J, I   )
               STEMP = CS*Z( J, I   ) - SN*Z( J, I-1 )
               Z( J, I-1 ) = CTEMP
               Z( J, I   ) = STEMP
            END DO
         END IF
      END IF
*
*     return to start of the main loop with new value of I.
*
      I = I - ND
      GO TO 20
*
  160 CONTINUE
      RETURN
*     End of SLAHQR
      END

      SUBROUTINE SLANVR2( AD, BC, CS, SN )
      IMPLICIT NONE
8     REAL          A, B, C, D, AD, BC, CS, SN, TAU
      REAL          SLAPY2, ZERO, HALF, ONE
      PARAMETER   ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
*
      CS = ONE
      SN = ZERO
*      IF( AD.NE.ZERO .AND. B.NE.ZERO .AND. C.NE.ZERO ) THEN
      IF( AD.NE.ZERO .AND. BC.NE.ZERO ) THEN
*        Make diagonal elements equal.
*        TAU = SLAPY2( SIGMA, TEMP )
         TAU = SQRT( AD*AD + BC*BC )
         CS = SQRT( ONE+ABS( BC ) / TAU )
         SN = -( AD / ( TAU*CS ) )*SIGN( ONE, BC )
         CS = CS*SQRT( HALF )
         SN = SN*SQRT( HALF )
      END IF
      RETURN
*     End of SLANV2
      END

*     =========================================================
*
*     Purpose
*     =======
*
*     SLAHQRE is an auxiliary routine called by SHSEQR to update the
*     eigenvalues already computed by SHSEQR, by dealing with
*     the Hessenberg submatrix in rows and columns ILO to IHI.
*
*     Arguments
*     =========
*
*     N       (input) INTEGER
*          The order of the matrix H.  N >= 0.
*
*     ILO     (input) INTEGER
*     IHI     (input) INTEGER
*          It is assumed that H is already upper quasi-triangular in
*          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
*          ILO = 1). SLAHQRE only works with the Hessenberg submatrix
*          in rows and columns ILO to IHI.
*          1 <= ILO <= max(1,IHI); IHI <= N.
*
*     H       (input/output) REAL array, dimension (LDH,N)
*          On entry, the upper Hessenberg matrix H.
*          On exit, if INFO is zero, the contents of H
*          are unspecified on exit.
*          The output state of H if INFO is nonzero is given below
*          under the description of INFO.
*
*     LDH     (input) INTEGER
*          The leading dimension of the array H. LDH >= max(1,N).
*
*     WR      (output) REAL array, dimension (N)
*     WI      (output) REAL array, dimension (N)
*          The real and imaginary parts, respectively, of the computed
*          eigenvalues ILO to IHI are stored in the corresponding
*          elements of WR and WI. If two eigenvalues are computed as a
*          complex conjugate pair, they are stored in consecutive
*          elements of WR and WI, say the i-th and (i+1)th, with
*          WI(i) > 0 and WI(i+1) < 0, i.e. WI(i+1) = -WI(i).
*
*     INFO    (output) INTEGER
*           =   0: successful exit
*          .GT. 0: If INFO = i, SLAHQR failed to compute all the
*                  eigenvalues ILO to IHI in a total of 30 iterations
*                  per eigenvalue; elements i+1:ihi of WR and WI
*                  contain those eigenvalues which have been
*                  successfully computed.
*
*                  The remaining unconverged eigenvalues are the
*                  eigenvalues of the upper Hessenberg matrix rows
*                  and columns ILO thorugh INFO of the final, output
*                  value of H.
*
*     Further Details
*     ===============
*
*     02-96 Based on modifications by
*     David Day, Sandia National Laboratory, USA
*
*     12-04 Further modifications by
*     Ralph Byers, University of Kansas, USA
*     This is a modified version of SLAHQR from LAPACK version 3.0.
*     It is (1) more robust against overflow and underflow and
*     (2) adopts the more conservative Ahues & Tisseur stopping
*     criterion (LAWN 122, 1997).
*
* https://www.youtube.com/watch?v=UUsYcKj-SS8
*     =========================================================
*
      SUBROUTINE SLAHQRE( N, ILO, IHI, H, LDH, WR, WI, INFO )
      IMPLICIT NONE
      INTEGER            IHI, ILO, INFO, LDH, N
      INTEGER            I, ITS, J, K, L, M, NH
      REAL               H( LDH, * ), WI( * ), WR( * ), V( 3 )
      INTEGER            ITMAX
      PARAMETER          ( ITMAX = 30 )
      REAL               AA, AB, BA, BB, DET, H11, H12, H21, H21S, H22,
     $                   H31S, RT1I, RT1R, RT2I, RT2R, RTDISC, S,
     $                   SMLNUM, CS, SN, SUM, TR, TST,
     $                   ULP, V1, V2, V3, BETA
      REAL               ZERO, ONE, TWO, DAT1, DAT2
      PARAMETER          ( ZERO = 0.0e0, ONE = 1.0e0, TWO = 2.0e0 )
      PARAMETER          ( DAT1 = 3.0e0 / 4.0e0, DAT2 = -0.4375e0 )
*
      INFO = 0
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
      IF( ILO.EQ.IHI ) THEN
         WR( ILO ) = H( ILO, ILO )
         WI( ILO ) = ZERO
         RETURN
      END IF
*
*     ==== clear out the trash ====
      DO 10 J = ILO, IHI - 3
         H( J+2, J ) = ZERO
         H( J+3, J ) = ZERO
   10 CONTINUE
      IF( ILO.LE.IHI-2 )
     $   H( IHI, IHI-2 ) = ZERO
*
*     Set machine-dependent constants for the stopping criterion.
*
      NH = IHI - ILO + 1
      ULP = TWO * EPSILON( ZERO )
      SMLNUM = TINY( ZERO )*( REAL( NH ) / ULP )
*
*     The main loop begins here. I is the loop index and decreases from
*     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
*     with the active submatrix in rows and columns L to I.
*     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
*     H(L,L-1) is negligible so that the matrix splits.
*
      I = IHI
   20 CONTINUE
      L = ILO
      IF( I.LT.ILO )
     $   GO TO 160
*
*     Perform QR iterations on rows and columns ILO to I until a
*     submatrix of order 1 or 2 splits off at the bottom because a
*     subdiagonal element has become negligible.
*
      DO 140 ITS = 0, ITMAX
*
*        Look for a single small subdiagonal element.
*
         DO 30 K = I, L + 1, -1
            IF( ABS( H( K, K-1 ) ).LE.SMLNUM )
     $         GO TO 40
            TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
            IF( TST.EQ.ZERO ) THEN
               IF( K-2.GE.ILO )
     $            TST = TST + ABS( H( K-1, K-2 ) )
               IF( K+1.LE.IHI )
     $            TST = TST + ABS( H( K+1, K ) )
*               TST = SLANHS1( I-L+1, H( L, L ), LDH )
            END IF
*           ==== The following is a conservative small subdiagonal
*           .    deflation  criterion due to Ahues & Tisseur (LAWN 122,
*           .    1997). It has better mathematical foundation and
*           .    improves accuracy in some cases.  ====
            IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN
               AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
               BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
               AA = MAX( ABS( H( K, K ) ),
     $              ABS( H( K-1, K-1 )-H( K, K ) ) )
               BB = MIN( ABS( H( K, K ) ),
     $              ABS( H( K-1, K-1 )-H( K, K ) ) )
               S = AA + AB
               IF( BA*( AB / S ).LE.MAX( SMLNUM,
     $             ULP*( BB*( AA / S ) ) ) )GO TO 40
            END IF
   30    CONTINUE
   40    CONTINUE
         L = K
         IF( L.GT.ILO ) THEN
*           H(L,L-1) is negligible
            H( L, L-1 ) = ZERO
         END IF
*
*        Exit from loop if a submatrix of order 1 or 2 has split off.
*
         IF( L.GE.I-1 )
     $      GO TO 150
*
*        Now the active submatrix is in rows and columns L to I.
*
         IF( ITS.EQ.10 ) THEN
*           Exceptional shift.
            S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
            RT1R = DAT1*S + H( L, L )
            RT1I = S*SQRT( -DAT2 )
            RT2R =  RT1R
         ELSE IF( ITS.EQ.20 ) THEN
*           Exceptional shift.
            S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
            RT1R = DAT1*S + H( I, I )
            RT1I = S*SQRT( -DAT2 )
            RT2R =  RT1R
         ELSE
*           Prepare to use Francis' double shift
*           (i.e. 2nd degree generalized Rayleigh quotient)
            CALL SLANVE2( H( I-1, I-1 ), LDH,
     $                    RT1R, RT1I, RT2R, RT2I )
         END IF
*
*        Look for two consecutive small subdiagonal elements.
*
         DO 50 M = I - 2, L, -1
*           Determine the effect of starting the double-shift QR
*           iteration at row M, and see if this would make H(M,M-1)
*           negligible.  (The following uses scaling to avoid
*           overflows and most underflows.)
*
            H21S = H( M+1, M )
            H31S = H( M+2, M )
            S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + 
     $          ABS( H21S ) + ABS( H31S )
            H21S = H( M+1, M ) / S
            H31S = H( M+2, M ) / S
            V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )*
     $               ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S )
            V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R )
     $             + H31S*H( M+1, M+2 )
            V( 3 ) = H31S*( H( M, M )+H( M+2, M+2 )-RT1R-RT2R )
     $             + H21S*H( M+2, M+1 )
*            V( 3 ) = H21S*H( M+2, M+1 )
            S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) )
            V( 1 ) = V( 1 ) / S
            V( 2 ) = V( 2 ) / S
            V( 3 ) = V( 3 ) / S
            IF( M.EQ.L )
     $         GO TO 60
            IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE.
     $          ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M,
     $          M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60
   50    CONTINUE
   60    CONTINUE
*
*        Double-shift QR step
*
         DO 130 K = M, I - 2
*
*           The first iteration of this loop determines a reflection G
*           from the vector V and applies it from left and right to H,
*           thus creating a nonzero bulge below the subdiagonal.
*
*           Each subsequent iteration determines a reflection G to
*           restore the Hessenberg form in the (K-1)th column, and thus
*           chases the bulge one step toward the bottom of the active submatrix.
*
            IF( K.GT.M ) THEN
               V( 1 ) = H( K,   K-1 )
               V( 2 ) = H( K+1, K-1 )
               V( 3 ) = H( K+2, K-1 )
            END IF
            BETA = SIGN(SQRT(V(1)*V(1) + V(2)*V(2) + V(3)*V(3)), V(1))
            V1 = V( 1 ) + BETA
            V2 = V( 2 ) / V1
            V3 = V( 3 ) / V1
            V1 = V1 / BETA 
            IF( K.GT.M ) THEN
               H( K,   K-1 ) = -BETA
               H( K+1, K-1 ) =  ZERO
               H( K+2, K-1 ) =  ZERO
            ELSE IF( M.GT.L ) THEN
*               ==== Use the following instead of H( K, K-1 ) = -H( K, K-1 )
*               .    to avoid a bug when v(2) and v(3) underflow. ====
               H( K, K-1 ) = H( K, K-1 )*( ONE-V1 )
            END IF
*           Apply G from the left to transform the rows of the matrix
*           in columns K to I.
            DO 70 J = K, I
               SUM = V1*( H(K, J) + V2*H(K+1, J) + V3*H(K+2, J) )
               H( K,   J ) = H( K,   J ) - SUM
               H( K+1, J ) = H( K+1, J ) - SUM*V2
               H( K+2, J ) = H( K+2, J ) - SUM*V3
   70       CONTINUE
*           Apply G from the right to transform the columns of the
*           matrix in rows L to min(K+3,I).
            DO 80 J = L, MIN( K+3, I )
               SUM = V1*( H(J, K) + V2*H(J, K+1) + V3*H(J, K+2)  )
               H( J, K   ) = H( J, K   ) - SUM
               H( J, K+1 ) = H( J, K+1 ) - SUM*V2
               H( J, K+2 ) = H( J, K+2 ) - SUM*V3
   80       CONTINUE
  130    CONTINUE
         IF( I-1.GE.M ) THEN
            V1 = H( I-1, I-2 )
            V2 = H( I,   I-2 )
            BETA = SIGN( SQRT( V1*V1 + V2*V2 ), V1 )
            V1 = V1 + BETA
            V2 = V2 / V1
            V1 = V1 / BETA 
            H( I-1, I-2 ) = -BETA
            H( I,   I-2 ) =  ZERO
*           Apply G from the left to transform the rows of the matrix
*           in columns I-1 to I.
            SUM = V1*( H( I-1, I-1 ) + V2*H( I, I-1 ) )
            H( I-1, I-1 ) = H( I-1, I-1 ) - SUM
            H( I,   I-1 ) = H( I,   I-1 ) - SUM*V2
            SUM = V1*( H( I-1, I ) + V2*H( I, I ) )
            H( I-1, I ) = H( I-1, I ) - SUM
            H( I,   I ) = H( I,   I ) - SUM*V2
*           Apply G from the right to transform the columns of the
*           matrix in rows L to I.
            DO 110 J = L, I
               SUM = V1*( H( J, I-1 ) + V2*H( J, I ) )
               H( J, I-1 ) = H( J, I-1 ) - SUM
               H( J, I   ) = H( J, I   ) - SUM*V2
  110       CONTINUE
         END IF
*
  140 CONTINUE
*
*     Failure to converge in remaining number of iterations
*
      INFO = I
      RETURN
*
  150 CONTINUE
      IF( H( I, I-1 ).LE.SMLNUM ) THEN
*        H(I,I-1) is negligible: one eigenvalue has converged.
         WR( I ) = H( I, I )
         WI( I ) = ZERO
         ND = 1
      ELSE IF( H( I-1, I-2 ).LE.SMLNUM ) THEN
*        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
         CALL SLANVE2( H( I-1, I-1 ), LDH,
     $                 WR( I-1 ), WI( I-1 ),
     $                 WR( I ), WI( I ) )
         ND = 2
      ELSE
         ND = 0 
      END IF
*     return to start of the main loop with new value of I.
      I = I - ND
      GO TO 20
*
  160 CONTINUE
      RETURN
*     End of SLAHQRE
      END

*  =====================================================================
*
*  Purpose
*  =======
*
*  SLANVE2 computes the eigenvalues of a real 2-by-2 nonsymmetric
*  matrix in standard form:
*
*       [ AA  BB ]
*       [ CC  DD ]
*
*  Arguments
*  =========
*
*  A       (input) REAL            
*  B       (input) REAL            
*  C       (input) REAL            
*  D       (input) REAL            
*          The elements of the input matrix.
*
*  RT1R    (output) REAL 
*  RT1I    (output) REAL            
*  RT2R    (output) REAL            
*  RT2I    (output) REAL            
*          The real and imaginary parts of the eigenvalues. If the
*          eigenvalues are a complex conjugate pair, RT1I > 0.
*
*  Further Details
*  ===============
*
*  Modified by V. Sima, Research Institute for Informatics, Bucharest,
*  Romania, to reduce the risk of cancellation errors,
*  when computing real eigenvalues
*
*  =====================================================================
*
      SUBROUTINE SLANVE2( AA, BB, CC, DD, RT1R, RT1I, RT2R, RT2I )
      IMPLICIT NONE
      REAL               AA, BB, CC, DD, RT1I, RT1R, RT2R, RT2I
      REAL               DELTA, P, Q, TEMP, Z
      REAL               ZERO, HALF, ONE
      PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
*     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
      IF( BB.EQ.ZERO .OR. CC.EQ.ZERO ) THEN
         RT1R = AA
         RT2R = DD
         RT1I = ZERO
      ELSE
         P    = HALF*(AA-DD)
         Q    = HALF*(AA+DD)
         TEMP = BB*CC
         DELTA = P*P + TEMP
         Z = SQRT( ABS( DELTA ) )
         IF( DELTA.GE.ZERO ) THEN
*           Real eigenvalues.
*            Q = P + SIGN( Z, P )
*            RT1R = DD + Q
*            RT2R = DD - ( TEMP / Q )
            RT1R = Q + Z
            RT2R = Q - Z
            RT1I = ZERO
         ELSE
*           Complex eigenvalues.
            RT1R = Q
            RT2R = Q
            RT1I = Z
         END IF
      END IF
      RT2I = -RT1I
      RETURN
*     End of SLANVE2
      END

*  =====================================================================
*
*  Purpose
*  =======
*
*  SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
*  matrix in standard form:
*
*       [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
*       [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]
*
*  where either
*  1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
*  2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
*  conjugate eigenvalues.
*
*  Arguments
*  =========
*
*  A       (input/output) REAL            
*  B       (input/output) REAL            
*  C       (input/output) REAL            
*  D       (input/output) REAL            
*          On entry, the elements of the input matrix.
*          On exit, they are overwritten by the elements of the
*          standardised Schur form.
*
*  RT1R    (output) REAL 
*  RT1I    (output) REAL            
*  RT2R    (output) REAL            
*  RT2I    (output) REAL            
*          The real and imaginary parts of the eigenvalues. If the
*          eigenvalues are a complex conjugate pair, RT1I > 0.
*
*  CS      (output) REAL            
*  SN      (output) REAL            
*          Parameters of the rotation matrix.
*
*  Further Details
*  ===============
*
*  Modified by V. Sima, Research Institute for Informatics, Bucharest,
*  Romania, to reduce the risk of cancellation errors,
*  when computing real eigenvalues, and to ensure, if possible, that
*  abs(RT1R) >= abs(RT2R).
*
*  =====================================================================
*
      SUBROUTINE SLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
      IMPLICIT NONE
      REAL               A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
      REAL               ZERO, HALF, ONE, MULTPL
      PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
      PARAMETER          ( MULTPL = 4.0E+0 )
      REAL               AA, BB, CC, CS1, DD, EPS, P, Q,
     $                   SAB, SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z
      REAL               SLAPY2
*
      IF( B.EQ.ZERO .OR. C.EQ.ZERO ) THEN
         CS = ONE
         SN = ZERO
         GO TO 10
      ELSE IF( (A-D).EQ.ZERO .AND. SIGN( ONE, B ).NE.
     $   SIGN( ONE, C ) ) THEN
         CS = ONE
         SN = ZERO
         GO TO 10
      ELSE
         P = HALF*(A - D)
         Q = HALF*(A + D)
         EPS = EPSILON( ZERO )
*         SCALE = MAX( ABS( P ), BCMAX )
*         Z = ( P / SCALE )*P + ( BCMAX / SCALE )*BCMIS
         SCALE = ONE
         Z = P*P + B*C
*        If Z is of the order of the machine accuracy, postpone the
*        decision on the nature of eigenvalues
         IF( Z.GE.MULTPL*EPS ) THEN
*           Real eigenvalues. Compute A and D.
            Z = P + SIGN( SQRT( SCALE )*SQRT( Z ), P )
            A = D + Z
            D = D - ( ( B*C ) / Z )
*           Compute B and the rotation matrix
*            TAU = ONE / SLAPY2( C, Z )
*            Z = Q + SQRT( Z )
            TAU = ONE / SQRT( C*C + Z*Z )
            CS = Z * TAU
            SN = C * TAU
            B = B - C
            C = ZERO
         ELSE
*           Complex eigenvalues, or real (almost) equal eigenvalues.
*           Make diagonal elements equal.
            TEMP  = A - D
            SIGMA = B + C
*            TAU = SLAPY2( SIGMA, TEMP )
            TAU = SQRT( SIGMA*SIGMA + TEMP*TEMP )
            CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
            SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA )
*
*           Compute [ AA  BB ] = [ A  B ] [ CS -SN ]
*                   [ CC  DD ]   [ C  D ] [ SN  CS ]
*
            AA = A*CS + B*SN
            BB = B*CS - A*SN
            CC = C*CS + D*SN
            DD = D*CS - C*SN
*
*           Compute [ A  B ] = [ CS  SN ] [ AA  BB ]
*                   [ C  D ]   [-SN  CS ] [ CC  DD ]
*
            A = AA*CS + CC*SN
            B = BB*CS + DD*SN
            C = CC*CS - AA*SN
            D = DD*CS - BB*SN
*
            TEMP = HALF*( A+D )
            A = TEMP
            D = TEMP
*
            IF( C.NE.ZERO ) THEN
               IF( B.NE.ZERO ) THEN
                  IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
*                    Real eigenvalues: reduce to upper triangular form
                     SAB = SQRT( ABS( B ) )
                     SAC = SQRT( ABS( C ) )
                     P = SIGN( SAB*SAC, C )
                     TAU = ONE / SQRT( ABS( B+C ) )
                     A = TEMP + P
                     D = TEMP - P
                     B = B - C
                     C = ZERO
                     CS1 = SAB*TAU
                     SN1 = SAC*TAU
                     TEMP = CS*CS1 - SN*SN1
                     SN =   CS*SN1 + SN*CS1
                     CS = TEMP
                  END IF
               END IF
            END IF
         END IF
      END IF
*
   10 CONTINUE
*     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
      RT1R = A
      RT2R = D
      IF( C.EQ.ZERO ) THEN
         RT1I = ZERO
      ELSE
         RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
      END IF
      RT2I = -RT1I
      RETURN
*     End of SLANV2
      END

