stdlib_lapack_eigv_gen2.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     module subroutine stdlib${ii}$_shseqr( job, compz, n, ilo, ihi, h, ldh, wr, wi, z,ldz, work, lwork, info )
     !! SHSEQR computes the eigenvalues of a Hessenberg matrix H
     !! and, optionally, the matrices T and Z from the Schur decomposition
     !! H = Z T Z**T, where T is an upper quasi-triangular matrix (the
     !! Schur form), and Z is the orthogonal matrix of Schur vectors.
     !! Optionally Z may be postmultiplied into an input orthogonal
     !! matrix Q so that this routine can give the Schur factorization
     !! of a matrix A which has been reduced to the Hessenberg form H
     !! by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
               
        ! -- lapack computational 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, ilo, ldh, ldz, lwork, n
           integer(${ik}$), intent(out) :: info
           character, intent(in) :: compz, job
           ! Array Arguments 
           real(sp), intent(inout) :: h(ldh,*), z(ldz,*)
           real(sp), intent(out) :: wi(*), work(*), wr(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: ntiny = 15_${ik}$
           integer(${ik}$), parameter :: nl = 49_${ik}$
           ! ==== matrices of order ntiny or smaller must be processed by
           ! .    stdlib${ii}$_slahqr because of insufficient subdiagonal scratch space.
           ! .    (this is a hard limit.) ====
           
           ! ==== nl allocates some local workspace to help small matrices
           ! .    through a rare stdlib${ii}$_slahqr failure.  nl > ntiny = 15 is
           ! .    required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom-
           ! .    mended.  (the default value of nmin is 75.)  using nl = 49
           ! .    allows up to six simultaneous shifts and a 16-by-16
           ! .    deflation window.  ====
           
           
           ! Local Arrays 
           real(sp) :: hl(nl,nl), workl(nl)
           ! Local Scalars 
           integer(${ik}$) :: i, kbot, nmin
           logical(lk) :: initz, lquery, wantt, wantz
           ! Intrinsic Functions 
           ! Executable Statements 
           ! ==== decode and check the input parameters. ====
           wantt = stdlib_lsame( job, 'S' )
           initz = stdlib_lsame( compz, 'I' )
           wantz = initz .or. stdlib_lsame( compz, 'V' )
           work( 1_${ik}$ ) = real( max( 1_${ik}$, n ),KIND=sp)
           lquery = lwork==-1_${ik}$
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ihi<min( ilo, n ) .or. ihi>n ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then
              info = -11_${ik}$
           else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              ! ==== quick return in case of invalid argument. ====
              call stdlib${ii}$_xerbla( 'SHSEQR', -info )
              return
           else if( n==0_${ik}$ ) then
              ! ==== quick return in case n = 0; nothing to do. ====
              return
           else if( lquery ) then
              ! ==== quick return in case of a workspace query ====
              call stdlib${ii}$_slaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                        work, lwork, info )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=sp), work( 1_${ik}$ ) )
              return
           else
              ! ==== copy eigenvalues isolated by stdlib${ii}$_sgebal ====
              do i = 1, ilo - 1
                 wr( i ) = h( i, i )
                 wi( i ) = zero
              end do
              do i = ihi + 1, n
                 wr( i ) = h( i, i )
                 wi( i ) = zero
              end do
              ! ==== initialize z, if requested ====
              if( initz )call stdlib${ii}$_slaset( 'A', n, n, zero, one, z, ldz )
              ! ==== quick return if possible ====
              if( ilo==ihi ) then
                 wr( ilo ) = h( ilo, ilo )
                 wi( ilo ) = zero
                 return
              end if
              ! ==== stdlib${ii}$_slahqr/stdlib${ii}$_slaqr0 crossover point ====
              nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'SHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork )
                        
              nmin = max( ntiny, nmin )
              ! ==== stdlib${ii}$_slaqr0 for big matrices; stdlib${ii}$_slahqr for small ones ====
              if( n>nmin ) then
                 call stdlib${ii}$_slaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                           work, lwork, info )
              else
                 ! ==== small matrix ====
                 call stdlib${ii}$_slahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                           info )
                 if( info>0_${ik}$ ) then
                    ! ==== a rare stdlib${ii}$_slahqr failure!  stdlib${ii}$_slaqr0 sometimes succeeds
                    ! .    when stdlib${ii}$_slahqr fails. ====
                    kbot = info
                    if( n>=nl ) then
                       ! ==== larger matrices have enough subdiagonal scratch
                       ! .    space to call stdlib${ii}$_slaqr0 directly. ====
                       call stdlib${ii}$_slaqr0( wantt, wantz, n, ilo, kbot, h, ldh, wr,wi, ilo, ihi, z,&
                                  ldz, work, lwork, info )
                    else
                       ! ==== tiny matrices don't have enough subdiagonal
                       ! .    scratch space to benefit from stdlib${ii}$_slaqr0.  hence,
                       ! .    tiny matrices must be copied into a larger
                       ! .    array before calling stdlib${ii}$_slaqr0. ====
                       call stdlib${ii}$_slacpy( 'A', n, n, h, ldh, hl, nl )
                       hl( n+1, n ) = zero
                       call stdlib${ii}$_slaset( 'A', nl, nl-n, zero, zero, hl( 1_${ik}$, n+1 ),nl )
                       call stdlib${ii}$_slaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, wr,wi, ilo, ihi, &
                                 z, ldz, workl, nl, info )
                       if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_slacpy( 'A', n, n, hl, nl, h, ldh )
                                 
                    end if
                 end if
              end if
              ! ==== clear out the trash, if necessary. ====
              if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_slaset( 'L', n-2, n-2, zero, zero,&
                         h( 3_${ik}$, 1_${ik}$ ), ldh )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=sp), work( 1_${ik}$ ) )
           end if
     end subroutine stdlib${ii}$_shseqr

     module subroutine stdlib${ii}$_dhseqr( job, compz, n, ilo, ihi, h, ldh, wr, wi, z,ldz, work, lwork, info )
     !! DHSEQR computes the eigenvalues of a Hessenberg matrix H
     !! and, optionally, the matrices T and Z from the Schur decomposition
     !! H = Z T Z**T, where T is an upper quasi-triangular matrix (the
     !! Schur form), and Z is the orthogonal matrix of Schur vectors.
     !! Optionally Z may be postmultiplied into an input orthogonal
     !! matrix Q so that this routine can give the Schur factorization
     !! of a matrix A which has been reduced to the Hessenberg form H
     !! by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
               
        ! -- lapack computational 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, ilo, ldh, ldz, lwork, n
           integer(${ik}$), intent(out) :: info
           character, intent(in) :: compz, job
           ! Array Arguments 
           real(dp), intent(inout) :: h(ldh,*), z(ldz,*)
           real(dp), intent(out) :: wi(*), work(*), wr(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: ntiny = 15_${ik}$
           integer(${ik}$), parameter :: nl = 49_${ik}$
           ! ==== matrices of order ntiny or smaller must be processed by
           ! .    stdlib${ii}$_dlahqr because of insufficient subdiagonal scratch space.
           ! .    (this is a hard limit.) ====
           
           ! ==== nl allocates some local workspace to help small matrices
           ! .    through a rare stdlib${ii}$_dlahqr failure.  nl > ntiny = 15 is
           ! .    required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom-
           ! .    mended.  (the default value of nmin is 75.)  using nl = 49
           ! .    allows up to six simultaneous shifts and a 16-by-16
           ! .    deflation window.  ====
           
           
           ! Local Arrays 
           real(dp) :: hl(nl,nl), workl(nl)
           ! Local Scalars 
           integer(${ik}$) :: i, kbot, nmin
           logical(lk) :: initz, lquery, wantt, wantz
           ! Intrinsic Functions 
           ! Executable Statements 
           ! ==== decode and check the input parameters. ====
           wantt = stdlib_lsame( job, 'S' )
           initz = stdlib_lsame( compz, 'I' )
           wantz = initz .or. stdlib_lsame( compz, 'V' )
           work( 1_${ik}$ ) = real( max( 1_${ik}$, n ),KIND=dp)
           lquery = lwork==-1_${ik}$
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ihi<min( ilo, n ) .or. ihi>n ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then
              info = -11_${ik}$
           else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              ! ==== quick return in case of invalid argument. ====
              call stdlib${ii}$_xerbla( 'DHSEQR', -info )
              return
           else if( n==0_${ik}$ ) then
              ! ==== quick return in case n = 0; nothing to do. ====
              return
           else if( lquery ) then
              ! ==== quick return in case of a workspace query ====
              call stdlib${ii}$_dlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                        work, lwork, info )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=dp), work( 1_${ik}$ ) )
              return
           else
              ! ==== copy eigenvalues isolated by stdlib${ii}$_dgebal ====
              do i = 1, ilo - 1
                 wr( i ) = h( i, i )
                 wi( i ) = zero
              end do
              do i = ihi + 1, n
                 wr( i ) = h( i, i )
                 wi( i ) = zero
              end do
              ! ==== initialize z, if requested ====
              if( initz )call stdlib${ii}$_dlaset( 'A', n, n, zero, one, z, ldz )
              ! ==== quick return if possible ====
              if( ilo==ihi ) then
                 wr( ilo ) = h( ilo, ilo )
                 wi( ilo ) = zero
                 return
              end if
              ! ==== stdlib${ii}$_dlahqr/stdlib${ii}$_dlaqr0 crossover point ====
              nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'DHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork )
                        
              nmin = max( ntiny, nmin )
              ! ==== stdlib${ii}$_dlaqr0 for big matrices; stdlib${ii}$_dlahqr for small ones ====
              if( n>nmin ) then
                 call stdlib${ii}$_dlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                           work, lwork, info )
              else
                 ! ==== small matrix ====
                 call stdlib${ii}$_dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                           info )
                 if( info>0_${ik}$ ) then
                    ! ==== a rare stdlib${ii}$_dlahqr failure!  stdlib${ii}$_dlaqr0 sometimes succeeds
                    ! .    when stdlib${ii}$_dlahqr fails. ====
                    kbot = info
                    if( n>=nl ) then
                       ! ==== larger matrices have enough subdiagonal scratch
                       ! .    space to call stdlib${ii}$_dlaqr0 directly. ====
                       call stdlib${ii}$_dlaqr0( wantt, wantz, n, ilo, kbot, h, ldh, wr,wi, ilo, ihi, z,&
                                  ldz, work, lwork, info )
                    else
                       ! ==== tiny matrices don't have enough subdiagonal
                       ! .    scratch space to benefit from stdlib${ii}$_dlaqr0.  hence,
                       ! .    tiny matrices must be copied into a larger
                       ! .    array before calling stdlib${ii}$_dlaqr0. ====
                       call stdlib${ii}$_dlacpy( 'A', n, n, h, ldh, hl, nl )
                       hl( n+1, n ) = zero
                       call stdlib${ii}$_dlaset( 'A', nl, nl-n, zero, zero, hl( 1_${ik}$, n+1 ),nl )
                       call stdlib${ii}$_dlaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, wr,wi, ilo, ihi, &
                                 z, ldz, workl, nl, info )
                       if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_dlacpy( 'A', n, n, hl, nl, h, ldh )
                                 
                    end if
                 end if
              end if
              ! ==== clear out the trash, if necessary. ====
              if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_dlaset( 'L', n-2, n-2, zero, zero,&
                         h( 3_${ik}$, 1_${ik}$ ), ldh )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=dp), work( 1_${ik}$ ) )
           end if
     end subroutine stdlib${ii}$_dhseqr

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$hseqr( job, compz, n, ilo, ihi, h, ldh, wr, wi, z,ldz, work, lwork, info )
     !! DHSEQR: computes the eigenvalues of a Hessenberg matrix H
     !! and, optionally, the matrices T and Z from the Schur decomposition
     !! H = Z T Z**T, where T is an upper quasi-triangular matrix (the
     !! Schur form), and Z is the orthogonal matrix of Schur vectors.
     !! Optionally Z may be postmultiplied into an input orthogonal
     !! matrix Q so that this routine can give the Schur factorization
     !! of a matrix A which has been reduced to the Hessenberg form H
     !! by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
               
        ! -- lapack computational 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, ilo, ldh, ldz, lwork, n
           integer(${ik}$), intent(out) :: info
           character, intent(in) :: compz, job
           ! Array Arguments 
           real(${rk}$), intent(inout) :: h(ldh,*), z(ldz,*)
           real(${rk}$), intent(out) :: wi(*), work(*), wr(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: ntiny = 15_${ik}$
           integer(${ik}$), parameter :: nl = 49_${ik}$
           ! ==== matrices of order ntiny or smaller must be processed by
           ! .    stdlib${ii}$_${ri}$lahqr because of insufficient subdiagonal scratch space.
           ! .    (this is a hard limit.) ====
           
           ! ==== nl allocates some local workspace to help small matrices
           ! .    through a rare stdlib${ii}$_${ri}$lahqr failure.  nl > ntiny = 15 is
           ! .    required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom-
           ! .    mended.  (the default value of nmin is 75.)  using nl = 49
           ! .    allows up to six simultaneous shifts and a 16-by-16
           ! .    deflation window.  ====
           
           
           ! Local Arrays 
           real(${rk}$) :: hl(nl,nl), workl(nl)
           ! Local Scalars 
           integer(${ik}$) :: i, kbot, nmin
           logical(lk) :: initz, lquery, wantt, wantz
           ! Intrinsic Functions 
           ! Executable Statements 
           ! ==== decode and check the input parameters. ====
           wantt = stdlib_lsame( job, 'S' )
           initz = stdlib_lsame( compz, 'I' )
           wantz = initz .or. stdlib_lsame( compz, 'V' )
           work( 1_${ik}$ ) = real( max( 1_${ik}$, n ),KIND=${rk}$)
           lquery = lwork==-1_${ik}$
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ihi<min( ilo, n ) .or. ihi>n ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then
              info = -11_${ik}$
           else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              ! ==== quick return in case of invalid argument. ====
              call stdlib${ii}$_xerbla( 'DHSEQR', -info )
              return
           else if( n==0_${ik}$ ) then
              ! ==== quick return in case n = 0; nothing to do. ====
              return
           else if( lquery ) then
              ! ==== quick return in case of a workspace query ====
              call stdlib${ii}$_${ri}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                        work, lwork, info )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=${rk}$), work( 1_${ik}$ ) )
              return
           else
              ! ==== copy eigenvalues isolated by stdlib${ii}$_${ri}$gebal ====
              do i = 1, ilo - 1
                 wr( i ) = h( i, i )
                 wi( i ) = zero
              end do
              do i = ihi + 1, n
                 wr( i ) = h( i, i )
                 wi( i ) = zero
              end do
              ! ==== initialize z, if requested ====
              if( initz )call stdlib${ii}$_${ri}$laset( 'A', n, n, zero, one, z, ldz )
              ! ==== quick return if possible ====
              if( ilo==ihi ) then
                 wr( ilo ) = h( ilo, ilo )
                 wi( ilo ) = zero
                 return
              end if
              ! ==== stdlib${ii}$_${ri}$lahqr/stdlib${ii}$_${ri}$laqr0 crossover point ====
              nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'DHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork )
                        
              nmin = max( ntiny, nmin )
              ! ==== stdlib${ii}$_${ri}$laqr0 for big matrices; stdlib${ii}$_${ri}$lahqr for small ones ====
              if( n>nmin ) then
                 call stdlib${ii}$_${ri}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                           work, lwork, info )
              else
                 ! ==== small matrix ====
                 call stdlib${ii}$_${ri}$lahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, &
                           info )
                 if( info>0_${ik}$ ) then
                    ! ==== a rare stdlib${ii}$_${ri}$lahqr failure!  stdlib${ii}$_${ri}$laqr0 sometimes succeeds
                    ! .    when stdlib${ii}$_${ri}$lahqr fails. ====
                    kbot = info
                    if( n>=nl ) then
                       ! ==== larger matrices have enough subdiagonal scratch
                       ! .    space to call stdlib${ii}$_${ri}$laqr0 directly. ====
                       call stdlib${ii}$_${ri}$laqr0( wantt, wantz, n, ilo, kbot, h, ldh, wr,wi, ilo, ihi, z,&
                                  ldz, work, lwork, info )
                    else
                       ! ==== tiny matrices don't have enough subdiagonal
                       ! .    scratch space to benefit from stdlib${ii}$_${ri}$laqr0.  hence,
                       ! .    tiny matrices must be copied into a larger
                       ! .    array before calling stdlib${ii}$_${ri}$laqr0. ====
                       call stdlib${ii}$_${ri}$lacpy( 'A', n, n, h, ldh, hl, nl )
                       hl( n+1, n ) = zero
                       call stdlib${ii}$_${ri}$laset( 'A', nl, nl-n, zero, zero, hl( 1_${ik}$, n+1 ),nl )
                       call stdlib${ii}$_${ri}$laqr0( wantt, wantz, nl, ilo, kbot, hl, nl, wr,wi, ilo, ihi, &
                                 z, ldz, workl, nl, info )
                       if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_${ri}$lacpy( 'A', n, n, hl, nl, h, ldh )
                                 
                    end if
                 end if
              end if
              ! ==== clear out the trash, if necessary. ====
              if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_${ri}$laset( 'L', n-2, n-2, zero, zero,&
                         h( 3_${ik}$, 1_${ik}$ ), ldh )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=${rk}$), work( 1_${ik}$ ) )
           end if
     end subroutine stdlib${ii}$_${ri}$hseqr

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_chseqr( job, compz, n, ilo, ihi, h, ldh, w, z, ldz,work, lwork, info )
     !! CHSEQR computes the eigenvalues of a Hessenberg matrix H
     !! and, optionally, the matrices T and Z from the Schur decomposition
     !! H = Z T Z**H, where T is an upper triangular matrix (the
     !! Schur form), and Z is the unitary matrix of Schur vectors.
     !! Optionally Z may be postmultiplied into an input unitary
     !! matrix Q so that this routine can give the Schur factorization
     !! of a matrix A which has been reduced to the Hessenberg form H
     !! by the unitary matrix Q:  A = Q*H*Q**H = (QZ)*T*(QZ)**H.
               
        ! -- lapack computational 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, ilo, ldh, ldz, lwork, n
           integer(${ik}$), intent(out) :: info
           character, intent(in) :: compz, job
           ! Array Arguments 
           complex(sp), intent(inout) :: h(ldh,*), z(ldz,*)
           complex(sp), intent(out) :: w(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: ntiny = 15_${ik}$
           integer(${ik}$), parameter :: nl = 49_${ik}$
           ! ==== matrices of order ntiny or smaller must be processed by
           ! .    stdlib${ii}$_clahqr because of insufficient subdiagonal scratch space.
           ! .    (this is a hard limit.) ====
           
           ! ==== nl allocates some local workspace to help small matrices
           ! .    through a rare stdlib${ii}$_clahqr failure.  nl > ntiny = 15 is
           ! .    required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom-
           ! .    mended.  (the default value of nmin is 75.)  using nl = 49
           ! .    allows up to six simultaneous shifts and a 16-by-16
           ! .    deflation window.  ====
           
           
           
           ! Local Arrays 
           complex(sp) :: hl(nl,nl), workl(nl)
           ! Local Scalars 
           integer(${ik}$) :: kbot, nmin
           logical(lk) :: initz, lquery, wantt, wantz
           ! Intrinsic Functions 
           ! Executable Statements 
           ! ==== decode and check the input parameters. ====
           wantt = stdlib_lsame( job, 'S' )
           initz = stdlib_lsame( compz, 'I' )
           wantz = initz .or. stdlib_lsame( compz, 'V' )
           work( 1_${ik}$ ) = cmplx( real( max( 1_${ik}$, n ),KIND=sp), zero,KIND=sp)
           lquery = lwork==-1_${ik}$
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ihi<min( ilo, n ) .or. ihi>n ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then
              info = -10_${ik}$
           else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then
              info = -12_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              ! ==== quick return in case of invalid argument. ====
              call stdlib${ii}$_xerbla( 'CHSEQR', -info )
              return
           else if( n==0_${ik}$ ) then
              ! ==== quick return in case n = 0; nothing to do. ====
              return
           else if( lquery ) then
              ! ==== quick return in case of a workspace query ====
              call stdlib${ii}$_claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,ldz, work, &
                        lwork, info )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = cmplx( max( real( work( 1_${ik}$ ),KIND=sp), real( max( 1_${ik}$,n ),KIND=sp) ), &
                        zero,KIND=sp)
              return
           else
              ! ==== copy eigenvalues isolated by stdlib${ii}$_cgebal ====
              if( ilo>1_${ik}$ )call stdlib${ii}$_ccopy( ilo-1, h, ldh+1, w, 1_${ik}$ )
              if( ihi<n )call stdlib${ii}$_ccopy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1_${ik}$ )
                        
              ! ==== initialize z, if requested ====
              if( initz )call stdlib${ii}$_claset( 'A', n, n, czero, cone, z, ldz )
              ! ==== quick return if possible ====
              if( ilo==ihi ) then
                 w( ilo ) = h( ilo, ilo )
                 return
              end if
              ! ==== stdlib${ii}$_clahqr/stdlib${ii}$_claqr0 crossover point ====
              nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'CHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork )
                        
              nmin = max( ntiny, nmin )
              ! ==== stdlib${ii}$_claqr0 for big matrices; stdlib${ii}$_clahqr for small ones ====
              if( n>nmin ) then
                 call stdlib${ii}$_claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, work, &
                           lwork, info )
              else
                 ! ==== small matrix ====
                 call stdlib${ii}$_clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, info )
                           
                 if( info>0_${ik}$ ) then
                    ! ==== a rare stdlib${ii}$_clahqr failure!  stdlib${ii}$_claqr0 sometimes succeeds
                    ! .    when stdlib${ii}$_clahqr fails. ====
                    kbot = info
                    if( n>=nl ) then
                       ! ==== larger matrices have enough subdiagonal scratch
                       ! .    space to call stdlib${ii}$_claqr0 directly. ====
                       call stdlib${ii}$_claqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,ilo, ihi, z, ldz,&
                                  work, lwork, info )
                    else
                       ! ==== tiny matrices don't have enough subdiagonal
                       ! .    scratch space to benefit from stdlib${ii}$_claqr0.  hence,
                       ! .    tiny matrices must be copied into a larger
                       ! .    array before calling stdlib${ii}$_claqr0. ====
                       call stdlib${ii}$_clacpy( 'A', n, n, h, ldh, hl, nl )
                       hl( n+1, n ) = czero
                       call stdlib${ii}$_claset( 'A', nl, nl-n, czero, czero, hl( 1_${ik}$, n+1 ),nl )
                       call stdlib${ii}$_claqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,ilo, ihi, z, &
                                 ldz, workl, nl, info )
                       if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_clacpy( 'A', n, n, hl, nl, h, ldh )
                                 
                    end if
                 end if
              end if
              ! ==== clear out the trash, if necessary. ====
              if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_claset( 'L', n-2, n-2, czero, &
                        czero, h( 3_${ik}$, 1_${ik}$ ), ldh )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = cmplx( max( real( max( 1_${ik}$, n ),KIND=sp),real( work( 1_${ik}$ ),KIND=sp) ), &
                        zero,KIND=sp)
           end if
     end subroutine stdlib${ii}$_chseqr

     pure module subroutine stdlib${ii}$_zhseqr( job, compz, n, ilo, ihi, h, ldh, w, z, ldz,work, lwork, info )
     !! ZHSEQR computes the eigenvalues of a Hessenberg matrix H
     !! and, optionally, the matrices T and Z from the Schur decomposition
     !! H = Z T Z**H, where T is an upper triangular matrix (the
     !! Schur form), and Z is the unitary matrix of Schur vectors.
     !! Optionally Z may be postmultiplied into an input unitary
     !! matrix Q so that this routine can give the Schur factorization
     !! of a matrix A which has been reduced to the Hessenberg form H
     !! by the unitary matrix Q:  A = Q*H*Q**H = (QZ)*T*(QZ)**H.
               
        ! -- lapack computational 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, ilo, ldh, ldz, lwork, n
           integer(${ik}$), intent(out) :: info
           character, intent(in) :: compz, job
           ! Array Arguments 
           complex(dp), intent(inout) :: h(ldh,*), z(ldz,*)
           complex(dp), intent(out) :: w(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: ntiny = 15_${ik}$
           integer(${ik}$), parameter :: nl = 49_${ik}$
           ! ==== matrices of order ntiny or smaller must be processed by
           ! .    stdlib${ii}$_zlahqr because of insufficient subdiagonal scratch space.
           ! .    (this is a hard limit.) ====
           
           ! ==== nl allocates some local workspace to help small matrices
           ! .    through a rare stdlib${ii}$_zlahqr failure.  nl > ntiny = 15 is
           ! .    required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom-
           ! .    mended.  (the default value of nmin is 75.)  using nl = 49
           ! .    allows up to six simultaneous shifts and a 16-by-16
           ! .    deflation window.  ====
           
           
           
           ! Local Arrays 
           complex(dp) :: hl(nl,nl), workl(nl)
           ! Local Scalars 
           integer(${ik}$) :: kbot, nmin
           logical(lk) :: initz, lquery, wantt, wantz
           ! Intrinsic Functions 
           ! Executable Statements 
           ! ==== decode and check the input parameters. ====
           wantt = stdlib_lsame( job, 'S' )
           initz = stdlib_lsame( compz, 'I' )
           wantz = initz .or. stdlib_lsame( compz, 'V' )
           work( 1_${ik}$ ) = cmplx( real( max( 1_${ik}$, n ),KIND=dp), zero,KIND=dp)
           lquery = lwork==-1_${ik}$
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ihi<min( ilo, n ) .or. ihi>n ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then
              info = -10_${ik}$
           else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then
              info = -12_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              ! ==== quick return in case of invalid argument. ====
              call stdlib${ii}$_xerbla( 'ZHSEQR', -info )
              return
           else if( n==0_${ik}$ ) then
              ! ==== quick return in case n = 0; nothing to do. ====
              return
           else if( lquery ) then
              ! ==== quick return in case of a workspace query ====
              call stdlib${ii}$_zlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,ldz, work, &
                        lwork, info )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = cmplx( max( real( work( 1_${ik}$ ),KIND=dp), real( max( 1_${ik}$,n ),KIND=dp) ), &
                        zero,KIND=dp)
              return
           else
              ! ==== copy eigenvalues isolated by stdlib${ii}$_zgebal ====
              if( ilo>1_${ik}$ )call stdlib${ii}$_zcopy( ilo-1, h, ldh+1, w, 1_${ik}$ )
              if( ihi<n )call stdlib${ii}$_zcopy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1_${ik}$ )
                        
              ! ==== initialize z, if requested ====
              if( initz )call stdlib${ii}$_zlaset( 'A', n, n, czero, cone, z, ldz )
              ! ==== quick return if possible ====
              if( ilo==ihi ) then
                 w( ilo ) = h( ilo, ilo )
                 return
              end if
              ! ==== stdlib${ii}$_zlahqr/stdlib${ii}$_zlaqr0 crossover point ====
              nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'ZHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork )
                        
              nmin = max( ntiny, nmin )
              ! ==== stdlib${ii}$_zlaqr0 for big matrices; stdlib${ii}$_zlahqr for small ones ====
              if( n>nmin ) then
                 call stdlib${ii}$_zlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, work, &
                           lwork, info )
              else
                 ! ==== small matrix ====
                 call stdlib${ii}$_zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, info )
                           
                 if( info>0_${ik}$ ) then
                    ! ==== a rare stdlib${ii}$_zlahqr failure!  stdlib${ii}$_zlaqr0 sometimes succeeds
                    ! .    when stdlib${ii}$_zlahqr fails. ====
                    kbot = info
                    if( n>=nl ) then
                       ! ==== larger matrices have enough subdiagonal scratch
                       ! .    space to call stdlib${ii}$_zlaqr0 directly. ====
                       call stdlib${ii}$_zlaqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,ilo, ihi, z, ldz,&
                                  work, lwork, info )
                    else
                       ! ==== tiny matrices don't have enough subdiagonal
                       ! .    scratch space to benefit from stdlib${ii}$_zlaqr0.  hence,
                       ! .    tiny matrices must be copied into a larger
                       ! .    array before calling stdlib${ii}$_zlaqr0. ====
                       call stdlib${ii}$_zlacpy( 'A', n, n, h, ldh, hl, nl )
                       hl( n+1, n ) = czero
                       call stdlib${ii}$_zlaset( 'A', nl, nl-n, czero, czero, hl( 1_${ik}$, n+1 ),nl )
                       call stdlib${ii}$_zlaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,ilo, ihi, z, &
                                 ldz, workl, nl, info )
                       if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_zlacpy( 'A', n, n, hl, nl, h, ldh )
                                 
                    end if
                 end if
              end if
              ! ==== clear out the trash, if necessary. ====
              if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_zlaset( 'L', n-2, n-2, czero, &
                        czero, h( 3_${ik}$, 1_${ik}$ ), ldh )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = cmplx( max( real( max( 1_${ik}$, n ),KIND=dp),real( work( 1_${ik}$ ),KIND=dp) ), &
                        zero,KIND=dp)
           end if
     end subroutine stdlib${ii}$_zhseqr

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$hseqr( job, compz, n, ilo, ihi, h, ldh, w, z, ldz,work, lwork, info )
     !! ZHSEQR: computes the eigenvalues of a Hessenberg matrix H
     !! and, optionally, the matrices T and Z from the Schur decomposition
     !! H = Z T Z**H, where T is an upper triangular matrix (the
     !! Schur form), and Z is the unitary matrix of Schur vectors.
     !! Optionally Z may be postmultiplied into an input unitary
     !! matrix Q so that this routine can give the Schur factorization
     !! of a matrix A which has been reduced to the Hessenberg form H
     !! by the unitary matrix Q:  A = Q*H*Q**H = (QZ)*T*(QZ)**H.
               
        ! -- lapack computational 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_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ihi, ilo, ldh, ldz, lwork, n
           integer(${ik}$), intent(out) :: info
           character, intent(in) :: compz, job
           ! Array Arguments 
           complex(${ck}$), intent(inout) :: h(ldh,*), z(ldz,*)
           complex(${ck}$), intent(out) :: w(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: ntiny = 15_${ik}$
           integer(${ik}$), parameter :: nl = 49_${ik}$
           ! ==== matrices of order ntiny or smaller must be processed by
           ! .    stdlib${ii}$_${ci}$lahqr because of insufficient subdiagonal scratch space.
           ! .    (this is a hard limit.) ====
           
           ! ==== nl allocates some local workspace to help small matrices
           ! .    through a rare stdlib${ii}$_${ci}$lahqr failure.  nl > ntiny = 15 is
           ! .    required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom-
           ! .    mended.  (the default value of nmin is 75.)  using nl = 49
           ! .    allows up to six simultaneous shifts and a 16-by-16
           ! .    deflation window.  ====
           
           
           
           ! Local Arrays 
           complex(${ck}$) :: hl(nl,nl), workl(nl)
           ! Local Scalars 
           integer(${ik}$) :: kbot, nmin
           logical(lk) :: initz, lquery, wantt, wantz
           ! Intrinsic Functions 
           ! Executable Statements 
           ! ==== decode and check the input parameters. ====
           wantt = stdlib_lsame( job, 'S' )
           initz = stdlib_lsame( compz, 'I' )
           wantz = initz .or. stdlib_lsame( compz, 'V' )
           work( 1_${ik}$ ) = cmplx( real( max( 1_${ik}$, n ),KIND=${ck}$), zero,KIND=${ck}$)
           lquery = lwork==-1_${ik}$
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ihi<min( ilo, n ) .or. ihi>n ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then
              info = -10_${ik}$
           else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then
              info = -12_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              ! ==== quick return in case of invalid argument. ====
              call stdlib${ii}$_xerbla( 'ZHSEQR', -info )
              return
           else if( n==0_${ik}$ ) then
              ! ==== quick return in case n = 0; nothing to do. ====
              return
           else if( lquery ) then
              ! ==== quick return in case of a workspace query ====
              call stdlib${ii}$_${ci}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,ldz, work, &
                        lwork, info )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = cmplx( max( real( work( 1_${ik}$ ),KIND=${ck}$), real( max( 1_${ik}$,n ),KIND=${ck}$) ), &
                        zero,KIND=${ck}$)
              return
           else
              ! ==== copy eigenvalues isolated by stdlib${ii}$_${ci}$gebal ====
              if( ilo>1_${ik}$ )call stdlib${ii}$_${ci}$copy( ilo-1, h, ldh+1, w, 1_${ik}$ )
              if( ihi<n )call stdlib${ii}$_${ci}$copy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1_${ik}$ )
                        
              ! ==== initialize z, if requested ====
              if( initz )call stdlib${ii}$_${ci}$laset( 'A', n, n, czero, cone, z, ldz )
              ! ==== quick return if possible ====
              if( ilo==ihi ) then
                 w( ilo ) = h( ilo, ilo )
                 return
              end if
              ! ==== stdlib${ii}$_${ci}$lahqr/stdlib${ii}$_${ci}$laqr0 crossover point ====
              nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'ZHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork )
                        
              nmin = max( ntiny, nmin )
              ! ==== stdlib${ii}$_${ci}$laqr0 for big matrices; stdlib${ii}$_${ci}$lahqr for small ones ====
              if( n>nmin ) then
                 call stdlib${ii}$_${ci}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, work, &
                           lwork, info )
              else
                 ! ==== small matrix ====
                 call stdlib${ii}$_${ci}$lahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, info )
                           
                 if( info>0_${ik}$ ) then
                    ! ==== a rare stdlib${ii}$_${ci}$lahqr failure!  stdlib${ii}$_${ci}$laqr0 sometimes succeeds
                    ! .    when stdlib${ii}$_${ci}$lahqr fails. ====
                    kbot = info
                    if( n>=nl ) then
                       ! ==== larger matrices have enough subdiagonal scratch
                       ! .    space to call stdlib${ii}$_${ci}$laqr0 directly. ====
                       call stdlib${ii}$_${ci}$laqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,ilo, ihi, z, ldz,&
                                  work, lwork, info )
                    else
                       ! ==== tiny matrices don't have enough subdiagonal
                       ! .    scratch space to benefit from stdlib${ii}$_${ci}$laqr0.  hence,
                       ! .    tiny matrices must be copied into a larger
                       ! .    array before calling stdlib${ii}$_${ci}$laqr0. ====
                       call stdlib${ii}$_${ci}$lacpy( 'A', n, n, h, ldh, hl, nl )
                       hl( n+1, n ) = czero
                       call stdlib${ii}$_${ci}$laset( 'A', nl, nl-n, czero, czero, hl( 1_${ik}$, n+1 ),nl )
                       call stdlib${ii}$_${ci}$laqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,ilo, ihi, z, &
                                 ldz, workl, nl, info )
                       if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_${ci}$lacpy( 'A', n, n, hl, nl, h, ldh )
                                 
                    end if
                 end if
              end if
              ! ==== clear out the trash, if necessary. ====
              if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_${ci}$laset( 'L', n-2, n-2, czero, &
                        czero, h( 3_${ik}$, 1_${ik}$ ), ldh )
              ! ==== ensure reported workspace size is backward-compatible with
              ! .    previous lapack versions. ====
              work( 1_${ik}$ ) = cmplx( max( real( max( 1_${ik}$, n ),KIND=${ck}$),real( work( 1_${ik}$ ),KIND=${ck}$) ), &
                        zero,KIND=${ck}$)
           end if
     end subroutine stdlib${ii}$_${ci}$hseqr

