stdlib_lapack_eigv_tridiag3.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_sstev( jobz, n, d, e, z, ldz, work, info )
     !! SSTEV computes all eigenvalues and, optionally, eigenvectors of a
     !! real symmetric tridiagonal matrix A.
        ! -- lapack driver 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) :: jobz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(sp), intent(inout) :: d(*), e(*)
           real(sp), intent(out) :: work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: wantz
           integer(${ik}$) :: imax, iscale
           real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSTEV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = sqrt( bignum )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           tnrm = stdlib${ii}$_slanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
           end if
           ! for eigenvalues only, call stdlib${ii}$_ssterf.  for eigenvalues and
           ! eigenvectors, call stdlib${ii}$_ssteqr.
           if( .not.wantz ) then
              call stdlib${ii}$_ssterf( n, d, e, info )
           else
              call stdlib${ii}$_ssteqr( 'I', n, d, e, z, ldz, work, info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = n
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_sscal( imax, one / sigma, d, 1_${ik}$ )
           end if
           return
     end subroutine stdlib${ii}$_sstev

     pure module subroutine stdlib${ii}$_dstev( jobz, n, d, e, z, ldz, work, info )
     !! DSTEV computes all eigenvalues and, optionally, eigenvectors of a
     !! real symmetric tridiagonal matrix A.
        ! -- lapack driver 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) :: jobz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(dp), intent(inout) :: d(*), e(*)
           real(dp), intent(out) :: work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: wantz
           integer(${ik}$) :: imax, iscale
           real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_dlamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = sqrt( bignum )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           tnrm = stdlib${ii}$_dlanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
           end if
           ! for eigenvalues only, call stdlib${ii}$_dsterf.  for eigenvalues and
           ! eigenvectors, call stdlib${ii}$_dsteqr.
           if( .not.wantz ) then
              call stdlib${ii}$_dsterf( n, d, e, info )
           else
              call stdlib${ii}$_dsteqr( 'I', n, d, e, z, ldz, work, info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = n
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_dscal( imax, one / sigma, d, 1_${ik}$ )
           end if
           return
     end subroutine stdlib${ii}$_dstev

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$stev( jobz, n, d, e, z, ldz, work, info )
     !! DSTEV: computes all eigenvalues and, optionally, eigenvectors of a
     !! real symmetric tridiagonal matrix A.
        ! -- lapack driver 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) :: jobz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: d(*), e(*)
           real(${rk}$), intent(out) :: work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: wantz
           integer(${ik}$) :: imax, iscale
           real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = sqrt( bignum )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
           end if
           ! for eigenvalues only, call stdlib${ii}$_${ri}$sterf.  for eigenvalues and
           ! eigenvectors, call stdlib${ii}$_${ri}$steqr.
           if( .not.wantz ) then
              call stdlib${ii}$_${ri}$sterf( n, d, e, info )
           else
              call stdlib${ii}$_${ri}$steqr( 'I', n, d, e, z, ldz, work, info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = n
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_${ri}$scal( imax, one / sigma, d, 1_${ik}$ )
           end if
           return
     end subroutine stdlib${ii}$_${ri}$stev

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sstevd( jobz, n, d, e, z, ldz, work, lwork, iwork,liwork, info )
     !! SSTEVD computes all eigenvalues and, optionally, eigenvectors of a
     !! real symmetric tridiagonal matrix. If eigenvectors are desired, it
     !! uses a divide and conquer algorithm.
     !! The divide and conquer algorithm makes very mild assumptions about
     !! floating point arithmetic. It will work on machines with a guard
     !! digit in add/subtract, or on those binary machines without guard
     !! digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
     !! Cray-2. It could conceivably fail on hexadecimal or decimal machines
     !! without guard digits, but we know of none.
               
        ! -- lapack driver 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) :: jobz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: d(*), e(*)
           real(sp), intent(out) :: work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, wantz
           integer(${ik}$) :: iscale, liwmin, lwmin
           real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           liwmin = 1_${ik}$
           lwmin = 1_${ik}$
           if( n>1_${ik}$ .and. wantz ) then
              lwmin = 1_${ik}$ + 4_${ik}$*n + n**2_${ik}$
              liwmin = 3_${ik}$ + 5_${ik}$*n
           end if
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -6_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -8_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -10_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSTEVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = sqrt( bignum )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           tnrm = stdlib${ii}$_slanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
           end if
           ! for eigenvalues only, call stdlib${ii}$_ssterf.  for eigenvalues and
           ! eigenvectors, call stdlib${ii}$_sstedc.
           if( .not.wantz ) then
              call stdlib${ii}$_ssterf( n, d, e, info )
           else
              call stdlib${ii}$_sstedc( 'I', n, d, e, z, ldz, work, lwork, iwork, liwork,info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           if( iscale==1_${ik}$ )call stdlib${ii}$_sscal( n, one / sigma, d, 1_${ik}$ )
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_sstevd

     pure module subroutine stdlib${ii}$_dstevd( jobz, n, d, e, z, ldz, work, lwork, iwork,liwork, info )
     !! DSTEVD computes all eigenvalues and, optionally, eigenvectors of a
     !! real symmetric tridiagonal matrix. If eigenvectors are desired, it
     !! uses a divide and conquer algorithm.
     !! The divide and conquer algorithm makes very mild assumptions about
     !! floating point arithmetic. It will work on machines with a guard
     !! digit in add/subtract, or on those binary machines without guard
     !! digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
     !! Cray-2. It could conceivably fail on hexadecimal or decimal machines
     !! without guard digits, but we know of none.
               
        ! -- lapack driver 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) :: jobz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: d(*), e(*)
           real(dp), intent(out) :: work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, wantz
           integer(${ik}$) :: iscale, liwmin, lwmin
           real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           liwmin = 1_${ik}$
           lwmin = 1_${ik}$
           if( n>1_${ik}$ .and. wantz ) then
              lwmin = 1_${ik}$ + 4_${ik}$*n + n**2_${ik}$
              liwmin = 3_${ik}$ + 5_${ik}$*n
           end if
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -6_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -8_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -10_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_dlamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = sqrt( bignum )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           tnrm = stdlib${ii}$_dlanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
           end if
           ! for eigenvalues only, call stdlib${ii}$_dsterf.  for eigenvalues and
           ! eigenvectors, call stdlib${ii}$_dstedc.
           if( .not.wantz ) then
              call stdlib${ii}$_dsterf( n, d, e, info )
           else
              call stdlib${ii}$_dstedc( 'I', n, d, e, z, ldz, work, lwork, iwork, liwork,info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           if( iscale==1_${ik}$ )call stdlib${ii}$_dscal( n, one / sigma, d, 1_${ik}$ )
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_dstevd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$stevd( jobz, n, d, e, z, ldz, work, lwork, iwork,liwork, info )
     !! DSTEVD: computes all eigenvalues and, optionally, eigenvectors of a
     !! real symmetric tridiagonal matrix. If eigenvectors are desired, it
     !! uses a divide and conquer algorithm.
     !! The divide and conquer algorithm makes very mild assumptions about
     !! floating point arithmetic. It will work on machines with a guard
     !! digit in add/subtract, or on those binary machines without guard
     !! digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
     !! Cray-2. It could conceivably fail on hexadecimal or decimal machines
     !! without guard digits, but we know of none.
               
        ! -- lapack driver 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) :: jobz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(inout) :: d(*), e(*)
           real(${rk}$), intent(out) :: work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, wantz
           integer(${ik}$) :: iscale, liwmin, lwmin
           real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           liwmin = 1_${ik}$
           lwmin = 1_${ik}$
           if( n>1_${ik}$ .and. wantz ) then
              lwmin = 1_${ik}$ + 4_${ik}$*n + n**2_${ik}$
              liwmin = 3_${ik}$ + 5_${ik}$*n
           end if
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -6_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -8_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -10_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = sqrt( bignum )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
           end if
           ! for eigenvalues only, call stdlib${ii}$_${ri}$sterf.  for eigenvalues and
           ! eigenvectors, call stdlib${ii}$_${ri}$stedc.
           if( .not.wantz ) then
              call stdlib${ii}$_${ri}$sterf( n, d, e, info )
           else
              call stdlib${ii}$_${ri}$stedc( 'I', n, d, e, z, ldz, work, lwork, iwork, liwork,info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           if( iscale==1_${ik}$ )call stdlib${ii}$_${ri}$scal( n, one / sigma, d, 1_${ik}$ )
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_${ri}$stevd

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sstevr( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, &
     !! SSTEVR computes selected eigenvalues and, optionally, eigenvectors
     !! of a real symmetric tridiagonal matrix T.  Eigenvalues and
     !! eigenvectors can be selected by specifying either a range of values
     !! or a range of indices for the desired eigenvalues.
     !! Whenever possible, SSTEVR calls SSTEMR to compute the
     !! eigenspectrum using Relatively Robust Representations.  SSTEMR
     !! computes eigenvalues by the dqds algorithm, while orthogonal
     !! eigenvectors are computed from various "good" L D L^T representations
     !! (also known as Relatively Robust Representations). Gram-Schmidt
     !! orthogonalization is avoided as far as possible. More specifically,
     !! the various steps of the algorithm are as follows. For the i-th
     !! unreduced block of T,
     !! (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
     !! is a relatively robust representation,
     !! (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
     !! relative accuracy by the dqds algorithm,
     !! (c) If there is a cluster of close eigenvalues, "choose" sigma_i
     !! close to the cluster, and go to step (a),
     !! (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
     !! compute the corresponding eigenvector by forming a
     !! rank-revealing twisted factorization.
     !! The desired accuracy of the output can be specified by the input
     !! parameter ABSTOL.
     !! For more details, see "A new O(n^2) algorithm for the symmetric
     !! tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
     !! Computer Science Division Technical Report No. UCB//CSD-97-971,
     !! UC Berkeley, May 1997.
     !! Note 1 : SSTEVR calls SSTEMR when the full spectrum is requested
     !! on machines which conform to the ieee-754 floating point standard.
     !! SSTEVR calls SSTEBZ and SSTEIN on non-ieee machines and
     !! when partial spectrum requests are made.
     !! Normal execution of SSTEMR may create NaNs and infinities and
     !! hence may abort due to a floating point exception in environments
     !! which do not handle NaNs and infinities in the ieee standard default
     !! manner.
               isuppz, work, lwork, iwork,liwork, info )
        ! -- lapack driver 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) :: jobz, range
           integer(${ik}$), intent(in) :: il, iu, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(sp), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: isuppz(*), iwork(*)
           real(sp), intent(inout) :: d(*), e(*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, lquery, valeig, wantz, tryrac
           character :: order
           integer(${ik}$) :: i, ieeeok, imax, indibl, indifl, indisp, indiwo, iscale, j, jj, liwmin,&
                      lwmin, nsplit
           real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, &
                     vuu
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           ieeeok = stdlib${ii}$_ilaenv( 10_${ik}$, 'SSTEVR', 'N', 1_${ik}$, 2_${ik}$, 3_${ik}$, 4_${ik}$ )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           lquery = ( ( lwork==-1_${ik}$ ) .or. ( liwork==-1_${ik}$ ) )
           lwmin = max( 1_${ik}$, 20_${ik}$*n )
           liwmin = max(1_${ik}$, 10_${ik}$*n )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -8_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -9_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -14_${ik}$
              end if
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -17_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -19_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSTEVR', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( alleig .or. indeig ) then
                 m = 1_${ik}$
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
              else
                 if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then
                    m = 1_${ik}$
                    w( 1_${ik}$ ) = d( 1_${ik}$ )
                 end if
              end if
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           if( valeig ) then
              vll = vl
              vuu = vu
           end if
           tnrm = stdlib${ii}$_slanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
              if( valeig ) then
                 vll = vl*sigma
                 vuu = vu*sigma
              end if
           end if
           ! initialize indices into workspaces.  note: these indices are used only
           ! if stdlib${ii}$_ssterf or stdlib${ii}$_sstemr fail.
           ! iwork(indibl:indibl+m-1) corresponds to iblock in stdlib${ii}$_sstebz and
           ! stores the block indices of each of the m<=n eigenvalues.
           indibl = 1_${ik}$
           ! iwork(indisp:indisp+nsplit-1) corresponds to isplit in stdlib${ii}$_sstebz and
           ! stores the starting and finishing indices of each block.
           indisp = indibl + n
           ! iwork(indifl:indifl+n-1) stores the indices of eigenvectors
           ! that corresponding to eigenvectors that fail to converge in
           ! stdlib${ii}$_sstein.  this information is discarded; if any fail, the driver
           ! returns info > 0.
           indifl = indisp + n
           ! indiwo is the offset of the remaining integer workspace.
           indiwo = indisp + n
           ! if all eigenvalues are desired, then
           ! call stdlib${ii}$_ssterf or stdlib${ii}$_sstemr.  if this fails for some eigenvalue, then
           ! try stdlib${ii}$_sstebz.
           test = .false.
           if( indeig ) then
              if( il==1_${ik}$ .and. iu==n ) then
                 test = .true.
              end if
           end if
           if( ( alleig .or. test ) .and. ieeeok==1_${ik}$ ) then
              call stdlib${ii}$_scopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ )
              if( .not.wantz ) then
                 call stdlib${ii}$_scopy( n, d, 1_${ik}$, w, 1_${ik}$ )
                 call stdlib${ii}$_ssterf( n, w, work, info )
              else
                 call stdlib${ii}$_scopy( n, d, 1_${ik}$, work( n+1 ), 1_${ik}$ )
                 if (abstol <= two*n*eps) then
                    tryrac = .true.
                 else
                    tryrac = .false.
                 end if
                 call stdlib${ii}$_sstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,iu, m, w, z, ldz,&
                            n, isuppz, tryrac,work( 2_${ik}$*n+1 ), lwork-2*n, iwork, liwork, info )
              end if
              if( info==0_${ik}$ ) then
                 m = n
                 go to 10
              end if
              info = 0_${ik}$
           end if
           ! otherwise, call stdlib${ii}$_sstebz and, if eigenvectors are desired, stdlib${ii}$_sstein.
           if( wantz ) then
              order = 'B'
           else
              order = 'E'
           end if
           call stdlib${ii}$_sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, &
                     iwork( indibl ), iwork( indisp ), work,iwork( indiwo ), info )
           if( wantz ) then
              call stdlib${ii}$_sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work, &
                        iwork( indiwo ), iwork( indifl ),info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           10 continue
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = m
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_sscal( imax, one / sigma, w, 1_${ik}$ )
           end if
           ! if eigenvalues are not in order, then sort them, along with
           ! eigenvectors.
           if( wantz ) then
              do j = 1, m - 1
                 i = 0_${ik}$
                 tmp1 = w( j )
                 do jj = j + 1, m
                    if( w( jj )<tmp1 ) then
                       i = jj
                       tmp1 = w( jj )
                    end if
                 end do
                 if( i/=0_${ik}$ ) then
                    w( i ) = w( j )
                    w( j ) = tmp1
                    call stdlib${ii}$_sswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ )
                 end if
              end do
           end if
            ! causes problems with tests 19
            ! if (wantz .and. indeig ) z( 1,1) = z(1,1) / 1.002_sp + .002
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_sstevr

     pure module subroutine stdlib${ii}$_dstevr( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, &
     !! DSTEVR computes selected eigenvalues and, optionally, eigenvectors
     !! of a real symmetric tridiagonal matrix T.  Eigenvalues and
     !! eigenvectors can be selected by specifying either a range of values
     !! or a range of indices for the desired eigenvalues.
     !! Whenever possible, DSTEVR calls DSTEMR to compute the
     !! eigenspectrum using Relatively Robust Representations.  DSTEMR
     !! computes eigenvalues by the dqds algorithm, while orthogonal
     !! eigenvectors are computed from various "good" L D L^T representations
     !! (also known as Relatively Robust Representations). Gram-Schmidt
     !! orthogonalization is avoided as far as possible. More specifically,
     !! the various steps of the algorithm are as follows. For the i-th
     !! unreduced block of T,
     !! (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
     !! is a relatively robust representation,
     !! (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
     !! relative accuracy by the dqds algorithm,
     !! (c) If there is a cluster of close eigenvalues, "choose" sigma_i
     !! close to the cluster, and go to step (a),
     !! (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
     !! compute the corresponding eigenvector by forming a
     !! rank-revealing twisted factorization.
     !! The desired accuracy of the output can be specified by the input
     !! parameter ABSTOL.
     !! For more details, see "A new O(n^2) algorithm for the symmetric
     !! tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
     !! Computer Science Division Technical Report No. UCB//CSD-97-971,
     !! UC Berkeley, May 1997.
     !! Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
     !! on machines which conform to the ieee-754 floating point standard.
     !! DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
     !! when partial spectrum requests are made.
     !! Normal execution of DSTEMR may create NaNs and infinities and
     !! hence may abort due to a floating point exception in environments
     !! which do not handle NaNs and infinities in the ieee standard default
     !! manner.
               isuppz, work, lwork, iwork,liwork, info )
        ! -- lapack driver 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) :: jobz, range
           integer(${ik}$), intent(in) :: il, iu, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(dp), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: isuppz(*), iwork(*)
           real(dp), intent(inout) :: d(*), e(*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, lquery, valeig, wantz, tryrac
           character :: order
           integer(${ik}$) :: i, ieeeok, imax, indibl, indifl, indisp, indiwo, iscale, itmp1, j, jj, &
                     liwmin, lwmin, nsplit
           real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, &
                     vuu
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           ieeeok = stdlib${ii}$_ilaenv( 10_${ik}$, 'DSTEVR', 'N', 1_${ik}$, 2_${ik}$, 3_${ik}$, 4_${ik}$ )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           lquery = ( ( lwork==-1_${ik}$ ) .or. ( liwork==-1_${ik}$ ) )
           lwmin = max( 1_${ik}$, 20_${ik}$*n )
           liwmin = max( 1_${ik}$, 10_${ik}$*n )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -8_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -9_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -14_${ik}$
              end if
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -17_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -19_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEVR', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( alleig .or. indeig ) then
                 m = 1_${ik}$
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
              else
                 if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then
                    m = 1_${ik}$
                    w( 1_${ik}$ ) = d( 1_${ik}$ )
                 end if
              end if
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_dlamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           if( valeig ) then
              vll = vl
              vuu = vu
           end if
           tnrm = stdlib${ii}$_dlanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
              if( valeig ) then
                 vll = vl*sigma
                 vuu = vu*sigma
              end if
           end if
           ! initialize indices into workspaces.  note: these indices are used only
           ! if stdlib${ii}$_dsterf or stdlib${ii}$_dstemr fail.
           ! iwork(indibl:indibl+m-1) corresponds to iblock in stdlib${ii}$_dstebz and
           ! stores the block indices of each of the m<=n eigenvalues.
           indibl = 1_${ik}$
           ! iwork(indisp:indisp+nsplit-1) corresponds to isplit in stdlib${ii}$_dstebz and
           ! stores the starting and finishing indices of each block.
           indisp = indibl + n
           ! iwork(indifl:indifl+n-1) stores the indices of eigenvectors
           ! that corresponding to eigenvectors that fail to converge in
           ! stdlib${ii}$_dstein.  this information is discarded; if any fail, the driver
           ! returns info > 0.
           indifl = indisp + n
           ! indiwo is the offset of the remaining integer workspace.
           indiwo = indisp + n
           ! if all eigenvalues are desired, then
           ! call stdlib${ii}$_dsterf or stdlib${ii}$_dstemr.  if this fails for some eigenvalue, then
           ! try stdlib${ii}$_dstebz.
           test = .false.
           if( indeig ) then
              if( il==1_${ik}$ .and. iu==n ) then
                 test = .true.
              end if
           end if
           if( ( alleig .or. test ) .and. ieeeok==1_${ik}$ ) then
              call stdlib${ii}$_dcopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ )
              if( .not.wantz ) then
                 call stdlib${ii}$_dcopy( n, d, 1_${ik}$, w, 1_${ik}$ )
                 call stdlib${ii}$_dsterf( n, w, work, info )
              else
                 call stdlib${ii}$_dcopy( n, d, 1_${ik}$, work( n+1 ), 1_${ik}$ )
                 if (abstol <= two*n*eps) then
                    tryrac = .true.
                 else
                    tryrac = .false.
                 end if
                 call stdlib${ii}$_dstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,iu, m, w, z, ldz,&
                            n, isuppz, tryrac,work( 2_${ik}$*n+1 ), lwork-2*n, iwork, liwork, info )
              end if
              if( info==0_${ik}$ ) then
                 m = n
                 go to 10
              end if
              info = 0_${ik}$
           end if
           ! otherwise, call stdlib${ii}$_dstebz and, if eigenvectors are desired, stdlib${ii}$_dstein.
           if( wantz ) then
              order = 'B'
           else
              order = 'E'
           end if
           call stdlib${ii}$_dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, &
                     iwork( indibl ), iwork( indisp ), work,iwork( indiwo ), info )
           if( wantz ) then
              call stdlib${ii}$_dstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work, &
                        iwork( indiwo ), iwork( indifl ),info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           10 continue
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = m
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_dscal( imax, one / sigma, w, 1_${ik}$ )
           end if
           ! if eigenvalues are not in order, then sort them, along with
           ! eigenvectors.
           if( wantz ) then
              do j = 1, m - 1
                 i = 0_${ik}$
                 tmp1 = w( j )
                 do jj = j + 1, m
                    if( w( jj )<tmp1 ) then
                       i = jj
                       tmp1 = w( jj )
                    end if
                 end do
                 if( i/=0_${ik}$ ) then
                    itmp1 = iwork( i )
                    w( i ) = w( j )
                    iwork( i ) = iwork( j )
                    w( j ) = tmp1
                    iwork( j ) = itmp1
                    call stdlib${ii}$_dswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ )
                 end if
              end do
           end if
            ! causes problems with tests 19
            ! if (wantz .and. indeig ) z( 1,1) = z(1,1) / 1.002_dp + .002
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_dstevr

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$stevr( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, &
     !! DSTEVR: computes selected eigenvalues and, optionally, eigenvectors
     !! of a real symmetric tridiagonal matrix T.  Eigenvalues and
     !! eigenvectors can be selected by specifying either a range of values
     !! or a range of indices for the desired eigenvalues.
     !! Whenever possible, DSTEVR calls DSTEMR to compute the
     !! eigenspectrum using Relatively Robust Representations.  DSTEMR
     !! computes eigenvalues by the dqds algorithm, while orthogonal
     !! eigenvectors are computed from various "good" L D L^T representations
     !! (also known as Relatively Robust Representations). Gram-Schmidt
     !! orthogonalization is avoided as far as possible. More specifically,
     !! the various steps of the algorithm are as follows. For the i-th
     !! unreduced block of T,
     !! (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
     !! is a relatively robust representation,
     !! (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
     !! relative accuracy by the dqds algorithm,
     !! (c) If there is a cluster of close eigenvalues, "choose" sigma_i
     !! close to the cluster, and go to step (a),
     !! (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
     !! compute the corresponding eigenvector by forming a
     !! rank-revealing twisted factorization.
     !! The desired accuracy of the output can be specified by the input
     !! parameter ABSTOL.
     !! For more details, see "A new O(n^2) algorithm for the symmetric
     !! tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
     !! Computer Science Division Technical Report No. UCB//CSD-97-971,
     !! UC Berkeley, May 1997.
     !! Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
     !! on machines which conform to the ieee-754 floating point standard.
     !! DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
     !! when partial spectrum requests are made.
     !! Normal execution of DSTEMR may create NaNs and infinities and
     !! hence may abort due to a floating point exception in environments
     !! which do not handle NaNs and infinities in the ieee standard default
     !! manner.
               isuppz, work, lwork, iwork,liwork, info )
        ! -- lapack driver 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) :: jobz, range
           integer(${ik}$), intent(in) :: il, iu, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(${rk}$), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: isuppz(*), iwork(*)
           real(${rk}$), intent(inout) :: d(*), e(*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, lquery, valeig, wantz, tryrac
           character :: order
           integer(${ik}$) :: i, ieeeok, imax, indibl, indifl, indisp, indiwo, iscale, itmp1, j, jj, &
                     liwmin, lwmin, nsplit
           real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, &
                     vuu
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           ieeeok = stdlib${ii}$_ilaenv( 10_${ik}$, 'DSTEVR', 'N', 1_${ik}$, 2_${ik}$, 3_${ik}$, 4_${ik}$ )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           lquery = ( ( lwork==-1_${ik}$ ) .or. ( liwork==-1_${ik}$ ) )
           lwmin = max( 1_${ik}$, 20_${ik}$*n )
           liwmin = max( 1_${ik}$, 10_${ik}$*n )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -8_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -9_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -14_${ik}$
              end if
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -17_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -19_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEVR', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( alleig .or. indeig ) then
                 m = 1_${ik}$
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
              else
                 if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then
                    m = 1_${ik}$
                    w( 1_${ik}$ ) = d( 1_${ik}$ )
                 end if
              end if
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           if( valeig ) then
              vll = vl
              vuu = vu
           end if
           tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
              if( valeig ) then
                 vll = vl*sigma
                 vuu = vu*sigma
              end if
           end if
           ! initialize indices into workspaces.  note: these indices are used only
           ! if stdlib${ii}$_${ri}$sterf or stdlib${ii}$_${ri}$stemr fail.
           ! iwork(indibl:indibl+m-1) corresponds to iblock in stdlib${ii}$_${ri}$stebz and
           ! stores the block indices of each of the m<=n eigenvalues.
           indibl = 1_${ik}$
           ! iwork(indisp:indisp+nsplit-1) corresponds to isplit in stdlib${ii}$_${ri}$stebz and
           ! stores the starting and finishing indices of each block.
           indisp = indibl + n
           ! iwork(indifl:indifl+n-1) stores the indices of eigenvectors
           ! that corresponding to eigenvectors that fail to converge in
           ! stdlib${ii}$_${ri}$stein.  this information is discarded; if any fail, the driver
           ! returns info > 0.
           indifl = indisp + n
           ! indiwo is the offset of the remaining integer workspace.
           indiwo = indisp + n
           ! if all eigenvalues are desired, then
           ! call stdlib${ii}$_${ri}$sterf or stdlib${ii}$_${ri}$stemr.  if this fails for some eigenvalue, then
           ! try stdlib${ii}$_${ri}$stebz.
           test = .false.
           if( indeig ) then
              if( il==1_${ik}$ .and. iu==n ) then
                 test = .true.
              end if
           end if
           if( ( alleig .or. test ) .and. ieeeok==1_${ik}$ ) then
              call stdlib${ii}$_${ri}$copy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ )
              if( .not.wantz ) then
                 call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, w, 1_${ik}$ )
                 call stdlib${ii}$_${ri}$sterf( n, w, work, info )
              else
                 call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, work( n+1 ), 1_${ik}$ )
                 if (abstol <= two*n*eps) then
                    tryrac = .true.
                 else
                    tryrac = .false.
                 end if
                 call stdlib${ii}$_${ri}$stemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,iu, m, w, z, ldz,&
                            n, isuppz, tryrac,work( 2_${ik}$*n+1 ), lwork-2*n, iwork, liwork, info )
              end if
              if( info==0_${ik}$ ) then
                 m = n
                 go to 10
              end if
              info = 0_${ik}$
           end if
           ! otherwise, call stdlib${ii}$_${ri}$stebz and, if eigenvectors are desired, stdlib${ii}$_${ri}$stein.
           if( wantz ) then
              order = 'B'
           else
              order = 'E'
           end if
           call stdlib${ii}$_${ri}$stebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, &
                     iwork( indibl ), iwork( indisp ), work,iwork( indiwo ), info )
           if( wantz ) then
              call stdlib${ii}$_${ri}$stein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work, &
                        iwork( indiwo ), iwork( indifl ),info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           10 continue
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = m
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_${ri}$scal( imax, one / sigma, w, 1_${ik}$ )
           end if
           ! if eigenvalues are not in order, then sort them, along with
           ! eigenvectors.
           if( wantz ) then
              do j = 1, m - 1
                 i = 0_${ik}$
                 tmp1 = w( j )
                 do jj = j + 1, m
                    if( w( jj )<tmp1 ) then
                       i = jj
                       tmp1 = w( jj )
                    end if
                 end do
                 if( i/=0_${ik}$ ) then
                    itmp1 = iwork( i )
                    w( i ) = w( j )
                    iwork( i ) = iwork( j )
                    w( j ) = tmp1
                    iwork( j ) = itmp1
                    call stdlib${ii}$_${ri}$swap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ )
                 end if
              end do
           end if
            ! causes problems with tests 19
            ! if (wantz .and. indeig ) z( 1,1) = z(1,1) / 1.002_${rk}$ + .002
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_${ri}$stevr

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sstevx( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, &
     !! SSTEVX computes selected eigenvalues and, optionally, eigenvectors
     !! of a real symmetric tridiagonal matrix A.  Eigenvalues and
     !! eigenvectors can be selected by specifying either a range of values
     !! or a range of indices for the desired eigenvalues.
               work, iwork, ifail, info )
        ! -- lapack driver 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) :: jobz, range
           integer(${ik}$), intent(in) :: il, iu, ldz, n
           integer(${ik}$), intent(out) :: info, m
           real(sp), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ifail(*), iwork(*)
           real(sp), intent(inout) :: d(*), e(*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, valeig, wantz
           character :: order
           integer(${ik}$) :: i, imax, indibl, indisp, indiwo, indwrk, iscale, itmp1, j, jj, &
                     nsplit
           real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, &
                     vuu
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -8_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -9_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) )info = -14_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSTEVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( alleig .or. indeig ) then
                 m = 1_${ik}$
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
              else
                 if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then
                    m = 1_${ik}$
                    w( 1_${ik}$ ) = d( 1_${ik}$ )
                 end if
              end if
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_slamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           if ( valeig ) then
              vll = vl
              vuu = vu
           else
              vll = zero
              vuu = zero
           endif
           tnrm = stdlib${ii}$_slanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
              if( valeig ) then
                 vll = vl*sigma
                 vuu = vu*sigma
              end if
           end if
           ! if all eigenvalues are desired and abstol is less than zero, then
           ! call stdlib${ii}$_ssterf or stdlib${ii}$_ssteqr.  if this fails for some eigenvalue, then
           ! try stdlib${ii}$_sstebz.
           test = .false.
           if( indeig ) then
              if( il==1_${ik}$ .and. iu==n ) then
                 test = .true.
              end if
           end if
           if( ( alleig .or. test ) .and. ( abstol<=zero ) ) then
              call stdlib${ii}$_scopy( n, d, 1_${ik}$, w, 1_${ik}$ )
              call stdlib${ii}$_scopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ )
              indwrk = n + 1_${ik}$
              if( .not.wantz ) then
                 call stdlib${ii}$_ssterf( n, w, work, info )
              else
                 call stdlib${ii}$_ssteqr( 'I', n, w, work, z, ldz, work( indwrk ), info )
                 if( info==0_${ik}$ ) then
                    do i = 1, n
                       ifail( i ) = 0_${ik}$
                    end do
                 end if
              end if
              if( info==0_${ik}$ ) then
                 m = n
                 go to 20
              end if
              info = 0_${ik}$
           end if
           ! otherwise, call stdlib${ii}$_sstebz and, if eigenvectors are desired, stdlib${ii}$_sstein.
           if( wantz ) then
              order = 'B'
           else
              order = 'E'
           end if
           indwrk = 1_${ik}$
           indibl = 1_${ik}$
           indisp = indibl + n
           indiwo = indisp + n
           call stdlib${ii}$_sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, &
                     iwork( indibl ), iwork( indisp ),work( indwrk ), iwork( indiwo ), info )
           if( wantz ) then
              call stdlib${ii}$_sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work( &
                        indwrk ), iwork( indiwo ), ifail,info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           20 continue
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = m
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_sscal( imax, one / sigma, w, 1_${ik}$ )
           end if
           ! if eigenvalues are not in order, then sort them, along with
           ! eigenvectors.
           if( wantz ) then
              do j = 1, m - 1
                 i = 0_${ik}$
                 tmp1 = w( j )
                 do jj = j + 1, m
                    if( w( jj )<tmp1 ) then
                       i = jj
                       tmp1 = w( jj )
                    end if
                 end do
                 if( i/=0_${ik}$ ) then
                    itmp1 = iwork( indibl+i-1 )
                    w( i ) = w( j )
                    iwork( indibl+i-1 ) = iwork( indibl+j-1 )
                    w( j ) = tmp1
                    iwork( indibl+j-1 ) = itmp1
                    call stdlib${ii}$_sswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ )
                    if( info/=0_${ik}$ ) then
                       itmp1 = ifail( i )
                       ifail( i ) = ifail( j )
                       ifail( j ) = itmp1
                    end if
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_sstevx

     pure module subroutine stdlib${ii}$_dstevx( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, &
     !! DSTEVX computes selected eigenvalues and, optionally, eigenvectors
     !! of a real symmetric tridiagonal matrix A.  Eigenvalues and
     !! eigenvectors can be selected by specifying either a range of values
     !! or a range of indices for the desired eigenvalues.
               work, iwork, ifail, info )
        ! -- lapack driver 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) :: jobz, range
           integer(${ik}$), intent(in) :: il, iu, ldz, n
           integer(${ik}$), intent(out) :: info, m
           real(dp), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ifail(*), iwork(*)
           real(dp), intent(inout) :: d(*), e(*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, valeig, wantz
           character :: order
           integer(${ik}$) :: i, imax, indibl, indisp, indiwo, indwrk, iscale, itmp1, j, jj, &
                     nsplit
           real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, &
                     vuu
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -8_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -9_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) )info = -14_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( alleig .or. indeig ) then
                 m = 1_${ik}$
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
              else
                 if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then
                    m = 1_${ik}$
                    w( 1_${ik}$ ) = d( 1_${ik}$ )
                 end if
              end if
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_dlamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           if( valeig ) then
              vll = vl
              vuu = vu
           else
              vll = zero
              vuu = zero
           end if
           tnrm = stdlib${ii}$_dlanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
              if( valeig ) then
                 vll = vl*sigma
                 vuu = vu*sigma
              end if
           end if
           ! if all eigenvalues are desired and abstol is less than zero, then
           ! call stdlib${ii}$_dsterf or stdlib${ii}$_ssteqr.  if this fails for some eigenvalue, then
           ! try stdlib${ii}$_dstebz.
           test = .false.
           if( indeig ) then
              if( il==1_${ik}$ .and. iu==n ) then
                 test = .true.
              end if
           end if
           if( ( alleig .or. test ) .and. ( abstol<=zero ) ) then
              call stdlib${ii}$_dcopy( n, d, 1_${ik}$, w, 1_${ik}$ )
              call stdlib${ii}$_dcopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ )
              indwrk = n + 1_${ik}$
              if( .not.wantz ) then
                 call stdlib${ii}$_dsterf( n, w, work, info )
              else
                 call stdlib${ii}$_dsteqr( 'I', n, w, work, z, ldz, work( indwrk ), info )
                 if( info==0_${ik}$ ) then
                    do i = 1, n
                       ifail( i ) = 0_${ik}$
                    end do
                 end if
              end if
              if( info==0_${ik}$ ) then
                 m = n
                 go to 20
              end if
              info = 0_${ik}$
           end if
           ! otherwise, call stdlib${ii}$_dstebz and, if eigenvectors are desired, stdlib${ii}$_sstein.
           if( wantz ) then
              order = 'B'
           else
              order = 'E'
           end if
           indwrk = 1_${ik}$
           indibl = 1_${ik}$
           indisp = indibl + n
           indiwo = indisp + n
           call stdlib${ii}$_dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, &
                     iwork( indibl ), iwork( indisp ),work( indwrk ), iwork( indiwo ), info )
           if( wantz ) then
              call stdlib${ii}$_dstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work( &
                        indwrk ), iwork( indiwo ), ifail,info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           20 continue
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = m
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_dscal( imax, one / sigma, w, 1_${ik}$ )
           end if
           ! if eigenvalues are not in order, then sort them, along with
           ! eigenvectors.
           if( wantz ) then
              do j = 1, m - 1
                 i = 0_${ik}$
                 tmp1 = w( j )
                 do jj = j + 1, m
                    if( w( jj )<tmp1 ) then
                       i = jj
                       tmp1 = w( jj )
                    end if
                 end do
                 if( i/=0_${ik}$ ) then
                    itmp1 = iwork( indibl+i-1 )
                    w( i ) = w( j )
                    iwork( indibl+i-1 ) = iwork( indibl+j-1 )
                    w( j ) = tmp1
                    iwork( indibl+j-1 ) = itmp1
                    call stdlib${ii}$_dswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ )
                    if( info/=0_${ik}$ ) then
                       itmp1 = ifail( i )
                       ifail( i ) = ifail( j )
                       ifail( j ) = itmp1
                    end if
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_dstevx

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$stevx( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, &
     !! DSTEVX: computes selected eigenvalues and, optionally, eigenvectors
     !! of a real symmetric tridiagonal matrix A.  Eigenvalues and
     !! eigenvectors can be selected by specifying either a range of values
     !! or a range of indices for the desired eigenvalues.
               work, iwork, ifail, info )
        ! -- lapack driver 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) :: jobz, range
           integer(${ik}$), intent(in) :: il, iu, ldz, n
           integer(${ik}$), intent(out) :: info, m
           real(${rk}$), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ifail(*), iwork(*)
           real(${rk}$), intent(inout) :: d(*), e(*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, valeig, wantz
           character :: order
           integer(${ik}$) :: i, imax, indibl, indisp, indiwo, indwrk, iscale, itmp1, j, jj, &
                     nsplit
           real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, &
                     vuu
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -8_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -9_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) )info = -14_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( alleig .or. indeig ) then
                 m = 1_${ik}$
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
              else
                 if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then
                    m = 1_${ik}$
                    w( 1_${ik}$ ) = d( 1_${ik}$ )
                 end if
              end if
              if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           ! get machine constants.
           safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
           eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' )
           smlnum = safmin / eps
           bignum = one / smlnum
           rmin = sqrt( smlnum )
           rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
           ! scale matrix to allowable range, if necessary.
           iscale = 0_${ik}$
           if( valeig ) then
              vll = vl
              vuu = vu
           else
              vll = zero
              vuu = zero
           end if
           tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e )
           if( tnrm>zero .and. tnrm<rmin ) then
              iscale = 1_${ik}$
              sigma = rmin / tnrm
           else if( tnrm>rmax ) then
              iscale = 1_${ik}$
              sigma = rmax / tnrm
           end if
           if( iscale==1_${ik}$ ) then
              call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ )
              call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ )
              if( valeig ) then
                 vll = vl*sigma
                 vuu = vu*sigma
              end if
           end if
           ! if all eigenvalues are desired and abstol is less than zero, then
           ! call stdlib${ii}$_${ri}$sterf or stdlib${ii}$_dsteqr.  if this fails for some eigenvalue, then
           ! try stdlib${ii}$_${ri}$stebz.
           test = .false.
           if( indeig ) then
              if( il==1_${ik}$ .and. iu==n ) then
                 test = .true.
              end if
           end if
           if( ( alleig .or. test ) .and. ( abstol<=zero ) ) then
              call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, w, 1_${ik}$ )
              call stdlib${ii}$_${ri}$copy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ )
              indwrk = n + 1_${ik}$
              if( .not.wantz ) then
                 call stdlib${ii}$_${ri}$sterf( n, w, work, info )
              else
                 call stdlib${ii}$_${ri}$steqr( 'I', n, w, work, z, ldz, work( indwrk ), info )
                 if( info==0_${ik}$ ) then
                    do i = 1, n
                       ifail( i ) = 0_${ik}$
                    end do
                 end if
              end if
              if( info==0_${ik}$ ) then
                 m = n
                 go to 20
              end if
              info = 0_${ik}$
           end if
           ! otherwise, call stdlib${ii}$_${ri}$stebz and, if eigenvectors are desired, stdlib${ii}$_dstein.
           if( wantz ) then
              order = 'B'
           else
              order = 'E'
           end if
           indwrk = 1_${ik}$
           indibl = 1_${ik}$
           indisp = indibl + n
           indiwo = indisp + n
           call stdlib${ii}$_${ri}$stebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, &
                     iwork( indibl ), iwork( indisp ),work( indwrk ), iwork( indiwo ), info )
           if( wantz ) then
              call stdlib${ii}$_${ri}$stein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work( &
                        indwrk ), iwork( indiwo ), ifail,info )
           end if
           ! if matrix was scaled, then rescale eigenvalues appropriately.
           20 continue
           if( iscale==1_${ik}$ ) then
              if( info==0_${ik}$ ) then
                 imax = m
              else
                 imax = info - 1_${ik}$
              end if
              call stdlib${ii}$_${ri}$scal( imax, one / sigma, w, 1_${ik}$ )
           end if
           ! if eigenvalues are not in order, then sort them, along with
           ! eigenvectors.
           if( wantz ) then
              do j = 1, m - 1
                 i = 0_${ik}$
                 tmp1 = w( j )
                 do jj = j + 1, m
                    if( w( jj )<tmp1 ) then
                       i = jj
                       tmp1 = w( jj )
                    end if
                 end do
                 if( i/=0_${ik}$ ) then
                    itmp1 = iwork( indibl+i-1 )
                    w( i ) = w( j )
                    iwork( indibl+i-1 ) = iwork( indibl+j-1 )
                    w( j ) = tmp1
                    iwork( indibl+j-1 ) = itmp1
                    call stdlib${ii}$_${ri}$swap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ )
                    if( info/=0_${ik}$ ) then
                       itmp1 = ifail( i )
                       ifail( i ) = ifail( j )
                       ifail( j ) = itmp1
                    end if
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_${ri}$stevx

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_spteqr( compz, n, d, e, z, ldz, work, info )
     !! SPTEQR computes all eigenvalues and, optionally, eigenvectors of a
     !! symmetric positive definite tridiagonal matrix by first factoring the
     !! matrix using SPTTRF, and then calling SBDSQR to compute the singular
     !! values of the bidiagonal factor.
     !! This routine computes the eigenvalues of the positive definite
     !! tridiagonal matrix to high relative accuracy.  This means that if the
     !! eigenvalues range over many orders of magnitude in size, then the
     !! small eigenvalues and corresponding eigenvectors will be computed
     !! more accurately than, for example, with the standard QR method.
     !! The eigenvectors of a full or band symmetric positive definite matrix
     !! can also be found if SSYTRD, SSPTRD, or SSBTRD has been used to
     !! reduce this matrix to tridiagonal form. (The reduction to tridiagonal
     !! form, however, may preclude the possibility of obtaining high
     !! relative accuracy in the small eigenvalues of the original matrix, if
     !! these eigenvalues range over many orders of magnitude.)
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: compz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(sp), intent(inout) :: d(*), e(*), z(ldz,*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Arrays 
           real(sp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$)
           ! Local Scalars 
           integer(${ik}$) :: i, icompz, nru
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( stdlib_lsame( compz, 'N' ) ) then
              icompz = 0_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              icompz = 2_${ik}$
           else
              icompz = -1_${ik}$
           end if
           if( icompz<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SPTEQR', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           if( icompz==2_${ik}$ )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, z, ldz )
           ! call stdlib${ii}$_spttrf to factor the matrix.
           call stdlib${ii}$_spttrf( n, d, e, info )
           if( info/=0 )return
           do i = 1, n
              d( i ) = sqrt( d( i ) )
           end do
           do i = 1, n - 1
              e( i ) = e( i )*d( i )
           end do
           ! call stdlib${ii}$_sbdsqr to compute the singular values/vectors of the
           ! bidiagonal factor.
           if( icompz>0_${ik}$ ) then
              nru = n
           else
              nru = 0_${ik}$
           end if
           call stdlib${ii}$_sbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info )
                     
           ! square the singular values.
           if( info==0_${ik}$ ) then
              do i = 1, n
                 d( i ) = d( i )*d( i )
              end do
           else
              info = n + info
           end if
           return
     end subroutine stdlib${ii}$_spteqr

     pure module subroutine stdlib${ii}$_dpteqr( compz, n, d, e, z, ldz, work, info )
     !! DPTEQR computes all eigenvalues and, optionally, eigenvectors of a
     !! symmetric positive definite tridiagonal matrix by first factoring the
     !! matrix using DPTTRF, and then calling DBDSQR to compute the singular
     !! values of the bidiagonal factor.
     !! This routine computes the eigenvalues of the positive definite
     !! tridiagonal matrix to high relative accuracy.  This means that if the
     !! eigenvalues range over many orders of magnitude in size, then the
     !! small eigenvalues and corresponding eigenvectors will be computed
     !! more accurately than, for example, with the standard QR method.
     !! The eigenvectors of a full or band symmetric positive definite matrix
     !! can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to
     !! reduce this matrix to tridiagonal form. (The reduction to tridiagonal
     !! form, however, may preclude the possibility of obtaining high
     !! relative accuracy in the small eigenvalues of the original matrix, if
     !! these eigenvalues range over many orders of magnitude.)
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: compz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(dp), intent(inout) :: d(*), e(*), z(ldz,*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Arrays 
           real(dp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$)
           ! Local Scalars 
           integer(${ik}$) :: i, icompz, nru
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( stdlib_lsame( compz, 'N' ) ) then
              icompz = 0_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              icompz = 2_${ik}$
           else
              icompz = -1_${ik}$
           end if
           if( icompz<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DPTEQR', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           if( icompz==2_${ik}$ )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, z, ldz )
           ! call stdlib${ii}$_dpttrf to factor the matrix.
           call stdlib${ii}$_dpttrf( n, d, e, info )
           if( info/=0 )return
           do i = 1, n
              d( i ) = sqrt( d( i ) )
           end do
           do i = 1, n - 1
              e( i ) = e( i )*d( i )
           end do
           ! call stdlib${ii}$_dbdsqr to compute the singular values/vectors of the
           ! bidiagonal factor.
           if( icompz>0_${ik}$ ) then
              nru = n
           else
              nru = 0_${ik}$
           end if
           call stdlib${ii}$_dbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info )
                     
           ! square the singular values.
           if( info==0_${ik}$ ) then
              do i = 1, n
                 d( i ) = d( i )*d( i )
              end do
           else
              info = n + info
           end if
           return
     end subroutine stdlib${ii}$_dpteqr

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$pteqr( compz, n, d, e, z, ldz, work, info )
     !! DPTEQR: computes all eigenvalues and, optionally, eigenvectors of a
     !! symmetric positive definite tridiagonal matrix by first factoring the
     !! matrix using DPTTRF, and then calling DBDSQR to compute the singular
     !! values of the bidiagonal factor.
     !! This routine computes the eigenvalues of the positive definite
     !! tridiagonal matrix to high relative accuracy.  This means that if the
     !! eigenvalues range over many orders of magnitude in size, then the
     !! small eigenvalues and corresponding eigenvectors will be computed
     !! more accurately than, for example, with the standard QR method.
     !! The eigenvectors of a full or band symmetric positive definite matrix
     !! can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to
     !! reduce this matrix to tridiagonal form. (The reduction to tridiagonal
     !! form, however, may preclude the possibility of obtaining high
     !! relative accuracy in the small eigenvalues of the original matrix, if
     !! these eigenvalues range over many orders of magnitude.)
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: compz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: d(*), e(*), z(ldz,*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Arrays 
           real(${rk}$) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$)
           ! Local Scalars 
           integer(${ik}$) :: i, icompz, nru
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( stdlib_lsame( compz, 'N' ) ) then
              icompz = 0_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              icompz = 2_${ik}$
           else
              icompz = -1_${ik}$
           end if
           if( icompz<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DPTEQR', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = one
              return
           end if
           if( icompz==2_${ik}$ )call stdlib${ii}$_${ri}$laset( 'FULL', n, n, zero, one, z, ldz )
           ! call stdlib${ii}$_${ri}$pttrf to factor the matrix.
           call stdlib${ii}$_${ri}$pttrf( n, d, e, info )
           if( info/=0 )return
           do i = 1, n
              d( i ) = sqrt( d( i ) )
           end do
           do i = 1, n - 1
              e( i ) = e( i )*d( i )
           end do
           ! call stdlib${ii}$_${ri}$bdsqr to compute the singular values/vectors of the
           ! bidiagonal factor.
           if( icompz>0_${ik}$ ) then
              nru = n
           else
              nru = 0_${ik}$
           end if
           call stdlib${ii}$_${ri}$bdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info )
                     
           ! square the singular values.
           if( info==0_${ik}$ ) then
              do i = 1, n
                 d( i ) = d( i )*d( i )
              end do
           else
              info = n + info
           end if
           return
     end subroutine stdlib${ii}$_${ri}$pteqr

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_cpteqr( compz, n, d, e, z, ldz, work, info )
     !! CPTEQR computes all eigenvalues and, optionally, eigenvectors of a
     !! symmetric positive definite tridiagonal matrix by first factoring the
     !! matrix using SPTTRF and then calling CBDSQR to compute the singular
     !! values of the bidiagonal factor.
     !! This routine computes the eigenvalues of the positive definite
     !! tridiagonal matrix to high relative accuracy.  This means that if the
     !! eigenvalues range over many orders of magnitude in size, then the
     !! small eigenvalues and corresponding eigenvectors will be computed
     !! more accurately than, for example, with the standard QR method.
     !! The eigenvectors of a full or band positive definite Hermitian matrix
     !! can also be found if CHETRD, CHPTRD, or CHBTRD has been used to
     !! reduce this matrix to tridiagonal form.  (The reduction to
     !! tridiagonal form, however, may preclude the possibility of obtaining
     !! high relative accuracy in the small eigenvalues of the original
     !! matrix, if these eigenvalues range over many orders of magnitude.)
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: compz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(sp), intent(inout) :: d(*), e(*)
           real(sp), intent(out) :: work(*)
           complex(sp), intent(inout) :: z(ldz,*)
        ! ====================================================================
           
           ! Local Arrays 
           complex(sp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$)
           ! Local Scalars 
           integer(${ik}$) :: i, icompz, nru
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( stdlib_lsame( compz, 'N' ) ) then
              icompz = 0_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              icompz = 2_${ik}$
           else
              icompz = -1_${ik}$
           end if
           if( icompz<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CPTEQR', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = cone
              return
           end if
           if( icompz==2_${ik}$ )call stdlib${ii}$_claset( 'FULL', n, n, czero, cone, z, ldz )
           ! call stdlib${ii}$_spttrf to factor the matrix.
           call stdlib${ii}$_spttrf( n, d, e, info )
           if( info/=0 )return
           do i = 1, n
              d( i ) = sqrt( d( i ) )
           end do
           do i = 1, n - 1
              e( i ) = e( i )*d( i )
           end do
           ! call stdlib${ii}$_cbdsqr to compute the singular values/vectors of the
           ! bidiagonal factor.
           if( icompz>0_${ik}$ ) then
              nru = n
           else
              nru = 0_${ik}$
           end if
           call stdlib${ii}$_cbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info )
                     
           ! square the singular values.
           if( info==0_${ik}$ ) then
              do i = 1, n
                 d( i ) = d( i )*d( i )
              end do
           else
              info = n + info
           end if
           return
     end subroutine stdlib${ii}$_cpteqr

     pure module subroutine stdlib${ii}$_zpteqr( compz, n, d, e, z, ldz, work, info )
     !! ZPTEQR computes all eigenvalues and, optionally, eigenvectors of a
     !! symmetric positive definite tridiagonal matrix by first factoring the
     !! matrix using DPTTRF and then calling ZBDSQR to compute the singular
     !! values of the bidiagonal factor.
     !! This routine computes the eigenvalues of the positive definite
     !! tridiagonal matrix to high relative accuracy.  This means that if the
     !! eigenvalues range over many orders of magnitude in size, then the
     !! small eigenvalues and corresponding eigenvectors will be computed
     !! more accurately than, for example, with the standard QR method.
     !! The eigenvectors of a full or band positive definite Hermitian matrix
     !! can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to
     !! reduce this matrix to tridiagonal form.  (The reduction to
     !! tridiagonal form, however, may preclude the possibility of obtaining
     !! high relative accuracy in the small eigenvalues of the original
     !! matrix, if these eigenvalues range over many orders of magnitude.)
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: compz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(dp), intent(inout) :: d(*), e(*)
           real(dp), intent(out) :: work(*)
           complex(dp), intent(inout) :: z(ldz,*)
        ! ====================================================================
           
           ! Local Arrays 
           complex(dp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$)
           ! Local Scalars 
           integer(${ik}$) :: i, icompz, nru
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( stdlib_lsame( compz, 'N' ) ) then
              icompz = 0_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              icompz = 2_${ik}$
           else
              icompz = -1_${ik}$
           end if
           if( icompz<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZPTEQR', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = cone
              return
           end if
           if( icompz==2_${ik}$ )call stdlib${ii}$_zlaset( 'FULL', n, n, czero, cone, z, ldz )
           ! call stdlib${ii}$_dpttrf to factor the matrix.
           call stdlib${ii}$_dpttrf( n, d, e, info )
           if( info/=0 )return
           do i = 1, n
              d( i ) = sqrt( d( i ) )
           end do
           do i = 1, n - 1
              e( i ) = e( i )*d( i )
           end do
           ! call stdlib${ii}$_zbdsqr to compute the singular values/vectors of the
           ! bidiagonal factor.
           if( icompz>0_${ik}$ ) then
              nru = n
           else
              nru = 0_${ik}$
           end if
           call stdlib${ii}$_zbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info )
                     
           ! square the singular values.
           if( info==0_${ik}$ ) then
              do i = 1, n
                 d( i ) = d( i )*d( i )
              end do
           else
              info = n + info
           end if
           return
     end subroutine stdlib${ii}$_zpteqr

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$pteqr( compz, n, d, e, z, ldz, work, info )
     !! ZPTEQR: computes all eigenvalues and, optionally, eigenvectors of a
     !! symmetric positive definite tridiagonal matrix by first factoring the
     !! matrix using DPTTRF and then calling ZBDSQR to compute the singular
     !! values of the bidiagonal factor.
     !! This routine computes the eigenvalues of the positive definite
     !! tridiagonal matrix to high relative accuracy.  This means that if the
     !! eigenvalues range over many orders of magnitude in size, then the
     !! small eigenvalues and corresponding eigenvectors will be computed
     !! more accurately than, for example, with the standard QR method.
     !! The eigenvectors of a full or band positive definite Hermitian matrix
     !! can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to
     !! reduce this matrix to tridiagonal form.  (The reduction to
     !! tridiagonal form, however, may preclude the possibility of obtaining
     !! high relative accuracy in the small eigenvalues of the original
     !! matrix, if these eigenvalues range over many orders of magnitude.)
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: compz
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldz, n
           ! Array Arguments 
           real(${ck}$), intent(inout) :: d(*), e(*)
           real(${ck}$), intent(out) :: work(*)
           complex(${ck}$), intent(inout) :: z(ldz,*)
        ! ====================================================================
           
           ! Local Arrays 
           complex(${ck}$) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$)
           ! Local Scalars 
           integer(${ik}$) :: i, icompz, nru
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( stdlib_lsame( compz, 'N' ) ) then
              icompz = 0_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              icompz = 2_${ik}$
           else
              icompz = -1_${ik}$
           end if
           if( icompz<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZPTEQR', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           if( n==1_${ik}$ ) then
              if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = cone
              return
           end if
           if( icompz==2_${ik}$ )call stdlib${ii}$_${ci}$laset( 'FULL', n, n, czero, cone, z, ldz )
           ! call stdlib${ii}$_${c2ri(ci)}$pttrf to factor the matrix.
           call stdlib${ii}$_${c2ri(ci)}$pttrf( n, d, e, info )
           if( info/=0 )return
           do i = 1, n
              d( i ) = sqrt( d( i ) )
           end do
           do i = 1, n - 1
              e( i ) = e( i )*d( i )
           end do
           ! call stdlib${ii}$_${ci}$bdsqr to compute the singular values/vectors of the
           ! bidiagonal factor.
           if( icompz>0_${ik}$ ) then
              nru = n
           else
              nru = 0_${ik}$
           end if
           call stdlib${ii}$_${ci}$bdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info )
                     
           ! square the singular values.
           if( info==0_${ik}$ ) then
              do i = 1, n
                 d( i ) = d( i )*d( i )
              end do
           else
              info = n + info
           end if
           return
     end subroutine stdlib${ii}$_${ci}$pteqr

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sstebz( range, order, n, vl, vu, il, iu, abstol, d, e,m, nsplit, w, &
     !! SSTEBZ computes the eigenvalues of a symmetric tridiagonal
     !! matrix T.  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.
               iblock, isplit, work, iwork,info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: order, range
           integer(${ik}$), intent(in) :: il, iu, n
           integer(${ik}$), intent(out) :: info, m, nsplit
           real(sp), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iblock(*), isplit(*), iwork(*)
           real(sp), intent(in) :: d(*), e(*)
           real(sp), intent(out) :: w(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: fudge = 2.1_sp
           real(sp), parameter :: relfac = two
           
           
           ! Local Scalars 
           logical(lk) :: ncnvrg, toofew
           integer(${ik}$) :: ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iorder, &
                     iout, irange, itmax, itmp1, iw, iwoff, j, jb, jdisc, je, nb, nwl, nwu
           real(sp) :: atoli, bnorm, gl, gu, pivmin, rtoli, safemn, tmp1, tmp2, tnorm, ulp, wkill,&
                      wl, wlu, wu, wul
           ! Local Arrays 
           integer(${ik}$) :: idumma(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! decode range
           if( stdlib_lsame( range, 'A' ) ) then
              irange = 1_${ik}$
           else if( stdlib_lsame( range, 'V' ) ) then
              irange = 2_${ik}$
           else if( stdlib_lsame( range, 'I' ) ) then
              irange = 3_${ik}$
           else
              irange = 0_${ik}$
           end if
           ! decode order
           if( stdlib_lsame( order, 'B' ) ) then
              iorder = 2_${ik}$
           else if( stdlib_lsame( order, 'E' ) ) then
              iorder = 1_${ik}$
           else
              iorder = 0_${ik}$
           end if
           ! check for errors
           if( irange<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( iorder<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( irange==2_${ik}$ ) then
              if( vl>=vu ) info = -5_${ik}$
           else if( irange==3_${ik}$ .and. ( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) )then
              info = -6_${ik}$
           else if( irange==3_${ik}$ .and. ( iu<min( n, il ) .or. iu>n ) )then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSTEBZ', -info )
              return
           end if
           ! initialize error flags
           info = 0_${ik}$
           ncnvrg = .false.
           toofew = .false.
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           ! simplifications:
           if( irange==3_${ik}$ .and. il==1_${ik}$ .and. iu==n )irange = 1_${ik}$
           ! get machine constants
           ! nb is the minimum vector length for vector bisection, or 0
           ! if only scalar is to be done.
           safemn = stdlib${ii}$_slamch( 'S' )
           ulp = stdlib${ii}$_slamch( 'P' )
           rtoli = ulp*relfac
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ )nb = 0_${ik}$
           ! special case when n=1
           if( n==1_${ik}$ ) then
              nsplit = 1_${ik}$
              isplit( 1_${ik}$ ) = 1_${ik}$
              if( irange==2_${ik}$ .and. ( vl>=d( 1_${ik}$ ) .or. vu<d( 1_${ik}$ ) ) ) then
                 m = 0_${ik}$
              else
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
                 iblock( 1_${ik}$ ) = 1_${ik}$
                 m = 1_${ik}$
              end if
              return
           end if
           ! compute splitting points
           nsplit = 1_${ik}$
           work( n ) = zero
           pivmin = one
           do j = 2, n
              tmp1 = e( j-1 )**2_${ik}$
              if( abs( d( j )*d( j-1 ) )*ulp**2_${ik}$+safemn>tmp1 ) then
                 isplit( nsplit ) = j - 1_${ik}$
                 nsplit = nsplit + 1_${ik}$
                 work( j-1 ) = zero
              else
                 work( j-1 ) = tmp1
                 pivmin = max( pivmin, tmp1 )
              end if
           end do
           isplit( nsplit ) = n
           pivmin = pivmin*safemn
           ! compute interval and atoli
           if( irange==3_${ik}$ ) then
              ! range='i': compute the interval containing eigenvalues
                         ! il through iu.
              ! compute gershgorin interval for entire (split) matrix
              ! and use it as the initial interval
              gu = d( 1_${ik}$ )
              gl = d( 1_${ik}$ )
              tmp1 = zero
              do j = 1, n - 1
                 tmp2 = sqrt( work( j ) )
                 gu = max( gu, d( j )+tmp1+tmp2 )
                 gl = min( gl, d( j )-tmp1-tmp2 )
                 tmp1 = tmp2
              end do
              gu = max( gu, d( n )+tmp1 )
              gl = min( gl, d( n )-tmp1 )
              tnorm = max( abs( gl ), abs( gu ) )
              gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
              gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
              ! compute iteration parameters
              itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
              if( abstol<=zero ) then
                 atoli = ulp*tnorm
              else
                 atoli = abstol
              end if
              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,work, iwork( &
                        5_${ik}$ ), work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo )
              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
              if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then
                 info = 4_${ik}$
                 return
              end if
           else
              ! range='a' or 'v' -- set atoli
              tnorm = max( abs( d( 1_${ik}$ ) )+abs( e( 1_${ik}$ ) ),abs( d( n ) )+abs( e( n-1 ) ) )
              do j = 2, n - 1
                 tnorm = max( tnorm, abs( d( j ) )+abs( e( j-1 ) )+abs( e( j ) ) )
              end do
              if( abstol<=zero ) then
                 atoli = ulp*tnorm
              else
                 atoli = abstol
              end if
              if( irange==2_${ik}$ ) then
                 wl = vl
                 wu = vu
              else
                 wl = zero
                 wu = zero
              end if
           end if
           ! 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 jb = 1, nsplit
              ioff = iend
              ibegin = ioff + 1_${ik}$
              iend = isplit( jb )
              in = iend - ioff
              if( in==1_${ik}$ ) then
                 ! special case -- in=1
                 if( irange==1_${ik}$ .or. wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$
                 if( irange==1_${ik}$ .or. wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$
                 if( irange==1_${ik}$ .or. ( wl<d( ibegin )-pivmin .and. wu>=d( ibegin )-pivmin ) ) &
                           then
                    m = m + 1_${ik}$
                    w( m ) = d( ibegin )
                    iblock( m ) = jb
                 end if
              else
                 ! general case -- in > 1
                 ! compute gershgorin interval
                 ! and use it as the initial interval
                 gu = d( ibegin )
                 gl = d( ibegin )
                 tmp1 = zero
                 do j = ibegin, iend - 1
                    tmp2 = abs( e( j ) )
                    gu = max( gu, d( j )+tmp1+tmp2 )
                    gl = min( gl, d( j )-tmp1-tmp2 )
                    tmp1 = tmp2
                 end do
                 gu = max( gu, d( iend )+tmp1 )
                 gl = min( gl, d( iend )-tmp1 )
                 bnorm = max( abs( gl ), abs( gu ) )
                 gl = gl - fudge*bnorm*ulp*in - fudge*pivmin
                 gu = gu + fudge*bnorm*ulp*in + fudge*pivmin
                 ! compute atoli for the current submatrix
                 if( abstol<=zero ) then
                    atoli = ulp*max( abs( gl ), abs( gu ) )
                 else
                    atoli = abstol
                 end if
                 if( irange>1_${ik}$ ) then
                    if( gu<wl ) then
                       nwl = nwl + in
                       nwu = nwu + in
                       cycle loop_70
                    end if
                    gl = max( gl, wl )
                    gu = min( gu, wu )
                    if( gl>=gu )cycle loop_70
                 end if
                 ! set up initial interval
                 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 ), work( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 &
                           ), iblock( m+1 ), iinfo )
                 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 ), work( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( &
                            m+1 ), iblock( m+1 ), iinfo )
                 ! copy eigenvalues into w and iblock
                 ! use -jb for block number for unconverged eigenvalues.
                 do j = 1, iout
                    tmp1 = half*( work( j+n )+work( j+in+n ) )
                    ! flag non-convergence.
                    if( j>iout-iinfo ) then
                       ncnvrg = .true.
                       ib = -jb
                    else
                       ib = jb
                    end if
                    do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff
                       w( je ) = tmp1
                       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==3_${ik}$ ) then
              im = 0_${ik}$
              idiscl = il - 1_${ik}$ - nwl
              idiscu = nwu - iu
              if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
                 do je = 1, m
                    if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then
                       idiscl = idiscl - 1_${ik}$
                    else if( w( je )>=wul .and. idiscu>0_${ik}$ ) then
                       idiscu = idiscu - 1_${ik}$
                    else
                       im = im + 1_${ik}$
                       w( im ) = w( je )
                       iblock( im ) = iblock( je )
                    end if
                 end do
                 m = im
              end if
              if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
                 ! code to deal with effects of bad arithmetic:
                 ! 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 simply finding the smallest/largest
                 ! eigenvalue(s).
                 ! (if n(w) is monotone non-decreasing, this should never
                     ! happen.)
                 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
                 im = 0_${ik}$
                 do je = 1, m
                    if( iblock( je )/=0_${ik}$ ) then
                       im = im + 1_${ik}$
                       w( im ) = w( 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 order='b', do nothing -- the eigenvalues are already sorted
              ! by block.
           ! if order='e', sort the eigenvalues from smallest to largest
           if( iorder==1_${ik}$ .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
                    itmp1 = iblock( ie )
                    w( ie ) = w( je )
                    iblock( ie ) = iblock( je )
                    w( je ) = tmp1
                    iblock( je ) = itmp1
                 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}$_sstebz

     pure module subroutine stdlib${ii}$_dstebz( range, order, n, vl, vu, il, iu, abstol, d, e,m, nsplit, w, &
     !! DSTEBZ computes the eigenvalues of a symmetric tridiagonal
     !! matrix T.  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.
               iblock, isplit, work, iwork,info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: order, range
           integer(${ik}$), intent(in) :: il, iu, n
           integer(${ik}$), intent(out) :: info, m, nsplit
           real(dp), intent(in) :: abstol, vl, vu
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iblock(*), isplit(*), iwork(*)
           real(dp), intent(in) :: d(*), e(*)
           real(dp), intent(out) :: w(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: fudge = 2.1_dp
           real(dp), parameter :: relfac = two
           
           
           ! Local Scalars 
           logical(lk) :: ncnvrg, toofew
           integer(${ik}$) :: ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iorder, &
                     iout, irange, itmax, itmp1, iw, iwoff, j, jb, jdisc, je, nb, nwl, nwu
           real(dp) :: atoli, bnorm, gl, gu, pivmin, rtoli, safemn, tmp1, tmp2, tnorm, ulp, wkill,&
                      wl, wlu, wu, wul
           ! Local Arrays 
           integer(${ik}$) :: idumma(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! decode range
           if( stdlib_lsame( range, 'A' ) ) then
              irange = 1_${ik}$
           else if( stdlib_lsame( range, 'V' ) ) then
              irange = 2_${ik}$
           else if( stdlib_lsame( range, 'I' ) ) then
              irange = 3_${ik}$
           else
              irange = 0_${ik}$
           end if
           ! decode order
           if( stdlib_lsame( order, 'B' ) ) then
              iorder = 2_${ik}$
           else if( stdlib_lsame( order, 'E' ) ) then
              iorder = 1_${ik}$
           else
              iorder = 0_${ik}$
           end if
           ! check for errors
           if( irange<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( iorder<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( irange==2_${ik}$ ) then
              if( vl>=vu )info = -5_${ik}$
           else if( irange==3_${ik}$ .and. ( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) )then
              info = -6_${ik}$
           else if( irange==3_${ik}$ .and. ( iu<min( n, il ) .or. iu>n ) )then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSTEBZ', -info )
              return
           end if
           ! initialize error flags
           info = 0_${ik}$
           ncnvrg = .false.
           toofew = .false.
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           ! simplifications:
           if( irange==3_${ik}$ .and. il==1_${ik}$ .and. iu==n )irange = 1_${ik}$
           ! get machine constants
           ! nb is the minimum vector length for vector bisection, or 0
           ! if only scalar is to be done.
           safemn = stdlib${ii}$_dlamch( 'S' )
           ulp = stdlib${ii}$_dlamch( 'P' )
           rtoli = ulp*relfac
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ )nb = 0_${ik}$
           ! special case when n=1
           if( n==1_${ik}$ ) then
              nsplit = 1_${ik}$
              isplit( 1_${ik}$ ) = 1_${ik}$
              if( irange==2_${ik}$ .and. ( vl>=d( 1_${ik}$ ) .or. vu<d( 1_${ik}$ ) ) ) then
                 m = 0_${ik}$
              else
                 w( 1_${ik}$ ) = d( 1_${ik}$ )
                 iblock( 1_${ik}$ ) = 1_${ik}$
                 m = 1_${ik}$
              end if
              return
           end if
           ! compute splitting points
           nsplit = 1_${ik}$
           work( n ) = zero
           pivmin = one
           do j = 2, n
              tmp1 = e( j-1 )**2_${ik}$
              if( abs( d( j )*d( j-1 ) )*ulp**2_${ik}$+safemn>tmp1 ) then
                 isplit( nsplit ) = j - 1_${ik}$
                 nsplit = nsplit + 1_${ik}$
                 work( j-1 ) = zero
              else
                 work( j-1 ) = tmp1
                 pivmin = max( pivmin, tmp1 )
              end if
           end do
           isplit( nsplit ) = n
           pivmin = pivmin*safemn
           ! compute interval and atoli
           if( irange==3_${ik}$ ) then
              ! range='i': compute the interval containing eigenvalues
                         ! il through iu.
              ! compute gershgorin interval for entire (split) matrix
              ! and use it as the initial interval
              gu = d( 1_${ik}$ )
              gl = d( 1_${ik}$ )
              tmp1 = zero
              do j = 1, n - 1
                 tmp2 = sqrt( work( j ) )
                 gu = max( gu, d( j )+tmp1+tmp2 )
                 gl = min( gl, d( j )-tmp1-tmp2 )
                 tmp1 = tmp2
              end do
              gu = max( gu, d( n )+tmp1 )
              gl = min( gl, d( n )-tmp1 )
              tnorm = max( abs( gl ), abs( gu ) )
              gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
              gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
              ! compute iteration parameters
              itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
              if( abstol<=zero ) then
                 atoli = ulp*tnorm
              else
                 atoli = abstol
              end if
              work( n+1 ) = gl
              work( n+2 ) = gl
              work( n+3 ) = gu
              work( n+</