stdlib_lapack_eigv_tridiag2.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_slamrg( n1, n2, a, strd1, strd2, index )
     !! SLAMRG will create a permutation list which will merge the elements
     !! of A (which is composed of two independently sorted sets) into a
     !! single set which is sorted in ascending order.
        ! -- 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) :: n1, n2, strd1, strd2
           ! Array Arguments 
           integer(${ik}$), intent(out) :: index(*)
           real(sp), intent(in) :: a(*)
        ! =====================================================================
           ! Local Scalars 
           integer(${ik}$) :: i, ind1, ind2, n1sv, n2sv
           ! Executable Statements 
           n1sv = n1
           n2sv = n2
           if( strd1>0_${ik}$ ) then
              ind1 = 1_${ik}$
           else
              ind1 = n1
           end if
           if( strd2>0_${ik}$ ) then
              ind2 = 1_${ik}$ + n1
           else
              ind2 = n1 + n2
           end if
           i = 1_${ik}$
           ! while ( (n1sv > 0)
           10 continue
           if( n1sv>0_${ik}$ .and. n2sv>0_${ik}$ ) then
              if( a( ind1 )<=a( ind2 ) ) then
                 index( i ) = ind1
                 i = i + 1_${ik}$
                 ind1 = ind1 + strd1
                 n1sv = n1sv - 1_${ik}$
              else
                 index( i ) = ind2
                 i = i + 1_${ik}$
                 ind2 = ind2 + strd2
                 n2sv = n2sv - 1_${ik}$
              end if
              go to 10
           end if
           ! end while
           if( n1sv==0_${ik}$ ) then
              do n1sv = 1, n2sv
                 index( i ) = ind2
                 i = i + 1_${ik}$
                 ind2 = ind2 + strd2
              end do
           else
           ! n2sv == 0
              do n2sv = 1, n1sv
                 index( i ) = ind1
                 i = i + 1_${ik}$
                 ind1 = ind1 + strd1
              end do
           end if
           return
     end subroutine stdlib${ii}$_slamrg

     pure module subroutine stdlib${ii}$_dlamrg( n1, n2, a, dtrd1, dtrd2, index )
     !! DLAMRG will create a permutation list which will merge the elements
     !! of A (which is composed of two independently sorted sets) into a
     !! single set which is sorted in ascending order.
        ! -- 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) :: dtrd1, dtrd2, n1, n2
           ! Array Arguments 
           integer(${ik}$), intent(out) :: index(*)
           real(dp), intent(in) :: a(*)
        ! =====================================================================
           ! Local Scalars 
           integer(${ik}$) :: i, ind1, ind2, n1sv, n2sv
           ! Executable Statements 
           n1sv = n1
           n2sv = n2
           if( dtrd1>0_${ik}$ ) then
              ind1 = 1_${ik}$
           else
              ind1 = n1
           end if
           if( dtrd2>0_${ik}$ ) then
              ind2 = 1_${ik}$ + n1
           else
              ind2 = n1 + n2
           end if
           i = 1_${ik}$
           ! while ( (n1sv > 0)
           10 continue
           if( n1sv>0_${ik}$ .and. n2sv>0_${ik}$ ) then
              if( a( ind1 )<=a( ind2 ) ) then
                 index( i ) = ind1
                 i = i + 1_${ik}$
                 ind1 = ind1 + dtrd1
                 n1sv = n1sv - 1_${ik}$
              else
                 index( i ) = ind2
                 i = i + 1_${ik}$
                 ind2 = ind2 + dtrd2
                 n2sv = n2sv - 1_${ik}$
              end if
              go to 10
           end if
           ! end while
           if( n1sv==0_${ik}$ ) then
              do n1sv = 1, n2sv
                 index( i ) = ind2
                 i = i + 1_${ik}$
                 ind2 = ind2 + dtrd2
              end do
           else
           ! n2sv == 0
              do n2sv = 1, n1sv
                 index( i ) = ind1
                 i = i + 1_${ik}$
                 ind1 = ind1 + dtrd1
              end do
           end if
           return
     end subroutine stdlib${ii}$_dlamrg

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$lamrg( n1, n2, a, dtrd1, dtrd2, index )
     !! DLAMRG: will create a permutation list which will merge the elements
     !! of A (which is composed of two independently sorted sets) into a
     !! single set which is sorted in ascending order.
        ! -- 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) :: dtrd1, dtrd2, n1, n2
           ! Array Arguments 
           integer(${ik}$), intent(out) :: index(*)
           real(${rk}$), intent(in) :: a(*)
        ! =====================================================================
           ! Local Scalars 
           integer(${ik}$) :: i, ind1, ind2, n1sv, n2sv
           ! Executable Statements 
           n1sv = n1
           n2sv = n2
           if( dtrd1>0_${ik}$ ) then
              ind1 = 1_${ik}$
           else
              ind1 = n1
           end if
           if( dtrd2>0_${ik}$ ) then
              ind2 = 1_${ik}$ + n1
           else
              ind2 = n1 + n2
           end if
           i = 1_${ik}$
           ! while ( (n1sv > 0)
           10 continue
           if( n1sv>0_${ik}$ .and. n2sv>0_${ik}$ ) then
              if( a( ind1 )<=a( ind2 ) ) then
                 index( i ) = ind1
                 i = i + 1_${ik}$
                 ind1 = ind1 + dtrd1
                 n1sv = n1sv - 1_${ik}$
              else
                 index( i ) = ind2
                 i = i + 1_${ik}$
                 ind2 = ind2 + dtrd2
                 n2sv = n2sv - 1_${ik}$
              end if
              go to 10
           end if
           ! end while
           if( n1sv==0_${ik}$ ) then
              do n1sv = 1, n2sv
                 index( i ) = ind2
                 i = i + 1_${ik}$
                 ind2 = ind2 + dtrd2
              end do
           else
           ! n2sv == 0
              do n2sv = 1, n1sv
                 index( i ) = ind1
                 i = i + 1_${ik}$
                 ind1 = ind1 + dtrd1
              end do
           end if
           return
     end subroutine stdlib${ii}$_${ri}$lamrg

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slaeda( n, tlvls, curlvl, curpbm, prmptr, perm, givptr,givcol, givnum,&
     !! SLAEDA computes the Z vector corresponding to the merge step in the
     !! CURLVLth step of the merge process with TLVLS steps for the CURPBMth
     !! problem.
                q, qptr, z, ztemp, 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 
           integer(${ik}$), intent(in) :: curlvl, curpbm, n, tlvls
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           integer(${ik}$), intent(in) :: givcol(2_${ik}$,*), givptr(*), perm(*), prmptr(*), qptr(*)
           real(sp), intent(in) :: givnum(2_${ik}$,*), q(*)
           real(sp), intent(out) :: z(*), ztemp(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: bsiz1, bsiz2, curr, i, k, mid, psiz1, psiz2, ptr, zptr1
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( n<0_${ik}$ ) then
              info = -1_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SLAEDA', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! determine location of first number in second half.
           mid = n / 2_${ik}$ + 1_${ik}$
           ! gather last/first rows of appropriate eigenblocks into center of z
           ptr = 1_${ik}$
           ! determine location of lowest level subproblem in the full storage
           ! scheme
           curr = ptr + curpbm*2_${ik}$**curlvl + 2_${ik}$**( curlvl-1 ) - 1_${ik}$
           ! determine size of these matrices.  we add half to the value of
           ! the sqrt in case the machine underestimates one of these square
           ! roots.
           bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=sp) ),KIND=${ik}$)
           bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=sp) ),KIND=${ik}$)
                     
           do k = 1, mid - bsiz1 - 1
              z( k ) = zero
           end do
           call stdlib${ii}$_scopy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,z( mid-bsiz1 ), 1_${ik}$ )
           call stdlib${ii}$_scopy( bsiz2, q( qptr( curr+1 ) ), bsiz2, z( mid ), 1_${ik}$ )
           do k = mid + bsiz2, n
              z( k ) = zero
           end do
           ! loop through remaining levels 1 -> curlvl applying the givens
           ! rotations and permutation and then multiplying the center matrices
           ! against the current z.
           ptr = 2_${ik}$**tlvls + 1_${ik}$
           loop_70: do k = 1, curlvl - 1
              curr = ptr + curpbm*2_${ik}$**( curlvl-k ) + 2_${ik}$**( curlvl-k-1 ) - 1_${ik}$
              psiz1 = prmptr( curr+1 ) - prmptr( curr )
              psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
              zptr1 = mid - psiz1
             ! apply givens at curr and curr+1
              do i = givptr( curr ), givptr( curr+1 ) - 1
                 call stdlib${ii}$_srot( 1_${ik}$, z( zptr1+givcol( 1_${ik}$, i )-1_${ik}$ ), 1_${ik}$,z( zptr1+givcol( 2_${ik}$, i )-1_${ik}$ ), &
                           1_${ik}$, givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
              end do
              do i = givptr( curr+1 ), givptr( curr+2 ) - 1
                 call stdlib${ii}$_srot( 1_${ik}$, z( mid-1+givcol( 1_${ik}$, i ) ), 1_${ik}$,z( mid-1+givcol( 2_${ik}$, i ) ), 1_${ik}$, &
                           givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
              end do
              psiz1 = prmptr( curr+1 ) - prmptr( curr )
              psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
              do i = 0, psiz1 - 1
                 ztemp( i+1 ) = z( zptr1+perm( prmptr( curr )+i )-1_${ik}$ )
              end do
              do i = 0, psiz2 - 1
                 ztemp( psiz1+i+1 ) = z( mid+perm( prmptr( curr+1 )+i )-1_${ik}$ )
              end do
              ! multiply blocks at curr and curr+1
              ! determine size of these matrices.  we add half to the value of
              ! the sqrt in case the machine underestimates one of these
              ! square roots.
              bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=sp) ),KIND=${ik}$)
                        
              bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=sp) ),KIND=${ik}$)
                        
              if( bsiz1>0_${ik}$ ) then
                 call stdlib${ii}$_sgemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),bsiz1, ztemp( 1_${ik}$ ), &
                           1_${ik}$, zero, z( zptr1 ), 1_${ik}$ )
              end if
              call stdlib${ii}$_scopy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1_${ik}$, z( zptr1+bsiz1 ),1_${ik}$ )
              if( bsiz2>0_${ik}$ ) then
                 call stdlib${ii}$_sgemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),bsiz2, ztemp( &
                           psiz1+1 ), 1_${ik}$, zero, z( mid ), 1_${ik}$ )
              end if
              call stdlib${ii}$_scopy( psiz2-bsiz2, ztemp( psiz1+bsiz2+1 ), 1_${ik}$,z( mid+bsiz2 ), 1_${ik}$ )
                        
              ptr = ptr + 2_${ik}$**( tlvls-k )
           end do loop_70
           return
     end subroutine stdlib${ii}$_slaeda

     pure module subroutine stdlib${ii}$_dlaeda( n, tlvls, curlvl, curpbm, prmptr, perm, givptr,givcol, givnum,&
     !! DLAEDA computes the Z vector corresponding to the merge step in the
     !! CURLVLth step of the merge process with TLVLS steps for the CURPBMth
     !! problem.
                q, qptr, z, ztemp, 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 
           integer(${ik}$), intent(in) :: curlvl, curpbm, n, tlvls
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           integer(${ik}$), intent(in) :: givcol(2_${ik}$,*), givptr(*), perm(*), prmptr(*), qptr(*)
           real(dp), intent(in) :: givnum(2_${ik}$,*), q(*)
           real(dp), intent(out) :: z(*), ztemp(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: bsiz1, bsiz2, curr, i, k, mid, psiz1, psiz2, ptr, zptr1
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( n<0_${ik}$ ) then
              info = -1_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DLAEDA', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! determine location of first number in second half.
           mid = n / 2_${ik}$ + 1_${ik}$
           ! gather last/first rows of appropriate eigenblocks into center of z
           ptr = 1_${ik}$
           ! determine location of lowest level subproblem in the full storage
           ! scheme
           curr = ptr + curpbm*2_${ik}$**curlvl + 2_${ik}$**( curlvl-1 ) - 1_${ik}$
           ! determine size of these matrices.  we add half to the value of
           ! the sqrt in case the machine underestimates one of these square
           ! roots.
           bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=dp) ),KIND=${ik}$)
           bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=dp) ),KIND=${ik}$)
                     
           do k = 1, mid - bsiz1 - 1
              z( k ) = zero
           end do
           call stdlib${ii}$_dcopy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,z( mid-bsiz1 ), 1_${ik}$ )
           call stdlib${ii}$_dcopy( bsiz2, q( qptr( curr+1 ) ), bsiz2, z( mid ), 1_${ik}$ )
           do k = mid + bsiz2, n
              z( k ) = zero
           end do
           ! loop through remaining levels 1 -> curlvl applying the givens
           ! rotations and permutation and then multiplying the center matrices
           ! against the current z.
           ptr = 2_${ik}$**tlvls + 1_${ik}$
           loop_70: do k = 1, curlvl - 1
              curr = ptr + curpbm*2_${ik}$**( curlvl-k ) + 2_${ik}$**( curlvl-k-1 ) - 1_${ik}$
              psiz1 = prmptr( curr+1 ) - prmptr( curr )
              psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
              zptr1 = mid - psiz1
             ! apply givens at curr and curr+1
              do i = givptr( curr ), givptr( curr+1 ) - 1
                 call stdlib${ii}$_drot( 1_${ik}$, z( zptr1+givcol( 1_${ik}$, i )-1_${ik}$ ), 1_${ik}$,z( zptr1+givcol( 2_${ik}$, i )-1_${ik}$ ), &
                           1_${ik}$, givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
              end do
              do i = givptr( curr+1 ), givptr( curr+2 ) - 1
                 call stdlib${ii}$_drot( 1_${ik}$, z( mid-1+givcol( 1_${ik}$, i ) ), 1_${ik}$,z( mid-1+givcol( 2_${ik}$, i ) ), 1_${ik}$, &
                           givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
              end do
              psiz1 = prmptr( curr+1 ) - prmptr( curr )
              psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
              do i = 0, psiz1 - 1
                 ztemp( i+1 ) = z( zptr1+perm( prmptr( curr )+i )-1_${ik}$ )
              end do
              do i = 0, psiz2 - 1
                 ztemp( psiz1+i+1 ) = z( mid+perm( prmptr( curr+1 )+i )-1_${ik}$ )
              end do
              ! multiply blocks at curr and curr+1
              ! determine size of these matrices.  we add half to the value of
              ! the sqrt in case the machine underestimates one of these
              ! square roots.
              bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=dp) ),KIND=${ik}$)
                        
              bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=dp) ),KIND=${ik}$)
                        
              if( bsiz1>0_${ik}$ ) then
                 call stdlib${ii}$_dgemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),bsiz1, ztemp( 1_${ik}$ ), &
                           1_${ik}$, zero, z( zptr1 ), 1_${ik}$ )
              end if
              call stdlib${ii}$_dcopy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1_${ik}$, z( zptr1+bsiz1 ),1_${ik}$ )
              if( bsiz2>0_${ik}$ ) then
                 call stdlib${ii}$_dgemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),bsiz2, ztemp( &
                           psiz1+1 ), 1_${ik}$, zero, z( mid ), 1_${ik}$ )
              end if
              call stdlib${ii}$_dcopy( psiz2-bsiz2, ztemp( psiz1+bsiz2+1 ), 1_${ik}$,z( mid+bsiz2 ), 1_${ik}$ )
                        
              ptr = ptr + 2_${ik}$**( tlvls-k )
           end do loop_70
           return
     end subroutine stdlib${ii}$_dlaeda

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$laeda( n, tlvls, curlvl, curpbm, prmptr, perm, givptr,givcol, givnum,&
     !! DLAEDA: computes the Z vector corresponding to the merge step in the
     !! CURLVLth step of the merge process with TLVLS steps for the CURPBMth
     !! problem.
                q, qptr, z, ztemp, 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 
           integer(${ik}$), intent(in) :: curlvl, curpbm, n, tlvls
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           integer(${ik}$), intent(in) :: givcol(2_${ik}$,*), givptr(*), perm(*), prmptr(*), qptr(*)
           real(${rk}$), intent(in) :: givnum(2_${ik}$,*), q(*)
           real(${rk}$), intent(out) :: z(*), ztemp(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: bsiz1, bsiz2, curr, i, k, mid, psiz1, psiz2, ptr, zptr1
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( n<0_${ik}$ ) then
              info = -1_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DLAEDA', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! determine location of first number in second half.
           mid = n / 2_${ik}$ + 1_${ik}$
           ! gather last/first rows of appropriate eigenblocks into center of z
           ptr = 1_${ik}$
           ! determine location of lowest level subproblem in the full storage
           ! scheme
           curr = ptr + curpbm*2_${ik}$**curlvl + 2_${ik}$**( curlvl-1 ) - 1_${ik}$
           ! determine size of these matrices.  we add half to the value of
           ! the sqrt in case the machine underestimates one of these square
           ! roots.
           bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=${rk}$) ),KIND=${ik}$)
           bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=${rk}$) ),KIND=${ik}$)
                     
           do k = 1, mid - bsiz1 - 1
              z( k ) = zero
           end do
           call stdlib${ii}$_${ri}$copy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,z( mid-bsiz1 ), 1_${ik}$ )
           call stdlib${ii}$_${ri}$copy( bsiz2, q( qptr( curr+1 ) ), bsiz2, z( mid ), 1_${ik}$ )
           do k = mid + bsiz2, n
              z( k ) = zero
           end do
           ! loop through remaining levels 1 -> curlvl applying the givens
           ! rotations and permutation and then multiplying the center matrices
           ! against the current z.
           ptr = 2_${ik}$**tlvls + 1_${ik}$
           loop_70: do k = 1, curlvl - 1
              curr = ptr + curpbm*2_${ik}$**( curlvl-k ) + 2_${ik}$**( curlvl-k-1 ) - 1_${ik}$
              psiz1 = prmptr( curr+1 ) - prmptr( curr )
              psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
              zptr1 = mid - psiz1
             ! apply givens at curr and curr+1
              do i = givptr( curr ), givptr( curr+1 ) - 1
                 call stdlib${ii}$_${ri}$rot( 1_${ik}$, z( zptr1+givcol( 1_${ik}$, i )-1_${ik}$ ), 1_${ik}$,z( zptr1+givcol( 2_${ik}$, i )-1_${ik}$ ), &
                           1_${ik}$, givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
              end do
              do i = givptr( curr+1 ), givptr( curr+2 ) - 1
                 call stdlib${ii}$_${ri}$rot( 1_${ik}$, z( mid-1+givcol( 1_${ik}$, i ) ), 1_${ik}$,z( mid-1+givcol( 2_${ik}$, i ) ), 1_${ik}$, &
                           givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
              end do
              psiz1 = prmptr( curr+1 ) - prmptr( curr )
              psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
              do i = 0, psiz1 - 1
                 ztemp( i+1 ) = z( zptr1+perm( prmptr( curr )+i )-1_${ik}$ )
              end do
              do i = 0, psiz2 - 1
                 ztemp( psiz1+i+1 ) = z( mid+perm( prmptr( curr+1 )+i )-1_${ik}$ )
              end do
              ! multiply blocks at curr and curr+1
              ! determine size of these matrices.  we add half to the value of
              ! the sqrt in case the machine underestimates one of these
              ! square roots.
              bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=${rk}$) ),KIND=${ik}$)
                        
              bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=${rk}$) ),KIND=${ik}$)
                        
              if( bsiz1>0_${ik}$ ) then
                 call stdlib${ii}$_${ri}$gemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),bsiz1, ztemp( 1_${ik}$ ), &
                           1_${ik}$, zero, z( zptr1 ), 1_${ik}$ )
              end if
              call stdlib${ii}$_${ri}$copy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1_${ik}$, z( zptr1+bsiz1 ),1_${ik}$ )
              if( bsiz2>0_${ik}$ ) then
                 call stdlib${ii}$_${ri}$gemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),bsiz2, ztemp( &
                           psiz1+1 ), 1_${ik}$, zero, z( mid ), 1_${ik}$ )
              end if
              call stdlib${ii}$_${ri}$copy( psiz2-bsiz2, ztemp( psiz1+bsiz2+1 ), 1_${ik}$,z( mid+bsiz2 ), 1_${ik}$ )
                        
              ptr = ptr + 2_${ik}$**( tlvls-k )
           end do loop_70
           return
     end subroutine stdlib${ii}$_${ri}$laeda

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slarra( n, d, e, e2, spltol, tnrm,nsplit, isplit, info )
     !! Compute the splitting points with threshold SPLTOL.
     !! SLARRA sets any "small" off-diagonal elements to zero.
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info, nsplit
           integer(${ik}$), intent(in) :: n
           real(sp), intent(in) :: spltol, tnrm
           ! Array Arguments 
           integer(${ik}$), intent(out) :: isplit(*)
           real(sp), intent(in) :: d(*)
           real(sp), intent(inout) :: e(*), e2(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           real(sp) :: eabs, tmp1
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           ! compute splitting points
           nsplit = 1_${ik}$
           if(spltol<zero) then
              ! criterion based on absolute off-diagonal value
              tmp1 = abs(spltol)* tnrm
              do i = 1, n-1
                 eabs = abs( e(i) )
                 if( eabs <= tmp1) then
                    e(i) = zero
                    e2(i) = zero
                    isplit( nsplit ) = i
                    nsplit = nsplit + 1_${ik}$
                 end if
              end do
           else
              ! criterion that guarantees relative accuracy
              do i = 1, n-1
                 eabs = abs( e(i) )
                 if( eabs <= spltol * sqrt(abs(d(i)))*sqrt(abs(d(i+1))) )then
                    e(i) = zero
                    e2(i) = zero
                    isplit( nsplit ) = i
                    nsplit = nsplit + 1_${ik}$
                 end if
              end do
           endif
           isplit( nsplit ) = n
           return
     end subroutine stdlib${ii}$_slarra

     pure module subroutine stdlib${ii}$_dlarra( n, d, e, e2, spltol, tnrm,nsplit, isplit, info )
     !! Compute the splitting points with threshold SPLTOL.
     !! DLARRA sets any "small" off-diagonal elements to zero.
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info, nsplit
           integer(${ik}$), intent(in) :: n
           real(dp), intent(in) :: spltol, tnrm
           ! Array Arguments 
           integer(${ik}$), intent(out) :: isplit(*)
           real(dp), intent(in) :: d(*)
           real(dp), intent(inout) :: e(*), e2(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           real(dp) :: eabs, tmp1
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           ! compute splitting points
           nsplit = 1_${ik}$
           if(spltol<zero) then
              ! criterion based on absolute off-diagonal value
              tmp1 = abs(spltol)* tnrm
              do i = 1, n-1
                 eabs = abs( e(i) )
                 if( eabs <= tmp1) then
                    e(i) = zero
                    e2(i) = zero
                    isplit( nsplit ) = i
                    nsplit = nsplit + 1_${ik}$
                 end if
              end do
           else
              ! criterion that guarantees relative accuracy
              do i = 1, n-1
                 eabs = abs( e(i) )
                 if( eabs <= spltol * sqrt(abs(d(i)))*sqrt(abs(d(i+1))) )then
                    e(i) = zero
                    e2(i) = zero
                    isplit( nsplit ) = i
                    nsplit = nsplit + 1_${ik}$
                 end if
              end do
           endif
           isplit( nsplit ) = n
           return
     end subroutine stdlib${ii}$_dlarra

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$larra( n, d, e, e2, spltol, tnrm,nsplit, isplit, info )
     !! Compute the splitting points with threshold SPLTOL.
     !! DLARRA: sets any "small" off-diagonal elements to zero.
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info, nsplit
           integer(${ik}$), intent(in) :: n
           real(${rk}$), intent(in) :: spltol, tnrm
           ! Array Arguments 
           integer(${ik}$), intent(out) :: isplit(*)
           real(${rk}$), intent(in) :: d(*)
           real(${rk}$), intent(inout) :: e(*), e2(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           real(${rk}$) :: eabs, tmp1
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           ! compute splitting points
           nsplit = 1_${ik}$
           if(spltol<zero) then
              ! criterion based on absolute off-diagonal value
              tmp1 = abs(spltol)* tnrm
              do i = 1, n-1
                 eabs = abs( e(i) )
                 if( eabs <= tmp1) then
                    e(i) = zero
                    e2(i) = zero
                    isplit( nsplit ) = i
                    nsplit = nsplit + 1_${ik}$
                 end if
              end do
           else
              ! criterion that guarantees relative accuracy
              do i = 1, n-1
                 eabs = abs( e(i) )
                 if( eabs <= spltol * sqrt(abs(d(i)))*sqrt(abs(d(i+1))) )then
                    e(i) = zero
                    e2(i) = zero
                    isplit( nsplit ) = i
                    nsplit = nsplit + 1_${ik}$
                 end if
              end do
           endif
           isplit( nsplit ) = n
           return
     end subroutine stdlib${ii}$_${ri}$larra

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slarrb( n, d, lld, ifirst, ilast, rtol1,rtol2, offset, w, wgap, werr, &
     !! Given the relatively robust representation(RRR) L D L^T, SLARRB:
     !! does "limited" bisection to refine the eigenvalues of L D L^T,
     !! W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
     !! guesses for these eigenvalues are input in W, the corresponding estimate
     !! of the error in these guesses and their gaps are input in WERR
     !! and WGAP, respectively. During bisection, intervals
     !! [left, right] are maintained by storing their mid-points and
     !! semi-widths in the arrays W and WERR respectively.
               work, iwork,pivmin, spdiam, twist, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ifirst, ilast, n, offset, twist
           integer(${ik}$), intent(out) :: info
           real(sp), intent(in) :: pivmin, rtol1, rtol2, spdiam
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(in) :: d(*), lld(*)
           real(sp), intent(inout) :: w(*), werr(*), wgap(*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           integer(${ik}$) :: maxitr
           ! Local Scalars 
           integer(${ik}$) :: i, i1, ii, ip, iter, k, negcnt, next, nint, olnint, prev, r
           real(sp) :: back, cvrgd, gap, left, lgap, mid, mnwdth, rgap, right, tmp, width
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
           mnwdth = two * pivmin
           r = twist
           if((r<1_${ik}$).or.(r>n)) r = n
           ! initialize unconverged intervals in [ work(2*i-1), work(2*i) ].
           ! the sturm count, count( work(2*i-1) ) is arranged to be i-1, while
           ! count( work(2*i) ) is stored in iwork( 2*i ). the integer iwork( 2*i-1 )
           ! for an unconverged interval is set to the index of the next unconverged
           ! interval, and is -1 or 0 for a converged interval. thus a linked
           ! list of unconverged intervals is set up.
           i1 = ifirst
           ! the number of unconverged intervals
           nint = 0_${ik}$
           ! the last unconverged interval found
           prev = 0_${ik}$
           rgap = wgap( i1-offset )
           loop_75: do i = i1, ilast
              k = 2_${ik}$*i
              ii = i - offset
              left = w( ii ) - werr( ii )
              right = w( ii ) + werr( ii )
              lgap = rgap
              rgap = wgap( ii )
              gap = min( lgap, rgap )
              ! make sure that [left,right] contains the desired eigenvalue
              ! compute negcount from dstqds facto l+d+l+^t = l d l^t - left
              ! do while( negcnt(left)>i-1 )
              back = werr( ii )
              20 continue
              negcnt = stdlib${ii}$_slaneg( n, d, lld, left, pivmin, r )
              if( negcnt>i-1 ) then
                 left = left - back
                 back = two*back
                 go to 20
              end if
              ! do while( negcnt(right)<i )
              ! compute negcount from dstqds facto l+d+l+^t = l d l^t - right
              back = werr( ii )
              50 continue
              negcnt = stdlib${ii}$_slaneg( n, d, lld, right, pivmin, r )
               if( negcnt<i ) then
                  right = right + back
                  back = two*back
                  go to 50
               end if
              width = half*abs( left - right )
              tmp = max( abs( left ), abs( right ) )
              cvrgd = max(rtol1*gap,rtol2*tmp)
              if( width<=cvrgd .or. width<=mnwdth ) then
                 ! this interval has already converged and does not need refinement.
                 ! (note that the gaps might change through refining the
                  ! eigenvalues, however, they can only get bigger.)
                 ! remove it from the list.
                 iwork( k-1 ) = -1_${ik}$
                 ! make sure that i1 always points to the first unconverged interval
                 if((i==i1).and.(i<ilast)) i1 = i + 1_${ik}$
                 if((prev>=i1).and.(i<=ilast)) iwork( 2_${ik}$*prev-1 ) = i + 1_${ik}$
              else
                 ! unconverged interval found
                 prev = i
                 nint = nint + 1_${ik}$
                 iwork( k-1 ) = i + 1_${ik}$
                 iwork( k ) = negcnt
              end if
              work( k-1 ) = left
              work( k ) = right
           end do loop_75
           ! do while( nint>0 ), i.e. there are still unconverged intervals
           ! and while (iter<maxitr)
           iter = 0_${ik}$
           80 continue
           prev = i1 - 1_${ik}$
           i = i1
           olnint = nint
           loop_100: do ip = 1, olnint
              k = 2_${ik}$*i
              ii = i - offset
              rgap = wgap( ii )
              lgap = rgap
              if(ii>1_${ik}$) lgap = wgap( ii-1 )
              gap = min( lgap, rgap )
              next = iwork( k-1 )
              left = work( k-1 )
              right = work( k )
              mid = half*( left + right )
              ! semiwidth of interval
              width = right - mid
              tmp = max( abs( left ), abs( right ) )
              cvrgd = max(rtol1*gap,rtol2*tmp)
              if( ( width<=cvrgd ) .or. ( width<=mnwdth ).or.( iter==maxitr ) )then
                 ! reduce number of unconverged intervals
                 nint = nint - 1_${ik}$
                 ! mark interval as converged.
                 iwork( k-1 ) = 0_${ik}$
                 if( i1==i ) then
                    i1 = next
                 else
                    ! prev holds the last unconverged interval previously examined
                    if(prev>=i1) iwork( 2_${ik}$*prev-1 ) = next
                 end if
                 i = next
                 cycle loop_100
              end if
              prev = i
              ! perform one bisection step
              negcnt = stdlib${ii}$_slaneg( n, d, lld, mid, pivmin, r )
              if( negcnt<=i-1 ) then
                 work( k-1 ) = mid
              else
                 work( k ) = mid
              end if
              i = next
           end do loop_100
           iter = iter + 1_${ik}$
           ! do another loop if there are still unconverged intervals
           ! however, in the last iteration, all intervals are accepted
           ! since this is the best we can do.
           if( ( nint>0 ).and.(iter<=maxitr) ) go to 80
           ! at this point, all the intervals have converged
           do i = ifirst, ilast
              k = 2_${ik}$*i
              ii = i - offset
              ! all intervals marked by '0' have been refined.
              if( iwork( k-1 )==0_${ik}$ ) then
                 w( ii ) = half*( work( k-1 )+work( k ) )
                 werr( ii ) = work( k ) - w( ii )
              end if
           end do
           do i = ifirst+1, ilast
              k = 2_${ik}$*i
              ii = i - offset
              wgap( ii-1 ) = max( zero,w(ii) - werr (ii) - w( ii-1 ) - werr( ii-1 ))
           end do
           return
     end subroutine stdlib${ii}$_slarrb

     pure module subroutine stdlib${ii}$_dlarrb( n, d, lld, ifirst, ilast, rtol1,rtol2, offset, w, wgap, werr, &
     !! Given the relatively robust representation(RRR) L D L^T, DLARRB:
     !! does "limited" bisection to refine the eigenvalues of L D L^T,
     !! W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
     !! guesses for these eigenvalues are input in W, the corresponding estimate
     !! of the error in these guesses and their gaps are input in WERR
     !! and WGAP, respectively. During bisection, intervals
     !! [left, right] are maintained by storing their mid-points and
     !! semi-widths in the arrays W and WERR respectively.
               work, iwork,pivmin, spdiam, twist, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ifirst, ilast, n, offset, twist
           integer(${ik}$), intent(out) :: info
           real(dp), intent(in) :: pivmin, rtol1, rtol2, spdiam
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(in) :: d(*), lld(*)
           real(dp), intent(inout) :: w(*), werr(*), wgap(*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           
           integer(${ik}$) :: maxitr
           ! Local Scalars 
           integer(${ik}$) :: i, i1, ii, ip, iter, k, negcnt, next, nint, olnint, prev, r
           real(dp) :: back, cvrgd, gap, left, lgap, mid, mnwdth, rgap, right, tmp, width
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
           mnwdth = two * pivmin
           r = twist
           if((r<1_${ik}$).or.(r>n)) r = n
           ! initialize unconverged intervals in [ work(2*i-1), work(2*i) ].
           ! the sturm count, count( work(2*i-1) ) is arranged to be i-1, while
           ! count( work(2*i) ) is stored in iwork( 2*i ). the integer iwork( 2*i-1 )
           ! for an unconverged interval is set to the index of the next unconverged
           ! interval, and is -1 or 0 for a converged interval. thus a linked
           ! list of unconverged intervals is set up.
           i1 = ifirst
           ! the number of unconverged intervals
           nint = 0_${ik}$
           ! the last unconverged interval found
           prev = 0_${ik}$
           rgap = wgap( i1-offset )
           loop_75: do i = i1, ilast
              k = 2_${ik}$*i
              ii = i - offset
              left = w( ii ) - werr( ii )
              right = w( ii ) + werr( ii )
              lgap = rgap
              rgap = wgap( ii )
              gap = min( lgap, rgap )
              ! make sure that [left,right] contains the desired eigenvalue
              ! compute negcount from dstqds facto l+d+l+^t = l d l^t - left
              ! do while( negcnt(left)>i-1 )
              back = werr( ii )
              20 continue
              negcnt = stdlib${ii}$_dlaneg( n, d, lld, left, pivmin, r )
              if( negcnt>i-1 ) then
                 left = left - back
                 back = two*back
                 go to 20
              end if
              ! do while( negcnt(right)<i )
              ! compute negcount from dstqds facto l+d+l+^t = l d l^t - right
              back = werr( ii )
              50 continue
              negcnt = stdlib${ii}$_dlaneg( n, d, lld, right, pivmin, r )
               if( negcnt<i ) then
                  right = right + back
                  back = two*back
                  go to 50
               end if
              width = half*abs( left - right )
              tmp = max( abs( left ), abs( right ) )
              cvrgd = max(rtol1*gap,rtol2*tmp)
              if( width<=cvrgd .or. width<=mnwdth ) then
                 ! this interval has already converged and does not need refinement.
                 ! (note that the gaps might change through refining the
                  ! eigenvalues, however, they can only get bigger.)
                 ! remove it from the list.
                 iwork( k-1 ) = -1_${ik}$
                 ! make sure that i1 always points to the first unconverged interval
                 if((i==i1).and.(i<ilast)) i1 = i + 1_${ik}$
                 if((prev>=i1).and.(i<=ilast)) iwork( 2_${ik}$*prev-1 ) = i + 1_${ik}$
              else
                 ! unconverged interval found
                 prev = i
                 nint = nint + 1_${ik}$
                 iwork( k-1 ) = i + 1_${ik}$
                 iwork( k ) = negcnt
              end if
              work( k-1 ) = left
              work( k ) = right
           end do loop_75
           ! do while( nint>0 ), i.e. there are still unconverged intervals
           ! and while (iter<maxitr)
           iter = 0_${ik}$
           80 continue
           prev = i1 - 1_${ik}$
           i = i1
           olnint = nint
           loop_100: do ip = 1, olnint
              k = 2_${ik}$*i
              ii = i - offset
              rgap = wgap( ii )
              lgap = rgap
              if(ii>1_${ik}$) lgap = wgap( ii-1 )
              gap = min( lgap, rgap )
              next = iwork( k-1 )
              left = work( k-1 )
              right = work( k )
              mid = half*( left + right )
              ! semiwidth of interval
              width = right - mid
              tmp = max( abs( left ), abs( right ) )
              cvrgd = max(rtol1*gap,rtol2*tmp)
              if( ( width<=cvrgd ) .or. ( width<=mnwdth ).or.( iter==maxitr ) )then
                 ! reduce number of unconverged intervals
                 nint = nint - 1_${ik}$
                 ! mark interval as converged.
                 iwork( k-1 ) = 0_${ik}$
                 if( i1==i ) then
                    i1 = next
                 else
                    ! prev holds the last unconverged interval previously examined
                    if(prev>=i1) iwork( 2_${ik}$*prev-1 ) = next
                 end if
                 i = next
                 cycle loop_100
              end if
              prev = i
              ! perform one bisection step
              negcnt = stdlib${ii}$_dlaneg( n, d, lld, mid, pivmin, r )
              if( negcnt<=i-1 ) then
                 work( k-1 ) = mid
              else
                 work( k ) = mid
              end if
              i = next
           end do loop_100
           iter = iter + 1_${ik}$
           ! do another loop if there are still unconverged intervals
           ! however, in the last iteration, all intervals are accepted
           ! since this is the best we can do.
           if( ( nint>0 ).and.(iter<=maxitr) ) go to 80
           ! at this point, all the intervals have converged
           do i = ifirst, ilast
              k = 2_${ik}$*i
              ii = i - offset
              ! all intervals marked by '0' have been refined.
              if( iwork( k-1 )==0_${ik}$ ) then
                 w( ii ) = half*( work( k-1 )+work( k ) )
                 werr( ii ) = work( k ) - w( ii )
              end if
           end do
           do i = ifirst+1, ilast
              k = 2_${ik}$*i
              ii = i - offset
              wgap( ii-1 ) = max( zero,w(ii) - werr (ii) - w( ii-1 ) - werr( ii-1 ))
           end do
           return
     end subroutine stdlib${ii}$_dlarrb

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$larrb( n, d, lld, ifirst, ilast, rtol1,rtol2, offset, w, wgap, werr, &
     !! Given the relatively robust representation(RRR) L D L^T, DLARRB:
     !! does "limited" bisection to refine the eigenvalues of L D L^T,
     !! W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
     !! guesses for these eigenvalues are input in W, the corresponding estimate
     !! of the error in these guesses and their gaps are input in WERR
     !! and WGAP, respectively. During bisection, intervals
     !! [left, right] are maintained by storing their mid-points and
     !! semi-widths in the arrays W and WERR respectively.
               work, iwork,pivmin, spdiam, twist, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ifirst, ilast, n, offset, twist
           integer(${ik}$), intent(out) :: info
           real(${rk}$), intent(in) :: pivmin, rtol1, rtol2, spdiam
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(in) :: d(*), lld(*)
           real(${rk}$), intent(inout) :: w(*), werr(*), wgap(*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           
           integer(${ik}$) :: maxitr
           ! Local Scalars 
           integer(${ik}$) :: i, i1, ii, ip, iter, k, negcnt, next, nint, olnint, prev, r
           real(${rk}$) :: back, cvrgd, gap, left, lgap, mid, mnwdth, rgap, right, tmp, width
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
           mnwdth = two * pivmin
           r = twist
           if((r<1_${ik}$).or.(r>n)) r = n
           ! initialize unconverged intervals in [ work(2*i-1), work(2*i) ].
           ! the sturm count, count( work(2*i-1) ) is arranged to be i-1, while
           ! count( work(2*i) ) is stored in iwork( 2*i ). the integer iwork( 2*i-1 )
           ! for an unconverged interval is set to the index of the next unconverged
           ! interval, and is -1 or 0 for a converged interval. thus a linked
           ! list of unconverged intervals is set up.
           i1 = ifirst
           ! the number of unconverged intervals
           nint = 0_${ik}$
           ! the last unconverged interval found
           prev = 0_${ik}$
           rgap = wgap( i1-offset )
           loop_75: do i = i1, ilast
              k = 2_${ik}$*i
              ii = i - offset
              left = w( ii ) - werr( ii )
              right = w( ii ) + werr( ii )
              lgap = rgap
              rgap = wgap( ii )
              gap = min( lgap, rgap )
              ! make sure that [left,right] contains the desired eigenvalue
              ! compute negcount from dstqds facto l+d+l+^t = l d l^t - left
              ! do while( negcnt(left)>i-1 )
              back = werr( ii )
              20 continue
              negcnt = stdlib${ii}$_${ri}$laneg( n, d, lld, left, pivmin, r )
              if( negcnt>i-1 ) then
                 left = left - back
                 back = two*back
                 go to 20
              end if
              ! do while( negcnt(right)<i )
              ! compute negcount from dstqds facto l+d+l+^t = l d l^t - right
              back = werr( ii )
              50 continue
              negcnt = stdlib${ii}$_${ri}$laneg( n, d, lld, right, pivmin, r )
               if( negcnt<i ) then
                  right = right + back
                  back = two*back
                  go to 50
               end if
              width = half*abs( left - right )
              tmp = max( abs( left ), abs( right ) )
              cvrgd = max(rtol1*gap,rtol2*tmp)
              if( width<=cvrgd .or. width<=mnwdth ) then
                 ! this interval has already converged and does not need refinement.
                 ! (note that the gaps might change through refining the
                  ! eigenvalues, however, they can only get bigger.)
                 ! remove it from the list.
                 iwork( k-1 ) = -1_${ik}$
                 ! make sure that i1 always points to the first unconverged interval
                 if((i==i1).and.(i<ilast)) i1 = i + 1_${ik}$
                 if((prev>=i1).and.(i<=ilast)) iwork( 2_${ik}$*prev-1 ) = i + 1_${ik}$
              else
                 ! unconverged interval found
                 prev = i
                 nint = nint + 1_${ik}$
                 iwork( k-1 ) = i + 1_${ik}$
                 iwork( k ) = negcnt
              end if
              work( k-1 ) = left
              work( k ) = right
           end do loop_75
           ! do while( nint>0 ), i.e. there are still unconverged intervals
           ! and while (iter<maxitr)
           iter = 0_${ik}$
           80 continue
           prev = i1 - 1_${ik}$
           i = i1
           olnint = nint
           loop_100: do ip = 1, olnint
              k = 2_${ik}$*i
              ii = i - offset
              rgap = wgap( ii )
              lgap = rgap
              if(ii>1_${ik}$) lgap = wgap( ii-1 )
              gap = min( lgap, rgap )
              next = iwork( k-1 )
              left = work( k-1 )
              right = work( k )
              mid = half*( left + right )
              ! semiwidth of interval
              width = right - mid
              tmp = max( abs( left ), abs( right ) )
              cvrgd = max(rtol1*gap,rtol2*tmp)
              if( ( width<=cvrgd ) .or. ( width<=mnwdth ).or.( iter==maxitr ) )then
                 ! reduce number of unconverged intervals
                 nint = nint - 1_${ik}$
                 ! mark interval as converged.
                 iwork( k-1 ) = 0_${ik}$
                 if( i1==i ) then
                    i1 = next
                 else
                    ! prev holds the last unconverged interval previously examined
                    if(prev>=i1) iwork( 2_${ik}$*prev-1 ) = next
                 end if
                 i = next
                 cycle loop_100
              end if
              prev = i
              ! perform one bisection step
              negcnt = stdlib${ii}$_${ri}$laneg( n, d, lld, mid, pivmin, r )
              if( negcnt<=i-1 ) then
                 work( k-1 ) = mid
              else
                 work( k ) = mid
              end if
              i = next
           end do loop_100
           iter = iter + 1_${ik}$
           ! do another loop if there are still unconverged intervals
           ! however, in the last iteration, all intervals are accepted
           ! since this is the best we can do.
           if( ( nint>0 ).and.(iter<=maxitr) ) go to 80
           ! at this point, all the intervals have converged
           do i = ifirst, ilast
              k = 2_${ik}$*i
              ii = i - offset
              ! all intervals marked by '0' have been refined.
              if( iwork( k-1 )==0_${ik}$ ) then
                 w( ii ) = half*( work( k-1 )+work( k ) )
                 werr( ii ) = work( k ) - w( ii )
              end if
           end do
           do i = ifirst+1, ilast
              k = 2_${ik}$*i
              ii = i - offset
              wgap( ii-1 ) = max( zero,w(ii) - werr (ii) - w( ii-1 ) - werr( ii-1 ))
           end do
           return
     end subroutine stdlib${ii}$_${ri}$larrb

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slarrc( jobt, n, vl, vu, d, e, pivmin,eigcnt, lcnt, rcnt, info )
     !! Find the number of eigenvalues of the symmetric tridiagonal matrix T
     !! that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T
     !! if JOBT = 'L'.
               
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: jobt
           integer(${ik}$), intent(out) :: eigcnt, info, lcnt, rcnt
           integer(${ik}$), intent(in) :: n
           real(sp), intent(in) :: pivmin, vl, vu
           ! Array Arguments 
           real(sp), intent(in) :: d(*), e(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           logical(lk) :: matt
           real(sp) :: lpivot, rpivot, sl, su, tmp, tmp2
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           lcnt = 0_${ik}$
           rcnt = 0_${ik}$
           eigcnt = 0_${ik}$
           matt = stdlib_lsame( jobt, 'T' )
           if (matt) then
              ! sturm sequence count on t
              lpivot = d( 1_${ik}$ ) - vl
              rpivot = d( 1_${ik}$ ) - vu
              if( lpivot<=zero ) then
                 lcnt = lcnt + 1_${ik}$
              endif
              if( rpivot<=zero ) then
                 rcnt = rcnt + 1_${ik}$
              endif
              do i = 1, n-1
                 tmp = e(i)**2_${ik}$
                 lpivot = ( d( i+1 )-vl ) - tmp/lpivot
                 rpivot = ( d( i+1 )-vu ) - tmp/rpivot
                 if( lpivot<=zero ) then
                    lcnt = lcnt + 1_${ik}$
                 endif
                 if( rpivot<=zero ) then
                    rcnt = rcnt + 1_${ik}$
                 endif
              end do
           else
              ! sturm sequence count on l d l^t
              sl = -vl
              su = -vu
              do i = 1, n - 1
                 lpivot = d( i ) + sl
                 rpivot = d( i ) + su
                 if( lpivot<=zero ) then
                    lcnt = lcnt + 1_${ik}$
                 endif
                 if( rpivot<=zero ) then
                    rcnt = rcnt + 1_${ik}$
                 endif
                 tmp = e(i) * d(i) * e(i)
                 tmp2 = tmp / lpivot
                 if( tmp2==zero ) then
                    sl =  tmp - vl
                 else
                    sl = sl*tmp2 - vl
                 end if
                 tmp2 = tmp / rpivot
                 if( tmp2==zero ) then
                    su =  tmp - vu
                 else
                    su = su*tmp2 - vu
                 end if
              end do
              lpivot = d( n ) + sl
              rpivot = d( n ) + su
              if( lpivot<=zero ) then
                 lcnt = lcnt + 1_${ik}$
              endif
              if( rpivot<=zero ) then
                 rcnt = rcnt + 1_${ik}$
              endif
           endif
           eigcnt = rcnt - lcnt
           return
     end subroutine stdlib${ii}$_slarrc

     pure module subroutine stdlib${ii}$_dlarrc( jobt, n, vl, vu, d, e, pivmin,eigcnt, lcnt, rcnt, info )
     !! Find the number of eigenvalues of the symmetric tridiagonal matrix T
     !! that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T
     !! if JOBT = 'L'.
               
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: jobt
           integer(${ik}$), intent(out) :: eigcnt, info, lcnt, rcnt
           integer(${ik}$), intent(in) :: n
           real(dp), intent(in) :: pivmin, vl, vu
           ! Array Arguments 
           real(dp), intent(in) :: d(*), e(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           logical(lk) :: matt
           real(dp) :: lpivot, rpivot, sl, su, tmp, tmp2
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           lcnt = 0_${ik}$
           rcnt = 0_${ik}$
           eigcnt = 0_${ik}$
           matt = stdlib_lsame( jobt, 'T' )
           if (matt) then
              ! sturm sequence count on t
              lpivot = d( 1_${ik}$ ) - vl
              rpivot = d( 1_${ik}$ ) - vu
              if( lpivot<=zero ) then
                 lcnt = lcnt + 1_${ik}$
              endif
              if( rpivot<=zero ) then
                 rcnt = rcnt + 1_${ik}$
              endif
              do i = 1, n-1
                 tmp = e(i)**2_${ik}$
                 lpivot = ( d( i+1 )-vl ) - tmp/lpivot
                 rpivot = ( d( i+1 )-vu ) - tmp/rpivot
                 if( lpivot<=zero ) then
                    lcnt = lcnt + 1_${ik}$
                 endif
                 if( rpivot<=zero ) then
                    rcnt = rcnt + 1_${ik}$
                 endif
              end do
           else
              ! sturm sequence count on l d l^t
              sl = -vl
              su = -vu
              do i = 1, n - 1
                 lpivot = d( i ) + sl
                 rpivot = d( i ) + su
                 if( lpivot<=zero ) then
                    lcnt = lcnt + 1_${ik}$
                 endif
                 if( rpivot<=zero ) then
                    rcnt = rcnt + 1_${ik}$
                 endif
                 tmp = e(i) * d(i) * e(i)
                 tmp2 = tmp / lpivot
                 if( tmp2==zero ) then
                    sl =  tmp - vl
                 else
                    sl = sl*tmp2 - vl
                 end if
                 tmp2 = tmp / rpivot
                 if( tmp2==zero ) then
                    su =  tmp - vu
                 else
                    su = su*tmp2 - vu
                 end if
              end do
              lpivot = d( n ) + sl
              rpivot = d( n ) + su
              if( lpivot<=zero ) then
                 lcnt = lcnt + 1_${ik}$
              endif
              if( rpivot<=zero ) then
                 rcnt = rcnt + 1_${ik}$
              endif
           endif
           eigcnt = rcnt - lcnt
           return
     end subroutine stdlib${ii}$_dlarrc

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$larrc( jobt, n, vl, vu, d, e, pivmin,eigcnt, lcnt, rcnt, info )
     !! Find the number of eigenvalues of the symmetric tridiagonal matrix T
     !! that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T
     !! if JOBT = 'L'.
               
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: jobt
           integer(${ik}$), intent(out) :: eigcnt, info, lcnt, rcnt
           integer(${ik}$), intent(in) :: n
           real(${rk}$), intent(in) :: pivmin, vl, vu
           ! Array Arguments 
           real(${rk}$), intent(in) :: d(*), e(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           logical(lk) :: matt
           real(${rk}$) :: lpivot, rpivot, sl, su, tmp, tmp2
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           lcnt = 0_${ik}$
           rcnt = 0_${ik}$
           eigcnt = 0_${ik}$
           matt = stdlib_lsame( jobt, 'T' )
           if (matt) then
              ! sturm sequence count on t
              lpivot = d( 1_${ik}$ ) - vl
              rpivot = d( 1_${ik}$ ) - vu
              if( lpivot<=zero ) then
                 lcnt = lcnt + 1_${ik}$
              endif
              if( rpivot<=zero ) then
                 rcnt = rcnt + 1_${ik}$
              endif
              do i = 1, n-1
                 tmp = e(i)**2_${ik}$
                 lpivot = ( d( i+1 )-vl ) - tmp/lpivot
                 rpivot = ( d( i+1 )-vu ) - tmp/rpivot
                 if( lpivot<=zero ) then
                    lcnt = lcnt + 1_${ik}$
                 endif
                 if( rpivot<=zero ) then
                    rcnt = rcnt + 1_${ik}$
                 endif
              end do
           else
              ! sturm sequence count on l d l^t
              sl = -vl
              su = -vu
              do i = 1, n - 1
                 lpivot = d( i ) + sl
                 rpivot = d( i ) + su
                 if( lpivot<=zero ) then
                    lcnt = lcnt + 1_${ik}$
                 endif
                 if( rpivot<=zero ) then
                    rcnt = rcnt + 1_${ik}$
                 endif
                 tmp = e(i) * d(i) * e(i)
                 tmp2 = tmp / lpivot
                 if( tmp2==zero ) then
                    sl =  tmp - vl
                 else
                    sl = sl*tmp2 - vl
                 end if
                 tmp2 = tmp / rpivot
                 if( tmp2==zero ) then
                    su =  tmp - vu
                 else
                    su = su*tmp2 - vu
                 end if
              end do
              lpivot = d( n ) + sl
              rpivot = d( n ) + su
              if( lpivot<=zero ) then
                 lcnt = lcnt + 1_${ik}$
              endif
              if( rpivot<=zero ) then
                 rcnt = rcnt + 1_${ik}$
              endif
           endif
           eigcnt = rcnt - lcnt
           return
     end subroutine stdlib${ii}$_${ri}$larrc

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slarrd( range, order, n, vl, vu, il, iu, gers,reltol, d, e, e2, &
     !! SLARRD computes the eigenvalues of a symmetric tridiagonal
     !! matrix T to suitable accuracy. This is an auxiliary code to be
     !! called from SSTEMR.
     !! The user may ask for all eigenvalues, all eigenvalues
     !! in the half-open interval (VL, VU], or the IL-th through IU-th
     !! eigenvalues.
     !! To avoid overflow, the matrix must be scaled so that its
     !! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
     !! accuracy, it should not be much smaller than that.
     !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
     !! Matrix", Report CS41, Computer Science Dept., Stanford
     !! University, July 21, 1966.
               pivmin, nsplit, isplit,m, w, werr, wl, wu, iblock, indexw,work, iwork, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: order, range
           integer(${ik}$), intent(in) :: il, iu, n, nsplit
           integer(${ik}$), intent(out) :: info, m
           real(sp), intent(in) :: pivmin, reltol, vl, vu
           real(sp), intent(out) :: wl, wu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iblock(*), indexw(*), iwork(*)
           integer(${ik}$), intent(in) :: isplit(*)
           real(sp), intent(in) :: d(*), e(*), e2(*), gers(*)
           real(sp), intent(out) :: w(*), werr(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: fudge = two
           integer(${ik}$), parameter :: allrng = 1_${ik}$
           integer(${ik}$), parameter :: valrng = 2_${ik}$
           integer(${ik}$), parameter :: indrng = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: ncnvrg, toofew
           integer(${ik}$) :: i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iout, &
                     irange, itmax, itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb, nwl, nwu
           real(sp) :: atoli, eps, gl, gu, rtoli, tmp1, tmp2, tnorm, uflow, wkill, wlu, &
                     wul
           ! Local Arrays 
           integer(${ik}$) :: idumma(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           ! decode range
           if( stdlib_lsame( range, 'A' ) ) then
              irange = allrng
           else if( stdlib_lsame( range, 'V' ) ) then
              irange = valrng
           else if( stdlib_lsame( range, 'I' ) ) then
              irange = indrng
           else
              irange = 0_${ik}$
           end if
           ! check for errors
           if( irange<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.(stdlib_lsame(order,'B').or.stdlib_lsame(order,'E')) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( irange==valrng ) then
              if( vl>=vu )info = -5_${ik}$
           else if( irange==indrng .and.( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) ) then
              info = -6_${ik}$
           else if( irange==indrng .and.( iu<min( n, il ) .or. iu>n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              return
           end if
           ! initialize error flags
           info = 0_${ik}$
           ncnvrg = .false.
           toofew = .false.
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 ) return
           ! simplification:
           if( irange==indrng .and. il==1_${ik}$ .and. iu==n ) irange = 1_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_slamch( 'P' )
           uflow = stdlib${ii}$_slamch( 'U' )
           ! special case when n=1
           ! treat case of 1x1 matrix for quick return
           if( n==1_${ik}$ ) then
              if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
                        irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
                 m = 1_${ik}$
                 w(1_${ik}$) = d(1_${ik}$)
                 ! the computation error of the eigenvalue is zero
                 werr(1_${ik}$) = zero
                 iblock( 1_${ik}$ ) = 1_${ik}$
                 indexw( 1_${ik}$ ) = 1_${ik}$
              endif
              return
           end if
           ! nb is the minimum vector length for vector bisection, or 0
           ! if only scalar is to be done.
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ ) nb = 0_${ik}$
           ! find global spectral radius
           gl = d(1_${ik}$)
           gu = d(1_${ik}$)
           do i = 1,n
              gl =  min( gl, gers( 2_${ik}$*i - 1_${ik}$))
              gu = max( gu, gers(2_${ik}$*i) )
           end do
           ! compute global gerschgorin bounds and spectral diameter
           tnorm = max( abs( gl ), abs( gu ) )
           gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
           gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
           ! [jan/28/2009] remove the line below since spdiam variable not use
           ! spdiam = gu - gl
           ! input arguments for stdlib${ii}$_slaebz:
           ! the relative tolerance.  an interval (a,b] lies within
           ! "relative tolerance" if  b-a < reltol*max(|a|,|b|),
           rtoli = reltol
           ! set the absolute tolerance for interval convergence to zero to force
           ! interval convergence based on relative size of the interval.
           ! this is dangerous because intervals might not converge when reltol is
           ! small. but at least a very small number should be selected so that for
           ! strongly graded matrices, the code can get relatively accurate
           ! eigenvalues.
           atoli = fudge*two*uflow + fudge*two*pivmin
           if( irange==indrng ) then
              ! range='i': compute an interval containing eigenvalues
              ! il through iu. the initial interval [gl,gu] from the global
              ! gerschgorin bounds gl and gu is refined by stdlib${ii}$_slaebz.
              itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
              work( n+1 ) = gl
              work( n+2 ) = gl
              work( n+3 ) = gu
              work( n+4 ) = gu
              work( n+5 ) = gl
              work( n+6 ) = gu
              iwork( 1_${ik}$ ) = -1_${ik}$
              iwork( 2_${ik}$ ) = -1_${ik}$
              iwork( 3_${ik}$ ) = n + 1_${ik}$
              iwork( 4_${ik}$ ) = n + 1_${ik}$
              iwork( 5_${ik}$ ) = il - 1_${ik}$
              iwork( 6_${ik}$ ) = iu
              call stdlib${ii}$_slaebz( 3_${ik}$, itmax, n, 2_${ik}$, 2_${ik}$, nb, atoli, rtoli, pivmin,d, e, e2, iwork( 5_${ik}$ )&
                        , work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo )
              if( iinfo /= 0_${ik}$ ) then
                 info = iinfo
                 return
              end if
              ! on exit, output intervals may not be ordered by ascending negcount
              if( iwork( 6_${ik}$ )==iu ) then
                 wl = work( n+1 )
                 wlu = work( n+3 )
                 nwl = iwork( 1_${ik}$ )
                 wu = work( n+4 )
                 wul = work( n+2 )
                 nwu = iwork( 4_${ik}$ )
              else
                 wl = work( n+2 )
                 wlu = work( n+4 )
                 nwl = iwork( 2_${ik}$ )
                 wu = work( n+3 )
                 wul = work( n+1 )
                 nwu = iwork( 3_${ik}$ )
              end if
              ! on exit, the interval [wl, wlu] contains a value with negcount nwl,
              ! and [wul, wu] contains a value with negcount nwu.
              if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then
                 info = 4_${ik}$
                 return
              end if
           elseif( irange==valrng ) then
              wl = vl
              wu = vu
           elseif( irange==allrng ) then
              wl = gl
              wu = gu
           endif
           ! find eigenvalues -- loop over blocks and recompute nwl and nwu.
           ! nwl accumulates the number of eigenvalues .le. wl,
           ! nwu accumulates the number of eigenvalues .le. wu
           m = 0_${ik}$
           iend = 0_${ik}$
           info = 0_${ik}$
           nwl = 0_${ik}$
           nwu = 0_${ik}$
           loop_70: do jblk = 1, nsplit
              ioff = iend
              ibegin = ioff + 1_${ik}$
              iend = isplit( jblk )
              in = iend - ioff
              if( in==1_${ik}$ ) then
                 ! 1x1 block
                 if( wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$
                 if( wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$
                 if( irange==allrng .or.( wl<d( ibegin )-pivmin.and. wu>= d( ibegin )-pivmin ) ) &
                           then
                    m = m + 1_${ik}$
                    w( m ) = d( ibegin )
                    werr(m) = zero
                    ! the gap for a single block doesn't matter for the later
                    ! algorithm and is assigned an arbitrary large value
                    iblock( m ) = jblk
                    indexw( m ) = 1_${ik}$
                 end if
              ! disabled 2x2 case because of a failure on the following matrix
              ! range = 'i', il = iu = 4
                ! original tridiagonal, d = [
                 ! -0.150102010615740e+00_sp
                 ! -0.849897989384260e+00_sp
                 ! -0.128208148052635e-15_sp
                  ! 0.128257718286320e-15_sp
                ! ];
                ! e = [
                 ! -0.357171383266986e+00_sp
                 ! -0.180411241501588e-15_sp
                 ! -0.175152352710251e-15_sp
                ! ];
               ! else if( in==2 ) then
      ! *           2x2 block
                  ! disc = sqrt( (half*(d(ibegin)-d(iend)))**2 + e(ibegin)**2 )
                  ! tmp1 = half*(d(ibegin)+d(iend))
                  ! l1 = tmp1 - disc
                  ! if( wl>= l1-pivmin )
           ! $         nwl = nwl + 1
                  ! if( wu>= l1-pivmin )
           ! $         nwu = nwu + 1
                  ! if( irange==allrng .or. ( wl<l1-pivmin .and. wu>=
           ! $          l1-pivmin ) ) then
                     ! m = m + 1
                     ! w( m ) = l1
      ! *              the uncertainty of eigenvalues of a 2x2 matrix is very small
                     ! werr( m ) = eps * abs( w( m ) ) * two
                     ! iblock( m ) = jblk
                     ! indexw( m ) = 1
                  ! endif
                  ! l2 = tmp1 + disc
                  ! if( wl>= l2-pivmin )
           ! $         nwl = nwl + 1
                  ! if( wu>= l2-pivmin )
           ! $         nwu = nwu + 1
                  ! if( irange==allrng .or. ( wl<l2-pivmin .and. wu>=
           ! $          l2-pivmin ) ) then
                     ! m = m + 1
                     ! w( m ) = l2
      ! *              the uncertainty of eigenvalues of a 2x2 matrix is very small
                     ! werr( m ) = eps * abs( w( m ) ) * two
                     ! iblock( m ) = jblk
                     ! indexw( m ) = 2
                  ! endif
              else
                 ! general case - block of size in >= 2
                 ! compute local gerschgorin interval and use it as the initial
                 ! interval for stdlib${ii}$_slaebz
                 gu = d( ibegin )
                 gl = d( ibegin )
                 tmp1 = zero
                 do j = ibegin, iend
                    gl =  min( gl, gers( 2_${ik}$*j - 1_${ik}$))
                    gu = max( gu, gers(2_${ik}$*j) )
                 end do
                 ! [jan/28/2009]
                 ! change spdiam by tnorm in lines 2 and 3 thereafter
                 ! line 1: remove computation of spdiam (not useful anymore)
                 ! spdiam = gu - gl
                 ! gl = gl - fudge*spdiam*eps*in - fudge*pivmin
                 ! gu = gu + fudge*spdiam*eps*in + fudge*pivmin
                 gl = gl - fudge*tnorm*eps*in - fudge*pivmin
                 gu = gu + fudge*tnorm*eps*in + fudge*pivmin
                 if( irange>1_${ik}$ ) then
                    if( gu<wl ) then
                       ! the local block contains none of the wanted eigenvalues
                       nwl = nwl + in
                       nwu = nwu + in
                       cycle loop_70
                    end if
                    ! refine search interval if possible, only range (wl,wu] matters
                    gl = max( gl, wl )
                    gu = min( gu, wu )
                    if( gl>=gu )cycle loop_70
                 end if
                 ! find negcount of initial interval boundaries gl and gu
                 work( n+1 ) = gl
                 work( n+in+1 ) = gu
                 call stdlib${ii}$_slaebz( 1_${ik}$, 0_${ik}$, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e( &
                 ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 ),&
                            iblock( m+1 ), iinfo )
                 if( iinfo /= 0_${ik}$ ) then
                    info = iinfo
                    return
                 end if
                 nwl = nwl + iwork( 1_${ik}$ )
                 nwu = nwu + iwork( in+1 )
                 iwoff = m - iwork( 1_${ik}$ )
                 ! compute eigenvalues
                 itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + &
                           2_${ik}$
                 call stdlib${ii}$_slaebz( 2_${ik}$, itmax, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e(&
                  ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( m+&
                            1_${ik}$ ), iblock( m+1 ), iinfo )
                 if( iinfo /= 0_${ik}$ ) then
                    info = iinfo
                    return
                 end if
                 ! copy eigenvalues into w and iblock
                 ! use -jblk for block number for unconverged eigenvalues.
                 ! loop over the number of output intervals from stdlib${ii}$_slaebz
                 do j = 1, iout
                    ! eigenvalue approximation is middle point of interval
                    tmp1 = half*( work( j+n )+work( j+in+n ) )
                    ! semi length of error interval
                    tmp2 = half*abs( work( j+n )-work( j+in+n ) )
                    if( j>iout-iinfo ) then
                       ! flag non-convergence.
                       ncnvrg = .true.
                       ib = -jblk
                    else
                       ib = jblk
                    end if
                    do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff
                       w( je ) = tmp1
                       werr( je ) = tmp2
                       indexw( je ) = je - iwoff
                       iblock( je ) = ib
                    end do
                 end do
                 m = m + im
              end if
           end do loop_70
           ! if range='i', then (wl,wu) contains eigenvalues nwl+1,...,nwu
           ! if nwl+1 < il or nwu > iu, discard extra eigenvalues.
           if( irange==indrng ) then
              idiscl = il - 1_${ik}$ - nwl
              idiscu = nwu - iu
              if( idiscl>0_${ik}$ ) then
                 im = 0_${ik}$
                 do je = 1, m
                    ! remove some of the smallest eigenvalues from the left so that
                    ! at the end idiscl =0. move all eigenvalues up to the left.
                    if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then
                       idiscl = idiscl - 1_${ik}$
                    else
                       im = im + 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 m = im
              end if
              if( idiscu>0_${ik}$ ) then
                 ! remove some of the largest eigenvalues from the right so that
                 ! at the end idiscu =0. move all eigenvalues up to the left.
                 im=m+1
                 do je = m, 1, -1
                    if( w( je )>=wul .and. idiscu>0_${ik}$ ) then
                       idiscu = idiscu - 1_${ik}$
                    else
                       im = im - 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 jee = 0_${ik}$
                 do je = im, m
                    jee = jee + 1_${ik}$
                    w( jee ) = w( je )
                    werr( jee ) = werr( je )
                    indexw( jee ) = indexw( je )
                    iblock( jee ) = iblock( je )
                 end do
                 m = m-im+1
              end if
              if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
                 ! code to deal with effects of bad arithmetic. (if n(w) is
                 ! monotone non-decreasing, this should never happen.)
                 ! some low eigenvalues to be discarded are not in (wl,wlu],
                 ! or high eigenvalues to be discarded are not in (wul,wu]
                 ! so just kill off the smallest idiscl/largest idiscu
                 ! eigenvalues, by marking the corresponding iblock = 0
                 if( idiscl>0_${ik}$ ) then
                    wkill = wu
                    do jdisc = 1, idiscl
                       iw = 0_${ik}$
                       do je = 1, m
                          if( iblock( je )/=0_${ik}$ .and.( w( je )<wkill .or. iw==0_${ik}$ ) ) then
                             iw = je
                             wkill = w( je )
                          end if
                       end do
                       iblock( iw ) = 0_${ik}$
                    end do
                 end if
                 if( idiscu>0_${ik}$ ) then
                    wkill = wl
                    do jdisc = 1, idiscu
                       iw = 0_${ik}$
                       do je = 1, m
                          if( iblock( je )/=0_${ik}$ .and.( w( je )>=wkill .or. iw==0_${ik}$ ) ) then
                             iw = je
                             wkill = w( je )
                          end if
                       end do
                       iblock( iw ) = 0_${ik}$
                    end do
                 end if
                 ! now erase all eigenvalues with iblock set to zero
                 im = 0_${ik}$
                 do je = 1, m
                    if( iblock( je )/=0_${ik}$ ) then
                       im = im + 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 m = im
              end if
              if( idiscl<0_${ik}$ .or. idiscu<0_${ik}$ ) then
                 toofew = .true.
              end if
           end if
           if(( irange==allrng .and. m/=n ).or.( irange==indrng .and. m/=iu-il+1 ) ) then
              toofew = .true.
           end if
           ! if order='b', do nothing the eigenvalues are already sorted by
              ! block.
           ! if order='e', sort the eigenvalues from smallest to largest
           if( stdlib_lsame(order,'E') .and. nsplit>1_${ik}$ ) then
              do je = 1, m - 1
                 ie = 0_${ik}$
                 tmp1 = w( je )
                 do j = je + 1, m
                    if( w( j )<tmp1 ) then
                       ie = j
                       tmp1 = w( j )
                    end if
                 end do
                 if( ie/=0_${ik}$ ) then
                    tmp2 = werr( ie )
                    itmp1 = iblock( ie )
                    itmp2 = indexw( ie )
                    w( ie ) = w( je )
                    werr( ie ) = werr( je )
                    iblock( ie ) = iblock( je )
                    indexw( ie ) = indexw( je )
                    w( je ) = tmp1
                    werr( je ) = tmp2
                    iblock( je ) = itmp1
                    indexw( je ) = itmp2
                 end if
              end do
           end if
           info = 0_${ik}$
           if( ncnvrg )info = info + 1_${ik}$
           if( toofew )info = info + 2_${ik}$
           return
     end subroutine stdlib${ii}$_slarrd

     pure module subroutine stdlib${ii}$_dlarrd( range, order, n, vl, vu, il, iu, gers,reltol, d, e, e2, &
     !! DLARRD computes the eigenvalues of a symmetric tridiagonal
     !! matrix T to suitable accuracy. This is an auxiliary code to be
     !! called from DSTEMR.
     !! The user may ask for all eigenvalues, all eigenvalues
     !! in the half-open interval (VL, VU], or the IL-th through IU-th
     !! eigenvalues.
     !! To avoid overflow, the matrix must be scaled so that its
     !! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
     !! accuracy, it should not be much smaller than that.
     !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
     !! Matrix", Report CS41, Computer Science Dept., Stanford
     !! University, July 21, 1966.
               pivmin, nsplit, isplit,m, w, werr, wl, wu, iblock, indexw,work, iwork, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: order, range
           integer(${ik}$), intent(in) :: il, iu, n, nsplit
           integer(${ik}$), intent(out) :: info, m
           real(dp), intent(in) :: pivmin, reltol, vl, vu
           real(dp), intent(out) :: wl, wu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iblock(*), indexw(*), iwork(*)
           integer(${ik}$), intent(in) :: isplit(*)
           real(dp), intent(in) :: d(*), e(*), e2(*), gers(*)
           real(dp), intent(out) :: w(*), werr(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: fudge = two
           integer(${ik}$), parameter :: allrng = 1_${ik}$
           integer(${ik}$), parameter :: valrng = 2_${ik}$
           integer(${ik}$), parameter :: indrng = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: ncnvrg, toofew
           integer(${ik}$) :: i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iout, &
                     irange, itmax, itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb, nwl, nwu
           real(dp) :: atoli, eps, gl, gu, rtoli, tmp1, tmp2, tnorm, uflow, wkill, wlu, &
                     wul
           ! Local Arrays 
           integer(${ik}$) :: idumma(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           ! decode range
           if( stdlib_lsame( range, 'A' ) ) then
              irange = allrng
           else if( stdlib_lsame( range, 'V' ) ) then
              irange = valrng
           else if( stdlib_lsame( range, 'I' ) ) then
              irange = indrng
           else
              irange = 0_${ik}$
           end if
           ! check for errors
           if( irange<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.(stdlib_lsame(order,'B').or.stdlib_lsame(order,'E')) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( irange==valrng ) then
              if( vl>=vu )info = -5_${ik}$
           else if( irange==indrng .and.( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) ) then
              info = -6_${ik}$
           else if( irange==indrng .and.( iu<min( n, il ) .or. iu>n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              return
           end if
           ! initialize error flags
           info = 0_${ik}$
           ncnvrg = .false.
           toofew = .false.
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 ) return
           ! simplification:
           if( irange==indrng .and. il==1_${ik}$ .and. iu==n ) irange = 1_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_dlamch( 'P' )
           uflow = stdlib${ii}$_dlamch( 'U' )
           ! special case when n=1
           ! treat case of 1x1 matrix for quick return
           if( n==1_${ik}$ ) then
              if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
                        irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
                 m = 1_${ik}$
                 w(1_${ik}$) = d(1_${ik}$)
                 ! the computation error of the eigenvalue is zero
                 werr(1_${ik}$) = zero
                 iblock( 1_${ik}$ ) = 1_${ik}$
                 indexw( 1_${ik}$ ) = 1_${ik}$
              endif
              return
           end if
           ! nb is the minimum vector length for vector bisection, or 0
           ! if only scalar is to be done.
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ ) nb = 0_${ik}$
           ! find global spectral radius
           gl = d(1_${ik}$)
           gu = d(1_${ik}$)
           do i = 1,n
              gl =  min( gl, gers( 2_${ik}$*i - 1_${ik}$))
              gu = max( gu, gers(2_${ik}$*i) )
           end do
           ! compute global gerschgorin bounds and spectral diameter
           tnorm = max( abs( gl ), abs( gu ) )
           gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
           gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
           ! [jan/28/2009] remove the line below since spdiam variable not use
           ! spdiam = gu - gl
           ! input arguments for stdlib${ii}$_dlaebz:
           ! the relative tolerance.  an interval (a,b] lies within
           ! "relative tolerance" if  b-a < reltol*max(|a|,|b|),
           rtoli = reltol
           ! set the absolute tolerance for interval convergence to zero to force
           ! interval convergence based on relative size of the interval.
           ! this is dangerous because intervals might not converge when reltol is
           ! small. but at least a very small number should be selected so that for
           ! strongly graded matrices, the code can get relatively accurate
           ! eigenvalues.
           atoli = fudge*two*uflow + fudge*two*pivmin
           if( irange==indrng ) then
              ! range='i': compute an interval containing eigenvalues
              ! il through iu. the initial interval [gl,gu] from the global
              ! gerschgorin bounds gl and gu is refined by stdlib${ii}$_dlaebz.
              itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
              work( n+1 ) = gl
              work( n+2 ) = gl
              work( n+3 ) = gu
              work( n+4 ) = gu
              work( n+5 ) = gl
              work( n+6 ) = gu
              iwork( 1_${ik}$ ) = -1_${ik}$
              iwork( 2_${ik}$ ) = -1_${ik}$
              iwork( 3_${ik}$ ) = n + 1_${ik}$
              iwork( 4_${ik}$ ) = n + 1_${ik}$
              iwork( 5_${ik}$ ) = il - 1_${ik}$
              iwork( 6_${ik}$ ) = iu
              call stdlib${ii}$_dlaebz( 3_${ik}$, itmax, n, 2_${ik}$, 2_${ik}$, nb, atoli, rtoli, pivmin,d, e, e2, iwork( 5_${ik}$ )&
                        , work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo )
              if( iinfo /= 0_${ik}$ ) then
                 info = iinfo
                 return
              end if
              ! on exit, output intervals may not be ordered by ascending negcount
              if( iwork( 6_${ik}$ )==iu ) then
                 wl = work( n+1 )
                 wlu = work( n+3 )
                 nwl = iwork( 1_${ik}$ )
                 wu = work( n+4 )
                 wul = work( n+2 )
                 nwu = iwork( 4_${ik}$ )
              else
                 wl = work( n+2 )
                 wlu = work( n+4 )
                 nwl = iwork( 2_${ik}$ )
                 wu = work( n+3 )
                 wul = work( n+1 )
                 nwu = iwork( 3_${ik}$ )
              end if
              ! on exit, the interval [wl, wlu] contains a value with negcount nwl,
              ! and [wul, wu] contains a value with negcount nwu.
              if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then
                 info = 4_${ik}$
                 return
              end if
           elseif( irange==valrng ) then
              wl = vl
              wu = vu
           elseif( irange==allrng ) then
              wl = gl
              wu = gu
           endif
           ! find eigenvalues -- loop over blocks and recompute nwl and nwu.
           ! nwl accumulates the number of eigenvalues .le. wl,
           ! nwu accumulates the number of eigenvalues .le. wu
           m = 0_${ik}$
           iend = 0_${ik}$
           info = 0_${ik}$
           nwl = 0_${ik}$
           nwu = 0_${ik}$
           loop_70: do jblk = 1, nsplit
              ioff = iend
              ibegin = ioff + 1_${ik}$
              iend = isplit( jblk )
              in = iend - ioff
              if( in==1_${ik}$ ) then
                 ! 1x1 block
                 if( wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$
                 if( wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$
                 if( irange==allrng .or.( wl<d( ibegin )-pivmin.and. wu>= d( ibegin )-pivmin ) ) &
                           then
                    m = m + 1_${ik}$
                    w( m ) = d( ibegin )
                    werr(m) = zero
                    ! the gap for a single block doesn't matter for the later
                    ! algorithm and is assigned an arbitrary large value
                    iblock( m ) = jblk
                    indexw( m ) = 1_${ik}$
                 end if
              ! disabled 2x2 case because of a failure on the following matrix
              ! range = 'i', il = iu = 4
                ! original tridiagonal, d = [
                 ! -0.150102010615740e+00_dp
                 ! -0.849897989384260e+00_dp
                 ! -0.128208148052635e-15_dp
                  ! 0.128257718286320e-15_dp
                ! ];
                ! e = [
                 ! -0.357171383266986e+00_dp
                 ! -0.180411241501588e-15_dp
                 ! -0.175152352710251e-15_dp
                ! ];
               ! else if( in==2 ) then
      ! *           2x2 block
                  ! disc = sqrt( (half*(d(ibegin)-d(iend)))**2 + e(ibegin)**2 )
                  ! tmp1 = half*(d(ibegin)+d(iend))
                  ! l1 = tmp1 - disc
                  ! if( wl>= l1-pivmin )
           ! $         nwl = nwl + 1
                  ! if( wu>= l1-pivmin )
           ! $         nwu = nwu + 1
                  ! if( irange==allrng .or. ( wl<l1-pivmin .and. wu>=
           ! $          l1-pivmin ) ) then
                     ! m = m + 1
                     ! w( m ) = l1
      ! *              the uncertainty of eigenvalues of a 2x2 matrix is very small
                     ! werr( m ) = eps * abs( w( m ) ) * two
                     ! iblock( m ) = jblk
                     ! indexw( m ) = 1
                  ! endif
                  ! l2 = tmp1 + disc
                  ! if( wl>= l2-pivmin )
           ! $         nwl = nwl + 1
                  ! if( wu>= l2-pivmin )
           ! $         nwu = nwu + 1
                  ! if( irange==allrng .or. ( wl<l2-pivmin .and. wu>=
           ! $          l2-pivmin ) ) then
                     ! m = m + 1
                     ! w( m ) = l2
      ! *              the uncertainty of eigenvalues of a 2x2 matrix is very small
                     ! werr( m ) = eps * abs( w( m ) ) * two
                     ! iblock( m ) = jblk
                     ! indexw( m ) = 2
                  ! endif
              else
                 ! general case - block of size in >= 2
                 ! compute local gerschgorin interval and use it as the initial
                 ! interval for stdlib${ii}$_dlaebz
                 gu = d( ibegin )
                 gl = d( ibegin )
                 tmp1 = zero
                 do j = ibegin, iend
                    gl =  min( gl, gers( 2_${ik}$*j - 1_${ik}$))
                    gu = max( gu, gers(2_${ik}$*j) )
                 end do
                 ! [jan/28/2009]
                 ! change spdiam by tnorm in lines 2 and 3 thereafter
                 ! line 1: remove computation of spdiam (not useful anymore)
                 ! spdiam = gu - gl
                 ! gl = gl - fudge*spdiam*eps*in - fudge*pivmin
                 ! gu = gu + fudge*spdiam*eps*in + fudge*pivmin
                 gl = gl - fudge*tnorm*eps*in - fudge*pivmin
                 gu = gu + fudge*tnorm*eps*in + fudge*pivmin
                 if( irange>1_${ik}$ ) then
                    if( gu<wl ) then
                       ! the local block contains none of the wanted eigenvalues
                       nwl = nwl + in
                       nwu = nwu + in
                       cycle loop_70
                    end if
                    ! refine search interval if possible, only range (wl,wu] matters
                    gl = max( gl, wl )
                    gu = min( gu, wu )
                    if( gl>=gu )cycle loop_70
                 end if
                 ! find negcount of initial interval boundaries gl and gu
                 work( n+1 ) = gl
                 work( n+in+1 ) = gu
                 call stdlib${ii}$_dlaebz( 1_${ik}$, 0_${ik}$, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e( &
                 ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 ),&
                            iblock( m+1 ), iinfo )
                 if( iinfo /= 0_${ik}$ ) then
                    info = iinfo
                    return
                 end if
                 nwl = nwl + iwork( 1_${ik}$ )
                 nwu = nwu + iwork( in+1 )
                 iwoff = m - iwork( 1_${ik}$ )
                 ! compute eigenvalues
                 itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + &
                           2_${ik}$
                 call stdlib${ii}$_dlaebz( 2_${ik}$, itmax, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e(&
                  ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( m+&
                            1_${ik}$ ), iblock( m+1 ), iinfo )
                 if( iinfo /= 0_${ik}$ ) then
                    info = iinfo
                    return
                 end if
                 ! copy eigenvalues into w and iblock
                 ! use -jblk for block number for unconverged eigenvalues.
                 ! loop over the number of output intervals from stdlib${ii}$_dlaebz
                 do j = 1, iout
                    ! eigenvalue approximation is middle point of interval
                    tmp1 = half*( work( j+n )+work( j+in+n ) )
                    ! semi length of error interval
                    tmp2 = half*abs( work( j+n )-work( j+in+n ) )
                    if( j>iout-iinfo ) then
                       ! flag non-convergence.
                       ncnvrg = .true.
                       ib = -jblk
                    else
                       ib = jblk
                    end if
                    do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff
                       w( je ) = tmp1
                       werr( je ) = tmp2
                       indexw( je ) = je - iwoff
                       iblock( je ) = ib
                    end do
                 end do
                 m = m + im
              end if
           end do loop_70
           ! if range='i', then (wl,wu) contains eigenvalues nwl+1,...,nwu
           ! if nwl+1 < il or nwu > iu, discard extra eigenvalues.
           if( irange==indrng ) then
              idiscl = il - 1_${ik}$ - nwl
              idiscu = nwu - iu
              if( idiscl>0_${ik}$ ) then
                 im = 0_${ik}$
                 do je = 1, m
                    ! remove some of the smallest eigenvalues from the left so that
                    ! at the end idiscl =0. move all eigenvalues up to the left.
                    if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then
                       idiscl = idiscl - 1_${ik}$
                    else
                       im = im + 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 m = im
              end if
              if( idiscu>0_${ik}$ ) then
                 ! remove some of the largest eigenvalues from the right so that
                 ! at the end idiscu =0. move all eigenvalues up to the left.
                 im=m+1
                 do je = m, 1, -1
                    if( w( je )>=wul .and. idiscu>0_${ik}$ ) then
                       idiscu = idiscu - 1_${ik}$
                    else
                       im = im - 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 jee = 0_${ik}$
                 do je = im, m
                    jee = jee + 1_${ik}$
                    w( jee ) = w( je )
                    werr( jee ) = werr( je )
                    indexw( jee ) = indexw( je )
                    iblock( jee ) = iblock( je )
                 end do
                 m = m-im+1
              end if
              if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
                 ! code to deal with effects of bad arithmetic. (if n(w) is
                 ! monotone non-decreasing, this should never happen.)
                 ! some low eigenvalues to be discarded are not in (wl,wlu],
                 ! or high eigenvalues to be discarded are not in (wul,wu]
                 ! so just kill off the smallest idiscl/largest idiscu
                 ! eigenvalues, by marking the corresponding iblock = 0
                 if( idiscl>0_${ik}$ ) then
                    wkill = wu
                    do jdisc = 1, idiscl
                       iw = 0_${ik}$
                       do je = 1, m
                          if( iblock( je )/=0_${ik}$ .and.( w( je )<wkill .or. iw==0_${ik}$ ) ) then
                             iw = je
                             wkill = w( je )
                          end if
                       end do
                       iblock( iw ) = 0_${ik}$
                    end do
                 end if
                 if( idiscu>0_${ik}$ ) then
                    wkill = wl
                    do jdisc = 1, idiscu
                       iw = 0_${ik}$
                       do je = 1, m
                          if( iblock( je )/=0_${ik}$ .and.( w( je )>=wkill .or. iw==0_${ik}$ ) ) then
                             iw = je
                             wkill = w( je )
                          end if
                       end do
                       iblock( iw ) = 0_${ik}$
                    end do
                 end if
                 ! now erase all eigenvalues with iblock set to zero
                 im = 0_${ik}$
                 do je = 1, m
                    if( iblock( je )/=0_${ik}$ ) then
                       im = im + 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 m = im
              end if
              if( idiscl<0_${ik}$ .or. idiscu<0_${ik}$ ) then
                 toofew = .true.
              end if
           end if
           if(( irange==allrng .and. m/=n ).or.( irange==indrng .and. m/=iu-il+1 ) ) then
              toofew = .true.
           end if
           ! if order='b', do nothing the eigenvalues are already sorted by
              ! block.
           ! if order='e', sort the eigenvalues from smallest to largest
           if( stdlib_lsame(order,'E') .and. nsplit>1_${ik}$ ) then
              do je = 1, m - 1
                 ie = 0_${ik}$
                 tmp1 = w( je )
                 do j = je + 1, m
                    if( w( j )<tmp1 ) then
                       ie = j
                       tmp1 = w( j )
                    end if
                 end do
                 if( ie/=0_${ik}$ ) then
                    tmp2 = werr( ie )
                    itmp1 = iblock( ie )
                    itmp2 = indexw( ie )
                    w( ie ) = w( je )
                    werr( ie ) = werr( je )
                    iblock( ie ) = iblock( je )
                    indexw( ie ) = indexw( je )
                    w( je ) = tmp1
                    werr( je ) = tmp2
                    iblock( je ) = itmp1
                    indexw( je ) = itmp2
                 end if
              end do
           end if
           info = 0_${ik}$
           if( ncnvrg )info = info + 1_${ik}$
           if( toofew )info = info + 2_${ik}$
           return
     end subroutine stdlib${ii}$_dlarrd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$larrd( range, order, n, vl, vu, il, iu, gers,reltol, d, e, e2, &
     !! DLARRD: computes the eigenvalues of a symmetric tridiagonal
     !! matrix T to suitable accuracy. This is an auxiliary code to be
     !! called from DSTEMR.
     !! The user may ask for all eigenvalues, all eigenvalues
     !! in the half-open interval (VL, VU], or the IL-th through IU-th
     !! eigenvalues.
     !! To avoid overflow, the matrix must be scaled so that its
     !! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
     !! accuracy, it should not be much smaller than that.
     !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
     !! Matrix", Report CS41, Computer Science Dept., Stanford
     !! University, July 21, 1966.
               pivmin, nsplit, isplit,m, w, werr, wl, wu, iblock, indexw,work, iwork, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: order, range
           integer(${ik}$), intent(in) :: il, iu, n, nsplit
           integer(${ik}$), intent(out) :: info, m
           real(${rk}$), intent(in) :: pivmin, reltol, vl, vu
           real(${rk}$), intent(out) :: wl, wu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iblock(*), indexw(*), iwork(*)
           integer(${ik}$), intent(in) :: isplit(*)
           real(${rk}$), intent(in) :: d(*), e(*), e2(*), gers(*)
           real(${rk}$), intent(out) :: w(*), werr(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(${rk}$), parameter :: fudge = two
           integer(${ik}$), parameter :: allrng = 1_${ik}$
           integer(${ik}$), parameter :: valrng = 2_${ik}$
           integer(${ik}$), parameter :: indrng = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: ncnvrg, toofew
           integer(${ik}$) :: i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iout, &
                     irange, itmax, itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb, nwl, nwu
           real(${rk}$) :: atoli, eps, gl, gu, rtoli, tmp1, tmp2, tnorm, uflow, wkill, wlu, &
                     wul
           ! Local Arrays 
           integer(${ik}$) :: idumma(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           ! decode range
           if( stdlib_lsame( range, 'A' ) ) then
              irange = allrng
           else if( stdlib_lsame( range, 'V' ) ) then
              irange = valrng
           else if( stdlib_lsame( range, 'I' ) ) then
              irange = indrng
           else
              irange = 0_${ik}$
           end if
           ! check for errors
           if( irange<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.(stdlib_lsame(order,'B').or.stdlib_lsame(order,'E')) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( irange==valrng ) then
              if( vl>=vu )info = -5_${ik}$
           else if( irange==indrng .and.( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) ) then
              info = -6_${ik}$
           else if( irange==indrng .and.( iu<min( n, il ) .or. iu>n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              return
           end if
           ! initialize error flags
           info = 0_${ik}$
           ncnvrg = .false.
           toofew = .false.
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 ) return
           ! simplification:
           if( irange==indrng .and. il==1_${ik}$ .and. iu==n ) irange = 1_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_${ri}$lamch( 'P' )
           uflow = stdlib${ii}$_${ri}$lamch( 'U' )
           ! special case when n=1
           ! treat case of 1x1 matrix for quick return
           if( n==1_${ik}$ ) then
              if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
                        irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
                 m = 1_${ik}$
                 w(1_${ik}$) = d(1_${ik}$)
                 ! the computation error of the eigenvalue is zero
                 werr(1_${ik}$) = zero
                 iblock( 1_${ik}$ ) = 1_${ik}$
                 indexw( 1_${ik}$ ) = 1_${ik}$
              endif
              return
           end if
           ! nb is the minimum vector length for vector bisection, or 0
           ! if only scalar is to be done.
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ ) nb = 0_${ik}$
           ! find global spectral radius
           gl = d(1_${ik}$)
           gu = d(1_${ik}$)
           do i = 1,n
              gl =  min( gl, gers( 2_${ik}$*i - 1_${ik}$))
              gu = max( gu, gers(2_${ik}$*i) )
           end do
           ! compute global gerschgorin bounds and spectral diameter
           tnorm = max( abs( gl ), abs( gu ) )
           gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
           gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
           ! [jan/28/2009] remove the line below since spdiam variable not use
           ! spdiam = gu - gl
           ! input arguments for stdlib${ii}$_${ri}$laebz:
           ! the relative tolerance.  an interval (a,b] lies within
           ! "relative tolerance" if  b-a < reltol*max(|a|,|b|),
           rtoli = reltol
           ! set the absolute tolerance for interval convergence to zero to force
           ! interval convergence based on relative size of the interval.
           ! this is dangerous because intervals might not converge when reltol is
           ! small. but at least a very small number should be selected so that for
           ! strongly graded matrices, the code can get relatively accurate
           ! eigenvalues.
           atoli = fudge*two*uflow + fudge*two*pivmin
           if( irange==indrng ) then
              ! range='i': compute an interval containing eigenvalues
              ! il through iu. the initial interval [gl,gu] from the global
              ! gerschgorin bounds gl and gu is refined by stdlib${ii}$_${ri}$laebz.
              itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
              work( n+1 ) = gl
              work( n+2 ) = gl
              work( n+3 ) = gu
              work( n+4 ) = gu
              work( n+5 ) = gl
              work( n+6 ) = gu
              iwork( 1_${ik}$ ) = -1_${ik}$
              iwork( 2_${ik}$ ) = -1_${ik}$
              iwork( 3_${ik}$ ) = n + 1_${ik}$
              iwork( 4_${ik}$ ) = n + 1_${ik}$
              iwork( 5_${ik}$ ) = il - 1_${ik}$
              iwork( 6_${ik}$ ) = iu
              call stdlib${ii}$_${ri}$laebz( 3_${ik}$, itmax, n, 2_${ik}$, 2_${ik}$, nb, atoli, rtoli, pivmin,d, e, e2, iwork( 5_${ik}$ )&
                        , work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo )
              if( iinfo /= 0_${ik}$ ) then
                 info = iinfo
                 return
              end if
              ! on exit, output intervals may not be ordered by ascending negcount
              if( iwork( 6_${ik}$ )==iu ) then
                 wl = work( n+1 )
                 wlu = work( n+3 )
                 nwl = iwork( 1_${ik}$ )
                 wu = work( n+4 )
                 wul = work( n+2 )
                 nwu = iwork( 4_${ik}$ )
              else
                 wl = work( n+2 )
                 wlu = work( n+4 )
                 nwl = iwork( 2_${ik}$ )
                 wu = work( n+3 )
                 wul = work( n+1 )
                 nwu = iwork( 3_${ik}$ )
              end if
              ! on exit, the interval [wl, wlu] contains a value with negcount nwl,
              ! and [wul, wu] contains a value with negcount nwu.
              if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then
                 info = 4_${ik}$
                 return
              end if
           elseif( irange==valrng ) then
              wl = vl
              wu = vu
           elseif( irange==allrng ) then
              wl = gl
              wu = gu
           endif
           ! find eigenvalues -- loop over blocks and recompute nwl and nwu.
           ! nwl accumulates the number of eigenvalues .le. wl,
           ! nwu accumulates the number of eigenvalues .le. wu
           m = 0_${ik}$
           iend = 0_${ik}$
           info = 0_${ik}$
           nwl = 0_${ik}$
           nwu = 0_${ik}$
           loop_70: do jblk = 1, nsplit
              ioff = iend
              ibegin = ioff + 1_${ik}$
              iend = isplit( jblk )
              in = iend - ioff
              if( in==1_${ik}$ ) then
                 ! 1x1 block
                 if( wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$
                 if( wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$
                 if( irange==allrng .or.( wl<d( ibegin )-pivmin.and. wu>= d( ibegin )-pivmin ) ) &
                           then
                    m = m + 1_${ik}$
                    w( m ) = d( ibegin )
                    werr(m) = zero
                    ! the gap for a single block doesn't matter for the later
                    ! algorithm and is assigned an arbitrary large value
                    iblock( m ) = jblk
                    indexw( m ) = 1_${ik}$
                 end if
              ! disabled 2x2 case because of a failure on the following matrix
              ! range = 'i', il = iu = 4
                ! original tridiagonal, d = [
                 ! -0.150102010615740e+00_${rk}$
                 ! -0.849897989384260e+00_${rk}$
                 ! -0.128208148052635e-15_${rk}$
                  ! 0.128257718286320e-15_${rk}$
                ! ];
                ! e = [
                 ! -0.357171383266986e+00_${rk}$
                 ! -0.180411241501588e-15_${rk}$
                 ! -0.175152352710251e-15_${rk}$
                ! ];
               ! else if( in==2 ) then
      ! *           2x2 block
                  ! disc = sqrt( (half*(d(ibegin)-d(iend)))**2 + e(ibegin)**2 )
                  ! tmp1 = half*(d(ibegin)+d(iend))
                  ! l1 = tmp1 - disc
                  ! if( wl>= l1-pivmin )
           ! $         nwl = nwl + 1
                  ! if( wu>= l1-pivmin )
           ! $         nwu = nwu + 1
                  ! if( irange==allrng .or. ( wl<l1-pivmin .and. wu>=
           ! $          l1-pivmin ) ) then
                     ! m = m + 1
                     ! w( m ) = l1
      ! *              the uncertainty of eigenvalues of a 2x2 matrix is very small
                     ! werr( m ) = eps * abs( w( m ) ) * two
                     ! iblock( m ) = jblk
                     ! indexw( m ) = 1
                  ! endif
                  ! l2 = tmp1 + disc
                  ! if( wl>= l2-pivmin )
           ! $         nwl = nwl + 1
                  ! if( wu>= l2-pivmin )
           ! $         nwu = nwu + 1
                  ! if( irange==allrng .or. ( wl<l2-pivmin .and. wu>=
           ! $          l2-pivmin ) ) then
                     ! m = m + 1
                     ! w( m ) = l2
      ! *              the uncertainty of eigenvalues of a 2x2 matrix is very small
                     ! werr( m ) = eps * abs( w( m ) ) * two
                     ! iblock( m ) = jblk
                     ! indexw( m ) = 2
                  ! endif
              else
                 ! general case - block of size in >= 2
                 ! compute local gerschgorin interval and use it as the initial
                 ! interval for stdlib${ii}$_${ri}$laebz
                 gu = d( ibegin )
                 gl = d( ibegin )
                 tmp1 = zero
                 do j = ibegin, iend
                    gl =  min( gl, gers( 2_${ik}$*j - 1_${ik}$))
                    gu = max( gu, gers(2_${ik}$*j) )
                 end do
                 ! [jan/28/2009]
                 ! change spdiam by tnorm in lines 2 and 3 thereafter
                 ! line 1: remove computation of spdiam (not useful anymore)
                 ! spdiam = gu - gl
                 ! gl = gl - fudge*spdiam*eps*in - fudge*pivmin
                 ! gu = gu + fudge*spdiam*eps*in + fudge*pivmin
                 gl = gl - fudge*tnorm*eps*in - fudge*pivmin
                 gu = gu + fudge*tnorm*eps*in + fudge*pivmin
                 if( irange>1_${ik}$ ) then
                    if( gu<wl ) then
                       ! the local block contains none of the wanted eigenvalues
                       nwl = nwl + in
                       nwu = nwu + in
                       cycle loop_70
                    end if
                    ! refine search interval if possible, only range (wl,wu] matters
                    gl = max( gl, wl )
                    gu = min( gu, wu )
                    if( gl>=gu )cycle loop_70
                 end if
                 ! find negcount of initial interval boundaries gl and gu
                 work( n+1 ) = gl
                 work( n+in+1 ) = gu
                 call stdlib${ii}$_${ri}$laebz( 1_${ik}$, 0_${ik}$, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e( &
                 ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 ),&
                            iblock( m+1 ), iinfo )
                 if( iinfo /= 0_${ik}$ ) then
                    info = iinfo
                    return
                 end if
                 nwl = nwl + iwork( 1_${ik}$ )
                 nwu = nwu + iwork( in+1 )
                 iwoff = m - iwork( 1_${ik}$ )
                 ! compute eigenvalues
                 itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + &
                           2_${ik}$
                 call stdlib${ii}$_${ri}$laebz( 2_${ik}$, itmax, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e(&
                  ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( m+&
                            1_${ik}$ ), iblock( m+1 ), iinfo )
                 if( iinfo /= 0_${ik}$ ) then
                    info = iinfo
                    return
                 end if
                 ! copy eigenvalues into w and iblock
                 ! use -jblk for block number for unconverged eigenvalues.
                 ! loop over the number of output intervals from stdlib${ii}$_${ri}$laebz
                 do j = 1, iout
                    ! eigenvalue approximation is middle point of interval
                    tmp1 = half*( work( j+n )+work( j+in+n ) )
                    ! semi length of error interval
                    tmp2 = half*abs( work( j+n )-work( j+in+n ) )
                    if( j>iout-iinfo ) then
                       ! flag non-convergence.
                       ncnvrg = .true.
                       ib = -jblk
                    else
                       ib = jblk
                    end if
                    do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff
                       w( je ) = tmp1
                       werr( je ) = tmp2
                       indexw( je ) = je - iwoff
                       iblock( je ) = ib
                    end do
                 end do
                 m = m + im
              end if
           end do loop_70
           ! if range='i', then (wl,wu) contains eigenvalues nwl+1,...,nwu
           ! if nwl+1 < il or nwu > iu, discard extra eigenvalues.
           if( irange==indrng ) then
              idiscl = il - 1_${ik}$ - nwl
              idiscu = nwu - iu
              if( idiscl>0_${ik}$ ) then
                 im = 0_${ik}$
                 do je = 1, m
                    ! remove some of the smallest eigenvalues from the left so that
                    ! at the end idiscl =0. move all eigenvalues up to the left.
                    if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then
                       idiscl = idiscl - 1_${ik}$
                    else
                       im = im + 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 m = im
              end if
              if( idiscu>0_${ik}$ ) then
                 ! remove some of the largest eigenvalues from the right so that
                 ! at the end idiscu =0. move all eigenvalues up to the left.
                 im=m+1
                 do je = m, 1, -1
                    if( w( je )>=wul .and. idiscu>0_${ik}$ ) then
                       idiscu = idiscu - 1_${ik}$
                    else
                       im = im - 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 jee = 0_${ik}$
                 do je = im, m
                    jee = jee + 1_${ik}$
                    w( jee ) = w( je )
                    werr( jee ) = werr( je )
                    indexw( jee ) = indexw( je )
                    iblock( jee ) = iblock( je )
                 end do
                 m = m-im+1
              end if
              if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
                 ! code to deal with effects of bad arithmetic. (if n(w) is
                 ! monotone non-decreasing, this should never happen.)
                 ! some low eigenvalues to be discarded are not in (wl,wlu],
                 ! or high eigenvalues to be discarded are not in (wul,wu]
                 ! so just kill off the smallest idiscl/largest idiscu
                 ! eigenvalues, by marking the corresponding iblock = 0
                 if( idiscl>0_${ik}$ ) then
                    wkill = wu
                    do jdisc = 1, idiscl
                       iw = 0_${ik}$
                       do je = 1, m
                          if( iblock( je )/=0_${ik}$ .and.( w( je )<wkill .or. iw==0_${ik}$ ) ) then
                             iw = je
                             wkill = w( je )
                          end if
                       end do
                       iblock( iw ) = 0_${ik}$
                    end do
                 end if
                 if( idiscu>0_${ik}$ ) then
                    wkill = wl
                    do jdisc = 1, idiscu
                       iw = 0_${ik}$
                       do je = 1, m
                          if( iblock( je )/=0_${ik}$ .and.( w( je )>=wkill .or. iw==0_${ik}$ ) ) then
                             iw = je
                             wkill = w( je )
                          end if
                       end do
                       iblock( iw ) = 0_${ik}$
                    end do
                 end if
                 ! now erase all eigenvalues with iblock set to zero
                 im = 0_${ik}$
                 do je = 1, m
                    if( iblock( je )/=0_${ik}$ ) then
                       im = im + 1_${ik}$
                       w( im ) = w( je )
                       werr( im ) = werr( je )
                       indexw( im ) = indexw( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 m = im
              end if
              if( idiscl<0_${ik}$ .or. idiscu<0_${ik}$ ) then
                 toofew = .true.
              end if
           end if
           if(( irange==allrng .and. m/=n ).or.( irange==indrng .and. m/=iu-il+1 ) ) then
              toofew = .true.
           end if
           ! if order='b', do nothing the eigenvalues are already sorted by
              ! block.
           ! if order='e', sort the eigenvalues from smallest to largest
           if( stdlib_lsame(order,'E') .and. nsplit>1_${ik}$ ) then
              do je = 1, m - 1
                 ie = 0_${ik}$
                 tmp1 = w( je )
                 do j = je + 1, m
                    if( w( j )<tmp1 ) then
                       ie = j
                       tmp1 = w( j )
                    end if
                 end do
                 if( ie/=0_${ik}$ ) then
                    tmp2 = werr( ie )
                    itmp1 = iblock( ie )
                    itmp2 = indexw( ie )
                    w( ie ) = w( je )
                    werr( ie ) = werr( je )
                    iblock( ie ) = iblock( je )
                    indexw( ie ) = indexw( je )
                    w( je ) = tmp1
                    werr( je ) = tmp2
                    iblock( je ) = itmp1
                    indexw( je ) = itmp2
                 end if
              end do
           end if
           info = 0_${ik}$
           if( ncnvrg )info = info + 1_${ik}$
           if( toofew )info = info + 2_${ik}$
           return
     end subroutine stdlib${ii}$_${ri}$larrd

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slarre( range, n, vl, vu, il, iu, d, e, e2,rtol1, rtol2, spltol, &
     !! To find the desired eigenvalues of a given real symmetric
     !! tridiagonal matrix T, SLARRE: sets any "small" off-diagonal
     !! elements to zero, and for each unreduced block T_i, it finds
     !! (a) a suitable shift at one end of the block's spectrum,
     !! (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
     !! (c) eigenvalues of each L_i D_i L_i^T.
     !! The representations and eigenvalues found are then used by
     !! SSTEMR to compute the eigenvectors of T.
     !! The accuracy varies depending on whether bisection is used to
     !! find a few eigenvalues or the dqds algorithm (subroutine SLASQ2) to
     !! conpute all and then discard any unwanted one.
     !! As an added benefit, SLARRE also outputs the n
     !! Gerschgorin intervals for the matrices L_i D_i L_i^T.
               nsplit, isplit, m,w, werr, wgap, iblock, indexw, gers, pivmin,work, iwork, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: range
           integer(${ik}$), intent(in) :: il, iu, n
           integer(${ik}$), intent(out) :: info, m, nsplit
           real(sp), intent(out) :: pivmin
           real(sp), intent(in) :: rtol1, rtol2, spltol
           real(sp), intent(inout) :: vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iblock(*), isplit(*), iwork(*), indexw(*)
           real(sp), intent(inout) :: d(*), e(*), e2(*)
           real(sp), intent(out) :: gers(*), w(*), werr(*), wgap(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: hndrd = 100.0_sp
           real(sp), parameter :: pert = 4.0_sp
           real(sp), parameter :: fourth = one/four
           real(sp), parameter :: fac = half
           real(sp), parameter :: maxgrowth = 64.0_sp
           real(sp), parameter :: fudge = two
           integer(${ik}$), parameter :: maxtry = 6_${ik}$
           integer(${ik}$), parameter :: allrng = 1_${ik}$
           integer(${ik}$), parameter :: indrng = 2_${ik}$
           integer(${ik}$), parameter :: valrng = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: forceb, norep, usedqd
           integer(${ik}$) :: cnt, cnt1, cnt2, i, ibegin, idum, iend, iinfo, in, indl, indu, irange, &
                     j, jblk, mb, mm, wbegin, wend
           real(sp) :: avgap, bsrtol, clwdth, dmax, dpivot, eabs, emax, eold, eps, gl, gu, isleft,&
                      isrght, rtl, rtol, s1, s2, safmin, sgndef, sigma, spdiam, tau, tmp, tmp1
           ! Local Arrays 
           integer(${ik}$) :: iseed(4_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! quick return if possible
           if( n<=0_${ik}$ ) then
              return
           end if
           ! decode range
           if( stdlib_lsame( range, 'A' ) ) then
              irange = allrng
           else if( stdlib_lsame( range, 'V' ) ) then
              irange = valrng
           else if( stdlib_lsame( range, 'I' ) ) then
              irange = indrng
           end if
           m = 0_${ik}$
           ! get machine constants
           safmin = stdlib${ii}$_slamch( 'S' )
           eps = stdlib${ii}$_slamch( 'P' )
           ! set parameters
           rtl = hndrd*eps
           ! if one were ever to ask for less initial precision in bsrtol,
           ! one should keep in mind that for the subset case, the extremal
           ! eigenvalues must be at least as accurate as the current setting
           ! (eigenvalues in the middle need not as much accuracy)
           bsrtol = sqrt(eps)*(0.5e-3_sp)
           ! treat case of 1x1 matrix for quick return
           if( n==1_${ik}$ ) then
              if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
                        irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
                 m = 1_${ik}$
                 w(1_${ik}$) = d(1_${ik}$)
                 ! the computation error of the eigenvalue is zero
                 werr(1_${ik}$) = zero
                 wgap(1_${ik}$) = zero
                 iblock( 1_${ik}$ ) = 1_${ik}$
                 indexw( 1_${ik}$ ) = 1_${ik}$
                 gers(1_${ik}$) = d( 1_${ik}$ )
                 gers(2_${ik}$) = d( 1_${ik}$ )
              endif
              ! store the shift for the initial rrr, which is zero in this case
              e(1_${ik}$) = zero
              return
           end if
           ! general case: tridiagonal matrix of order > 1
           ! init werr, wgap. compute gerschgorin intervals and spectral diameter.
           ! compute maximum off-diagonal entry and pivmin.
           gl = d(1_${ik}$)
           gu = d(1_${ik}$)
           eold = zero
           emax = zero
           e(n) = zero
           do i = 1,n
              werr(i) = zero
              wgap(i) = zero
              eabs = abs( e(i) )
              if( eabs >= emax ) then
                 emax = eabs
              end if
              tmp1 = eabs + eold
              gers( 2_${ik}$*i-1) = d(i) - tmp1
              gl =  min( gl, gers( 2_${ik}$*i - 1_${ik}$))
              gers( 2_${ik}$*i ) = d(i) + tmp1
              gu = max( gu, gers(2_${ik}$*i) )
              eold  = eabs
           end do
           ! the minimum pivot allowed in the sturm sequence for t
           pivmin = safmin * max( one, emax**2_${ik}$ )
           ! compute spectral diameter. the gerschgorin bounds give an
           ! estimate that is wrong by at most a factor of sqrt(2)
           spdiam = gu - gl
           ! compute splitting points
           call stdlib${ii}$_slarra( n, d, e, e2, spltol, spdiam,nsplit, isplit, iinfo )
           ! can force use of bisection instead of faster dqds.
           ! option left in the code for future multisection work.
           forceb = .false.
           ! initialize usedqd, dqds should be used for allrng unless someone
           ! explicitly wants bisection.
           usedqd = (( irange==allrng ) .and. (.not.forceb))
           if( (irange==allrng) .and. (.not. forceb) ) then
              ! set interval [vl,vu] that contains all eigenvalues
              vl = gl
              vu = gu
           else
              ! we call stdlib${ii}$_slarrd to find crude approximations to the eigenvalues
              ! in the desired range. in case irange = indrng, we also obtain the
              ! interval (vl,vu] that contains all the wanted eigenvalues.
              ! an interval [left,right] has converged if
              ! right-left<rtol*max(abs(left),abs(right))
              ! stdlib${ii}$_slarrd needs a work of size 4*n, iwork of size 3*n
              call stdlib${ii}$_slarrd( range, 'B', n, vl, vu, il, iu, gers,bsrtol, d, e, e2, pivmin, &
                        nsplit, isplit,mm, w, werr, vl, vu, iblock, indexw,work, iwork, iinfo )
              if( iinfo/=0_${ik}$ ) then
                 info = -1_${ik}$
                 return
              endif
              ! make sure that the entries m+1 to n in w, werr, iblock, indexw are 0
              do i = mm+1,n
                 w( i ) = zero
                 werr( i ) = zero
                 iblock( i ) = 0_${ik}$
                 indexw( i ) = 0_${ik}$
              end do
           end if
      ! **
           ! loop over unreduced blocks
           ibegin = 1_${ik}$
           wbegin = 1_${ik}$
           loop_170: do jblk = 1, nsplit
              iend = isplit( jblk )
              in = iend - ibegin + 1