#:endif
#:endfor



     module subroutine stdlib${ii}$_shsein( side, eigsrc, initv, select, n, h, ldh, wr, wi,vl, ldvl, vr, ldvr, &
     !! SHSEIN uses inverse iteration to find specified right and/or left
     !! eigenvectors of a real upper Hessenberg matrix H.
     !! The right eigenvector x and the left eigenvector y of the matrix H
     !! corresponding to an eigenvalue w are defined by:
     !! H * x = w * x,     y**h * H = w * y**h
     !! where y**h denotes the conjugate transpose of the vector y.
               mm, m, work, ifaill,ifailr, info )
        ! -- lapack computational 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 
           character, intent(in) :: eigsrc, initv, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(inout) :: select(*)
           integer(${ik}$), intent(out) :: ifaill(*), ifailr(*)
           real(sp), intent(in) :: h(ldh,*), wi(*)
           real(sp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), wr(*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: bothv, fromqr, leftv, noinit, pair, rightv
           integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ksi, ksr, ldwork
           real(sp) :: bignum, eps3, hnorm, smlnum, ulp, unfl, wki, wkr
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters.
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           fromqr = stdlib_lsame( eigsrc, 'Q' )
           noinit = stdlib_lsame( initv, 'N' )
           ! set m to the number of columns required to store the selected
           ! eigenvectors, and standardize the array select.
           m = 0_${ik}$
           pair = .false.
           do k = 1, n
              if( pair ) then
                 pair = .false.
                 select( k ) = .false.
              else
                 if( wi( k )==zero ) then
                    if( select( k ) )m = m + 1_${ik}$
                 else
                    pair = .true.
                    if( select( k ) .or. select( k+1 ) ) then
                       select( k ) = .true.
                       m = m + 2_${ik}$
                    end if
                 end if
              end if
           end do
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then
              info = -2_${ik}$
           else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -11_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info = -13_${ik}$
           else if( mm<m ) then
              info = -14_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SHSEIN', -info )
              return
           end if
           ! quick return if possible.
           if( n==0 )return
           ! set machine-dependent constants.
           unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           ulp = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = unfl*( n / ulp )
           bignum = ( one-ulp ) / smlnum
           ldwork = n + 1_${ik}$
           kl = 1_${ik}$
           kln = 0_${ik}$
           if( fromqr ) then
              kr = 0_${ik}$
           else
              kr = n
           end if
           ksr = 1_${ik}$
           loop_120: do k = 1, n
              if( select( k ) ) then
                 ! compute eigenvector(s) corresponding to w(k).
                 if( fromqr ) then
                    ! if affiliation of eigenvalues is known, check whether
                    ! the matrix splits.
                    ! determine kl and kr such that 1 <= kl <= k <= kr <= n
                    ! and h(kl,kl-1) and h(kr+1,kr) are zero (or kl = 1 or
                    ! kr = n).
                    ! then inverse iteration can be performed with the
                    ! submatrix h(kl:n,kl:n) for a left eigenvector, and with
                    ! the submatrix h(1:kr,1:kr) for a right eigenvector.
                    do i = k, kl + 1, -1
                       if( h( i, i-1 )==zero )go to 30
                    end do
                    30 continue
                    kl = i
                    if( k>kr ) then
                       do i = k, n - 1
                          if( h( i+1, i )==zero )go to 50
                       end do
                       50 continue
                       kr = i
                    end if
                 end if
                 if( kl/=kln ) then
                    kln = kl
                    ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it
                    ! has not ben computed before.
                    hnorm = stdlib${ii}$_slanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
                    if( stdlib${ii}$_sisnan( hnorm ) ) then
                       info = -6_${ik}$
                       return
                    else if( hnorm>zero ) then
                       eps3 = hnorm*ulp
                    else
                       eps3 = smlnum
                    end if
                 end if
                 ! perturb eigenvalue if it is close to any previous
                 ! selected eigenvalues affiliated to the submatrix
                 ! h(kl:kr,kl:kr). close roots are modified by eps3.
                 wkr = wr( k )
                 wki = wi( k )
                 60 continue
                 do i = k - 1, kl, -1
                    if( select( i ) .and. abs( wr( i )-wkr )+abs( wi( i )-wki )<eps3 ) &
                              then
                       wkr = wkr + eps3
                       go to 60
                    end if
                 end do
                 wr( k ) = wkr
                 pair = wki/=zero
                 if( pair ) then
                    ksi = ksr + 1_${ik}$
                 else
                    ksi = ksr
                 end if
                 if( leftv ) then
                    ! compute left eigenvector.
                    call stdlib${ii}$_slaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wkr, wki, vl( &
                    kl, ksr ), vl( kl, ksi ),work, ldwork, work( n*n+n+1 ), eps3, smlnum,bignum, &
                              iinfo )
                    if( iinfo>0_${ik}$ ) then
                       if( pair ) then
                          info = info + 2_${ik}$
                       else
                          info = info + 1_${ik}$
                       end if
                       ifaill( ksr ) = k
                       ifaill( ksi ) = k
                    else
                       ifaill( ksr ) = 0_${ik}$
                       ifaill( ksi ) = 0_${ik}$
                    end if
                    do i = 1, kl - 1
                       vl( i, ksr ) = zero
                    end do
                    if( pair ) then
                       do i = 1, kl - 1
                          vl( i, ksi ) = zero
                       end do
                    end if
                 end if
                 if( rightv ) then
                    ! compute right eigenvector.
                    call stdlib${ii}$_slaein( .true., noinit, kr, h, ldh, wkr, wki,vr( 1_${ik}$, ksr ), vr( 1_${ik}$, &
                              ksi ), work, ldwork,work( n*n+n+1 ), eps3, smlnum, bignum,iinfo )
                    if( iinfo>0_${ik}$ ) then
                       if( pair ) then
                          info = info + 2_${ik}$
                       else
                          info = info + 1_${ik}$
                       end if
                       ifailr( ksr ) = k
                       ifailr( ksi ) = k
                    else
                       ifailr( ksr ) = 0_${ik}$
                       ifailr( ksi ) = 0_${ik}$
                    end if
                    do i = kr + 1, n
                       vr( i, ksr ) = zero
                    end do
                    if( pair ) then
                       do i = kr + 1, n
                          vr( i, ksi ) = zero
                       end do
                    end if
                 end if
                 if( pair ) then
                    ksr = ksr + 2_${ik}$
                 else
                    ksr = ksr + 1_${ik}$
                 end if
              end if
           end do loop_120
           return
     end subroutine stdlib${ii}$_shsein

     module subroutine stdlib${ii}$_dhsein( side, eigsrc, initv, select, n, h, ldh, wr, wi,vl, ldvl, vr, ldvr, &
     !! DHSEIN uses inverse iteration to find specified right and/or left
     !! eigenvectors of a real upper Hessenberg matrix H.
     !! The right eigenvector x and the left eigenvector y of the matrix H
     !! corresponding to an eigenvalue w are defined by:
     !! H * x = w * x,     y**h * H = w * y**h
     !! where y**h denotes the conjugate transpose of the vector y.
               mm, m, work, ifaill,ifailr, info )
        ! -- lapack computational 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 
           character, intent(in) :: eigsrc, initv, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(inout) :: select(*)
           integer(${ik}$), intent(out) :: ifaill(*), ifailr(*)
           real(dp), intent(in) :: h(ldh,*), wi(*)
           real(dp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), wr(*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: bothv, fromqr, leftv, noinit, pair, rightv
           integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ksi, ksr, ldwork
           real(dp) :: bignum, eps3, hnorm, smlnum, ulp, unfl, wki, wkr
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters.
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           fromqr = stdlib_lsame( eigsrc, 'Q' )
           noinit = stdlib_lsame( initv, 'N' )
           ! set m to the number of columns required to store the selected
           ! eigenvectors, and standardize the array select.
           m = 0_${ik}$
           pair = .false.
           do k = 1, n
              if( pair ) then
                 pair = .false.
                 select( k ) = .false.
              else
                 if( wi( k )==zero ) then
                    if( select( k ) )m = m + 1_${ik}$
                 else
                    pair = .true.
                    if( select( k ) .or. select( k+1 ) ) then
                       select( k ) = .true.
                       m = m + 2_${ik}$
                    end if
                 end if
              end if
           end do
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then
              info = -2_${ik}$
           else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -11_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info = -13_${ik}$
           else if( mm<m ) then
              info = -14_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DHSEIN', -info )
              return
           end if
           ! quick return if possible.
           if( n==0 )return
           ! set machine-dependent constants.
           unfl = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           ulp = stdlib${ii}$_dlamch( 'PRECISION' )
           smlnum = unfl*( n / ulp )
           bignum = ( one-ulp ) / smlnum
           ldwork = n + 1_${ik}$
           kl = 1_${ik}$
           kln = 0_${ik}$
           if( fromqr ) then
              kr = 0_${ik}$
           else
              kr = n
           end if
           ksr = 1_${ik}$
           loop_120: do k = 1, n
              if( select( k ) ) then
                 ! compute eigenvector(s) corresponding to w(k).
                 if( fromqr ) then
                    ! if affiliation of eigenvalues is known, check whether
                    ! the matrix splits.
                    ! determine kl and kr such that 1 <= kl <= k <= kr <= n
                    ! and h(kl,kl-1) and h(kr+1,kr) are zero (or kl = 1 or
                    ! kr = n).
                    ! then inverse iteration can be performed with the
                    ! submatrix h(kl:n,kl:n) for a left eigenvector, and with
                    ! the submatrix h(1:kr,1:kr) for a right eigenvector.
                    do i = k, kl + 1, -1
                       if( h( i, i-1 )==zero )go to 30
                    end do
                    30 continue
                    kl = i
                    if( k>kr ) then
                       do i = k, n - 1
                          if( h( i+1, i )==zero )go to 50
                       end do
                       50 continue
                       kr = i
                    end if
                 end if
                 if( kl/=kln ) then
                    kln = kl
                    ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it
                    ! has not ben computed before.
                    hnorm = stdlib${ii}$_dlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
                    if( stdlib${ii}$_disnan( hnorm ) ) then
                       info = -6_${ik}$
                       return
                    else if( hnorm>zero ) then
                       eps3 = hnorm*ulp
                    else
                       eps3 = smlnum
                    end if
                 end if
                 ! perturb eigenvalue if it is close to any previous
                 ! selected eigenvalues affiliated to the submatrix
                 ! h(kl:kr,kl:kr). close roots are modified by eps3.
                 wkr = wr( k )
                 wki = wi( k )
                 60 continue
                 do i = k - 1, kl, -1
                    if( select( i ) .and. abs( wr( i )-wkr )+abs( wi( i )-wki )<eps3 ) &
                              then
                       wkr = wkr + eps3
                       go to 60
                    end if
                 end do
                 wr( k ) = wkr
                 pair = wki/=zero
                 if( pair ) then
                    ksi = ksr + 1_${ik}$
                 else
                    ksi = ksr
                 end if
                 if( leftv ) then
                    ! compute left eigenvector.
                    call stdlib${ii}$_dlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wkr, wki, vl( &
                    kl, ksr ), vl( kl, ksi ),work, ldwork, work( n*n+n+1 ), eps3, smlnum,bignum, &
                              iinfo )
                    if( iinfo>0_${ik}$ ) then
                       if( pair ) then
                          info = info + 2_${ik}$
                       else
                          info = info + 1_${ik}$
                       end if
                       ifaill( ksr ) = k
                       ifaill( ksi ) = k
                    else
                       ifaill( ksr ) = 0_${ik}$
                       ifaill( ksi ) = 0_${ik}$
                    end if
                    do i = 1, kl - 1
                       vl( i, ksr ) = zero
                    end do
                    if( pair ) then
                       do i = 1, kl - 1
                          vl( i, ksi ) = zero
                       end do
                    end if
                 end if
                 if( rightv ) then
                    ! compute right eigenvector.
                    call stdlib${ii}$_dlaein( .true., noinit, kr, h, ldh, wkr, wki,vr( 1_${ik}$, ksr ), vr( 1_${ik}$, &
                              ksi ), work, ldwork,work( n*n+n+1 ), eps3, smlnum, bignum,iinfo )
                    if( iinfo>0_${ik}$ ) then
                       if( pair ) then
                          info = info + 2_${ik}$
                       else
                          info = info + 1_${ik}$
                       end if
                       ifailr( ksr ) = k
                       ifailr( ksi ) = k
                    else
                       ifailr( ksr ) = 0_${ik}$
                       ifailr( ksi ) = 0_${ik}$
                    end if
                    do i = kr + 1, n
                       vr( i, ksr ) = zero
                    end do
                    if( pair ) then
                       do i = kr + 1, n
                          vr( i, ksi ) = zero
                       end do
                    end if
                 end if
                 if( pair ) then
                    ksr = ksr + 2_${ik}$
                 else
                    ksr = ksr + 1_${ik}$
                 end if
              end if
           end do loop_120
           return
     end subroutine stdlib${ii}$_dhsein

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$hsein( side, eigsrc, initv, select, n, h, ldh, wr, wi,vl, ldvl, vr, ldvr, &
     !! DHSEIN: uses inverse iteration to find specified right and/or left
     !! eigenvectors of a real upper Hessenberg matrix H.
     !! The right eigenvector x and the left eigenvector y of the matrix H
     !! corresponding to an eigenvalue w are defined by:
     !! H * x = w * x,     y**h * H = w * y**h
     !! where y**h denotes the conjugate transpose of the vector y.
               mm, m, work, ifaill,ifailr, info )
        ! -- lapack computational 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 
           character, intent(in) :: eigsrc, initv, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(inout) :: select(*)
           integer(${ik}$), intent(out) :: ifaill(*), ifailr(*)
           real(${rk}$), intent(in) :: h(ldh,*), wi(*)
           real(${rk}$), intent(inout) :: vl(ldvl,*), vr(ldvr,*), wr(*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: bothv, fromqr, leftv, noinit, pair, rightv
           integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ksi, ksr, ldwork
           real(${rk}$) :: bignum, eps3, hnorm, smlnum, ulp, unfl, wki, wkr
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters.
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           fromqr = stdlib_lsame( eigsrc, 'Q' )
           noinit = stdlib_lsame( initv, 'N' )
           ! set m to the number of columns required to store the selected
           ! eigenvectors, and standardize the array select.
           m = 0_${ik}$
           pair = .false.
           do k = 1, n
              if( pair ) then
                 pair = .false.
                 select( k ) = .false.
              else
                 if( wi( k )==zero ) then
                    if( select( k ) )m = m + 1_${ik}$
                 else
                    pair = .true.
                    if( select( k ) .or. select( k+1 ) ) then
                       select( k ) = .true.
                       m = m + 2_${ik}$
                    end if
                 end if
              end if
           end do
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then
              info = -2_${ik}$
           else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -11_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info = -13_${ik}$
           else if( mm<m ) then
              info = -14_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DHSEIN', -info )
              return
           end if
           ! quick return if possible.
           if( n==0 )return
           ! set machine-dependent constants.
           unfl = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
           ulp = stdlib${ii}$_${ri}$lamch( 'PRECISION' )
           smlnum = unfl*( n / ulp )
           bignum = ( one-ulp ) / smlnum
           ldwork = n + 1_${ik}$
           kl = 1_${ik}$
           kln = 0_${ik}$
           if( fromqr ) then
              kr = 0_${ik}$
           else
              kr = n
           end if
           ksr = 1_${ik}$
           loop_120: do k = 1, n
              if( select( k ) ) then
                 ! compute eigenvector(s) corresponding to w(k).
                 if( fromqr ) then
                    ! if affiliation of eigenvalues is known, check whether
                    ! the matrix splits.
                    ! determine kl and kr such that 1 <= kl <= k <= kr <= n
                    ! and h(kl,kl-1) and h(kr+1,kr) are zero (or kl = 1 or
                    ! kr = n).
                    ! then inverse iteration can be performed with the
                    ! submatrix h(kl:n,kl:n) for a left eigenvector, and with
                    ! the submatrix h(1:kr,1:kr) for a right eigenvector.
                    do i = k, kl + 1, -1
                       if( h( i, i-1 )==zero )go to 30
                    end do
                    30 continue
                    kl = i
                    if( k>kr ) then
                       do i = k, n - 1
                          if( h( i+1, i )==zero )go to 50
                       end do
                       50 continue
                       kr = i
                    end if
                 end if
                 if( kl/=kln ) then
                    kln = kl
                    ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it
                    ! has not ben computed before.
                    hnorm = stdlib${ii}$_${ri}$lanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
                    if( stdlib${ii}$_${ri}$isnan( hnorm ) ) then
                       info = -6_${ik}$
                       return
                    else if( hnorm>zero ) then
                       eps3 = hnorm*ulp
                    else
                       eps3 = smlnum
                    end if
                 end if
                 ! perturb eigenvalue if it is close to any previous
                 ! selected eigenvalues affiliated to the submatrix
                 ! h(kl:kr,kl:kr). close roots are modified by eps3.
                 wkr = wr( k )
                 wki = wi( k )
                 60 continue
                 do i = k - 1, kl, -1
                    if( select( i ) .and. abs( wr( i )-wkr )+abs( wi( i )-wki )<eps3 ) &
                              then
                       wkr = wkr + eps3
                       go to 60
                    end if
                 end do
                 wr( k ) = wkr
                 pair = wki/=zero
                 if( pair ) then
                    ksi = ksr + 1_${ik}$
                 else
                    ksi = ksr
                 end if
                 if( leftv ) then
                    ! compute left eigenvector.
                    call stdlib${ii}$_${ri}$laein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wkr, wki, vl( &
                    kl, ksr ), vl( kl, ksi ),work, ldwork, work( n*n+n+1 ), eps3, smlnum,bignum, &
                              iinfo )
                    if( iinfo>0_${ik}$ ) then
                       if( pair ) then
                          info = info + 2_${ik}$
                       else
                          info = info + 1_${ik}$
                       end if
                       ifaill( ksr ) = k
                       ifaill( ksi ) = k
                    else
                       ifaill( ksr ) = 0_${ik}$
                       ifaill( ksi ) = 0_${ik}$
                    end if
                    do i = 1, kl - 1
                       vl( i, ksr ) = zero
                    end do
                    if( pair ) then
                       do i = 1, kl - 1
                          vl( i, ksi ) = zero
                       end do
                    end if
                 end if
                 if( rightv ) then
                    ! compute right eigenvector.
                    call stdlib${ii}$_${ri}$laein( .true., noinit, kr, h, ldh, wkr, wki,vr( 1_${ik}$, ksr ), vr( 1_${ik}$, &
                              ksi ), work, ldwork,work( n*n+n+1 ), eps3, smlnum, bignum,iinfo )
                    if( iinfo>0_${ik}$ ) then
                       if( pair ) then
                          info = info + 2_${ik}$
                       else
                          info = info + 1_${ik}$
                       end if
                       ifailr( ksr ) = k
                       ifailr( ksi ) = k
                    else
                       ifailr( ksr ) = 0_${ik}$
                       ifailr( ksi ) = 0_${ik}$
                    end if
                    do i = kr + 1, n
                       vr( i, ksr ) = zero
                    end do
                    if( pair ) then
                       do i = kr + 1, n
                          vr( i, ksi ) = zero
                       end do
                    end if
                 end if
                 if( pair ) then
                    ksr = ksr + 2_${ik}$
                 else
                    ksr = ksr + 1_${ik}$
                 end if
              end if
           end do loop_120
           return
     end subroutine stdlib${ii}$_${ri}$hsein

