stdlib_lapack_eigv_gen3.fypp Source File


Source Code

#:include "common.fypp" 
submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_gen3
  implicit none


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     module subroutine stdlib${ii}$_slaqtr( ltran, lreal, n, t, ldt, b, w, scale, x, work,info )
     !! SLAQTR solves the real quasi-triangular system
     !! op(T)*p = scale*c,               if LREAL = .TRUE.
     !! or the complex quasi-triangular systems
     !! op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.
     !! in real arithmetic, where T is upper quasi-triangular.
     !! If LREAL = .FALSE., then the first diagonal block of T must be
     !! 1 by 1, B is the specially structured matrix
     !! B = [ b(1) b(2) ... b(n) ]
     !! [       w            ]
     !! [           w        ]
     !! [              .     ]
     !! [                 w  ]
     !! op(A) = A or A**T, A**T denotes the transpose of
     !! matrix A.
     !! On input, X = [ c ].  On output, X = [ p ].
     !! [ d ]                  [ q ]
     !! This subroutine is designed for the condition number estimation
     !! in routine STRSNA.
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           logical(lk), intent(in) :: lreal, ltran
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldt, n
           real(sp), intent(out) :: scale
           real(sp), intent(in) :: w
           ! Array Arguments 
           real(sp), intent(in) :: b(*), t(ldt,*)
           real(sp), intent(out) :: work(*)
           real(sp), intent(inout) :: x(*)
       ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: notran
           integer(${ik}$) :: i, ierr, j, j1, j2, jnext, k, n1, n2
           real(sp) :: bignum, eps, rec, scaloc, si, smin, sminw, smlnum, sr, tjj, tmp, xj, xmax, &
                     xnorm, z
           ! Local Arrays 
           real(sp) :: d(2_${ik}$,2_${ik}$), v(2_${ik}$,2_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! do not test the input parameters for errors
           notran = .not.ltran
           info = 0_${ik}$
           ! quick return if possible
           if( n==0 )return
           ! set constants to control overflow
           eps = stdlib${ii}$_slamch( 'P' )
           smlnum = stdlib${ii}$_slamch( 'S' ) / eps
           bignum = one / smlnum
           xnorm = stdlib${ii}$_slange( 'M', n, n, t, ldt, d )
           if( .not.lreal )xnorm = max( xnorm, abs( w ), stdlib${ii}$_slange( 'M', n, 1_${ik}$, b, n, d ) )
                     
           smin = max( smlnum, eps*xnorm )
           ! compute 1-norm of each column of strictly upper triangular
           ! part of t to control overflow in triangular solver.
           work( 1_${ik}$ ) = zero
           do j = 2, n
              work( j ) = stdlib${ii}$_sasum( j-1, t( 1_${ik}$, j ), 1_${ik}$ )
           end do
           if( .not.lreal ) then
              do i = 2, n
                 work( i ) = work( i ) + abs( b( i ) )
              end do
           end if
           n2 = 2_${ik}$*n
           n1 = n
           if( .not.lreal )n1 = n2
           k = stdlib${ii}$_isamax( n1, x, 1_${ik}$ )
           xmax = abs( x( k ) )
           scale = one
           if( xmax>bignum ) then
              scale = bignum / xmax
              call stdlib${ii}$_sscal( n1, scale, x, 1_${ik}$ )
              xmax = bignum
           end if
           if( lreal ) then
              if( notran ) then
                 ! solve t*p = scale*c
                 jnext = n
                 loop_30: do j = n, 1, -1
                    if( j>jnext )cycle loop_30
                    j1 = j
                    j2 = j
                    jnext = j - 1_${ik}$
                    if( j>1_${ik}$ ) then
                       if( t( j, j-1 )/=zero ) then
                          j1 = j - 1_${ik}$
                          jnext = j - 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! meet 1 by 1 diagonal block
                       ! scale to avoid overflow when computing
                           ! x(j) = b(j)/t(j,j)
                       xj = abs( x( j1 ) )
                       tjj = abs( t( j1, j1 ) )
                       tmp = t( j1, j1 )
                       if( tjj<smin ) then
                          tmp = smin
                          tjj = smin
                          info = 1_${ik}$
                       end if
                       if( xj==zero )cycle loop_30
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) / tmp
                       xj = abs( x( j1 ) )
                       ! scale x if necessary to avoid overflow when adding a
                       ! multiple of column j1 of t.
                       if( xj>one ) then
                          rec = one / xj
                          if( work( j1 )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          k = stdlib${ii}$_isamax( j1-1, x, 1_${ik}$ )
                          xmax = abs( x( k ) )
                       end if
                    else
                       ! meet 2 by 2 diagonal block
                       ! call 2 by 2 linear system solve, to take
                       ! care of possible overflow by scaling factor.
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 )
                       call stdlib${ii}$_slaln2( .false., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d,&
                                  2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_sscal( n, scaloc, x, 1_${ik}$ )
                          scale = scale*scaloc
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       ! scale v(1,1) (= x(j1)) and/or v(2,1) (=x(j2))
                       ! to avoid overflow in updating right-hand side.
                       xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) ), abs( v( 2_${ik}$, 1_${ik}$ ) ) )
                       if( xj>one ) then
                          rec = one / xj
                          if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       ! update right-hand side
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_saxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ )
                          k = stdlib${ii}$_isamax( j1-1, x, 1_${ik}$ )
                          xmax = abs( x( k ) )
                       end if
                    end if
                 end do loop_30
              else
                 ! solve t**t*p = scale*c
                 jnext = 1_${ik}$
                 loop_40: do j = 1, n
                    if( j<jnext )cycle loop_40
                    j1 = j
                    j2 = j
                    jnext = j + 1_${ik}$
                    if( j<n ) then
                       if( t( j+1, j )/=zero ) then
                          j2 = j + 1_${ik}$
                          jnext = j + 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = abs( x( j1 ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( work( j1 )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                       xj = abs( x( j1 ) )
                       tjj = abs( t( j1, j1 ) )
                       tmp = t( j1, j1 )
                       if( tjj<smin ) then
                          tmp = smin
                          tjj = smin
                          info = 1_${ik}$
                       end if
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) / tmp
                       xmax = max( xmax, abs( x( j1 ) ) )
                    else
                       ! 2 by 2 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side elements by inner product.
                       xj = max( abs( x( j1 ) ), abs( x( j2 ) ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( max( work( j2 ), work( j1 ) )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ )
                       call stdlib${ii}$_slaln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d, &
                                 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_sscal( n, scaloc, x, 1_${ik}$ )
                          scale = scale*scaloc
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax )
                    end if
                 end do loop_40
              end if
           else
              sminw = max( eps*abs( w ), smin )
              if( notran ) then
                 ! solve (t + ib)*(p+iq) = c+id
                 jnext = n
                 loop_70: do j = n, 1, -1
                    if( j>jnext )cycle loop_70
                    j1 = j
                    j2 = j
                    jnext = j - 1_${ik}$
                    if( j>1_${ik}$ ) then
                       if( t( j, j-1 )/=zero ) then
                          j1 = j - 1_${ik}$
                          jnext = j - 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in division
                       z = w
                       if( j1==1_${ik}$ )z = b( 1_${ik}$ )
                       xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
                       tjj = abs( t( j1, j1 ) ) + abs( z )
                       tmp = t( j1, j1 )
                       if( tjj<sminw ) then
                          tmp = sminw
                          tjj = sminw
                          info = 1_${ik}$
                       end if
                       if( xj==zero )cycle loop_70
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       call stdlib${ii}$_sladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si )
                       x( j1 ) = sr
                       x( n+j1 ) = si
                       xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
                       ! scale x if necessary to avoid overflow when adding a
                       ! multiple of column j1 of t.
                       if( xj>one ) then
                          rec = one / xj
                          if( work( j1 )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_saxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 )
                          x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 )
                          xmax = zero
                          do k = 1, j1 - 1
                             xmax = max( xmax, abs( x( k ) )+abs( x( k+n ) ) )
                          end do
                       end if
                    else
                       ! meet 2 by 2 diagonal block
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 )
                       d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 )
                       d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 )
                       call stdlib${ii}$_slaln2( .false., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, &
                                 d, 2_${ik}$, zero, -w, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_sscal( 2_${ik}$*n, scaloc, x, 1_${ik}$ )
                          scale = scaloc*scale
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ )
                       x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ )
                       ! scale x(j1), .... to avoid overflow in
                       ! updating right hand side.
                       xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) )+abs( v( 1_${ik}$, 2_${ik}$ ) ),abs( v( 2_${ik}$, 1_${ik}$ ) )+abs( v( 2_${ik}$, 2_${ik}$ )&
                                  ) )
                       if( xj>one ) then
                          rec = one / xj
                          if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       ! update the right-hand side.
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_saxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_saxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          call stdlib${ii}$_saxpy( j1-1, -x( n+j2 ), t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) +b( j2 )*x( n+j2 )
                          x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -b( j2 )*x( j2 )
                          xmax = zero
                          do k = 1, j1 - 1
                             xmax = max( abs( x( k ) )+abs( x( k+n ) ),xmax )
                          end do
                       end if
                    end if
                 end do loop_70
              else
                 ! solve (t + ib)**t*(p+iq) = c+id
                 jnext = 1_${ik}$
                 loop_80: do j = 1, n
                    if( j<jnext )cycle loop_80
                    j1 = j
                    j2 = j
                    jnext = j + 1_${ik}$
                    if( j<n ) then
                       if( t( j+1, j )/=zero ) then
                          j2 = j + 1_${ik}$
                          jnext = j + 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = abs( x( j1 ) ) + abs( x( j1+n ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( work( j1 )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                       x( n+j1 ) = x( n+j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       if( j1>1_${ik}$ ) then
                          x( j1 ) = x( j1 ) - b( j1 )*x( n+1 )
                          x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1_${ik}$ )
                       end if
                       xj = abs( x( j1 ) ) + abs( x( j1+n ) )
                       z = w
                       if( j1==1_${ik}$ )z = b( 1_${ik}$ )
                       ! scale if necessary to avoid overflow in
                       ! complex division
                       tjj = abs( t( j1, j1 ) ) + abs( z )
                       tmp = t( j1, j1 )
                       if( tjj<sminw ) then
                          tmp = sminw
                          tjj = sminw
                          info = 1_${ik}$
                       end if
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       call stdlib${ii}$_sladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si )
                       x( j1 ) = sr
                       x( j1+n ) = si
                       xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax )
                    else
                       ! 2 by 2 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) ) )
                                 
                       if( xmax>one ) then
                          rec = one / xmax
                          if( max( work( j1 ), work( j2 ) )>( bignum-xj ) / xmax ) then
                             call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ )
                       d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       d( 1_${ik}$, 1_${ik}$ ) = d( 1_${ik}$, 1_${ik}$ ) - b( j1 )*x( n+1 )
                       d( 2_${ik}$, 1_${ik}$ ) = d( 2_${ik}$, 1_${ik}$ ) - b( j2 )*x( n+1 )
                       d( 1_${ik}$, 2_${ik}$ ) = d( 1_${ik}$, 2_${ik}$ ) + b( j1 )*x( 1_${ik}$ )
                       d( 2_${ik}$, 2_${ik}$ ) = d( 2_${ik}$, 2_${ik}$ ) + b( j2 )*x( 1_${ik}$ )
                       call stdlib${ii}$_slaln2( .true., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, d,&
                                  2_${ik}$, zero, w, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_sscal( n2, scaloc, x, 1_${ik}$ )
                          scale = scaloc*scale
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ )
                       x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ )
                       xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) )&
                                 , xmax )
                    end if
                 end do loop_80
              end if
           end if
           return
     end subroutine stdlib${ii}$_slaqtr

     module subroutine stdlib${ii}$_dlaqtr( ltran, lreal, n, t, ldt, b, w, scale, x, work,info )
     !! DLAQTR solves the real quasi-triangular system
     !! op(T)*p = scale*c,               if LREAL = .TRUE.
     !! or the complex quasi-triangular systems
     !! op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.
     !! in real arithmetic, where T is upper quasi-triangular.
     !! If LREAL = .FALSE., then the first diagonal block of T must be
     !! 1 by 1, B is the specially structured matrix
     !! B = [ b(1) b(2) ... b(n) ]
     !! [       w            ]
     !! [           w        ]
     !! [              .     ]
     !! [                 w  ]
     !! op(A) = A or A**T, A**T denotes the transpose of
     !! matrix A.
     !! On input, X = [ c ].  On output, X = [ p ].
     !! [ d ]                  [ q ]
     !! This subroutine is designed for the condition number estimation
     !! in routine DTRSNA.
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           logical(lk), intent(in) :: lreal, ltran
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldt, n
           real(dp), intent(out) :: scale
           real(dp), intent(in) :: w
           ! Array Arguments 
           real(dp), intent(in) :: b(*), t(ldt,*)
           real(dp), intent(out) :: work(*)
           real(dp), intent(inout) :: x(*)
       ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: notran
           integer(${ik}$) :: i, ierr, j, j1, j2, jnext, k, n1, n2
           real(dp) :: bignum, eps, rec, scaloc, si, smin, sminw, smlnum, sr, tjj, tmp, xj, xmax, &
                     xnorm, z
           ! Local Arrays 
           real(dp) :: d(2_${ik}$,2_${ik}$), v(2_${ik}$,2_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! do not test the input parameters for errors
           notran = .not.ltran
           info = 0_${ik}$
           ! quick return if possible
           if( n==0 )return
           ! set constants to control overflow
           eps = stdlib${ii}$_dlamch( 'P' )
           smlnum = stdlib${ii}$_dlamch( 'S' ) / eps
           bignum = one / smlnum
           xnorm = stdlib${ii}$_dlange( 'M', n, n, t, ldt, d )
           if( .not.lreal )xnorm = max( xnorm, abs( w ), stdlib${ii}$_dlange( 'M', n, 1_${ik}$, b, n, d ) )
                     
           smin = max( smlnum, eps*xnorm )
           ! compute 1-norm of each column of strictly upper triangular
           ! part of t to control overflow in triangular solver.
           work( 1_${ik}$ ) = zero
           do j = 2, n
              work( j ) = stdlib${ii}$_dasum( j-1, t( 1_${ik}$, j ), 1_${ik}$ )
           end do
           if( .not.lreal ) then
              do i = 2, n
                 work( i ) = work( i ) + abs( b( i ) )
              end do
           end if
           n2 = 2_${ik}$*n
           n1 = n
           if( .not.lreal )n1 = n2
           k = stdlib${ii}$_idamax( n1, x, 1_${ik}$ )
           xmax = abs( x( k ) )
           scale = one
           if( xmax>bignum ) then
              scale = bignum / xmax
              call stdlib${ii}$_dscal( n1, scale, x, 1_${ik}$ )
              xmax = bignum
           end if
           if( lreal ) then
              if( notran ) then
                 ! solve t*p = scale*c
                 jnext = n
                 loop_30: do j = n, 1, -1
                    if( j>jnext )cycle loop_30
                    j1 = j
                    j2 = j
                    jnext = j - 1_${ik}$
                    if( j>1_${ik}$ ) then
                       if( t( j, j-1 )/=zero ) then
                          j1 = j - 1_${ik}$
                          jnext = j - 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! meet 1 by 1 diagonal block
                       ! scale to avoid overflow when computing
                           ! x(j) = b(j)/t(j,j)
                       xj = abs( x( j1 ) )
                       tjj = abs( t( j1, j1 ) )
                       tmp = t( j1, j1 )
                       if( tjj<smin ) then
                          tmp = smin
                          tjj = smin
                          info = 1_${ik}$
                       end if
                       if( xj==zero )cycle loop_30
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) / tmp
                       xj = abs( x( j1 ) )
                       ! scale x if necessary to avoid overflow when adding a
                       ! multiple of column j1 of t.
                       if( xj>one ) then
                          rec = one / xj
                          if( work( j1 )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          k = stdlib${ii}$_idamax( j1-1, x, 1_${ik}$ )
                          xmax = abs( x( k ) )
                       end if
                    else
                       ! meet 2 by 2 diagonal block
                       ! call 2 by 2 linear system solve, to take
                       ! care of possible overflow by scaling factor.
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 )
                       call stdlib${ii}$_dlaln2( .false., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d,&
                                  2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_dscal( n, scaloc, x, 1_${ik}$ )
                          scale = scale*scaloc
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       ! scale v(1,1) (= x(j1)) and/or v(2,1) (=x(j2))
                       ! to avoid overflow in updating right-hand side.
                       xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) ), abs( v( 2_${ik}$, 1_${ik}$ ) ) )
                       if( xj>one ) then
                          rec = one / xj
                          if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       ! update right-hand side
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_daxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ )
                          k = stdlib${ii}$_idamax( j1-1, x, 1_${ik}$ )
                          xmax = abs( x( k ) )
                       end if
                    end if
                 end do loop_30
              else
                 ! solve t**t*p = scale*c
                 jnext = 1_${ik}$
                 loop_40: do j = 1, n
                    if( j<jnext )cycle loop_40
                    j1 = j
                    j2 = j
                    jnext = j + 1_${ik}$
                    if( j<n ) then
                       if( t( j+1, j )/=zero ) then
                          j2 = j + 1_${ik}$
                          jnext = j + 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = abs( x( j1 ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( work( j1 )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                       xj = abs( x( j1 ) )
                       tjj = abs( t( j1, j1 ) )
                       tmp = t( j1, j1 )
                       if( tjj<smin ) then
                          tmp = smin
                          tjj = smin
                          info = 1_${ik}$
                       end if
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) / tmp
                       xmax = max( xmax, abs( x( j1 ) ) )
                    else
                       ! 2 by 2 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side elements by inner product.
                       xj = max( abs( x( j1 ) ), abs( x( j2 ) ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( max( work( j2 ), work( j1 ) )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ )
                       call stdlib${ii}$_dlaln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d, &
                                 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_dscal( n, scaloc, x, 1_${ik}$ )
                          scale = scale*scaloc
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax )
                    end if
                 end do loop_40
              end if
           else
              sminw = max( eps*abs( w ), smin )
              if( notran ) then
                 ! solve (t + ib)*(p+iq) = c+id
                 jnext = n
                 loop_70: do j = n, 1, -1
                    if( j>jnext )cycle loop_70
                    j1 = j
                    j2 = j
                    jnext = j - 1_${ik}$
                    if( j>1_${ik}$ ) then
                       if( t( j, j-1 )/=zero ) then
                          j1 = j - 1_${ik}$
                          jnext = j - 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in division
                       z = w
                       if( j1==1_${ik}$ )z = b( 1_${ik}$ )
                       xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
                       tjj = abs( t( j1, j1 ) ) + abs( z )
                       tmp = t( j1, j1 )
                       if( tjj<sminw ) then
                          tmp = sminw
                          tjj = sminw
                          info = 1_${ik}$
                       end if
                       if( xj==zero )cycle loop_70
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       call stdlib${ii}$_dladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si )
                       x( j1 ) = sr
                       x( n+j1 ) = si
                       xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
                       ! scale x if necessary to avoid overflow when adding a
                       ! multiple of column j1 of t.
                       if( xj>one ) then
                          rec = one / xj
                          if( work( j1 )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_daxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 )
                          x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 )
                          xmax = zero
                          do k = 1, j1 - 1
                             xmax = max( xmax, abs( x( k ) )+abs( x( k+n ) ) )
                          end do
                       end if
                    else
                       ! meet 2 by 2 diagonal block
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 )
                       d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 )
                       d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 )
                       call stdlib${ii}$_dlaln2( .false., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, &
                                 d, 2_${ik}$, zero, -w, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_dscal( 2_${ik}$*n, scaloc, x, 1_${ik}$ )
                          scale = scaloc*scale
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ )
                       x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ )
                       ! scale x(j1), .... to avoid overflow in
                       ! updating right hand side.
                       xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) )+abs( v( 1_${ik}$, 2_${ik}$ ) ),abs( v( 2_${ik}$, 1_${ik}$ ) )+abs( v( 2_${ik}$, 2_${ik}$ )&
                                  ) )
                       if( xj>one ) then
                          rec = one / xj
                          if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       ! update the right-hand side.
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_daxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_daxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          call stdlib${ii}$_daxpy( j1-1, -x( n+j2 ), t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) +b( j2 )*x( n+j2 )
                          x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -b( j2 )*x( j2 )
                          xmax = zero
                          do k = 1, j1 - 1
                             xmax = max( abs( x( k ) )+abs( x( k+n ) ),xmax )
                          end do
                       end if
                    end if
                 end do loop_70
              else
                 ! solve (t + ib)**t*(p+iq) = c+id
                 jnext = 1_${ik}$
                 loop_80: do j = 1, n
                    if( j<jnext )cycle loop_80
                    j1 = j
                    j2 = j
                    jnext = j + 1_${ik}$
                    if( j<n ) then
                       if( t( j+1, j )/=zero ) then
                          j2 = j + 1_${ik}$
                          jnext = j + 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = abs( x( j1 ) ) + abs( x( j1+n ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( work( j1 )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                       x( n+j1 ) = x( n+j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       if( j1>1_${ik}$ ) then
                          x( j1 ) = x( j1 ) - b( j1 )*x( n+1 )
                          x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1_${ik}$ )
                       end if
                       xj = abs( x( j1 ) ) + abs( x( j1+n ) )
                       z = w
                       if( j1==1_${ik}$ )z = b( 1_${ik}$ )
                       ! scale if necessary to avoid overflow in
                       ! complex division
                       tjj = abs( t( j1, j1 ) ) + abs( z )
                       tmp = t( j1, j1 )
                       if( tjj<sminw ) then
                          tmp = sminw
                          tjj = sminw
                          info = 1_${ik}$
                       end if
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       call stdlib${ii}$_dladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si )
                       x( j1 ) = sr
                       x( j1+n ) = si
                       xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax )
                    else
                       ! 2 by 2 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) ) )
                                 
                       if( xmax>one ) then
                          rec = one / xmax
                          if( max( work( j1 ), work( j2 ) )>( bignum-xj ) / xmax ) then
                             call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ )
                       d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       d( 1_${ik}$, 1_${ik}$ ) = d( 1_${ik}$, 1_${ik}$ ) - b( j1 )*x( n+1 )
                       d( 2_${ik}$, 1_${ik}$ ) = d( 2_${ik}$, 1_${ik}$ ) - b( j2 )*x( n+1 )
                       d( 1_${ik}$, 2_${ik}$ ) = d( 1_${ik}$, 2_${ik}$ ) + b( j1 )*x( 1_${ik}$ )
                       d( 2_${ik}$, 2_${ik}$ ) = d( 2_${ik}$, 2_${ik}$ ) + b( j2 )*x( 1_${ik}$ )
                       call stdlib${ii}$_dlaln2( .true., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, d,&
                                  2_${ik}$, zero, w, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_dscal( n2, scaloc, x, 1_${ik}$ )
                          scale = scaloc*scale
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ )
                       x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ )
                       xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) )&
                                 , xmax )
                    end if
                 end do loop_80
              end if
           end if
           return
     end subroutine stdlib${ii}$_dlaqtr

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$laqtr( ltran, lreal, n, t, ldt, b, w, scale, x, work,info )
     !! DLAQTR: solves the real quasi-triangular system
     !! op(T)*p = scale*c,               if LREAL = .TRUE.
     !! or the complex quasi-triangular systems
     !! op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.
     !! in real arithmetic, where T is upper quasi-triangular.
     !! If LREAL = .FALSE., then the first diagonal block of T must be
     !! 1 by 1, B is the specially structured matrix
     !! B = [ b(1) b(2) ... b(n) ]
     !! [       w            ]
     !! [           w        ]
     !! [              .     ]
     !! [                 w  ]
     !! op(A) = A or A**T, A**T denotes the transpose of
     !! matrix A.
     !! On input, X = [ c ].  On output, X = [ p ].
     !! [ d ]                  [ q ]
     !! This subroutine is designed for the condition number estimation
     !! in routine DTRSNA.
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           logical(lk), intent(in) :: lreal, ltran
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldt, n
           real(${rk}$), intent(out) :: scale
           real(${rk}$), intent(in) :: w
           ! Array Arguments 
           real(${rk}$), intent(in) :: b(*), t(ldt,*)
           real(${rk}$), intent(out) :: work(*)
           real(${rk}$), intent(inout) :: x(*)
       ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: notran
           integer(${ik}$) :: i, ierr, j, j1, j2, jnext, k, n1, n2
           real(${rk}$) :: bignum, eps, rec, scaloc, si, smin, sminw, smlnum, sr, tjj, tmp, xj, xmax, &
                     xnorm, z
           ! Local Arrays 
           real(${rk}$) :: d(2_${ik}$,2_${ik}$), v(2_${ik}$,2_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! do not test the input parameters for errors
           notran = .not.ltran
           info = 0_${ik}$
           ! quick return if possible
           if( n==0 )return
           ! set constants to control overflow
           eps = stdlib${ii}$_${ri}$lamch( 'P' )
           smlnum = stdlib${ii}$_${ri}$lamch( 'S' ) / eps
           bignum = one / smlnum
           xnorm = stdlib${ii}$_${ri}$lange( 'M', n, n, t, ldt, d )
           if( .not.lreal )xnorm = max( xnorm, abs( w ), stdlib${ii}$_${ri}$lange( 'M', n, 1_${ik}$, b, n, d ) )
                     
           smin = max( smlnum, eps*xnorm )
           ! compute 1-norm of each column of strictly upper triangular
           ! part of t to control overflow in triangular solver.
           work( 1_${ik}$ ) = zero
           do j = 2, n
              work( j ) = stdlib${ii}$_${ri}$asum( j-1, t( 1_${ik}$, j ), 1_${ik}$ )
           end do
           if( .not.lreal ) then
              do i = 2, n
                 work( i ) = work( i ) + abs( b( i ) )
              end do
           end if
           n2 = 2_${ik}$*n
           n1 = n
           if( .not.lreal )n1 = n2
           k = stdlib${ii}$_i${ri}$amax( n1, x, 1_${ik}$ )
           xmax = abs( x( k ) )
           scale = one
           if( xmax>bignum ) then
              scale = bignum / xmax
              call stdlib${ii}$_${ri}$scal( n1, scale, x, 1_${ik}$ )
              xmax = bignum
           end if
           if( lreal ) then
              if( notran ) then
                 ! solve t*p = scale*c
                 jnext = n
                 loop_30: do j = n, 1, -1
                    if( j>jnext )cycle loop_30
                    j1 = j
                    j2 = j
                    jnext = j - 1_${ik}$
                    if( j>1_${ik}$ ) then
                       if( t( j, j-1 )/=zero ) then
                          j1 = j - 1_${ik}$
                          jnext = j - 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! meet 1 by 1 diagonal block
                       ! scale to avoid overflow when computing
                           ! x(j) = b(j)/t(j,j)
                       xj = abs( x( j1 ) )
                       tjj = abs( t( j1, j1 ) )
                       tmp = t( j1, j1 )
                       if( tjj<smin ) then
                          tmp = smin
                          tjj = smin
                          info = 1_${ik}$
                       end if
                       if( xj==zero )cycle loop_30
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) / tmp
                       xj = abs( x( j1 ) )
                       ! scale x if necessary to avoid overflow when adding a
                       ! multiple of column j1 of t.
                       if( xj>one ) then
                          rec = one / xj
                          if( work( j1 )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          k = stdlib${ii}$_i${ri}$amax( j1-1, x, 1_${ik}$ )
                          xmax = abs( x( k ) )
                       end if
                    else
                       ! meet 2 by 2 diagonal block
                       ! call 2 by 2 linear system solve, to take
                       ! care of possible overflow by scaling factor.
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 )
                       call stdlib${ii}$_${ri}$laln2( .false., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d,&
                                  2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_${ri}$scal( n, scaloc, x, 1_${ik}$ )
                          scale = scale*scaloc
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       ! scale v(1,1) (= x(j1)) and/or v(2,1) (=x(j2))
                       ! to avoid overflow in updating right-hand side.
                       xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) ), abs( v( 2_${ik}$, 1_${ik}$ ) ) )
                       if( xj>one ) then
                          rec = one / xj
                          if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       ! update right-hand side
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ )
                          k = stdlib${ii}$_i${ri}$amax( j1-1, x, 1_${ik}$ )
                          xmax = abs( x( k ) )
                       end if
                    end if
                 end do loop_30
              else
                 ! solve t**t*p = scale*c
                 jnext = 1_${ik}$
                 loop_40: do j = 1, n
                    if( j<jnext )cycle loop_40
                    j1 = j
                    j2 = j
                    jnext = j + 1_${ik}$
                    if( j<n ) then
                       if( t( j+1, j )/=zero ) then
                          j2 = j + 1_${ik}$
                          jnext = j + 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = abs( x( j1 ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( work( j1 )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                       xj = abs( x( j1 ) )
                       tjj = abs( t( j1, j1 ) )
                       tmp = t( j1, j1 )
                       if( tjj<smin ) then
                          tmp = smin
                          tjj = smin
                          info = 1_${ik}$
                       end if
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) / tmp
                       xmax = max( xmax, abs( x( j1 ) ) )
                    else
                       ! 2 by 2 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side elements by inner product.
                       xj = max( abs( x( j1 ) ), abs( x( j2 ) ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( max( work( j2 ), work( j1 ) )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ )
                       call stdlib${ii}$_${ri}$laln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d, &
                                 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_${ri}$scal( n, scaloc, x, 1_${ik}$ )
                          scale = scale*scaloc
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax )
                    end if
                 end do loop_40
              end if
           else
              sminw = max( eps*abs( w ), smin )
              if( notran ) then
                 ! solve (t + ib)*(p+iq) = c+id
                 jnext = n
                 loop_70: do j = n, 1, -1
                    if( j>jnext )cycle loop_70
                    j1 = j
                    j2 = j
                    jnext = j - 1_${ik}$
                    if( j>1_${ik}$ ) then
                       if( t( j, j-1 )/=zero ) then
                          j1 = j - 1_${ik}$
                          jnext = j - 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in division
                       z = w
                       if( j1==1_${ik}$ )z = b( 1_${ik}$ )
                       xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
                       tjj = abs( t( j1, j1 ) ) + abs( z )
                       tmp = t( j1, j1 )
                       if( tjj<sminw ) then
                          tmp = sminw
                          tjj = sminw
                          info = 1_${ik}$
                       end if
                       if( xj==zero )cycle loop_70
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       call stdlib${ii}$_${ri}$ladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si )
                       x( j1 ) = sr
                       x( n+j1 ) = si
                       xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
                       ! scale x if necessary to avoid overflow when adding a
                       ! multiple of column j1 of t.
                       if( xj>one ) then
                          rec = one / xj
                          if( work( j1 )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 )
                          x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 )
                          xmax = zero
                          do k = 1, j1 - 1
                             xmax = max( xmax, abs( x( k ) )+abs( x( k+n ) ) )
                          end do
                       end if
                    else
                       ! meet 2 by 2 diagonal block
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 )
                       d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 )
                       d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 )
                       call stdlib${ii}$_${ri}$laln2( .false., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, &
                                 d, 2_${ik}$, zero, -w, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_${ri}$scal( 2_${ik}$*n, scaloc, x, 1_${ik}$ )
                          scale = scaloc*scale
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ )
                       x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ )
                       ! scale x(j1), .... to avoid overflow in
                       ! updating right hand side.
                       xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) )+abs( v( 1_${ik}$, 2_${ik}$ ) ),abs( v( 2_${ik}$, 1_${ik}$ ) )+abs( v( 2_${ik}$, 2_${ik}$ )&
                                  ) )
                       if( xj>one ) then
                          rec = one / xj
                          if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then
                             call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                          end if
                       end if
                       ! update the right-hand side.
                       if( j1>1_${ik}$ ) then
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ )
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          call stdlib${ii}$_${ri}$axpy( j1-1, -x( n+j2 ), t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                          x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) +b( j2 )*x( n+j2 )
                          x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -b( j2 )*x( j2 )
                          xmax = zero
                          do k = 1, j1 - 1
                             xmax = max( abs( x( k ) )+abs( x( k+n ) ),xmax )
                          end do
                       end if
                    end if
                 end do loop_70
              else
                 ! solve (t + ib)**t*(p+iq) = c+id
                 jnext = 1_${ik}$
                 loop_80: do j = 1, n
                    if( j<jnext )cycle loop_80
                    j1 = j
                    j2 = j
                    jnext = j + 1_${ik}$
                    if( j<n ) then
                       if( t( j+1, j )/=zero ) then
                          j2 = j + 1_${ik}$
                          jnext = j + 2_${ik}$
                       end if
                    end if
                    if( j1==j2 ) then
                       ! 1 by 1 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = abs( x( j1 ) ) + abs( x( j1+n ) )
                       if( xmax>one ) then
                          rec = one / xmax
                          if( work( j1 )>( bignum-xj )*rec ) then
                             call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       x( j1 ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ )
                       x( n+j1 ) = x( n+j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       if( j1>1_${ik}$ ) then
                          x( j1 ) = x( j1 ) - b( j1 )*x( n+1 )
                          x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1_${ik}$ )
                       end if
                       xj = abs( x( j1 ) ) + abs( x( j1+n ) )
                       z = w
                       if( j1==1_${ik}$ )z = b( 1_${ik}$ )
                       ! scale if necessary to avoid overflow in
                       ! complex division
                       tjj = abs( t( j1, j1 ) ) + abs( z )
                       tmp = t( j1, j1 )
                       if( tjj<sminw ) then
                          tmp = sminw
                          tjj = sminw
                          info = 1_${ik}$
                       end if
                       if( tjj<one ) then
                          if( xj>bignum*tjj ) then
                             rec = one / xj
                             call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       call stdlib${ii}$_${ri}$ladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si )
                       x( j1 ) = sr
                       x( j1+n ) = si
                       xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax )
                    else
                       ! 2 by 2 diagonal block
                       ! scale if necessary to avoid overflow in forming the
                       ! right-hand side element by inner product.
                       xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) ) )
                                 
                       if( xmax>one ) then
                          rec = one / xmax
                          if( max( work( j1 ), work( j2 ) )>( bignum-xj ) / xmax ) then
                             call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ )
                             scale = scale*rec
                             xmax = xmax*rec
                          end if
                       end if
                       d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ )
                       d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ )
                       d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ )
                                 
                       d( 1_${ik}$, 1_${ik}$ ) = d( 1_${ik}$, 1_${ik}$ ) - b( j1 )*x( n+1 )
                       d( 2_${ik}$, 1_${ik}$ ) = d( 2_${ik}$, 1_${ik}$ ) - b( j2 )*x( n+1 )
                       d( 1_${ik}$, 2_${ik}$ ) = d( 1_${ik}$, 2_${ik}$ ) + b( j1 )*x( 1_${ik}$ )
                       d( 2_${ik}$, 2_${ik}$ ) = d( 2_${ik}$, 2_${ik}$ ) + b( j2 )*x( 1_${ik}$ )
                       call stdlib${ii}$_${ri}$laln2( .true., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, d,&
                                  2_${ik}$, zero, w, v, 2_${ik}$,scaloc, xnorm, ierr )
                       if( ierr/=0_${ik}$ )info = 2_${ik}$
                       if( scaloc/=one ) then
                          call stdlib${ii}$_${ri}$scal( n2, scaloc, x, 1_${ik}$ )
                          scale = scaloc*scale
                       end if
                       x( j1 ) = v( 1_${ik}$, 1_${ik}$ )
                       x( j2 ) = v( 2_${ik}$, 1_${ik}$ )
                       x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ )
                       x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ )
                       xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) )&
                                 , xmax )
                    end if
                 end do loop_80
              end if
           end if
           return
     end subroutine stdlib${ii}$_${ri}$laqtr

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,iloz, ihiz, z, ldz, &
     !! 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.
               info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n
           integer(${ik}$), intent(out) :: info
           logical(lk), intent(in) :: wantt, wantz
           ! Array Arguments 
           real(sp), intent(inout) :: h(ldh,*), z(ldz,*)
           real(sp), intent(out) :: wi(*), wr(*)
        ! =========================================================
           ! Parameters 
           real(sp), parameter :: dat1 = 3.0_sp/4.0_sp
           real(sp), parameter :: dat2 = -0.4375_sp
           integer(${ik}$), parameter :: kexsh = 10_${ik}$
           
           
           
           ! Local Scalars 
           real(sp) :: aa, ab, ba, bb, cs, det, h11, h12, h21, h21s, h22, rt1i, rt1r, rt2i, rt2r, &
                     rtdisc, s, safmax, safmin, smlnum, sn, sum, t1, t2, t3, tr, tst, ulp, v2, v3
           integer(${ik}$) :: i, i1, i2, its, itmax, j, k, l, m, nh, nr, nz, kdefl
           ! Local Arrays 
           real(sp) :: v(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n==0 )return
           if( ilo==ihi ) then
              wr( ilo ) = h( ilo, ilo )
              wi( ilo ) = zero
              return
           end if
           ! ==== clear out the trash ====
           do j = ilo, ihi - 3
              h( j+2, j ) = zero
              h( j+3, j ) = zero
           end do
           if( ilo<=ihi-2 )h( ihi, ihi-2 ) = zero
           nh = ihi - ilo + 1_${ik}$
           nz = ihiz - iloz + 1_${ik}$
           ! set machine-dependent constants for the stopping criterion.
           safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           safmax = one / safmin
           call stdlib${ii}$_slabad( safmin, safmax )
           ulp = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = safmin*( real( nh,KIND=sp) / 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.
           if( wantt ) then
              i1 = 1_${ik}$
              i2 = n
           end if
           ! itmax is the total number of qr iterations allowed.
           itmax = 30_${ik}$ * max( 10_${ik}$, nh )
           ! kdefl counts the number of iterations since a deflation
           kdefl = 0_${ik}$
           ! 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<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.
           loop_140: do its = 0, itmax
              ! look for a single small subdiagonal element.
              do k = i, l + 1, -1
                 if( abs( h( k, k-1 ) )<=smlnum )go to 40
                 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
                 if( tst==zero ) then
                    if( k-2>=ilo )tst = tst + abs( h( k-1, k-2 ) )
                    if( k+1<=ihi )tst = tst + abs( h( k+1, k ) )
                 end if
                 ! ==== the following is a conservative small subdiagonal
                 ! .    deflation  criterion due to ahues
                 ! .    1997). it has better mathematical foundation and
                 ! .    improves accuracy in some cases.  ====
                 if( abs( h( k, k-1 ) )<=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 )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 40
                 end if
              end do
              40 continue
              l = k
              if( l>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>=i-1 )go to 150
              kdefl = kdefl + 1_${ik}$
              ! now the active submatrix is in rows and columns l to i. if
              ! eigenvalues only are being computed, only the active submatrix
              ! need be transformed.
              if( .not.wantt ) then
                 i1 = l
                 i2 = i
              end if
              if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
                 h11 = dat1*s + h( i, i )
                 h12 = dat2*s
                 h21 = s
                 h22 = h11
              else if( mod(kdefl,kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
                 h11 = dat1*s + h( l, l )
                 h12 = dat2*s
                 h21 = s
                 h22 = h11
              else
                 ! prepare to use francis' double shift
                 ! (i.e. 2nd degree generalized rayleigh quotient)
                 h11 = h( i-1, i-1 )
                 h21 = h( i, i-1 )
                 h12 = h( i-1, i )
                 h22 = h( i, i )
              end if
              s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
              if( s==zero ) then
                 rt1r = zero
                 rt1i = zero
                 rt2r = zero
                 rt2i = zero
              else
                 h11 = h11 / s
                 h21 = h21 / s
                 h12 = h12 / s
                 h22 = h22 / s
                 tr = ( h11+h22 ) / two
                 det = ( h11-tr )*( h22-tr ) - h12*h21
                 rtdisc = sqrt( abs( det ) )
                 if( det>=zero ) then
                    ! ==== complex conjugate shifts ====
                    rt1r = tr*s
                    rt2r = rt1r
                    rt1i = rtdisc*s
                    rt2i = -rt1i
                 else
                    ! ==== realshifts (use only one of them,KIND=sp)  ====
                    rt1r = tr + rtdisc
                    rt2r = tr - rtdisc
                    if( abs( rt1r-h22 )<=abs( rt2r-h22 ) ) then
                       rt1r = rt1r*s
                       rt2r = rt1r
                    else
                       rt2r = rt2r*s
                       rt1r = rt2r
                    end if
                    rt1i = zero
                    rt2i = zero
                 end if
              end if
              ! look for two consecutive small subdiagonal elements.
              do 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 )
                 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
                 h21s = h( m+1, m ) / s
                 v( 1_${ik}$ ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*( ( h( m, m )-rt2r ) / s ) - &
                           rt1i*( rt2i / s )
                 v( 2_${ik}$ ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
                 v( 3_${ik}$ ) = h21s*h( m+2, m+1 )
                 s = abs( v( 1_${ik}$ ) ) + abs( v( 2_${ik}$ ) ) + abs( v( 3_${ik}$ ) )
                 v( 1_${ik}$ ) = v( 1_${ik}$ ) / s
                 v( 2_${ik}$ ) = v( 2_${ik}$ ) / s
                 v( 3_${ik}$ ) = v( 3_${ik}$ ) / s
                 if( m==l )go to 60
                 if( abs( h( m, m-1 ) )*( abs( v( 2_${ik}$ ) )+abs( v( 3_${ik}$ ) ) )<=ulp*abs( v( 1_${ik}$ ) )*( abs( &
                           h( m-1, m-1 ) )+abs( h( m,m ) )+abs( h( m+1, m+1 ) ) ) )go to 60
              end do
              60 continue
              ! double-shift qr step
              loop_130: do k = m, i - 1
                 ! 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. nr is the order of g.
                 nr = min( 3_${ik}$, i-k+1 )
                 if( k>m )call stdlib${ii}$_scopy( nr, h( k, k-1 ), 1_${ik}$, v, 1_${ik}$ )
                 call stdlib${ii}$_slarfg( nr, v( 1_${ik}$ ), v( 2_${ik}$ ), 1_${ik}$, t1 )
                 if( k>m ) then
                    h( k, k-1 ) = v( 1_${ik}$ )
                    h( k+1, k-1 ) = zero
                    if( k<i-1 )h( k+2, k-1 ) = zero
                 else if( m>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-t1 )
                 end if
                 v2 = v( 2_${ik}$ )
                 t2 = t1*v2
                 if( nr==3_${ik}$ ) then
                    v3 = v( 3_${ik}$ )
                    t3 = t1*v3
                    ! apply g from the left to transform the rows of the matrix
                    ! in columns k to i2.
                    do j = k, i2
                       sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
                       h( k, j ) = h( k, j ) - sum*t1
                       h( k+1, j ) = h( k+1, j ) - sum*t2
                       h( k+2, j ) = h( k+2, j ) - sum*t3
                    end do
                    ! apply g from the right to transform the columns of the
                    ! matrix in rows i1 to min(k+3,i).
                    do j = i1, min( k+3, i )
                       sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
                       h( j, k ) = h( j, k ) - sum*t1
                       h( j, k+1 ) = h( j, k+1 ) - sum*t2
                       h( j, k+2 ) = h( j, k+2 ) - sum*t3
                    end do
                    if( wantz ) then
                       ! accumulate transformations in the matrix z
                       do j = iloz, ihiz
                          sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
                          z( j, k ) = z( j, k ) - sum*t1
                          z( j, k+1 ) = z( j, k+1 ) - sum*t2
                          z( j, k+2 ) = z( j, k+2 ) - sum*t3
                       end do
                    end if
                 else if( nr==2_${ik}$ ) then
                    ! apply g from the left to transform the rows of the matrix
                    ! in columns k to i2.
                    do j = k, i2
                       sum = h( k, j ) + v2*h( k+1, j )
                       h( k, j ) = h( k, j ) - sum*t1
                       h( k+1, j ) = h( k+1, j ) - sum*t2
                    end do
                    ! apply g from the right to transform the columns of the
                    ! matrix in rows i1 to min(k+3,i).
                    do j = i1, i
                       sum = h( j, k ) + v2*h( j, k+1 )
                       h( j, k ) = h( j, k ) - sum*t1
                       h( j, k+1 ) = h( j, k+1 ) - sum*t2
                    end do
                    if( wantz ) then
                       ! accumulate transformations in the matrix z
                       do j = iloz, ihiz
                          sum = z( j, k ) + v2*z( j, k+1 )
                          z( j, k ) = z( j, k ) - sum*t1
                          z( j, k+1 ) = z( j, k+1 ) - sum*t2
                       end do
                    end if
                 end if
              end do loop_130
           end do loop_140
           ! failure to converge in remaining number of iterations
           info = i
           return
           150 continue
           if( l==i ) then
              ! h(i,i-1) is negligible: one eigenvalue has converged.
              wr( i ) = h( i, i )
              wi( i ) = zero
           else if( l==i-1 ) then
              ! h(i-1,i-2) is negligible: a pair of eigenvalues have converged.
              ! transform the 2-by-2 submatrix to standard schur form,
              ! and compute and store the eigenvalues.
              call stdlib${ii}$_slanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),h( i, i ), wr( i-1 ), &
                        wi( i-1 ), wr( i ), wi( i ),cs, sn )
              if( wantt ) then
                 ! apply the transformation to the rest of h.
                 if( i2>i )call stdlib${ii}$_srot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,cs, sn )
                           
                 call stdlib${ii}$_srot( i-i1-1, h( i1, i-1 ), 1_${ik}$, h( i1, i ), 1_${ik}$, cs, sn )
              end if
              if( wantz ) then
                 ! apply the transformation to z.
                 call stdlib${ii}$_srot( nz, z( iloz, i-1 ), 1_${ik}$, z( iloz, i ), 1_${ik}$, cs, sn )
              end if
           end if
           ! reset deflation counter
           kdefl = 0_${ik}$
           ! return to start of the main loop with new value of i.
           i = l - 1_${ik}$
           go to 20
           160 continue
           return
     end subroutine stdlib${ii}$_slahqr

     pure module subroutine stdlib${ii}$_dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,iloz, ihiz, z, ldz, &
     !! DLAHQR is an auxiliary routine called by DHSEQR to update the
     !! eigenvalues and Schur decomposition already computed by DHSEQR, by
     !! dealing with the Hessenberg submatrix in rows and columns ILO to
     !! IHI.
               info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n
           integer(${ik}$), intent(out) :: info
           logical(lk), intent(in) :: wantt, wantz
           ! Array Arguments 
           real(dp), intent(inout) :: h(ldh,*), z(ldz,*)
           real(dp), intent(out) :: wi(*), wr(*)
        ! =========================================================
           ! Parameters 
           real(dp), parameter :: dat1 = 3.0_dp/4.0_dp
           real(dp), parameter :: dat2 = -0.4375_dp
           integer(${ik}$), parameter :: kexsh = 10_${ik}$
           
           
           
           ! Local Scalars 
           real(dp) :: aa, ab, ba, bb, cs, det, h11, h12, h21, h21s, h22, rt1i, rt1r, rt2i, rt2r, &
                     rtdisc, s, safmax, safmin, smlnum, sn, sum, t1, t2, t3, tr, tst, ulp, v2, v3
           integer(${ik}$) :: i, i1, i2, its, itmax, j, k, l, m, nh, nr, nz, kdefl
           ! Local Arrays 
           real(dp) :: v(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n==0 )return
           if( ilo==ihi ) then
              wr( ilo ) = h( ilo, ilo )
              wi( ilo ) = zero
              return
           end if
           ! ==== clear out the trash ====
           do j = ilo, ihi - 3
              h( j+2, j ) = zero
              h( j+3, j ) = zero
           end do
           if( ilo<=ihi-2 )h( ihi, ihi-2 ) = zero
           nh = ihi - ilo + 1_${ik}$
           nz = ihiz - iloz + 1_${ik}$
           ! set machine-dependent constants for the stopping criterion.
           safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           safmax = one / safmin
           call stdlib${ii}$_dlabad( safmin, safmax )
           ulp = stdlib${ii}$_dlamch( 'PRECISION' )
           smlnum = safmin*( real( nh,KIND=dp) / 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.
           if( wantt ) then
              i1 = 1_${ik}$
              i2 = n
           end if
           ! itmax is the total number of qr iterations allowed.
           itmax = 30_${ik}$ * max( 10_${ik}$, nh )
           ! kdefl counts the number of iterations since a deflation
           kdefl = 0_${ik}$
           ! 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<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.
           loop_140: do its = 0, itmax
              ! look for a single small subdiagonal element.
              do k = i, l + 1, -1
                 if( abs( h( k, k-1 ) )<=smlnum )go to 40
                 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
                 if( tst==zero ) then
                    if( k-2>=ilo )tst = tst + abs( h( k-1, k-2 ) )
                    if( k+1<=ihi )tst = tst + abs( h( k+1, k ) )
                 end if
                 ! ==== the following is a conservative small subdiagonal
                 ! .    deflation  criterion due to ahues
                 ! .    1997). it has better mathematical foundation and
                 ! .    improves accuracy in some cases.  ====
                 if( abs( h( k, k-1 ) )<=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 )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 40
                 end if
              end do
              40 continue
              l = k
              if( l>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>=i-1 )go to 150
              kdefl = kdefl + 1_${ik}$
              ! now the active submatrix is in rows and columns l to i. if
              ! eigenvalues only are being computed, only the active submatrix
              ! need be transformed.
              if( .not.wantt ) then
                 i1 = l
                 i2 = i
              end if
              if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
                 h11 = dat1*s + h( i, i )
                 h12 = dat2*s
                 h21 = s
                 h22 = h11
              else if( mod(kdefl,kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
                 h11 = dat1*s + h( l, l )
                 h12 = dat2*s
                 h21 = s
                 h22 = h11
              else
                 ! prepare to use francis' double shift
                 ! (i.e. 2nd degree generalized rayleigh quotient)
                 h11 = h( i-1, i-1 )
                 h21 = h( i, i-1 )
                 h12 = h( i-1, i )
                 h22 = h( i, i )
              end if
              s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
              if( s==zero ) then
                 rt1r = zero
                 rt1i = zero
                 rt2r = zero
                 rt2i = zero
              else
                 h11 = h11 / s
                 h21 = h21 / s
                 h12 = h12 / s
                 h22 = h22 / s
                 tr = ( h11+h22 ) / two
                 det = ( h11-tr )*( h22-tr ) - h12*h21
                 rtdisc = sqrt( abs( det ) )
                 if( det>=zero ) then
                    ! ==== complex conjugate shifts ====
                    rt1r = tr*s
                    rt2r = rt1r
                    rt1i = rtdisc*s
                    rt2i = -rt1i
                 else
                    ! ==== realshifts (use only one of them,KIND=dp)  ====
                    rt1r = tr + rtdisc
                    rt2r = tr - rtdisc
                    if( abs( rt1r-h22 )<=abs( rt2r-h22 ) ) then
                       rt1r = rt1r*s
                       rt2r = rt1r
                    else
                       rt2r = rt2r*s
                       rt1r = rt2r
                    end if
                    rt1i = zero
                    rt2i = zero
                 end if
              end if
              ! look for two consecutive small subdiagonal elements.
              do 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 )
                 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
                 h21s = h( m+1, m ) / s
                 v( 1_${ik}$ ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*( ( h( m, m )-rt2r ) / s ) - &
                           rt1i*( rt2i / s )
                 v( 2_${ik}$ ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
                 v( 3_${ik}$ ) = h21s*h( m+2, m+1 )
                 s = abs( v( 1_${ik}$ ) ) + abs( v( 2_${ik}$ ) ) + abs( v( 3_${ik}$ ) )
                 v( 1_${ik}$ ) = v( 1_${ik}$ ) / s
                 v( 2_${ik}$ ) = v( 2_${ik}$ ) / s
                 v( 3_${ik}$ ) = v( 3_${ik}$ ) / s
                 if( m==l )go to 60
                 if( abs( h( m, m-1 ) )*( abs( v( 2_${ik}$ ) )+abs( v( 3_${ik}$ ) ) )<=ulp*abs( v( 1_${ik}$ ) )*( abs( &
                           h( m-1, m-1 ) )+abs( h( m,m ) )+abs( h( m+1, m+1 ) ) ) )go to 60
              end do
              60 continue
              ! double-shift qr step
              loop_130: do k = m, i - 1
                 ! 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. nr is the order of g.
                 nr = min( 3_${ik}$, i-k+1 )
                 if( k>m )call stdlib${ii}$_dcopy( nr, h( k, k-1 ), 1_${ik}$, v, 1_${ik}$ )
                 call stdlib${ii}$_dlarfg( nr, v( 1_${ik}$ ), v( 2_${ik}$ ), 1_${ik}$, t1 )
                 if( k>m ) then
                    h( k, k-1 ) = v( 1_${ik}$ )
                    h( k+1, k-1 ) = zero
                    if( k<i-1 )h( k+2, k-1 ) = zero
                 else if( m>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-t1 )
                 end if
                 v2 = v( 2_${ik}$ )
                 t2 = t1*v2
                 if( nr==3_${ik}$ ) then
                    v3 = v( 3_${ik}$ )
                    t3 = t1*v3
                    ! apply g from the left to transform the rows of the matrix
                    ! in columns k to i2.
                    do j = k, i2
                       sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
                       h( k, j ) = h( k, j ) - sum*t1
                       h( k+1, j ) = h( k+1, j ) - sum*t2
                       h( k+2, j ) = h( k+2, j ) - sum*t3
                    end do
                    ! apply g from the right to transform the columns of the
                    ! matrix in rows i1 to min(k+3,i).
                    do j = i1, min( k+3, i )
                       sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
                       h( j, k ) = h( j, k ) - sum*t1
                       h( j, k+1 ) = h( j, k+1 ) - sum*t2
                       h( j, k+2 ) = h( j, k+2 ) - sum*t3
                    end do
                    if( wantz ) then
                       ! accumulate transformations in the matrix z
                       do j = iloz, ihiz
                          sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
                          z( j, k ) = z( j, k ) - sum*t1
                          z( j, k+1 ) = z( j, k+1 ) - sum*t2
                          z( j, k+2 ) = z( j, k+2 ) - sum*t3
                       end do
                    end if
                 else if( nr==2_${ik}$ ) then
                    ! apply g from the left to transform the rows of the matrix
                    ! in columns k to i2.
                    do j = k, i2
                       sum = h( k, j ) + v2*h( k+1, j )
                       h( k, j ) = h( k, j ) - sum*t1
                       h( k+1, j ) = h( k+1, j ) - sum*t2
                    end do
                    ! apply g from the right to transform the columns of the
                    ! matrix in rows i1 to min(k+3,i).
                    do j = i1, i
                       sum = h( j, k ) + v2*h( j, k+1 )
                       h( j, k ) = h( j, k ) - sum*t1
                       h( j, k+1 ) = h( j, k+1 ) - sum*t2
                    end do
                    if( wantz ) then
                       ! accumulate transformations in the matrix z
                       do j = iloz, ihiz
                          sum = z( j, k ) + v2*z( j, k+1 )
                          z( j, k ) = z( j, k ) - sum*t1
                          z( j, k+1 ) = z( j, k+1 ) - sum*t2
                       end do
                    end if
                 end if
              end do loop_130
           end do loop_140
           ! failure to converge in remaining number of iterations
           info = i
           return
           150 continue
           if( l==i ) then
              ! h(i,i-1) is negligible: one eigenvalue has converged.
              wr( i ) = h( i, i )
              wi( i ) = zero
           else if( l==i-1 ) then
              ! h(i-1,i-2) is negligible: a pair of eigenvalues have converged.
              ! transform the 2-by-2 submatrix to standard schur form,
              ! and compute and store the eigenvalues.
              call stdlib${ii}$_dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),h( i, i ), wr( i-1 ), &
                        wi( i-1 ), wr( i ), wi( i ),cs, sn )
              if( wantt ) then
                 ! apply the transformation to the rest of h.
                 if( i2>i )call stdlib${ii}$_drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,cs, sn )
                           
                 call stdlib${ii}$_drot( i-i1-1, h( i1, i-1 ), 1_${ik}$, h( i1, i ), 1_${ik}$, cs, sn )
              end if
              if( wantz ) then
                 ! apply the transformation to z.
                 call stdlib${ii}$_drot( nz, z( iloz, i-1 ), 1_${ik}$, z( iloz, i ), 1_${ik}$, cs, sn )
              end if
           end if
           ! reset deflation counter
           kdefl = 0_${ik}$
           ! return to start of the main loop with new value of i.
           i = l - 1_${ik}$
           go to 20
           160 continue
           return
     end subroutine stdlib${ii}$_dlahqr

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$lahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,iloz, ihiz, z, ldz, &
     !! DLAHQR: is an auxiliary routine called by DHSEQR to update the
     !! eigenvalues and Schur decomposition already computed by DHSEQR, by
     !! dealing with the Hessenberg submatrix in rows and columns ILO to
     !! IHI.
               info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n
           integer(${ik}$), intent(out) :: info
           logical(lk), intent(in) :: wantt, wantz
           ! Array Arguments 
           real(${rk}$), intent(inout) :: h(ldh,*), z(ldz,*)
           real(${rk}$), intent(out) :: wi(*), wr(*)
        ! =========================================================
           ! Parameters 
           real(${rk}$), parameter :: dat1 = 3.0_${rk}$/4.0_${rk}$
           real(${rk}$), parameter :: dat2 = -0.4375_${rk}$
           integer(${ik}$), parameter :: kexsh = 10_${ik}$
           
           
           
           ! Local Scalars 
           real(${rk}$) :: aa, ab, ba, bb, cs, det, h11, h12, h21, h21s, h22, rt1i, rt1r, rt2i, rt2r, &
                     rtdisc, s, safmax, safmin, smlnum, sn, sum, t1, t2, t3, tr, tst, ulp, v2, v3
           integer(${ik}$) :: i, i1, i2, its, itmax, j, k, l, m, nh, nr, nz, kdefl
           ! Local Arrays 
           real(${rk}$) :: v(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n==0 )return
           if( ilo==ihi ) then
              wr( ilo ) = h( ilo, ilo )
              wi( ilo ) = zero
              return
           end if
           ! ==== clear out the trash ====
           do j = ilo, ihi - 3
              h( j+2, j ) = zero
              h( j+3, j ) = zero
           end do
           if( ilo<=ihi-2 )h( ihi, ihi-2 ) = zero
           nh = ihi - ilo + 1_${ik}$
           nz = ihiz - iloz + 1_${ik}$
           ! set machine-dependent constants for the stopping criterion.
           safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
           safmax = one / safmin
           call stdlib${ii}$_${ri}$labad( safmin, safmax )
           ulp = stdlib${ii}$_${ri}$lamch( 'PRECISION' )
           smlnum = safmin*( real( nh,KIND=${rk}$) / 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.
           if( wantt ) then
              i1 = 1_${ik}$
              i2 = n
           end if
           ! itmax is the total number of qr iterations allowed.
           itmax = 30_${ik}$ * max( 10_${ik}$, nh )
           ! kdefl counts the number of iterations since a deflation
           kdefl = 0_${ik}$
           ! 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<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.
           loop_140: do its = 0, itmax
              ! look for a single small subdiagonal element.
              do k = i, l + 1, -1
                 if( abs( h( k, k-1 ) )<=smlnum )go to 40
                 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
                 if( tst==zero ) then
                    if( k-2>=ilo )tst = tst + abs( h( k-1, k-2 ) )
                    if( k+1<=ihi )tst = tst + abs( h( k+1, k ) )
                 end if
                 ! ==== the following is a conservative small subdiagonal
                 ! .    deflation  criterion due to ahues
                 ! .    1997). it has better mathematical foundation and
                 ! .    improves accuracy in some cases.  ====
                 if( abs( h( k, k-1 ) )<=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 )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 40
                 end if
              end do
              40 continue
              l = k
              if( l>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>=i-1 )go to 150
              kdefl = kdefl + 1_${ik}$
              ! now the active submatrix is in rows and columns l to i. if
              ! eigenvalues only are being computed, only the active submatrix
              ! need be transformed.
              if( .not.wantt ) then
                 i1 = l
                 i2 = i
              end if
              if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
                 h11 = dat1*s + h( i, i )
                 h12 = dat2*s
                 h21 = s
                 h22 = h11
              else if( mod(kdefl,kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
                 h11 = dat1*s + h( l, l )
                 h12 = dat2*s
                 h21 = s
                 h22 = h11
              else
                 ! prepare to use francis' double shift
                 ! (i.e. 2nd degree generalized rayleigh quotient)
                 h11 = h( i-1, i-1 )
                 h21 = h( i, i-1 )
                 h12 = h( i-1, i )
                 h22 = h( i, i )
              end if
              s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
              if( s==zero ) then
                 rt1r = zero
                 rt1i = zero
                 rt2r = zero
                 rt2i = zero
              else
                 h11 = h11 / s
                 h21 = h21 / s
                 h12 = h12 / s
                 h22 = h22 / s
                 tr = ( h11+h22 ) / two
                 det = ( h11-tr )*( h22-tr ) - h12*h21
                 rtdisc = sqrt( abs( det ) )
                 if( det>=zero ) then
                    ! ==== complex conjugate shifts ====
                    rt1r = tr*s
                    rt2r = rt1r
                    rt1i = rtdisc*s
                    rt2i = -rt1i
                 else
                    ! ==== realshifts (use only one of them,KIND=${rk}$)  ====
                    rt1r = tr + rtdisc
                    rt2r = tr - rtdisc
                    if( abs( rt1r-h22 )<=abs( rt2r-h22 ) ) then
                       rt1r = rt1r*s
                       rt2r = rt1r
                    else
                       rt2r = rt2r*s
                       rt1r = rt2r
                    end if
                    rt1i = zero
                    rt2i = zero
                 end if
              end if
              ! look for two consecutive small subdiagonal elements.
              do 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 )
                 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
                 h21s = h( m+1, m ) / s
                 v( 1_${ik}$ ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*( ( h( m, m )-rt2r ) / s ) - &
                           rt1i*( rt2i / s )
                 v( 2_${ik}$ ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
                 v( 3_${ik}$ ) = h21s*h( m+2, m+1 )
                 s = abs( v( 1_${ik}$ ) ) + abs( v( 2_${ik}$ ) ) + abs( v( 3_${ik}$ ) )
                 v( 1_${ik}$ ) = v( 1_${ik}$ ) / s
                 v( 2_${ik}$ ) = v( 2_${ik}$ ) / s
                 v( 3_${ik}$ ) = v( 3_${ik}$ ) / s
                 if( m==l )go to 60
                 if( abs( h( m, m-1 ) )*( abs( v( 2_${ik}$ ) )+abs( v( 3_${ik}$ ) ) )<=ulp*abs( v( 1_${ik}$ ) )*( abs( &
                           h( m-1, m-1 ) )+abs( h( m,m ) )+abs( h( m+1, m+1 ) ) ) )go to 60
              end do
              60 continue
              ! double-shift qr step
              loop_130: do k = m, i - 1
                 ! 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. nr is the order of g.
                 nr = min( 3_${ik}$, i-k+1 )
                 if( k>m )call stdlib${ii}$_${ri}$copy( nr, h( k, k-1 ), 1_${ik}$, v, 1_${ik}$ )
                 call stdlib${ii}$_${ri}$larfg( nr, v( 1_${ik}$ ), v( 2_${ik}$ ), 1_${ik}$, t1 )
                 if( k>m ) then
                    h( k, k-1 ) = v( 1_${ik}$ )
                    h( k+1, k-1 ) = zero
                    if( k<i-1 )h( k+2, k-1 ) = zero
                 else if( m>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-t1 )
                 end if
                 v2 = v( 2_${ik}$ )
                 t2 = t1*v2
                 if( nr==3_${ik}$ ) then
                    v3 = v( 3_${ik}$ )
                    t3 = t1*v3
                    ! apply g from the left to transform the rows of the matrix
                    ! in columns k to i2.
                    do j = k, i2
                       sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
                       h( k, j ) = h( k, j ) - sum*t1
                       h( k+1, j ) = h( k+1, j ) - sum*t2
                       h( k+2, j ) = h( k+2, j ) - sum*t3
                    end do
                    ! apply g from the right to transform the columns of the
                    ! matrix in rows i1 to min(k+3,i).
                    do j = i1, min( k+3, i )
                       sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
                       h( j, k ) = h( j, k ) - sum*t1
                       h( j, k+1 ) = h( j, k+1 ) - sum*t2
                       h( j, k+2 ) = h( j, k+2 ) - sum*t3
                    end do
                    if( wantz ) then
                       ! accumulate transformations in the matrix z
                       do j = iloz, ihiz
                          sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
                          z( j, k ) = z( j, k ) - sum*t1
                          z( j, k+1 ) = z( j, k+1 ) - sum*t2
                          z( j, k+2 ) = z( j, k+2 ) - sum*t3
                       end do
                    end if
                 else if( nr==2_${ik}$ ) then
                    ! apply g from the left to transform the rows of the matrix
                    ! in columns k to i2.
                    do j = k, i2
                       sum = h( k, j ) + v2*h( k+1, j )
                       h( k, j ) = h( k, j ) - sum*t1
                       h( k+1, j ) = h( k+1, j ) - sum*t2
                    end do
                    ! apply g from the right to transform the columns of the
                    ! matrix in rows i1 to min(k+3,i).
                    do j = i1, i
                       sum = h( j, k ) + v2*h( j, k+1 )
                       h( j, k ) = h( j, k ) - sum*t1
                       h( j, k+1 ) = h( j, k+1 ) - sum*t2
                    end do
                    if( wantz ) then
                       ! accumulate transformations in the matrix z
                       do j = iloz, ihiz
                          sum = z( j, k ) + v2*z( j, k+1 )
                          z( j, k ) = z( j, k ) - sum*t1
                          z( j, k+1 ) = z( j, k+1 ) - sum*t2
                       end do
                    end if
                 end if
              end do loop_130
           end do loop_140
           ! failure to converge in remaining number of iterations
           info = i
           return
           150 continue
           if( l==i ) then
              ! h(i,i-1) is negligible: one eigenvalue has converged.
              wr( i ) = h( i, i )
              wi( i ) = zero
           else if( l==i-1 ) then
              ! h(i-1,i-2) is negligible: a pair of eigenvalues have converged.
              ! transform the 2-by-2 submatrix to standard schur form,
              ! and compute and store the eigenvalues.
              call stdlib${ii}$_${ri}$lanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),h( i, i ), wr( i-1 ), &
                        wi( i-1 ), wr( i ), wi( i ),cs, sn )
              if( wantt ) then
                 ! apply the transformation to the rest of h.
                 if( i2>i )call stdlib${ii}$_${ri}$rot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,cs, sn )
                           
                 call stdlib${ii}$_${ri}$rot( i-i1-1, h( i1, i-1 ), 1_${ik}$, h( i1, i ), 1_${ik}$, cs, sn )
              end if
              if( wantz ) then
                 ! apply the transformation to z.
                 call stdlib${ii}$_${ri}$rot( nz, z( iloz, i-1 ), 1_${ik}$, z( iloz, i ), 1_${ik}$, cs, sn )
              end if
           end if
           ! reset deflation counter
           kdefl = 0_${ik}$
           ! return to start of the main loop with new value of i.
           i = l - 1_${ik}$
           go to 20
           160 continue
           return
     end subroutine stdlib${ii}$_${ri}$lahqr

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,ihiz, z, ldz, info &
     !! CLAHQR is an auxiliary routine called by CHSEQR to update the
     !! eigenvalues and Schur decomposition already computed by CHSEQR, by
     !! dealing with the Hessenberg submatrix in rows and columns ILO to
     !! IHI.
               )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n
           integer(${ik}$), intent(out) :: info
           logical(lk), intent(in) :: wantt, wantz
           ! Array Arguments 
           complex(sp), intent(inout) :: h(ldh,*), z(ldz,*)
           complex(sp), intent(out) :: w(*)
        ! =========================================================
           ! Parameters 
           real(sp), parameter :: dat1 = 3.0_sp/4.0_sp
           integer(${ik}$), parameter :: kexsh = 10_${ik}$
           
           
           
           
           ! Local Scalars 
           complex(sp) :: cdum, h11, h11s, h22, sc, sum, t, t1, temp, u, v2, x, y
           real(sp) :: aa, ab, ba, bb, h10, h21, rtemp, s, safmax, safmin, smlnum, sx, t2, tst, &
                     ulp
           integer(${ik}$) :: i, i1, i2, its, itmax, j, jhi, jlo, k, l, m, nh, nz, kdefl
           ! Local Arrays 
           complex(sp) :: v(2_${ik}$)
           ! Statement Functions 
           real(sp) :: cabs1
           ! Intrinsic Functions 
           ! Statement Function Definitions 
           cabs1( cdum ) = abs( real( cdum,KIND=sp) ) + abs( aimag( cdum ) )
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n==0 )return
           if( ilo==ihi ) then
              w( ilo ) = h( ilo, ilo )
              return
           end if
           ! ==== clear out the trash ====
           do j = ilo, ihi - 3
              h( j+2, j ) = czero
              h( j+3, j ) = czero
           end do
           if( ilo<=ihi-2 )h( ihi, ihi-2 ) = czero
           ! ==== ensure that subdiagonal entries are real ====
           if( wantt ) then
              jlo = 1_${ik}$
              jhi = n
           else
              jlo = ilo
              jhi = ihi
           end if
           do i = ilo + 1, ihi
              if( aimag( h( i, i-1 ) )/=zero ) then
                 ! ==== the following redundant normalization
                 ! .    avoids problems with both gradual and
                 ! .    sudden underflow in abs(h(i,i-1)) ====
                 sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
                 sc = conjg( sc ) / abs( sc )
                 h( i, i-1 ) = abs( h( i, i-1 ) )
                 call stdlib${ii}$_cscal( jhi-i+1, sc, h( i, i ), ldh )
                 call stdlib${ii}$_cscal( min( jhi, i+1 )-jlo+1, conjg( sc ), h( jlo, i ),1_${ik}$ )
                 if( wantz )call stdlib${ii}$_cscal( ihiz-iloz+1, conjg( sc ), z( iloz, i ), 1_${ik}$ )
              end if
           end do
           nh = ihi - ilo + 1_${ik}$
           nz = ihiz - iloz + 1_${ik}$
           ! set machine-dependent constants for the stopping criterion.
           safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           safmax = one / safmin
           call stdlib${ii}$_slabad( safmin, safmax )
           ulp = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = safmin*( real( nh,KIND=sp) / 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.
           if( wantt ) then
              i1 = 1_${ik}$
              i2 = n
           end if
           ! itmax is the total number of qr iterations allowed.
           itmax = 30_${ik}$ * max( 10_${ik}$, nh )
           ! kdefl counts the number of iterations since a deflation
           kdefl = 0_${ik}$
           ! the main loop begins here. i is the loop index and decreases from
           ! ihi to ilo in steps of 1. 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
           30 continue
           if( i<ilo )go to 150
           ! perform qr iterations on rows and columns ilo to i until a
           ! submatrix of order 1 splits off at the bottom because a
           ! subdiagonal element has become negligible.
           l = ilo
           loop_130: do its = 0, itmax
              ! look for a single small subdiagonal element.
              do k = i, l + 1, -1
                 if( cabs1( h( k, k-1 ) )<=smlnum )go to 50
                 tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
                 if( tst==czero ) then
                    if( k-2>=ilo )tst = tst + abs( real( h( k-1, k-2 ),KIND=sp) )
                    if( k+1<=ihi )tst = tst + abs( real( h( k+1, k ),KIND=sp) )
                 end if
                 ! ==== the following is a conservative small subdiagonal
                 ! .    deflation criterion due to ahues
                 ! .    1997). it has better mathematical foundation and
                 ! .    improves accuracy in some examples.  ====
                 if( abs( real( h( k, k-1 ),KIND=sp) )<=ulp*tst ) then
                    ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
                    ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
                    aa = max( cabs1( h( k, k ) ),cabs1( h( k-1, k-1 )-h( k, k ) ) )
                    bb = min( cabs1( h( k, k ) ),cabs1( h( k-1, k-1 )-h( k, k ) ) )
                    s = aa + ab
                    if( ba*( ab / s )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 50
                 end if
              end do
              50 continue
              l = k
              if( l>ilo ) then
                 ! h(l,l-1) is negligible
                 h( l, l-1 ) = czero
              end if
              ! exit from loop if a submatrix of order 1 has split off.
              if( l>=i )go to 140
              kdefl = kdefl + 1_${ik}$
              ! now the active submatrix is in rows and columns l to i. if
              ! eigenvalues only are being computed, only the active submatrix
              ! need be transformed.
              if( .not.wantt ) then
                 i1 = l
                 i2 = i
              end if
              if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = dat1*abs( real( h( i, i-1 ),KIND=sp) )
                 t = s + h( i, i )
              else if( mod(kdefl,kexsh)==0_${ik}$ ) then
                 ! exceptional shift.
                 s = dat1*abs( real( h( l+1, l ),KIND=sp) )
                 t = s + h( l, l )
              else
                 ! wilkinson's shift.
                 t = h( i, i )
                 u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
                 s = cabs1( u )
                 if( s/=zero ) then
                    x = half*( h( i-1, i-1 )-t )
                    sx = cabs1( x )
                    s = max( s, cabs1( x ) )
                    y = s*sqrt( ( x / s )**2_${ik}$+( u / s )**2_${ik}$ )
                    if( sx>zero ) then
                       if( real( x / sx,KIND=sp)*real( y,KIND=sp)+aimag( x / sx )*aimag( y )&
                                 <zero )y = -y
                    end if
                    t = t - u*stdlib${ii}$_cladiv( u, ( x+y ) )
                 end if
              end if
              ! look for two consecutive small subdiagonal elements.
              do m = i - 1, l + 1, -1
                 ! determine the effect of starting the single-shift qr
                 ! iteration at row m, and see if this would make h(m,m-1)
                 ! negligible.
                 h11 = h( m, m )
                 h22 = h( m+1, m+1 )
                 h11s = h11 - t
                 h21 = real( h( m+1, m ),KIND=sp)
                 s = cabs1( h11s ) + abs( h21 )
                 h11s = h11s / s
                 h21