#:endif
#:endfor

     module subroutine stdlib${ii}$_chsein( side, eigsrc, initv, select, n, h, ldh, w, vl,ldvl, vr, ldvr, mm, &
     !! CHSEIN uses inverse iteration to find specified right and/or left
     !! eigenvectors of a complex upper Hessenberg matrix H.
     !! The right eigenvector x and the left eigenvector y of the matrix H
     !! corresponding to an eigenvalue w are defined by:
     !! H * x = w * x,     y**h * H = w * y**h
     !! where y**h denotes the conjugate transpose of the vector y.
               m, work, rwork, ifaill,ifailr, info )
        ! -- lapack computational 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 
           character, intent(in) :: eigsrc, initv, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: ifaill(*), ifailr(*)
           real(sp), intent(out) :: rwork(*)
           complex(sp), intent(in) :: h(ldh,*)
           complex(sp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), w(*)
           complex(sp), intent(out) :: work(*)
        ! =====================================================================
           ! Parameters 
           
           ! Local Scalars 
           logical(lk) :: bothv, fromqr, leftv, noinit, rightv
           integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ks, ldwork
           real(sp) :: eps3, hnorm, smlnum, ulp, unfl
           complex(sp) :: cdum, wk
           ! Intrinsic Functions 
           ! Statement Functions 
           real(sp) :: cabs1
           ! Statement Function Definitions 
           cabs1( cdum ) = abs( real( cdum,KIND=sp) ) + abs( aimag( cdum ) )
           ! Executable Statements 
           ! decode and test the input parameters.
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           fromqr = stdlib_lsame( eigsrc, 'Q' )
           noinit = stdlib_lsame( initv, 'N' )
           ! set m to the number of columns required to store the selected
           ! eigenvectors.
           m = 0_${ik}$
           do k = 1, n
              if( select( k ) )m = m + 1_${ik}$
           end do
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then
              info = -2_${ik}$
           else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -10_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info = -12_${ik}$
           else if( mm<m ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CHSEIN', -info )
              return
           end if
           ! quick return if possible.
           if( n==0 )return
           ! set machine-dependent constants.
           unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           ulp = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = unfl*( n / ulp )
           ldwork = n
           kl = 1_${ik}$
           kln = 0_${ik}$
           if( fromqr ) then
              kr = 0_${ik}$
           else
              kr = n
           end if
           ks = 1_${ik}$
           loop_100: do k = 1, n
              if( select( k ) ) then
                 ! compute eigenvector(s) corresponding to w(k).
                 if( fromqr ) then
                    ! if affiliation of eigenvalues is known, check whether
                    ! the matrix splits.
                    ! determine kl and kr such that 1 <= kl <= k <= kr <= n
                    ! and h(kl,kl-1) and h(kr+1,kr) are czero (or kl = 1 or
                    ! kr = n).
                    ! then inverse iteration can be performed with the
                    ! submatrix h(kl:n,kl:n) for a left eigenvector, and with
                    ! the submatrix h(1:kr,1:kr) for a right eigenvector.
                    do i = k, kl + 1, -1
                       if( h( i, i-1 )==czero )go to 30
                    end do
                    30 continue
                    kl = i
                    if( k>kr ) then
                       do i = k, n - 1
                          if( h( i+1, i )==czero )go to 50
                       end do
                       50 continue
                       kr = i
                    end if
                 end if
                 if( kl/=kln ) then
                    kln = kl
                    ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it
                    ! has not ben computed before.
                    hnorm = stdlib${ii}$_clanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork )
                    if( stdlib${ii}$_sisnan( hnorm ) ) then
                       info = -6_${ik}$
                       return
                    else if( (hnorm>zero) ) then
                       eps3 = hnorm*ulp
                    else
                       eps3 = smlnum
                    end if
                 end if
                 ! perturb eigenvalue if it is close to any previous
                 ! selected eigenvalues affiliated to the submatrix
                 ! h(kl:kr,kl:kr). close roots are modified by eps3.
                 wk = w( k )
                 60 continue
                 do i = k - 1, kl, -1
                    if( select( i ) .and. cabs1( w( i )-wk )<eps3 ) then
                       wk = wk + eps3
                       go to 60
                    end if
                 end do
                 w( k ) = wk
                 if( leftv ) then
                    ! compute left eigenvector.
                    call stdlib${ii}$_claein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wk, vl( kl, ks )&
                              , work, ldwork, rwork, eps3,smlnum, iinfo )
                    if( iinfo>0_${ik}$ ) then
                       info = info + 1_${ik}$
                       ifaill( ks ) = k
                    else
                       ifaill( ks ) = 0_${ik}$
                    end if
                    do i = 1, kl - 1
                       vl( i, ks ) = czero
                    end do
                 end if
                 if( rightv ) then
                    ! compute right eigenvector.
                    call stdlib${ii}$_claein( .true., noinit, kr, h, ldh, wk, vr( 1_${ik}$, ks ),work, ldwork, &
                              rwork, eps3, smlnum, iinfo )
                    if( iinfo>0_${ik}$ ) then
                       info = info + 1_${ik}$
                       ifailr( ks ) = k
                    else
                       ifailr( ks ) = 0_${ik}$
                    end if
                    do i = kr + 1, n
                       vr( i, ks ) = czero
                    end do
                 end if
                 ks = ks + 1_${ik}$
              end if
           end do loop_100
           return
     end subroutine stdlib${ii}$_chsein

     module subroutine stdlib${ii}$_zhsein( side, eigsrc, initv, select, n, h, ldh, w, vl,ldvl, vr, ldvr, mm, &
     !! ZHSEIN uses inverse iteration to find specified right and/or left
     !! eigenvectors of a complex upper Hessenberg matrix H.
     !! The right eigenvector x and the left eigenvector y of the matrix H
     !! corresponding to an eigenvalue w are defined by:
     !! H * x = w * x,     y**h * H = w * y**h
     !! where y**h denotes the conjugate transpose of the vector y.
               m, work, rwork, ifaill,ifailr, info )
        ! -- lapack computational 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 
           character, intent(in) :: eigsrc, initv, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: ifaill(*), ifailr(*)
           real(dp), intent(out) :: rwork(*)
           complex(dp), intent(in) :: h(ldh,*)
           complex(dp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), w(*)
           complex(dp), intent(out) :: work(*)
        ! =====================================================================
           ! Parameters 
           
           ! Local Scalars 
           logical(lk) :: bothv, fromqr, leftv, noinit, rightv
           integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ks, ldwork
           real(dp) :: eps3, hnorm, smlnum, ulp, unfl
           complex(dp) :: cdum, wk
           ! Intrinsic Functions 
           ! Statement Functions 
           real(dp) :: cabs1
           ! Statement Function Definitions 
           cabs1( cdum ) = abs( real( cdum,KIND=dp) ) + abs( aimag( cdum ) )
           ! Executable Statements 
           ! decode and test the input parameters.
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           fromqr = stdlib_lsame( eigsrc, 'Q' )
           noinit = stdlib_lsame( initv, 'N' )
           ! set m to the number of columns required to store the selected
           ! eigenvectors.
           m = 0_${ik}$
           do k = 1, n
              if( select( k ) )m = m + 1_${ik}$
           end do
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then
              info = -2_${ik}$
           else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -10_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info = -12_${ik}$
           else if( mm<m ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZHSEIN', -info )
              return
           end if
           ! quick return if possible.
           if( n==0 )return
           ! set machine-dependent constants.
           unfl = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           ulp = stdlib${ii}$_dlamch( 'PRECISION' )
           smlnum = unfl*( n / ulp )
           ldwork = n
           kl = 1_${ik}$
           kln = 0_${ik}$
           if( fromqr ) then
              kr = 0_${ik}$
           else
              kr = n
           end if
           ks = 1_${ik}$
           loop_100: do k = 1, n
              if( select( k ) ) then
                 ! compute eigenvector(s) corresponding to w(k).
                 if( fromqr ) then
                    ! if affiliation of eigenvalues is known, check whether
                    ! the matrix splits.
                    ! determine kl and kr such that 1 <= kl <= k <= kr <= n
                    ! and h(kl,kl-1) and h(kr+1,kr) are czero (or kl = 1 or
                    ! kr = n).
                    ! then inverse iteration can be performed with the
                    ! submatrix h(kl:n,kl:n) for a left eigenvector, and with
                    ! the submatrix h(1:kr,1:kr) for a right eigenvector.
                    do i = k, kl + 1, -1
                       if( h( i, i-1 )==czero )go to 30
                    end do
                    30 continue
                    kl = i
                    if( k>kr ) then
                       do i = k, n - 1
                          if( h( i+1, i )==czero )go to 50
                       end do
                       50 continue
                       kr = i
                    end if
                 end if
                 if( kl/=kln ) then
                    kln = kl
                    ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it
                    ! has not ben computed before.
                    hnorm = stdlib${ii}$_zlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork )
                    if( stdlib${ii}$_disnan( hnorm ) ) then
                       info = -6_${ik}$
                       return
                    else if( hnorm>zero ) then
                       eps3 = hnorm*ulp
                    else
                       eps3 = smlnum
                    end if
                 end if
                 ! perturb eigenvalue if it is close to any previous
                 ! selected eigenvalues affiliated to the submatrix
                 ! h(kl:kr,kl:kr). close roots are modified by eps3.
                 wk = w( k )
                 60 continue
                 do i = k - 1, kl, -1
                    if( select( i ) .and. cabs1( w( i )-wk )<eps3 ) then
                       wk = wk + eps3
                       go to 60
                    end if
                 end do
                 w( k ) = wk
                 if( leftv ) then
                    ! compute left eigenvector.
                    call stdlib${ii}$_zlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wk, vl( kl, ks )&
                              , work, ldwork, rwork, eps3,smlnum, iinfo )
                    if( iinfo>0_${ik}$ ) then
                       info = info + 1_${ik}$
                       ifaill( ks ) = k
                    else
                       ifaill( ks ) = 0_${ik}$
                    end if
                    do i = 1, kl - 1
                       vl( i, ks ) = czero
                    end do
                 end if
                 if( rightv ) then
                    ! compute right eigenvector.
                    call stdlib${ii}$_zlaein( .true., noinit, kr, h, ldh, wk, vr( 1_${ik}$, ks ),work, ldwork, &
                              rwork, eps3, smlnum, iinfo )
                    if( iinfo>0_${ik}$ ) then
                       info = info + 1_${ik}$
                       ifailr( ks ) = k
                    else
                       ifailr( ks ) = 0_${ik}$
                    end if
                    do i = kr + 1, n
                       vr( i, ks ) = czero
                    end do
                 end if
                 ks = ks + 1_${ik}$
              end if
           end do loop_100
           return
     end subroutine stdlib${ii}$_zhsein

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     module subroutine stdlib${ii}$_${ci}$hsein( side, eigsrc, initv, select, n, h, ldh, w, vl,ldvl, vr, ldvr, mm, &
     !! ZHSEIN: uses inverse iteration to find specified right and/or left
     !! eigenvectors of a complex upper Hessenberg matrix H.
     !! The right eigenvector x and the left eigenvector y of the matrix H
     !! corresponding to an eigenvalue w are defined by:
     !! H * x = w * x,     y**h * H = w * y**h
     !! where y**h denotes the conjugate transpose of the vector y.
               m, work, rwork, ifaill,ifailr, info )
        ! -- lapack computational 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_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: eigsrc, initv, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: ifaill(*), ifailr(*)
           real(${ck}$), intent(out) :: rwork(*)
           complex(${ck}$), intent(in) :: h(ldh,*)
           complex(${ck}$), intent(inout) :: vl(ldvl,*), vr(ldvr,*), w(*)
           complex(${ck}$), intent(out) :: work(*)
        ! =====================================================================
           ! Parameters 
           
           ! Local Scalars 
           logical(lk) :: bothv, fromqr, leftv, noinit, rightv
           integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ks, ldwork
           real(${ck}$) :: eps3, hnorm, smlnum, ulp, unfl
           complex(${ck}$) :: cdum, wk
           ! Intrinsic Functions 
           ! Statement Functions 
           real(${ck}$) :: cabs1
           ! Statement Function Definitions 
           cabs1( cdum ) = abs( real( cdum,KIND=${ck}$) ) + abs( aimag( cdum ) )
           ! Executable Statements 
           ! decode and test the input parameters.
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           fromqr = stdlib_lsame( eigsrc, 'Q' )
           noinit = stdlib_lsame( initv, 'N' )
           ! set m to the number of columns required to store the selected
           ! eigenvectors.
           m = 0_${ik}$
           do k = 1, n
              if( select( k ) )m = m + 1_${ik}$
           end do
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then
              info = -2_${ik}$
           else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ldh<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -10_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info = -12_${ik}$
           else if( mm<m ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZHSEIN', -info )
              return
           end if
           ! quick return if possible.
           if( n==0 )return
           ! set machine-dependent constants.
           unfl = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' )
           ulp = stdlib${ii}$_${c2ri(ci)}$lamch( 'PRECISION' )
           smlnum = unfl*( n / ulp )
           ldwork = n
           kl = 1_${ik}$
           kln = 0_${ik}$
           if( fromqr ) then
              kr = 0_${ik}$
           else
              kr = n
           end if
           ks = 1_${ik}$
           loop_100: do k = 1, n
              if( select( k ) ) then
                 ! compute eigenvector(s) corresponding to w(k).
                 if( fromqr ) then
                    ! if affiliation of eigenvalues is known, check whether
                    ! the matrix splits.
                    ! determine kl and kr such that 1 <= kl <= k <= kr <= n
                    ! and h(kl,kl-1) and h(kr+1,kr) are czero (or kl = 1 or
                    ! kr = n).
                    ! then inverse iteration can be performed with the
                    ! submatrix h(kl:n,kl:n) for a left eigenvector, and with
                    ! the submatrix h(1:kr,1:kr) for a right eigenvector.
                    do i = k, kl + 1, -1
                       if( h( i, i-1 )==czero )go to 30
                    end do
                    30 continue
                    kl = i
                    if( k>kr ) then
                       do i = k, n - 1
                          if( h( i+1, i )==czero )go to 50
                       end do
                       50 continue
                       kr = i
                    end if
                 end if
                 if( kl/=kln ) then
                    kln = kl
                    ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it
                    ! has not ben computed before.
                    hnorm = stdlib${ii}$_${ci}$lanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork )
                    if( stdlib${ii}$_${c2ri(ci)}$isnan( hnorm ) ) then
                       info = -6_${ik}$
                       return
                    else if( hnorm>zero ) then
                       eps3 = hnorm*ulp
                    else
                       eps3 = smlnum
                    end if
                 end if
                 ! perturb eigenvalue if it is close to any previous
                 ! selected eigenvalues affiliated to the submatrix
                 ! h(kl:kr,kl:kr). close roots are modified by eps3.
                 wk = w( k )
                 60 continue
                 do i = k - 1, kl, -1
                    if( select( i ) .and. cabs1( w( i )-wk )<eps3 ) then
                       wk = wk + eps3
                       go to 60
                    end if
                 end do
                 w( k ) = wk
                 if( leftv ) then
                    ! compute left eigenvector.
                    call stdlib${ii}$_${ci}$laein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wk, vl( kl, ks )&
                              , work, ldwork, rwork, eps3,smlnum, iinfo )
                    if( iinfo>0_${ik}$ ) then
                       info = info + 1_${ik}$
                       ifaill( ks ) = k
                    else
                       ifaill( ks ) = 0_${ik}$
                    end if
                    do i = 1, kl - 1
                       vl( i, ks ) = czero
                    end do
                 end if
                 if( rightv ) then
                    ! compute right eigenvector.
                    call stdlib${ii}$_${ci}$laein( .true., noinit, kr, h, ldh, wk, vr( 1_${ik}$, ks ),work, ldwork, &
                              rwork, eps3, smlnum, iinfo )
                    if( iinfo>0_${ik}$ ) then
                       info = info + 1_${ik}$
                       ifailr( ks ) = k
                    else
                       ifailr( ks ) = 0_${ik}$
                    end if
                    do i = kr + 1, n
                       vr( i, ks ) = czero
                    end do
                 end if
                 ks = ks + 1_${ik}$
              end if
           end do loop_100
           return
     end subroutine stdlib${ii}$_${ci}$hsein

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_strevc( side, howmny, select, n, t, ldt, vl, ldvl, vr,ldvr, mm, m, &
     !! STREVC computes some or all of the right and/or left eigenvectors of
     !! a real upper quasi-triangular matrix T.
     !! Matrices of this type are produced by the Schur factorization of
     !! a real general matrix:  A = Q*T*Q**T, as computed by SHSEQR.
     !! The right eigenvector x and the left eigenvector y of T corresponding
     !! to an eigenvalue w are defined by:
     !! T*x = w*x,     (y**H)*T = w*(y**H)
     !! where y**H denotes the conjugate transpose of y.
     !! The eigenvalues are not input to this routine, but are read directly
     !! from the diagonal blocks of T.
     !! This routine returns the matrices X and/or Y of right and left
     !! eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
     !! input matrix.  If Q is the orthogonal factor that reduces a matrix
     !! A to Schur form T, then Q*X and Q*Y are the matrices of right and
     !! left eigenvectors of A.
               work, info )
        ! -- lapack computational 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 
           character, intent(in) :: howmny, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldt, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(inout) :: select(*)
           real(sp), intent(in) :: t(ldt,*)
           real(sp), intent(inout) :: vl(ldvl,*), vr(ldvr,*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: allv, bothv, leftv, over, pair, rightv, somev
           integer(${ik}$) :: i, ierr, ii, ip, is, j, j1, j2, jnxt, k, ki, n2
           real(sp) :: beta, bignum, emax, ovfl, rec, remax, scale, smin, smlnum, ulp, unfl, &
                     vcrit, vmax, wi, wr, xnorm
           ! Intrinsic Functions 
           ! Local Arrays 
           real(sp) :: x(2_${ik}$,2_${ik}$)
           ! Executable Statements 
           ! decode and test the input parameters
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           allv = stdlib_lsame( howmny, 'A' )
           over = stdlib_lsame( howmny, 'B' )
           somev = stdlib_lsame( howmny, 'S' )
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.allv .and. .not.over .and. .not.somev ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldt<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -8_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info = -10_${ik}$
           else
              ! set m to the number of columns required to store the selected
              ! eigenvectors, standardize the array select if necessary, and
              ! test mm.
              if( somev ) then
                 m = 0_${ik}$
                 pair = .false.
                 do j = 1, n
                    if( pair ) then
                       pair = .false.
                       select( j ) = .false.
                    else
                       if( j<n ) then
                          if( t( j+1, j )==zero ) then
                             if( select( j ) )m = m + 1_${ik}$
                          else
                             pair = .true.
                             if( select( j ) .or. select( j+1 ) ) then
                                select( j ) = .true.
                                m = m + 2_${ik}$
                             end if
                          end if
                       else
                          if( select( n ) )m = m + 1_${ik}$
                       end if
                    end if
                 end do
              else
                 m = n
              end if
              if( mm<m ) then
                 info = -11_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'STREVC', -info )
              return
           end if
           ! quick return if possible.
           if( n==0 )return
           ! set the constants to control overflow.
           unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           ovfl = one / unfl
           call stdlib${ii}$_slabad( unfl, ovfl )
           ulp = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = unfl*( n / ulp )
           bignum = ( one-ulp ) / smlnum
           ! 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 ) = zero
              do i = 1, j - 1
                 work( j ) = work( j ) + abs( t( i, j ) )
              end do
           end do
           ! index ip is used to specify the real or complex eigenvalue:
             ! ip = 0, real eigenvalue,
                  ! 1, first of conjugate complex pair: (wr,wi)
                 ! -1, second of conjugate complex pair: (wr,wi)
           n2 = 2_${ik}$*n
           if( rightv ) then
              ! compute right eigenvectors.
              ip = 0_${ik}$
              is = m
              loop_140: do ki = n, 1, -1
                 if( ip==1 )go to 130
                 if( ki==1 )go to 40
                 if( t( ki, ki-1 )==zero )go to 40
                 ip = -1_${ik}$
                 40 continue
                 if( somev ) then
                    if( ip==0_${ik}$ ) then
                       if( .not.select( ki ) )go to 130
                    else
                       if( .not.select( ki-1 ) )go to 130
                    end if
                 end if
                 ! compute the ki-th eigenvalue (wr,wi).
                 wr = t( ki, ki )
                 wi = zero
                 if( ip/=0_${ik}$ )wi = sqrt( abs( t( ki, ki-1 ) ) )*sqrt( abs( t( ki-1, ki ) ) )
                 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
                 if( ip==0_${ik}$ ) then
                    ! real right eigenvector
                    work( ki+n ) = one
                    ! form right-hand side
                    do k = 1, ki - 1
                       work( k+n ) = -t( k, ki )
                    end do
                    ! solve the upper quasi-triangular system:
                       ! (t(1:ki-1,1:ki-1) - wr)*x = scale*work.
                    jnxt = ki - 1_${ik}$
                    loop_60: do j = ki - 1, 1, -1
                       if( j>jnxt )cycle loop_60
                       j1 = j
                       j2 = j
                       jnxt = j - 1_${ik}$
                       if( j>1_${ik}$ ) then
                          if( t( j, j-1 )/=zero ) then
                             j1 = j - 1_${ik}$
                             jnxt = j - 2_${ik}$
                          end if
                       end if
                       if( j1==j2 ) then
                          ! 1-by-1 diagonal block
                          call stdlib${ii}$_slaln2( .false., 1_${ik}$, 1_${ik}$, smin, one, t( j, j ),ldt, one, one, &
                                    work( j+n ), n, wr,zero, x, 2_${ik}$, scale, xnorm, ierr )
                          ! scale x(1,1) to avoid overflow when updating
                          ! the right-hand side.
                          if( xnorm>one ) then
                             if( work( j )>bignum / xnorm ) then
                                x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ ) / xnorm
                                scale = scale / xnorm
                             end if
                          end if
                          ! scale if necessary
                          if( scale/=one )call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ )
                          work( j+n ) = x( 1_${ik}$, 1_${ik}$ )
                          ! update right-hand side
                          call stdlib${ii}$_saxpy( j-1, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ )
                                    
                       else
                          ! 2-by-2 diagonal block
                          call stdlib${ii}$_slaln2( .false., 2_${ik}$, 1_${ik}$, smin, one,t( j-1, j-1 ), ldt, one, &
                                    one,work( j-1+n ), n, wr, zero, x, 2_${ik}$,scale, xnorm, ierr )
                          ! scale x(1,1) and x(2,1) to avoid overflow when
                          ! updating the right-hand side.
                          if( xnorm>one ) then
                             beta = max( work( j-1 ), work( j ) )
                             if( beta>bignum / xnorm ) then
                                x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ ) / xnorm
                                x( 2_${ik}$, 1_${ik}$ ) = x( 2_${ik}$, 1_${ik}$ ) / xnorm
                                scale = scale / xnorm
                             end if
                          end if
                          ! scale if necessary
                          if( scale/=one )call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ )
                          work( j-1+n ) = x( 1_${ik}$, 1_${ik}$ )
                          work( j+n ) = x( 2_${ik}$, 1_${ik}$ )
                          ! update right-hand side
                          call stdlib${ii}$_saxpy( j-2, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j-1 ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ )
                                    
                          call stdlib${ii}$_saxpy( j-2, -x( 2_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ )
                                    
                       end if
                    end do loop_60
                    ! copy the vector x or q*x to vr and normalize.
                    if( .not.over ) then
                       call stdlib${ii}$_scopy( ki, work( 1_${ik}$+n ), 1_${ik}$, vr( 1_${ik}$, is ), 1_${ik}$ )
                       ii = stdlib${ii}$_isamax( ki, vr( 1_${ik}$, is ), 1_${ik}$ )
                       remax = one / abs( vr( ii, is ) )
                       call stdlib${ii}$_sscal( ki, remax, vr( 1_${ik}$, is ), 1_${ik}$ )
                       do k = ki + 1, n
                          vr( k, is ) = zero
                       end do
                    else
                       if( ki>1_${ik}$ )call stdlib${ii}$_sgemv( 'N', n, ki-1, one, vr, ldvr,work( 1_${ik}$+n ), 1_${ik}$, &
                                 work( ki+n ),vr( 1_${ik}$, ki ), 1_${ik}$ )
                       ii = stdlib${ii}$_isamax( n, vr( 1_${ik}$, ki ), 1_${ik}$ )
                       remax = one / abs( vr( ii, ki ) )
                       call stdlib${ii}$_sscal( n, remax, vr( 1_${ik}$, ki ), 1_${ik}$ )
                    end if
                 else
                    ! complex right eigenvector.
                    ! initial solve
                      ! [ (t(ki-1,ki-1) t(ki-1,ki) ) - (wr + i* wi)]*x = 0.
                      ! [ (t(ki,ki-1)   t(ki,ki)   )               ]
                    if( abs( t( ki-1, ki ) )>=abs( t( ki, ki-1 ) ) ) then
                       work( ki-1+n ) = one
                       work( ki+n2 ) = wi / t( ki-1, ki )
                    else
                       work( ki-1+n ) = -wi / t( ki, ki-1 )
                       work( ki+n2 ) = one
                    end if
                    work( ki+n ) = zero
                    work( ki-1+n2 ) = zero
                    ! form right-hand side
                    do k = 1, ki - 2
                       work( k+n ) = -work( ki-1+n )*t( k, ki-1 )
                       work( k+n2 ) = -work( ki+n2 )*t( k, ki )
                    end do
                    ! solve upper quasi-triangular system:
                    ! (t(1:ki-2,1:ki-2) - (wr+i*wi))*x = scale*(work+i*work2)
                    jnxt = ki - 2_${ik}$
                    loop_90: do j = ki - 2, 1, -1
                       if( j>jnxt )cycle loop_90
                       j1 = j
                       j2 = j
                       jnxt = j - 1_${ik}$
                       if( j>1_${ik}$ ) then
                          if( t( j, j-1 )/=zero ) then
                             j1 = j - 1_${ik}$
                             jnxt = j - 2_${ik}$
                          end if
                       end if
                       if( j1==j2 ) then
                          ! 1-by-1 diagonal block
                          call stdlib${ii}$_slaln2( .false., 1_${ik}$, 2_${ik}$, smin, one, t( j, j ),ldt, one, one, &
                                    work( j+n ), n, wr, wi,x, 2_${ik}$, scale, xnorm, ierr )
                          ! scale x(1,1) and x(1,2) to avoid overflow when
                          ! updating the right-hand side.
                          if( xnorm>one ) then
                             if( work( j )>bignum / xnorm ) then
                                x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ ) / xnorm
                                x( 1_${ik}$, 2_${ik}$ ) = x( 1_${ik}$, 2_${ik}$ ) / xnorm
                                scale = scale / xnorm
                             end if
                          end if
                          ! scale if necessary
                          if( scale/=one ) then
                             call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ )
                             call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n2 ), 1_${ik}$ )
                          end if
                          work( j+n ) = x( 1_${ik}$, 1_${ik}$ )
                          work( j+n2 ) = x( 1_${ik}$, 2_${ik}$ )
                          ! update the right-hand side
                          call stdlib${ii}$_saxpy( j-1, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ )
                                    
                          call stdlib${ii}$_saxpy( j-1, -x( 1_${ik}$, 2_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n2 ), 1_${ik}$ )
                                    
                       else
                          ! 2-by-2 diagonal block
                          call stdlib${ii}$_slaln2( .false., 2_${ik}$, 2_${ik}$, smin, one,t( j-1, j-1 ), ldt, one, &
                                    one,work( j-1+n ), n, wr, wi, x, 2_${ik}$, scale,xnorm, ierr )
                          ! scale x to avoid overflow when updating
                          ! the right-hand side.
                          if( xnorm>one ) then
                             beta = max( work( j-1 ), work( j ) )
                             if( beta>bignum / xnorm ) then
                                rec = one / xnorm
                                x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ )*rec
                                x( 1_${ik}$, 2_${ik}$ ) = x( 1_${ik}$, 2_${ik}$ )*rec
                                x( 2_${ik}$, 1_${ik}$ ) = x( 2_${ik}$, 1_${ik}$ )*rec
                                x( 2_${ik}$, 2_${ik}$ ) = x( 2_${ik}$, 2_${ik}$ )*rec
                                scale = scale*rec
                             end if
                          end if
                          ! scale if necessary
                          if( scale/=one ) then
                             call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ )
                             call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n2 ), 1_${ik}$ )
                          end if
                          work( j-1+n ) = x( 1_${ik}$, 1_${ik}$ )
                          work( j+n ) = x( 2_${ik}$, 1_${ik}$ )
                          work( j-1+n2 ) = x( 1_${ik}$, 2_${ik}$ )
                          work( j+n2 ) = x( 2_${ik}$, 2_${ik}$ )
                          ! update the right-hand side
                          call stdlib${ii}$_saxpy( j-2, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j-1 ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ )
                                    
                          call stdlib${ii}$_saxpy( j-2, -x( 2_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ )
                                    
                          call stdlib${ii}$_saxpy( j-2, -x( 1_${ik}$, 2_${ik}$ ), t( 1_${ik}$, j-1 ), 1_${ik}$,work( 1_${ik}$+n2 ), 1_${ik}$ )
                                    
                          call stdlib${ii}$_saxpy( j-2, -x( 2_${ik}$, 2_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n2 ), 1_${ik}$ )
                                    
                       end if
                    end do loop_90
                    ! copy the vector x or q*x to vr and normalize.
                    if( .not.over ) then
                       call stdlib${ii}$_scopy( ki, work( 1_${ik}$+n ), 1_${ik}$, vr( 1_${ik}$, is-1 ), 1_${ik}$ )
                       call stdlib${ii}$_scopy( ki, work( 1_${ik}$+n2 ), 1_${ik}$, vr( 1_${ik}$, is ), 1_${ik}$ )
                       emax = zero
                       do k = 1, ki
                          emax = max( emax, abs( vr( k, is-1 ) )+abs( vr( k, is ) ) )
                       end do
                       remax = one / emax
                       call stdlib${ii}$_sscal( ki, remax, vr( 1_${ik}$, is-1 ), 1_${ik}$ )
                       call stdlib${ii}$_sscal( ki, remax, vr( 1_${ik}$, is ), 1_${ik}$ )
                       do k = ki + 1, n
                          vr( k, is-1 ) = zero
                          vr( k, is ) = zero
                       end do
                    else
                       if( ki>2_${ik}$ ) then
                          call stdlib${ii}$_sgemv( 'N', n, ki-2, one, vr, ldvr,work( 1_${ik}$+n ), 1_${ik}$, work( ki-&
                                    1_${ik}$+n ),vr( 1_${ik}$, ki-1 ), 1_${ik}$ )
                          call stdlib${ii}$_sgemv( 'N', n, ki-2, one, vr, ldvr,work( 1_${ik}$+n2 ), 1_${ik}$, work( &
                                    ki+n2 ),vr( 1_${ik}$, ki ), 1_${ik}$ )
                       else
                          call stdlib${ii}$_sscal( n, work( ki-1+n ), vr( 1_${ik}$, ki-1 ), 1_${ik}$ )
                          call stdlib${ii}$_sscal( n, work( ki+n2 ), vr( 1_${ik}$, ki ), 1_${ik}$ )
                       end if
                       emax = zero
                       do k = 1, n
                          emax = max( emax, abs( vr( k, ki-1 ) )+abs( vr( k, ki ) ) )
                       end do
                       remax = one / emax
                       call stdlib${ii}$_sscal( n, remax, vr( 1_${ik}$, ki-1 ), 1_${ik}$ )
                       call stdlib${ii}$_sscal( n, remax, vr( 1_${ik}$, ki ), 1_${ik}$ )
                    end if
                 end if
                 is = is - 1_${ik}$
                 if( ip/=0_${ik}$ )is = is - 1_${ik}$
                 130 continue
                 if( ip==1_${ik}$ )ip = 0_${ik}$
                 if( ip==-1_${ik}$ )ip = 1_${ik}$
              end do loop_140
           end if
           if( leftv ) then
              ! compute left eigenvectors.
              ip = 0_${ik}$
              is = 1_${ik}$
              loop_260: do ki = 1, n
                 if( ip==-1 )go to 250
                 if( ki==n )go to 150
                 if( t( ki+1, ki )==zero )go to 150
                 ip = 1_${ik}$
                 150 continue
                 if( somev ) then
                    if( .not.select( ki ) )go to 250
                 end if
                 ! compute the ki-th eigenvalue (wr,wi).
                 wr = t( ki, ki )
                 wi = zero
                 if( ip/=0_${ik}$ )wi = sqrt( abs( t( ki, ki+1 ) ) )*sqrt( abs( t( ki+1, ki ) ) )
                 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
                 if( ip==0_${ik}$ ) then
                    ! real left eigenvector.
                    work( ki+n ) = one
                    ! form right-hand side
                    do k = ki + 1, n
                       work( k+n ) = -t( ki, k )
                    end do
                    ! solve the quasi-triangular system:
                       ! (t(ki+1:n,ki+1:n) - wr)**t*x = scale*work
                    vmax = one
                    vcrit = bignum
                    jnxt = ki + 1_${ik}$
                    loop_170: do j = ki + 1, n
                       if( j<jnxt )cycle loop_170
                       j1 = j
                       j2 = j
                       jnxt = j + 1_${ik}$
                       if( j<n ) then
                          if( t( j+1, j )/=zero ) then
                             j2 = j + 1_${ik}$
                             jnxt = j + 2_${ik}$
                          end if
                       end if
                       if( j1==j2 ) then
                          ! 1-by-1 diagonal block
                          ! scale if necessary to avoid overflow when forming
                          ! the right-hand side.
                          if( work( j )>vcrit ) then
                             rec = one / vmax
                             call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ )
                             vmax = one
                             vcrit = bignum
                          end if
                          work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-1, t( ki+1, j ), 1_${ik}$,work( &
                                    ki+1+n ), 1_${ik}$ )
                          ! solve (t(j,j)-wr)**t*x = work
                          call stdlib${ii}$_slaln2( .false., 1_${ik}$, 1_${ik}$, smin, one, t( j, j ),ldt, one, one, &
                                    work( j+n ), n, wr,zero, x, 2_${ik}$, scale, xnorm, ierr )
                          ! scale if necessary
                          if( scale/=one )call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ )
                                    
                          work( j+n ) = x( 1_${ik}$, 1_${ik}$ )
                          vmax = max( abs( work( j+n ) ), vmax )
                          vcrit = bignum / vmax
                       else
                          ! 2-by-2 diagonal block
                          ! scale if necessary to avoid overflow when forming
                          ! the right-hand side.
                          beta = max( work( j ), work( j+1 ) )
                          if( beta>vcrit ) then
                             rec = one / vmax
                             call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ )
                             vmax = one
                             vcrit = bignum
                          end if
                          work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-1, t( ki+1, j ), 1_${ik}$,work( &
                                    ki+1+n ), 1_${ik}$ )
                          work( j+1+n ) = work( j+1+n ) -stdlib${ii}$_sdot( j-ki-1, t( ki+1, j+1 ), 1_${ik}$,&
                                    work( ki+1+n ), 1_${ik}$ )
                          ! solve
                            ! [t(j,j)-wr   t(j,j+1)     ]**t* x = scale*( work1 )
                            ! [t(j+1,j)    t(j+1,j+1)-wr]               ( work2 )
                          call stdlib${ii}$_slaln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j, j ),ldt, one, one, &
                                    work( j+n ), n, wr,zero, x, 2_${ik}$, scale, xnorm, ierr )
                          ! scale if necessary
                          if( scale/=one )call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ )
                                    
                          work( j+n ) = x( 1_${ik}$, 1_${ik}$ )
                          work( j+1+n ) = x( 2_${ik}$, 1_${ik}$ )
                          vmax = max( abs( work( j+n ) ),abs( work( j+1+n ) ), vmax )
                          vcrit = bignum / vmax
                       end if
                    end do loop_170
                    ! copy the vector x or q*x to vl and normalize.
                    if( .not.over ) then
                       call stdlib${ii}$_scopy( n-ki+1, work( ki+n ), 1_${ik}$, vl( ki, is ), 1_${ik}$ )
                       ii = stdlib${ii}$_isamax( n-ki+1, vl( ki, is ), 1_${ik}$ ) + ki - 1_${ik}$
                       remax = one / abs( vl( ii, is ) )
                       call stdlib${ii}$_sscal( n-ki+1, remax, vl( ki, is ), 1_${ik}$ )
                       do k = 1, ki - 1
                          vl( k, is ) = zero
                       end do
                    else
                       if( ki<n )call stdlib${ii}$_sgemv( 'N', n, n-ki, one, vl( 1_${ik}$, ki+1 ), ldvl,work( &
                                 ki+1+n ), 1_${ik}$, work( ki+n ),vl( 1_${ik}$, ki ), 1_${ik}$ )
                       ii = stdlib${ii}$_isamax( n, vl( 1_${ik}$, ki ), 1_${ik}$ )
                       remax = one / abs( vl( ii, ki ) )
                       call stdlib${ii}$_sscal( n, remax, vl( 1_${ik}$, ki ), 1_${ik}$ )
                    end if
                 else
                    ! complex left eigenvector.
                     ! initial solve:
                       ! ((t(ki,ki)    t(ki,ki+1) )**t - (wr - i* wi))*x = 0.
                       ! ((t(ki+1,ki) t(ki+1,ki+1))                )
                    if( abs( t( ki, ki+1 ) )>=abs( t( ki+1, ki ) ) ) then
                       work( ki+n ) = wi / t( ki, ki+1 )
                       work( ki+1+n2 ) = one
                    else
                       work( ki+n ) = one
                       work( ki+1+n2 ) = -wi / t( ki+1, ki )
                    end if
                    work( ki+1+n ) = zero
                    work( ki+n2 ) = zero
                    ! form right-hand side
                    do k = ki + 2, n
                       work( k+n ) = -work( ki+n )*t( ki, k )
                       work( k+n2 ) = -work( ki+1+n2 )*t( ki+1, k )
                    end do
                    ! solve complex quasi-triangular system:
                    ! ( t(ki+2,n:ki+2,n) - (wr-i*wi) )*x = work1+i*work2
                    vmax = one
                    vcrit = bignum
                    jnxt = ki + 2_${ik}$
                    loop_200: do j = ki + 2, n
                       if( j<jnxt )cycle loop_200
                       j1 = j
                       j2 = j
                       jnxt = j + 1_${ik}$
                       if( j<n ) then
                          if( t( j+1, j )/=zero ) then
                             j2 = j + 1_${ik}$
                             jnxt = j + 2_${ik}$
                          end if
                       end if
                       if( j1==j2 ) then
                          ! 1-by-1 diagonal block
                          ! scale if necessary to avoid overflow when
                          ! forming the right-hand side elements.
                          if( work( j )>vcrit ) then
                             rec = one / vmax
                             call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ )
                             call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n2 ), 1_${ik}$ )
                             vmax = one
                             vcrit = bignum
                          end if
                          work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( &
                                    ki+2+n ), 1_${ik}$ )
                          work( j+n2 ) = work( j+n2 ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( &
                                    ki+2+n2 ), 1_${ik}$ )
                          ! solve (t(j,j)-(wr-i*wi))*(x11+i*x12)= wk+i*wk2
                          call stdlib${ii}$_slaln2( .false., 1_${ik}$, 2_${ik}$, smin, one, t( j, j ),ldt, one, one, &
                                    work( j+n ), n, wr,-wi, x, 2_${ik}$, scale, xnorm, ierr )
                          ! scale if necessary
                          if( scale/=one ) then
                             call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ )
                             call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n2 ), 1_${ik}$ )
                          end if
                          work( j+n ) = x( 1_${ik}$, 1_${ik}$ )
                          work( j+n2 ) = x( 1_${ik}$, 2_${ik}$ )
                          vmax = max( abs( work( j+n ) ),abs( work( j+n2 ) ), vmax )
                          vcrit = bignum / vmax
                       else
                          ! 2-by-2 diagonal block
                          ! scale if necessary to avoid overflow when forming
                          ! the right-hand side elements.
                          beta = max( work( j ), work( j+1 ) )
                          if( beta>vcrit ) then
                             rec = one / vmax
                             call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ )
                             call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n2 ), 1_${ik}$ )
                             vmax = one
                             vcrit = bignum
                          end if
                          work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( &
                                    ki+2+n ), 1_${ik}$ )
                          work( j+n2 ) = work( j+n2 ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( &
                                    ki+2+n2 ), 1_${ik}$ )
                          work( j+1+n ) = work( j+1+n ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j+1 ), 1_${ik}$,&
                                    work( ki+2+n ), 1_${ik}$ )
                          work( j+1+n2 ) = work( j+1+n2 ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j+1 ), 1_${ik}$,&
                                    work( ki+2+n2 ), 1_${ik}$ )
                          ! solve 2-by-2 complex linear equation
                            ! ([t(j,j)   t(j,j+1)  ]**t-(wr-i*wi)*i)*x = scale*b
                            ! ([t(j+1,j) t(j+1,j+1)]               )
                          call stdlib${ii}$_slaln2( .true., 2_${ik}$, 2_${ik}$, smin, one, t( j, j ),ldt, one, one, &
                                    work( j+n ), n, wr,-wi, x, 2_${ik}$, scale, xnorm, ierr )
                          ! scale if necessary
                          if( scale/=one ) then
                             call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ )
                             call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n2 ), 1_${ik}$ )
                          end if
                          work( j+n ) = x( 1_${ik}$, 1_${ik}$ )
                          work( j+n2 ) = x( 1_${ik}$, 2_${ik}$ )
                          work( j+1+n ) = x( 2_${ik}$, 1_${ik}$ )
                          work( j+1+n2 ) = x( 2_${ik}$, 2_${ik}$ )
                          vmax = max( abs( x( 1_${ik}$, 1_${ik}$ ) ), abs( x( 1_${ik}$, 2_${ik}$ ) ),abs( x( 2_${ik}$, 1_${ik}$ ) ), abs( x(&
                                     2_${ik}$, 2_${ik}$ ) ), vmax )
                          vcrit = bignum / vmax
                       end if
                    end do loop_200
                    ! copy the vector x or q*x to vl and normalize.
                    if( .not.over ) then
                       call stdlib${ii}$_scopy( n-ki+1, work( ki+n ), 1_${ik}$, vl( ki, is ), 1_${ik}$ )
                       call stdlib${ii}$_scopy( n-ki+1, work( ki+n2 ), 1_${ik}$, vl( ki, is+1 ),1_${ik}$ )
                       emax = zero
                       do k = ki, n
                          emax = max( emax, abs( vl( k, is ) )+abs( vl( k, is+1 ) ) )
                       end do
                       remax = one / emax
                       call stdlib${ii}$_sscal( n-ki+1, remax, vl( ki, is ), 1_${ik}$ )
                       call stdlib${ii}$_sscal( n-ki+1, remax, vl( ki, is+1 ), 1_${ik}$ )
                       do k = 1, ki - 1
                          vl( k, is ) = zero
                          vl( k, is+1 ) = zero
                       end do
                    else
                       if( ki<n-1 ) then
                          call stdlib${ii}$_sgemv( 'N', n, n-ki-1, one, vl( 1_${ik}$, ki+2 ),ldvl, work( ki+2+&
                                    n ), 1_${ik}$, work( ki+n ),vl( 1_${ik}$, ki ), 1_${ik}$ )
                          call stdlib${ii}$_sgemv( 'N', n, n-ki-1, one, vl( 1_${ik}$, ki+2 ),ldvl, work( ki+2+&
                                    n2 ), 1_${ik}$,work( ki+1+n2 ), vl( 1_${ik}$, ki+1 ), 1_${ik}$ )
                       else
                          call stdlib${ii}$_sscal( n, work( ki+n ), vl( 1_${ik}$, ki ), 1_${ik}$ )
                          call stdlib${ii}$_sscal( n, work( ki+1+n2 ), vl( 1_${ik}$, ki+1 ), 1_${ik}$ )
                       end if
                       emax = zero
                       do k = 1, n
                          emax = max( emax, abs( vl( k, ki ) )+abs( vl( k, ki+1 ) ) )
                       end do
                       remax = one / emax
                       call stdlib${ii}$_sscal( n, remax, vl( 1_${ik}$, ki ), 1_${ik}$ )
                       call stdlib${ii}$_sscal( n, remax, vl( 1_${ik}$, ki+1 ), 1_${ik}$ )
                    end if
                 end if
                 is = is + 1_${ik}$
                 if( ip/=0_${ik}$ )is = is + 1_${ik}$
                 250 continue
                 if( ip==-1_${ik}$ )ip = 0_${ik}$
                 if( ip==1_${ik}$ )ip = -1_${ik}$
              end do loop_260
           end if
           return
     end subroutine stdlib${ii}$_strevc

     pure module subroutine stdlib${ii}$_dtrevc( side, howmny, select, n, t, ldt, vl, ldvl, vr,ldvr, mm, m, &
     !! DTREVC computes some or all of the right and/or left eigenvectors of
     !! a real upper quasi-triangular matrix T.
     !! Matrices of this type are produced by the Schur factorization of
     !! a real general matrix:  A = Q*T*Q**T, as computed by DHSEQR.
     !! The right eigenvector x and the left eigenvector y of T corresponding
     !! to an eigenvalue w are defined by:
     !! T*x = w*x,     (y**H)*T = w*(y**H)
     !! where y**H denotes the conjugate transpose of y.
     !! The eigenvalues are not input to this routine, but are read directly
     !! from the diagonal blocks of T.
     !! This routine returns the matrices X and/or Y of right and left
     !! eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
     !! input matrix.  If Q is the orthogonal factor that reduces a matrix
     !! A to Schur form T, then Q*X and Q*Y are the matrices of right and
     !! left eigenvectors of A.
               work, info )
        ! -- lapack computational 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 
           character, intent(in) :: howmny, side
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: ldt, ldvl, ldvr, mm, n
           ! Array Arguments 
           logical(lk), intent(inout) :: select(*)
           real(dp), intent(in) :: t(ldt,*)
           real(dp), intent(inout) :: vl(ldvl,*), vr(ldvr,*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: allv, bothv, leftv, over, pair, rightv, somev
           integer(${ik}$) :: i, ierr, ii, ip, is, j, j1, j2, jnxt, k, ki, n2
           real(dp) :: beta, bignum, emax, ovfl, rec, remax, scale, smin, smlnum, ulp, unfl, &
                     vcrit, vmax, wi, wr, xnorm
           ! Intrinsic Functions 
           ! Local Arrays 
           real(dp) :: x(2_${ik}$,2_${ik}$)
           ! Executable Statements 
           ! decode and test the input parameters
           bothv = stdlib_lsame( side, 'B' )
           rightv = stdlib_lsame( side, 'R' ) .or. bothv
           leftv = stdlib_lsame( side, 'L' ) .or. bothv
           allv = stdlib_lsame( howmny, 'A' )
           over = stdlib_lsame( howmny, 'B' )
           somev = stdlib_lsame( howmny, 'S' )
           info = 0_${ik}$
           if( .not.rightv .and. .not.leftv ) then
              info = -1_${ik}$
           else if( .not.allv .and. .not.over .and. .not.somev ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldt<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then
              info = -8_${ik}$
           else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then
              info =