stdlib_lapack_eigv_tridiag.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_slaebz( ijob, nitmax, n, mmax, minp, nbmin, abstol,reltol, pivmin, d, &
     !! SLAEBZ contains the iteration loops which compute and use the
     !! function N(w), which is the count of eigenvalues of a symmetric
     !! tridiagonal matrix T less than or equal to its argument  w.  It
     !! performs a choice of two types of loops:
     !! IJOB=1, followed by
     !! IJOB=2: It takes as input a list of intervals and returns a list of
     !! sufficiently small intervals whose union contains the same
     !! eigenvalues as the union of the original intervals.
     !! The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
     !! The output interval (AB(j,1),AB(j,2)] will contain
     !! eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
     !! IJOB=3: It performs a binary search in each input interval
     !! (AB(j,1),AB(j,2)] for a point  w(j)  such that
     !! N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
     !! the search.  If such a w(j) is found, then on output
     !! AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
     !! (AB(j,1),AB(j,2)] will be a small interval containing the
     !! point where N(w) jumps through NVAL(j), unless that point
     !! lies outside the initial interval.
     !! Note that the intervals are in all cases half-open intervals,
     !! i.e., of the form  (a,b] , which includes  b  but not  a .
     !! To avoid underflow, the matrix should be scaled so that its largest
     !! element is no greater than  overflow**(1/2) * underflow**(1/4)
     !! in absolute value.  To assure the most accurate computation
     !! of small eigenvalues, the matrix should be scaled to be
     !! not much smaller than that, either.
     !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
     !! Matrix", Report CS41, Computer Science Dept., Stanford
     !! University, July 21, 1966
     !! Note: the arguments are, in general, *not* checked for unreasonable
     !! values.
               e, e2, nval, ab, c, mout,nab, work, iwork, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ijob, minp, mmax, n, nbmin, nitmax
           integer(${ik}$), intent(out) :: info, mout
           real(sp), intent(in) :: abstol, pivmin, reltol
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           integer(${ik}$), intent(inout) :: nab(mmax,*), nval(*)
           real(sp), intent(inout) :: ab(mmax,*), c(*)
           real(sp), intent(in) :: d(*), e(*), e2(*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl, klnew
           real(sp) :: tmp1, tmp2
           ! Intrinsic Functions 
           ! Executable Statements 
           ! check for errors
           info = 0_${ik}$
           if( ijob<1_${ik}$ .or. ijob>3_${ik}$ ) then
              info = -1_${ik}$
              return
           end if
           ! initialize nab
           if( ijob==1_${ik}$ ) then
              ! compute the number of eigenvalues in the initial intervals.
              mout = 0_${ik}$
              do ji = 1, minp
                 do jp = 1, 2
                    tmp1 = d( 1_${ik}$ ) - ab( ji, jp )
                    if( abs( tmp1 )<pivmin )tmp1 = -pivmin
                    nab( ji, jp ) = 0_${ik}$
                    if( tmp1<=zero )nab( ji, jp ) = 1_${ik}$
                    do j = 2, n
                       tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp )
                       if( abs( tmp1 )<pivmin )tmp1 = -pivmin
                       if( tmp1<=zero )nab( ji, jp ) = nab( ji, jp ) + 1_${ik}$
                    end do
                 end do
                 mout = mout + nab( ji, 2_${ik}$ ) - nab( ji, 1_${ik}$ )
              end do
              return
           end if
           ! initialize for loop
           ! kf and kl have the following meaning:
              ! intervals 1,...,kf-1 have converged.
              ! intervals kf,...,kl  still need to be refined.
           kf = 1_${ik}$
           kl = minp
           ! if ijob=2, initialize c.
           ! if ijob=3, use the user-supplied starting point.
           if( ijob==2_${ik}$ ) then
              do ji = 1, minp
                 c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) )
              end do
           end if
           ! iteration loop
           loop_130: do jit = 1, nitmax
              ! loop over intervals
              if( kl-kf+1>=nbmin .and. nbmin>0_${ik}$ ) then
                 ! begin of parallel version of the loop
                 do ji = kf, kl
                    ! compute n(c), the number of eigenvalues less than c
                    work( ji ) = d( 1_${ik}$ ) - c( ji )
                    iwork( ji ) = 0_${ik}$
                    if( work( ji )<=pivmin ) then
                       iwork( ji ) = 1_${ik}$
                       work( ji ) = min( work( ji ), -pivmin )
                    end if
                    do j = 2, n
                       work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji )
                       if( work( ji )<=pivmin ) then
                          iwork( ji ) = iwork( ji ) + 1_${ik}$
                          work( ji ) = min( work( ji ), -pivmin )
                       end if
                    end do
                 end do
                 if( ijob<=2_${ik}$ ) then
                    ! ijob=2: choose all intervals containing eigenvalues.
                    klnew = kl
                    loop_70: do ji = kf, kl
                       ! insure that n(w) is monotone
                       iwork( ji ) = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), iwork( ji ) ) )
                       ! update the queue -- add intervals if both halves
                       ! contain eigenvalues.
                       if( iwork( ji )==nab( ji, 2_${ik}$ ) ) then
                          ! no eigenvalue in the upper interval:
                          ! just use the lower interval.
                          ab( ji, 2_${ik}$ ) = c( ji )
                       else if( iwork( ji )==nab( ji, 1_${ik}$ ) ) then
                          ! no eigenvalue in the lower interval:
                          ! just use the upper interval.
                          ab( ji, 1_${ik}$ ) = c( ji )
                       else
                          klnew = klnew + 1_${ik}$
                          if( klnew<=mmax ) then
                             ! eigenvalue in both intervals -- add upper to
                             ! queue.
                             ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ )
                             nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ )
                             ab( klnew, 1_${ik}$ ) = c( ji )
                             nab( klnew, 1_${ik}$ ) = iwork( ji )
                             ab( ji, 2_${ik}$ ) = c( ji )
                             nab( ji, 2_${ik}$ ) = iwork( ji )
                          else
                             info = mmax + 1_${ik}$
                          end if
                       end if
                    end do loop_70
                    if( info/=0 )return
                    kl = klnew
                 else
                    ! ijob=3: binary search.  keep only the interval containing
                            ! w   s.t. n(w) = nval
                    do ji = kf, kl
                       if( iwork( ji )<=nval( ji ) ) then
                          ab( ji, 1_${ik}$ ) = c( ji )
                          nab( ji, 1_${ik}$ ) = iwork( ji )
                       end if
                       if( iwork( ji )>=nval( ji ) ) then
                          ab( ji, 2_${ik}$ ) = c( ji )
                          nab( ji, 2_${ik}$ ) = iwork( ji )
                       end if
                    end do
                 end if
              else
                 ! end of parallel version of the loop
                 ! begin of serial version of the loop
                 klnew = kl
                 loop_100: do ji = kf, kl
                    ! compute n(w), the number of eigenvalues less than w
                    tmp1 = c( ji )
                    tmp2 = d( 1_${ik}$ ) - tmp1
                    itmp1 = 0_${ik}$
                    if( tmp2<=pivmin ) then
                       itmp1 = 1_${ik}$
                       tmp2 = min( tmp2, -pivmin )
                    end if
                    do j = 2, n
                       tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1
                       if( tmp2<=pivmin ) then
                          itmp1 = itmp1 + 1_${ik}$
                          tmp2 = min( tmp2, -pivmin )
                       end if
                    end do
                    if( ijob<=2_${ik}$ ) then
                       ! ijob=2: choose all intervals containing eigenvalues.
                       ! insure that n(w) is monotone
                       itmp1 = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), itmp1 ) )
                       ! update the queue -- add intervals if both halves
                       ! contain eigenvalues.
                       if( itmp1==nab( ji, 2_${ik}$ ) ) then
                          ! no eigenvalue in the upper interval:
                          ! just use the lower interval.
                          ab( ji, 2_${ik}$ ) = tmp1
                       else if( itmp1==nab( ji, 1_${ik}$ ) ) then
                          ! no eigenvalue in the lower interval:
                          ! just use the upper interval.
                          ab( ji, 1_${ik}$ ) = tmp1
                       else if( klnew<mmax ) then
                          ! eigenvalue in both intervals -- add upper to queue.
                          klnew = klnew + 1_${ik}$
                          ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ )
                          nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ )
                          ab( klnew, 1_${ik}$ ) = tmp1
                          nab( klnew, 1_${ik}$ ) = itmp1
                          ab( ji, 2_${ik}$ ) = tmp1
                          nab( ji, 2_${ik}$ ) = itmp1
                       else
                          info = mmax + 1_${ik}$
                          return
                       end if
                    else
                       ! ijob=3: binary search.  keep only the interval
                               ! containing  w  s.t. n(w) = nval
                       if( itmp1<=nval( ji ) ) then
                          ab( ji, 1_${ik}$ ) = tmp1
                          nab( ji, 1_${ik}$ ) = itmp1
                       end if
                       if( itmp1>=nval( ji ) ) then
                          ab( ji, 2_${ik}$ ) = tmp1
                          nab( ji, 2_${ik}$ ) = itmp1
                       end if
                    end if
                 end do loop_100
                 kl = klnew
              end if
              ! check for convergence
              kfnew = kf
              loop_110: do ji = kf, kl
                 tmp1 = abs( ab( ji, 2_${ik}$ )-ab( ji, 1_${ik}$ ) )
                 tmp2 = max( abs( ab( ji, 2_${ik}$ ) ), abs( ab( ji, 1_${ik}$ ) ) )
                 if( tmp1<max( abstol, pivmin, reltol*tmp2 ) .or.nab( ji, 1_${ik}$ )>=nab( ji, 2_${ik}$ ) ) &
                           then
                    ! converged -- swap with position kfnew,
                                 ! then increment kfnew
                    if( ji>kfnew ) then
                       tmp1 = ab( ji, 1_${ik}$ )
                       tmp2 = ab( ji, 2_${ik}$ )
                       itmp1 = nab( ji, 1_${ik}$ )
                       itmp2 = nab( ji, 2_${ik}$ )
                       ab( ji, 1_${ik}$ ) = ab( kfnew, 1_${ik}$ )
                       ab( ji, 2_${ik}$ ) = ab( kfnew, 2_${ik}$ )
                       nab( ji, 1_${ik}$ ) = nab( kfnew, 1_${ik}$ )
                       nab( ji, 2_${ik}$ ) = nab( kfnew, 2_${ik}$ )
                       ab( kfnew, 1_${ik}$ ) = tmp1
                       ab( kfnew, 2_${ik}$ ) = tmp2
                       nab( kfnew, 1_${ik}$ ) = itmp1
                       nab( kfnew, 2_${ik}$ ) = itmp2
                       if( ijob==3_${ik}$ ) then
                          itmp1 = nval( ji )
                          nval( ji ) = nval( kfnew )
                          nval( kfnew ) = itmp1
                       end if
                    end if
                    kfnew = kfnew + 1_${ik}$
                 end if
              end do loop_110
              kf = kfnew
              ! choose midpoints
              do ji = kf, kl
                 c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) )
              end do
              ! if no more intervals to refine, quit.
              if( kf>kl )go to 140
           end do loop_130
           ! converged
           140 continue
           info = max( kl+1-kf, 0_${ik}$ )
           mout = kl
           return
     end subroutine stdlib${ii}$_slaebz

     pure module subroutine stdlib${ii}$_dlaebz( ijob, nitmax, n, mmax, minp, nbmin, abstol,reltol, pivmin, d, &
     !! DLAEBZ contains the iteration loops which compute and use the
     !! function N(w), which is the count of eigenvalues of a symmetric
     !! tridiagonal matrix T less than or equal to its argument  w.  It
     !! performs a choice of two types of loops:
     !! IJOB=1, followed by
     !! IJOB=2: It takes as input a list of intervals and returns a list of
     !! sufficiently small intervals whose union contains the same
     !! eigenvalues as the union of the original intervals.
     !! The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
     !! The output interval (AB(j,1),AB(j,2)] will contain
     !! eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
     !! IJOB=3: It performs a binary search in each input interval
     !! (AB(j,1),AB(j,2)] for a point  w(j)  such that
     !! N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
     !! the search.  If such a w(j) is found, then on output
     !! AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
     !! (AB(j,1),AB(j,2)] will be a small interval containing the
     !! point where N(w) jumps through NVAL(j), unless that point
     !! lies outside the initial interval.
     !! Note that the intervals are in all cases half-open intervals,
     !! i.e., of the form  (a,b] , which includes  b  but not  a .
     !! To avoid underflow, the matrix should be scaled so that its largest
     !! element is no greater than  overflow**(1/2) * underflow**(1/4)
     !! in absolute value.  To assure the most accurate computation
     !! of small eigenvalues, the matrix should be scaled to be
     !! not much smaller than that, either.
     !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
     !! Matrix", Report CS41, Computer Science Dept., Stanford
     !! University, July 21, 1966
     !! Note: the arguments are, in general, *not* checked for unreasonable
     !! values.
               e, e2, nval, ab, c, mout,nab, work, iwork, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ijob, minp, mmax, n, nbmin, nitmax
           integer(${ik}$), intent(out) :: info, mout
           real(dp), intent(in) :: abstol, pivmin, reltol
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           integer(${ik}$), intent(inout) :: nab(mmax,*), nval(*)
           real(dp), intent(inout) :: ab(mmax,*), c(*)
           real(dp), intent(in) :: d(*), e(*), e2(*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl, klnew
           real(dp) :: tmp1, tmp2
           ! Intrinsic Functions 
           ! Executable Statements 
           ! check for errors
           info = 0_${ik}$
           if( ijob<1_${ik}$ .or. ijob>3_${ik}$ ) then
              info = -1_${ik}$
              return
           end if
           ! initialize nab
           if( ijob==1_${ik}$ ) then
              ! compute the number of eigenvalues in the initial intervals.
              mout = 0_${ik}$
              do ji = 1, minp
                 do jp = 1, 2
                    tmp1 = d( 1_${ik}$ ) - ab( ji, jp )
                    if( abs( tmp1 )<pivmin )tmp1 = -pivmin
                    nab( ji, jp ) = 0_${ik}$
                    if( tmp1<=zero )nab( ji, jp ) = 1_${ik}$
                    do j = 2, n
                       tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp )
                       if( abs( tmp1 )<pivmin )tmp1 = -pivmin
                       if( tmp1<=zero )nab( ji, jp ) = nab( ji, jp ) + 1_${ik}$
                    end do
                 end do
                 mout = mout + nab( ji, 2_${ik}$ ) - nab( ji, 1_${ik}$ )
              end do
              return
           end if
           ! initialize for loop
           ! kf and kl have the following meaning:
              ! intervals 1,...,kf-1 have converged.
              ! intervals kf,...,kl  still need to be refined.
           kf = 1_${ik}$
           kl = minp
           ! if ijob=2, initialize c.
           ! if ijob=3, use the user-supplied starting point.
           if( ijob==2_${ik}$ ) then
              do ji = 1, minp
                 c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) )
              end do
           end if
           ! iteration loop
           loop_130: do jit = 1, nitmax
              ! loop over intervals
              if( kl-kf+1>=nbmin .and. nbmin>0_${ik}$ ) then
                 ! begin of parallel version of the loop
                 do ji = kf, kl
                    ! compute n(c), the number of eigenvalues less than c
                    work( ji ) = d( 1_${ik}$ ) - c( ji )
                    iwork( ji ) = 0_${ik}$
                    if( work( ji )<=pivmin ) then
                       iwork( ji ) = 1_${ik}$
                       work( ji ) = min( work( ji ), -pivmin )
                    end if
                    do j = 2, n
                       work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji )
                       if( work( ji )<=pivmin ) then
                          iwork( ji ) = iwork( ji ) + 1_${ik}$
                          work( ji ) = min( work( ji ), -pivmin )
                       end if
                    end do
                 end do
                 if( ijob<=2_${ik}$ ) then
                    ! ijob=2: choose all intervals containing eigenvalues.
                    klnew = kl
                    loop_70: do ji = kf, kl
                       ! insure that n(w) is monotone
                       iwork( ji ) = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), iwork( ji ) ) )
                       ! update the queue -- add intervals if both halves
                       ! contain eigenvalues.
                       if( iwork( ji )==nab( ji, 2_${ik}$ ) ) then
                          ! no eigenvalue in the upper interval:
                          ! just use the lower interval.
                          ab( ji, 2_${ik}$ ) = c( ji )
                       else if( iwork( ji )==nab( ji, 1_${ik}$ ) ) then
                          ! no eigenvalue in the lower interval:
                          ! just use the upper interval.
                          ab( ji, 1_${ik}$ ) = c( ji )
                       else
                          klnew = klnew + 1_${ik}$
                          if( klnew<=mmax ) then
                             ! eigenvalue in both intervals -- add upper to
                             ! queue.
                             ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ )
                             nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ )
                             ab( klnew, 1_${ik}$ ) = c( ji )
                             nab( klnew, 1_${ik}$ ) = iwork( ji )
                             ab( ji, 2_${ik}$ ) = c( ji )
                             nab( ji, 2_${ik}$ ) = iwork( ji )
                          else
                             info = mmax + 1_${ik}$
                          end if
                       end if
                    end do loop_70
                    if( info/=0 )return
                    kl = klnew
                 else
                    ! ijob=3: binary search.  keep only the interval containing
                            ! w   s.t. n(w) = nval
                    do ji = kf, kl
                       if( iwork( ji )<=nval( ji ) ) then
                          ab( ji, 1_${ik}$ ) = c( ji )
                          nab( ji, 1_${ik}$ ) = iwork( ji )
                       end if
                       if( iwork( ji )>=nval( ji ) ) then
                          ab( ji, 2_${ik}$ ) = c( ji )
                          nab( ji, 2_${ik}$ ) = iwork( ji )
                       end if
                    end do
                 end if
              else
                 ! end of parallel version of the loop
                 ! begin of serial version of the loop
                 klnew = kl
                 loop_100: do ji = kf, kl
                    ! compute n(w), the number of eigenvalues less than w
                    tmp1 = c( ji )
                    tmp2 = d( 1_${ik}$ ) - tmp1
                    itmp1 = 0_${ik}$
                    if( tmp2<=pivmin ) then
                       itmp1 = 1_${ik}$
                       tmp2 = min( tmp2, -pivmin )
                    end if
                    do j = 2, n
                       tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1
                       if( tmp2<=pivmin ) then
                          itmp1 = itmp1 + 1_${ik}$
                          tmp2 = min( tmp2, -pivmin )
                       end if
                    end do
                    if( ijob<=2_${ik}$ ) then
                       ! ijob=2: choose all intervals containing eigenvalues.
                       ! insure that n(w) is monotone
                       itmp1 = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), itmp1 ) )
                       ! update the queue -- add intervals if both halves
                       ! contain eigenvalues.
                       if( itmp1==nab( ji, 2_${ik}$ ) ) then
                          ! no eigenvalue in the upper interval:
                          ! just use the lower interval.
                          ab( ji, 2_${ik}$ ) = tmp1
                       else if( itmp1==nab( ji, 1_${ik}$ ) ) then
                          ! no eigenvalue in the lower interval:
                          ! just use the upper interval.
                          ab( ji, 1_${ik}$ ) = tmp1
                       else if( klnew<mmax ) then
                          ! eigenvalue in both intervals -- add upper to queue.
                          klnew = klnew + 1_${ik}$
                          ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ )
                          nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ )
                          ab( klnew, 1_${ik}$ ) = tmp1
                          nab( klnew, 1_${ik}$ ) = itmp1
                          ab( ji, 2_${ik}$ ) = tmp1
                          nab( ji, 2_${ik}$ ) = itmp1
                       else
                          info = mmax + 1_${ik}$
                          return
                       end if
                    else
                       ! ijob=3: binary search.  keep only the interval
                               ! containing  w  s.t. n(w) = nval
                       if( itmp1<=nval( ji ) ) then
                          ab( ji, 1_${ik}$ ) = tmp1
                          nab( ji, 1_${ik}$ ) = itmp1
                       end if
                       if( itmp1>=nval( ji ) ) then
                          ab( ji, 2_${ik}$ ) = tmp1
                          nab( ji, 2_${ik}$ ) = itmp1
                       end if
                    end if
                 end do loop_100
                 kl = klnew
              end if
              ! check for convergence
              kfnew = kf
              loop_110: do ji = kf, kl
                 tmp1 = abs( ab( ji, 2_${ik}$ )-ab( ji, 1_${ik}$ ) )
                 tmp2 = max( abs( ab( ji, 2_${ik}$ ) ), abs( ab( ji, 1_${ik}$ ) ) )
                 if( tmp1<max( abstol, pivmin, reltol*tmp2 ) .or.nab( ji, 1_${ik}$ )>=nab( ji, 2_${ik}$ ) ) &
                           then
                    ! converged -- swap with position kfnew,
                                 ! then increment kfnew
                    if( ji>kfnew ) then
                       tmp1 = ab( ji, 1_${ik}$ )
                       tmp2 = ab( ji, 2_${ik}$ )
                       itmp1 = nab( ji, 1_${ik}$ )
                       itmp2 = nab( ji, 2_${ik}$ )
                       ab( ji, 1_${ik}$ ) = ab( kfnew, 1_${ik}$ )
                       ab( ji, 2_${ik}$ ) = ab( kfnew, 2_${ik}$ )
                       nab( ji, 1_${ik}$ ) = nab( kfnew, 1_${ik}$ )
                       nab( ji, 2_${ik}$ ) = nab( kfnew, 2_${ik}$ )
                       ab( kfnew, 1_${ik}$ ) = tmp1
                       ab( kfnew, 2_${ik}$ ) = tmp2
                       nab( kfnew, 1_${ik}$ ) = itmp1
                       nab( kfnew, 2_${ik}$ ) = itmp2
                       if( ijob==3_${ik}$ ) then
                          itmp1 = nval( ji )
                          nval( ji ) = nval( kfnew )
                          nval( kfnew ) = itmp1
                       end if
                    end if
                    kfnew = kfnew + 1_${ik}$
                 end if
              end do loop_110
              kf = kfnew
              ! choose midpoints
              do ji = kf, kl
                 c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) )
              end do
              ! if no more intervals to refine, quit.
              if( kf>kl )go to 140
           end do loop_130
           ! converged
           140 continue
           info = max( kl+1-kf, 0_${ik}$ )
           mout = kl
           return
     end subroutine stdlib${ii}$_dlaebz

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$laebz( ijob, nitmax, n, mmax, minp, nbmin, abstol,reltol, pivmin, d, &
     !! DLAEBZ: contains the iteration loops which compute and use the
     !! function N(w), which is the count of eigenvalues of a symmetric
     !! tridiagonal matrix T less than or equal to its argument  w.  It
     !! performs a choice of two types of loops:
     !! IJOB=1, followed by
     !! IJOB=2: It takes as input a list of intervals and returns a list of
     !! sufficiently small intervals whose union contains the same
     !! eigenvalues as the union of the original intervals.
     !! The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
     !! The output interval (AB(j,1),AB(j,2)] will contain
     !! eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
     !! IJOB=3: It performs a binary search in each input interval
     !! (AB(j,1),AB(j,2)] for a point  w(j)  such that
     !! N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
     !! the search.  If such a w(j) is found, then on output
     !! AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
     !! (AB(j,1),AB(j,2)] will be a small interval containing the
     !! point where N(w) jumps through NVAL(j), unless that point
     !! lies outside the initial interval.
     !! Note that the intervals are in all cases half-open intervals,
     !! i.e., of the form  (a,b] , which includes  b  but not  a .
     !! To avoid underflow, the matrix should be scaled so that its largest
     !! element is no greater than  overflow**(1/2) * underflow**(1/4)
     !! in absolute value.  To assure the most accurate computation
     !! of small eigenvalues, the matrix should be scaled to be
     !! not much smaller than that, either.
     !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
     !! Matrix", Report CS41, Computer Science Dept., Stanford
     !! University, July 21, 1966
     !! Note: the arguments are, in general, *not* checked for unreasonable
     !! values.
               e, e2, nval, ab, c, mout,nab, work, iwork, info )
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: ijob, minp, mmax, n, nbmin, nitmax
           integer(${ik}$), intent(out) :: info, mout
           real(${rk}$), intent(in) :: abstol, pivmin, reltol
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           integer(${ik}$), intent(inout) :: nab(mmax,*), nval(*)
           real(${rk}$), intent(inout) :: ab(mmax,*), c(*)
           real(${rk}$), intent(in) :: d(*), e(*), e2(*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl, klnew
           real(${rk}$) :: tmp1, tmp2
           ! Intrinsic Functions 
           ! Executable Statements 
           ! check for errors
           info = 0_${ik}$
           if( ijob<1_${ik}$ .or. ijob>3_${ik}$ ) then
              info = -1_${ik}$
              return
           end if
           ! initialize nab
           if( ijob==1_${ik}$ ) then
              ! compute the number of eigenvalues in the initial intervals.
              mout = 0_${ik}$
              do ji = 1, minp
                 do jp = 1, 2
                    tmp1 = d( 1_${ik}$ ) - ab( ji, jp )
                    if( abs( tmp1 )<pivmin )tmp1 = -pivmin
                    nab( ji, jp ) = 0_${ik}$
                    if( tmp1<=zero )nab( ji, jp ) = 1_${ik}$
                    do j = 2, n
                       tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp )
                       if( abs( tmp1 )<pivmin )tmp1 = -pivmin
                       if( tmp1<=zero )nab( ji, jp ) = nab( ji, jp ) + 1_${ik}$
                    end do
                 end do
                 mout = mout + nab( ji, 2_${ik}$ ) - nab( ji, 1_${ik}$ )
              end do
              return
           end if
           ! initialize for loop
           ! kf and kl have the following meaning:
              ! intervals 1,...,kf-1 have converged.
              ! intervals kf,...,kl  still need to be refined.
           kf = 1_${ik}$
           kl = minp
           ! if ijob=2, initialize c.
           ! if ijob=3, use the user-supplied starting point.
           if( ijob==2_${ik}$ ) then
              do ji = 1, minp
                 c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) )
              end do
           end if
           ! iteration loop
           loop_130: do jit = 1, nitmax
              ! loop over intervals
              if( kl-kf+1>=nbmin .and. nbmin>0_${ik}$ ) then
                 ! begin of parallel version of the loop
                 do ji = kf, kl
                    ! compute n(c), the number of eigenvalues less than c
                    work( ji ) = d( 1_${ik}$ ) - c( ji )
                    iwork( ji ) = 0_${ik}$
                    if( work( ji )<=pivmin ) then
                       iwork( ji ) = 1_${ik}$
                       work( ji ) = min( work( ji ), -pivmin )
                    end if
                    do j = 2, n
                       work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji )
                       if( work( ji )<=pivmin ) then
                          iwork( ji ) = iwork( ji ) + 1_${ik}$
                          work( ji ) = min( work( ji ), -pivmin )
                       end if
                    end do
                 end do
                 if( ijob<=2_${ik}$ ) then
                    ! ijob=2: choose all intervals containing eigenvalues.
                    klnew = kl
                    loop_70: do ji = kf, kl
                       ! insure that n(w) is monotone
                       iwork( ji ) = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), iwork( ji ) ) )
                       ! update the queue -- add intervals if both halves
                       ! contain eigenvalues.
                       if( iwork( ji )==nab( ji, 2_${ik}$ ) ) then
                          ! no eigenvalue in the upper interval:
                          ! just use the lower interval.
                          ab( ji, 2_${ik}$ ) = c( ji )
                       else if( iwork( ji )==nab( ji, 1_${ik}$ ) ) then
                          ! no eigenvalue in the lower interval:
                          ! just use the upper interval.
                          ab( ji, 1_${ik}$ ) = c( ji )
                       else
                          klnew = klnew + 1_${ik}$
                          if( klnew<=mmax ) then
                             ! eigenvalue in both intervals -- add upper to
                             ! queue.
                             ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ )
                             nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ )
                             ab( klnew, 1_${ik}$ ) = c( ji )
                             nab( klnew, 1_${ik}$ ) = iwork( ji )
                             ab( ji, 2_${ik}$ ) = c( ji )
                             nab( ji, 2_${ik}$ ) = iwork( ji )
                          else
                             info = mmax + 1_${ik}$
                          end if
                       end if
                    end do loop_70
                    if( info/=0 )return
                    kl = klnew
                 else
                    ! ijob=3: binary search.  keep only the interval containing
                            ! w   s.t. n(w) = nval
                    do ji = kf, kl
                       if( iwork( ji )<=nval( ji ) ) then
                          ab( ji, 1_${ik}$ ) = c( ji )
                          nab( ji, 1_${ik}$ ) = iwork( ji )
                       end if
                       if( iwork( ji )>=nval( ji ) ) then
                          ab( ji, 2_${ik}$ ) = c( ji )
                          nab( ji, 2_${ik}$ ) = iwork( ji )
                       end if
                    end do
                 end if
              else
                 ! end of parallel version of the loop
                 ! begin of serial version of the loop
                 klnew = kl
                 loop_100: do ji = kf, kl
                    ! compute n(w), the number of eigenvalues less than w
                    tmp1 = c( ji )
                    tmp2 = d( 1_${ik}$ ) - tmp1
                    itmp1 = 0_${ik}$
                    if( tmp2<=pivmin ) then
                       itmp1 = 1_${ik}$
                       tmp2 = min( tmp2, -pivmin )
                    end if
                    do j = 2, n
                       tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1
                       if( tmp2<=pivmin ) then
                          itmp1 = itmp1 + 1_${ik}$
                          tmp2 = min( tmp2, -pivmin )
                       end if
                    end do
                    if( ijob<=2_${ik}$ ) then
                       ! ijob=2: choose all intervals containing eigenvalues.
                       ! insure that n(w) is monotone
                       itmp1 = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), itmp1 ) )
                       ! update the queue -- add intervals if both halves
                       ! contain eigenvalues.
                       if( itmp1==nab( ji, 2_${ik}$ ) ) then
                          ! no eigenvalue in the upper interval:
                          ! just use the lower interval.
                          ab( ji, 2_${ik}$ ) = tmp1
                       else if( itmp1==nab( ji, 1_${ik}$ ) ) then
                          ! no eigenvalue in the lower interval:
                          ! just use the upper interval.
                          ab( ji, 1_${ik}$ ) = tmp1
                       else if( klnew<mmax ) then
                          ! eigenvalue in both intervals -- add upper to queue.
                          klnew = klnew + 1_${ik}$
                          ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ )
                          nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ )
                          ab( klnew, 1_${ik}$ ) = tmp1
                          nab( klnew, 1_${ik}$ ) = itmp1
                          ab( ji, 2_${ik}$ ) = tmp1
                          nab( ji, 2_${ik}$ ) = itmp1
                       else
                          info = mmax + 1_${ik}$
                          return
                       end if
                    else
                       ! ijob=3: binary search.  keep only the interval
                               ! containing  w  s.t. n(w) = nval
                       if( itmp1<=nval( ji ) ) then
                          ab( ji, 1_${ik}$ ) = tmp1
                          nab( ji, 1_${ik}$ ) = itmp1
                       end if
                       if( itmp1>=nval( ji ) ) then
                          ab( ji, 2_${ik}$ ) = tmp1
                          nab( ji, 2_${ik}$ ) = itmp1
                       end if
                    end if
                 end do loop_100
                 kl = klnew
              end if
              ! check for convergence
              kfnew = kf
              loop_110: do ji = kf, kl
                 tmp1 = abs( ab( ji, 2_${ik}$ )-ab( ji, 1_${ik}$ ) )
                 tmp2 = max( abs( ab( ji, 2_${ik}$ ) ), abs( ab( ji, 1_${ik}$ ) ) )
                 if( tmp1<max( abstol, pivmin, reltol*tmp2 ) .or.nab( ji, 1_${ik}$ )>=nab( ji, 2_${ik}$ ) ) &
                           then
                    ! converged -- swap with position kfnew,
                                 ! then increment kfnew
                    if( ji>kfnew ) then
                       tmp1 = ab( ji, 1_${ik}$ )
                       tmp2 = ab( ji, 2_${ik}$ )
                       itmp1 = nab( ji, 1_${ik}$ )
                       itmp2 = nab( ji, 2_${ik}$ )
                       ab( ji, 1_${ik}$ ) = ab( kfnew, 1_${ik}$ )
                       ab( ji, 2_${ik}$ ) = ab( kfnew, 2_${ik}$ )
                       nab( ji, 1_${ik}$ ) = nab( kfnew, 1_${ik}$ )
                       nab( ji, 2_${ik}$ ) = nab( kfnew, 2_${ik}$ )
                       ab( kfnew, 1_${ik}$ ) = tmp1
                       ab( kfnew, 2_${ik}$ ) = tmp2
                       nab( kfnew, 1_${ik}$ ) = itmp1
                       nab( kfnew, 2_${ik}$ ) = itmp2
                       if( ijob==3_${ik}$ ) then
                          itmp1 = nval( ji )
                          nval( ji ) = nval( kfnew )
                          nval( kfnew ) = itmp1
                       end if
                    end if
                    kfnew = kfnew + 1_${ik}$
                 end if
              end do loop_110
              kf = kfnew
              ! choose midpoints
              do ji = kf, kl
                 c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) )
              end do
              ! if no more intervals to refine, quit.
              if( kf>kl )go to 140
           end do loop_130
           ! converged
           140 continue
           info = max( kl+1-kf, 0_${ik}$ )
           mout = kl
           return
     end subroutine stdlib${ii}$_${ri}$laebz

#:endif
#:endfor



     pure integer(${ik}$) module function stdlib${ii}$_slaneg( n, d, lld, sigma, pivmin, r )
     !! SLANEG computes the Sturm count, the number of negative pivots
     !! encountered while factoring tridiagonal T - sigma I = L D L^T.
     !! This implementation works directly on the factors without forming
     !! the tridiagonal matrix T.  The Sturm count is also the number of
     !! eigenvalues of T less than sigma.
     !! This routine is called from SLARRB.
     !! The current routine does not use the PIVMIN parameter but rather
     !! requires IEEE-754 propagation of Infinities and NaNs.  This
     !! routine also has no input range restrictions but does require
     !! default exception handling such that x/0 produces Inf when x is
     !! non-zero, and Inf/Inf produces NaN.  For more information, see:
     !! Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
     !! Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
     !! Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
     !! (Tech report version in LAWN 172 with the same title.)
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: n, r
           real(sp), intent(in) :: pivmin, sigma
           ! Array Arguments 
           real(sp), intent(in) :: d(*), lld(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: blklen = 128_${ik}$
           
           ! some architectures propagate infinities and nans very slowly, so
           ! the code computes counts in blklen chunks.  then a nan can
           ! propagate at most blklen columns before being detected.  this is
           ! not a general tuning parameter; it needs only to be just large
           ! enough that the overhead is tiny in common cases.
           
           ! Local Scalars 
           integer(${ik}$) :: bj, j, neg1, neg2, negcnt
           real(sp) :: bsav, dminus, dplus, gamma, p, t, tmp
           logical(lk) :: sawnan
           ! Intrinsic Functions 
           ! Executable Statements 
           negcnt = 0_${ik}$
           ! i) upper part: l d l^t - sigma i = l+ d+ l+^t
           t = -sigma
           loop_210: do bj = 1, r-1, blklen
              neg1 = 0_${ik}$
              bsav = t
              do j = bj, min(bj+blklen-1, r-1)
                 dplus = d( j ) + t
                 if( dplus<zero ) neg1 = neg1 + 1_${ik}$
                 tmp = t / dplus
                 t = tmp * lld( j ) - sigma
              end do
              sawnan = stdlib${ii}$_sisnan( t )
           ! run a slower version of the above loop if a nan is detected.
           ! a nan should occur only with a zero pivot after an infinite
           ! pivot.  in that case, substituting 1 for t/dplus is the
           ! correct limit.
              if( sawnan ) then
                 neg1 = 0_${ik}$
                 t = bsav
                 do j = bj, min(bj+blklen-1, r-1)
                    dplus = d( j ) + t
                    if( dplus<zero ) neg1 = neg1 + 1_${ik}$
                    tmp = t / dplus
                    if (stdlib${ii}$_sisnan(tmp)) tmp = one
                    t = tmp * lld(j) - sigma
                 end do
              end if
              negcnt = negcnt + neg1
           end do loop_210
           ! ii) lower part: l d l^t - sigma i = u- d- u-^t
           p = d( n ) - sigma
           do bj = n-1, r, -blklen
              neg2 = 0_${ik}$
              bsav = p
              do j = bj, max(bj-blklen+1, r), -1
                 dminus = lld( j ) + p
                 if( dminus<zero ) neg2 = neg2 + 1_${ik}$
                 tmp = p / dminus
                 p = tmp * d( j ) - sigma
              end do
              sawnan = stdlib${ii}$_sisnan( p )
           ! as above, run a slower version that substitutes 1 for inf/inf.
              if( sawnan ) then
                 neg2 = 0_${ik}$
                 p = bsav
                 do j = bj, max(bj-blklen+1, r), -1
                    dminus = lld( j ) + p
                    if( dminus<zero ) neg2 = neg2 + 1_${ik}$
                    tmp = p / dminus
                    if (stdlib${ii}$_sisnan(tmp)) tmp = one
                    p = tmp * d(j) - sigma
                 end do
              end if
              negcnt = negcnt + neg2
           end do
           ! iii) twist index
             ! t was shifted by sigma initially.
           gamma = (t + sigma) + p
           if( gamma<zero ) negcnt = negcnt+1
           stdlib${ii}$_slaneg = negcnt
     end function stdlib${ii}$_slaneg

     pure integer(${ik}$) module function stdlib${ii}$_dlaneg( n, d, lld, sigma, pivmin, r )
     !! DLANEG computes the Sturm count, the number of negative pivots
     !! encountered while factoring tridiagonal T - sigma I = L D L^T.
     !! This implementation works directly on the factors without forming
     !! the tridiagonal matrix T.  The Sturm count is also the number of
     !! eigenvalues of T less than sigma.
     !! This routine is called from DLARRB.
     !! The current routine does not use the PIVMIN parameter but rather
     !! requires IEEE-754 propagation of Infinities and NaNs.  This
     !! routine also has no input range restrictions but does require
     !! default exception handling such that x/0 produces Inf when x is
     !! non-zero, and Inf/Inf produces NaN.  For more information, see:
     !! Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
     !! Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
     !! Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
     !! (Tech report version in LAWN 172 with the same title.)
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: n, r
           real(dp), intent(in) :: pivmin, sigma
           ! Array Arguments 
           real(dp), intent(in) :: d(*), lld(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: blklen = 128_${ik}$
           
           ! some architectures propagate infinities and nans very slowly, so
           ! the code computes counts in blklen chunks.  then a nan can
           ! propagate at most blklen columns before being detected.  this is
           ! not a general tuning parameter; it needs only to be just large
           ! enough that the overhead is tiny in common cases.
           
           ! Local Scalars 
           integer(${ik}$) :: bj, j, neg1, neg2, negcnt
           real(dp) :: bsav, dminus, dplus, gamma, p, t, tmp
           logical(lk) :: sawnan
           ! Intrinsic Functions 
           ! Executable Statements 
           negcnt = 0_${ik}$
           ! i) upper part: l d l^t - sigma i = l+ d+ l+^t
           t = -sigma
           loop_210: do bj = 1, r-1, blklen
              neg1 = 0_${ik}$
              bsav = t
              do j = bj, min(bj+blklen-1, r-1)
                 dplus = d( j ) + t
                 if( dplus<zero ) neg1 = neg1 + 1_${ik}$
                 tmp = t / dplus
                 t = tmp * lld( j ) - sigma
              end do
              sawnan = stdlib${ii}$_disnan( t )
           ! run a slower version of the above loop if a nan is detected.
           ! a nan should occur only with a zero pivot after an infinite
           ! pivot.  in that case, substituting 1 for t/dplus is the
           ! correct limit.
              if( sawnan ) then
                 neg1 = 0_${ik}$
                 t = bsav
                 do j = bj, min(bj+blklen-1, r-1)
                    dplus = d( j ) + t
                    if( dplus<zero ) neg1 = neg1 + 1_${ik}$
                    tmp = t / dplus
                    if (stdlib${ii}$_disnan(tmp)) tmp = one
                    t = tmp * lld(j) - sigma
                 end do
              end if
              negcnt = negcnt + neg1
           end do loop_210
           ! ii) lower part: l d l^t - sigma i = u- d- u-^t
           p = d( n ) - sigma
           do bj = n-1, r, -blklen
              neg2 = 0_${ik}$
              bsav = p
              do j = bj, max(bj-blklen+1, r), -1
                 dminus = lld( j ) + p
                 if( dminus<zero ) neg2 = neg2 + 1_${ik}$
                 tmp = p / dminus
                 p = tmp * d( j ) - sigma
              end do
              sawnan = stdlib${ii}$_disnan( p )
           ! as above, run a slower version that substitutes 1 for inf/inf.
              if( sawnan ) then
                 neg2 = 0_${ik}$
                 p = bsav
                 do j = bj, max(bj-blklen+1, r), -1
                    dminus = lld( j ) + p
                    if( dminus<zero ) neg2 = neg2 + 1_${ik}$
                    tmp = p / dminus
                    if (stdlib${ii}$_disnan(tmp)) tmp = one
                    p = tmp * d(j) - sigma
                 end do
              end if
              negcnt = negcnt + neg2
           end do
           ! iii) twist index
             ! t was shifted by sigma initially.
           gamma = (t + sigma) + p
           if( gamma<zero ) negcnt = negcnt+1
           stdlib${ii}$_dlaneg = negcnt
     end function stdlib${ii}$_dlaneg

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure integer(${ik}$) module function stdlib${ii}$_${ri}$laneg( n, d, lld, sigma, pivmin, r )
     !! DLANEG: computes the Sturm count, the number of negative pivots
     !! encountered while factoring tridiagonal T - sigma I = L D L^T.
     !! This implementation works directly on the factors without forming
     !! the tridiagonal matrix T.  The Sturm count is also the number of
     !! eigenvalues of T less than sigma.
     !! This routine is called from DLARRB.
     !! The current routine does not use the PIVMIN parameter but rather
     !! requires IEEE-754 propagation of Infinities and NaNs.  This
     !! routine also has no input range restrictions but does require
     !! default exception handling such that x/0 produces Inf when x is
     !! non-zero, and Inf/Inf produces NaN.  For more information, see:
     !! Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
     !! Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
     !! Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
     !! (Tech report version in LAWN 172 with the same title.)
        ! -- lapack auxiliary routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: n, r
           real(${rk}$), intent(in) :: pivmin, sigma
           ! Array Arguments 
           real(${rk}$), intent(in) :: d(*), lld(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: blklen = 128_${ik}$
           
           ! some architectures propagate infinities and nans very slowly, so
           ! the code computes counts in blklen chunks.  then a nan can
           ! propagate at most blklen columns before being detected.  this is
           ! not a general tuning parameter; it needs only to be just large
           ! enough that the overhead is tiny in common cases.
           
           ! Local Scalars 
           integer(${ik}$) :: bj, j, neg1, neg2, negcnt
           real(${rk}$) :: bsav, dminus, dplus, gamma, p, t, tmp
           logical(lk) :: sawnan
           ! Intrinsic Functions 
           ! Executable Statements 
           negcnt = 0_${ik}$
           ! i) upper part: l d l^t - sigma i = l+ d+ l+^t
           t = -sigma
           loop_210: do bj = 1, r-1, blklen
              neg1 = 0_${ik}$
              bsav = t
              do j = bj, min(bj+blklen-1, r-1)
                 dplus = d( j ) + t
                 if( dplus<zero ) neg1 = neg1 + 1_${ik}$
                 tmp = t / dplus
                 t = tmp * lld( j ) - sigma
              end do
              sawnan = stdlib${ii}$_${ri}$isnan( t )
           ! run a slower version of the above loop if a nan is detected.
           ! a nan should occur only with a zero pivot after an infinite
           ! pivot.  in that case, substituting 1 for t/dplus is the
           ! correct limit.
              if( sawnan ) then
                 neg1 = 0_${ik}$
                 t = bsav
                 do j = bj, min(bj+blklen-1, r-1)
                    dplus = d( j ) + t
                    if( dplus<zero ) neg1 = neg1 + 1_${ik}$
                    tmp = t / dplus
                    if (stdlib${ii}$_${ri}$isnan(tmp)) tmp = one
                    t = tmp * lld(j) - sigma
                 end do
              end if
              negcnt = negcnt + neg1
           end do loop_210
           ! ii) lower part: l d l^t - sigma i = u- d- u-^t
           p = d( n ) - sigma
           do bj = n-1, r, -blklen
              neg2 = 0_${ik}$
              bsav = p
              do j = bj, max(bj-blklen+1, r), -1
                 dminus = lld( j ) + p
                 if( dminus<zero ) neg2 = neg2 + 1_${ik}$
                 tmp = p / dminus
                 p = tmp * d( j ) - sigma
              end do
              sawnan = stdlib${ii}$_${ri}$isnan( p )
           ! as above, run a slower version that substitutes 1 for inf/inf.
              if( sawnan ) then
                 neg2 = 0_${ik}$
                 p = bsav
                 do j = bj, max(bj-blklen+1, r), -1
                    dminus = lld( j ) + p
                    if( dminus<zero ) neg2 = neg2 + 1_${ik}$
                    tmp = p / dminus
                    if (stdlib${ii}$_${ri}$isnan(tmp)) tmp = one
                    p = tmp * d(j) - sigma
                 end do
              end if
              negcnt = negcnt + neg2
           end do
           ! iii) twist index
             ! t was shifted by sigma initially.
           gamma = (t + sigma) + p
           if( gamma<zero ) negcnt = negcnt+1
           stdlib${ii}$_${ri}$laneg = negcnt
     end function stdlib${ii}$_${ri}$laneg

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slaed0( icompq, qsiz, n, d, e, q, ldq, qstore, ldqs,work, iwork, info &
     !! SLAED0 computes all eigenvalues and corresponding eigenvectors of a
     !! symmetric tridiagonal matrix using the divide and conquer method.
               )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: icompq, ldq, ldqs, n, qsiz
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: d(*), e(*), q(ldq,*)
           real(sp), intent(out) :: qstore(ldqs,*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, &
           iq, iqptr, iwrem, j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, subpbs, &
                     tlvls
           real(sp) :: temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<0_${ik}$ .or. icompq>2_${ik}$ ) then
              info = -1_${ik}$
           else if( ( icompq==1_${ik}$ ) .and. ( qsiz<max( 0_${ik}$, n ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldqs<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SLAED0', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'SLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ )
           ! determine the size and placement of the submatrices, and save in
           ! the leading elements of iwork.
           iwork( 1_${ik}$ ) = n
           subpbs = 1_${ik}$
           tlvls = 0_${ik}$
           10 continue
           if( iwork( subpbs )>smlsiz ) then
              do j = subpbs, 1, -1
                 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$
                 iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$
              end do
              tlvls = tlvls + 1_${ik}$
              subpbs = 2_${ik}$*subpbs
              go to 10
           end if
           do j = 2, subpbs
              iwork( j ) = iwork( j ) + iwork( j-1 )
           end do
           ! divide the matrix into subpbs submatrices of size at most smlsiz+1
           ! using rank-1 modifications (cuts).
           spm1 = subpbs - 1_${ik}$
           do i = 1, spm1
              submat = iwork( i ) + 1_${ik}$
              smm1 = submat - 1_${ik}$
              d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
              d( submat ) = d( submat ) - abs( e( smm1 ) )
           end do
           indxq = 4_${ik}$*n + 3_${ik}$
           if( icompq/=2_${ik}$ ) then
              ! set up workspaces for eigenvalues only/accumulate new vectors
              ! routine
              temp = log( real( n,KIND=sp) ) / log( two )
              lgn = int( temp,KIND=${ik}$)
              if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
              if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
              iprmpt = indxq + n + 1_${ik}$
              iperm = iprmpt + n*lgn
              iqptr = iperm + n*lgn
              igivpt = iqptr + n + 2_${ik}$
              igivcl = igivpt + n*lgn
              igivnm = 1_${ik}$
              iq = igivnm + 2_${ik}$*n*lgn
              iwrem = iq + n**2_${ik}$ + 1_${ik}$
              ! initialize pointers
              do i = 0, subpbs
                 iwork( iprmpt+i ) = 1_${ik}$
                 iwork( igivpt+i ) = 1_${ik}$
              end do
              iwork( iqptr ) = 1_${ik}$
           end if
           ! solve each submatrix eigenproblem at the bottom of the divide and
           ! conquer tree.
           curr = 0_${ik}$
           loop_70: do i = 0, spm1
              if( i==0_${ik}$ ) then
                 submat = 1_${ik}$
                 matsiz = iwork( 1_${ik}$ )
              else
                 submat = iwork( i ) + 1_${ik}$
                 matsiz = iwork( i+1 ) - iwork( i )
              end if
              if( icompq==2_${ik}$ ) then
                 call stdlib${ii}$_ssteqr( 'I', matsiz, d( submat ), e( submat ),q( submat, submat ), &
                           ldq, work, info )
                 if( info/=0 )go to 130
              else
                 call stdlib${ii}$_ssteqr( 'I', matsiz, d( submat ), e( submat ),work( iq-1+iwork( &
                           iqptr+curr ) ), matsiz, work,info )
                 if( info/=0 )go to 130
                 if( icompq==1_${ik}$ ) then
                    call stdlib${ii}$_sgemm( 'N', 'N', qsiz, matsiz, matsiz, one,q( 1_${ik}$, submat ), ldq, &
                    work( iq-1+iwork( iqptr+curr ) ), matsiz, zero, qstore( 1_${ik}$, submat ),ldqs )
                              
                 end if
                 iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$
                 curr = curr + 1_${ik}$
              end if
              k = 1_${ik}$
              do j = submat, iwork( i+1 )
                 iwork( indxq+j ) = k
                 k = k + 1_${ik}$
              end do
           end do loop_70
           ! successively merge eigensystems of adjacent submatrices
           ! into eigensystem for the corresponding larger matrix.
           ! while ( subpbs > 1 )
           curlvl = 1_${ik}$
           80 continue
           if( subpbs>1_${ik}$ ) then
              spm2 = subpbs - 2_${ik}$
              loop_90: do i = 0, spm2, 2
                 if( i==0_${ik}$ ) then
                    submat = 1_${ik}$
                    matsiz = iwork( 2_${ik}$ )
                    msd2 = iwork( 1_${ik}$ )
                    curprb = 0_${ik}$
                 else
                    submat = iwork( i ) + 1_${ik}$
                    matsiz = iwork( i+2 ) - iwork( i )
                    msd2 = matsiz / 2_${ik}$
                    curprb = curprb + 1_${ik}$
                 end if
           ! merge lower order eigensystems (of size msd2 and matsiz - msd2)
           ! into an eigensystem of size matsiz.
           ! stdlib${ii}$_slaed1 is used only for the full eigensystem of a tridiagonal
           ! matrix.
           ! stdlib${ii}$_slaed7 handles the cases in which eigenvalues only or eigenvalues
           ! and eigenvectors of a full symmetric matrix (which was reduced to
           ! tridiagonal form) are desired.
                 if( icompq==2_${ik}$ ) then
                    call stdlib${ii}$_slaed1( matsiz, d( submat ), q( submat, submat ),ldq, iwork( &
                    indxq+submat ),e( submat+msd2-1 ), msd2, work,iwork( subpbs+1 ), info )
                              
                 else
                    call stdlib${ii}$_slaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,d( submat ), &
                    qstore( 1_${ik}$, submat ), ldqs,iwork( indxq+submat ), e( submat+msd2-1 ),msd2, &
                    work( iq ), iwork( iqptr ),iwork( iprmpt ), iwork( iperm ),iwork( igivpt ), &
                    iwork( igivcl ),work( igivnm ), work( iwrem ),iwork( subpbs+1 ), info )
                              
                 end if
                 if( info/=0 )go to 130
                 iwork( i / 2_${ik}$+1 ) = iwork( i+2 )
              end do loop_90
              subpbs = subpbs / 2_${ik}$
              curlvl = curlvl + 1_${ik}$
              go to 80
           end if
           ! end while
           ! re-merge the eigenvalues/vectors which were deflated at the final
           ! merge step.
           if( icompq==1_${ik}$ ) then
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
                 call stdlib${ii}$_scopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ )
              end do
              call stdlib${ii}$_scopy( n, work, 1_${ik}$, d, 1_${ik}$ )
           else if( icompq==2_${ik}$ ) then
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
                 call stdlib${ii}$_scopy( n, q( 1_${ik}$, j ), 1_${ik}$, work( n*i+1 ), 1_${ik}$ )
              end do
              call stdlib${ii}$_scopy( n, work, 1_${ik}$, d, 1_${ik}$ )
              call stdlib${ii}$_slacpy( 'A', n, n, work( n+1 ), n, q, ldq )
           else
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
              end do
              call stdlib${ii}$_scopy( n, work, 1_${ik}$, d, 1_${ik}$ )
           end if
           go to 140
           130 continue
           info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
           140 continue
           return
     end subroutine stdlib${ii}$_slaed0

     pure module subroutine stdlib${ii}$_dlaed0( icompq, qsiz, n, d, e, q, ldq, qstore, ldqs,work, iwork, info &
     !! DLAED0 computes all eigenvalues and corresponding eigenvectors of a
     !! symmetric tridiagonal matrix using the divide and conquer method.
               )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: icompq, ldq, ldqs, n, qsiz
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: d(*), e(*), q(ldq,*)
           real(dp), intent(out) :: qstore(ldqs,*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, &
           iq, iqptr, iwrem, j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, subpbs, &
                     tlvls
           real(dp) :: temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<0_${ik}$ .or. icompq>2_${ik}$ ) then
              info = -1_${ik}$
           else if( ( icompq==1_${ik}$ ) .and. ( qsiz<max( 0_${ik}$, n ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldqs<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DLAED0', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'DLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ )
           ! determine the size and placement of the submatrices, and save in
           ! the leading elements of iwork.
           iwork( 1_${ik}$ ) = n
           subpbs = 1_${ik}$
           tlvls = 0_${ik}$
           10 continue
           if( iwork( subpbs )>smlsiz ) then
              do j = subpbs, 1, -1
                 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$
                 iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$
              end do
              tlvls = tlvls + 1_${ik}$
              subpbs = 2_${ik}$*subpbs
              go to 10
           end if
           do j = 2, subpbs
              iwork( j ) = iwork( j ) + iwork( j-1 )
           end do
           ! divide the matrix into subpbs submatrices of size at most smlsiz+1
           ! using rank-1 modifications (cuts).
           spm1 = subpbs - 1_${ik}$
           do i = 1, spm1
              submat = iwork( i ) + 1_${ik}$
              smm1 = submat - 1_${ik}$
              d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
              d( submat ) = d( submat ) - abs( e( smm1 ) )
           end do
           indxq = 4_${ik}$*n + 3_${ik}$
           if( icompq/=2_${ik}$ ) then
              ! set up workspaces for eigenvalues only/accumulate new vectors
              ! routine
              temp = log( real( n,KIND=dp) ) / log( two )
              lgn = int( temp,KIND=${ik}$)
              if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
              if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
              iprmpt = indxq + n + 1_${ik}$
              iperm = iprmpt + n*lgn
              iqptr = iperm + n*lgn
              igivpt = iqptr + n + 2_${ik}$
              igivcl = igivpt + n*lgn
              igivnm = 1_${ik}$
              iq = igivnm + 2_${ik}$*n*lgn
              iwrem = iq + n**2_${ik}$ + 1_${ik}$
              ! initialize pointers
              do i = 0, subpbs
                 iwork( iprmpt+i ) = 1_${ik}$
                 iwork( igivpt+i ) = 1_${ik}$
              end do
              iwork( iqptr ) = 1_${ik}$
           end if
           ! solve each submatrix eigenproblem at the bottom of the divide and
           ! conquer tree.
           curr = 0_${ik}$
           loop_70: do i = 0, spm1
              if( i==0_${ik}$ ) then
                 submat = 1_${ik}$
                 matsiz = iwork( 1_${ik}$ )
              else
                 submat = iwork( i ) + 1_${ik}$
                 matsiz = iwork( i+1 ) - iwork( i )
              end if
              if( icompq==2_${ik}$ ) then
                 call stdlib${ii}$_dsteqr( 'I', matsiz, d( submat ), e( submat ),q( submat, submat ), &
                           ldq, work, info )
                 if( info/=0 )go to 130
              else
                 call stdlib${ii}$_dsteqr( 'I', matsiz, d( submat ), e( submat ),work( iq-1+iwork( &
                           iqptr+curr ) ), matsiz, work,info )
                 if( info/=0 )go to 130
                 if( icompq==1_${ik}$ ) then
                    call stdlib${ii}$_dgemm( 'N', 'N', qsiz, matsiz, matsiz, one,q( 1_${ik}$, submat ), ldq, &
                    work( iq-1+iwork( iqptr+curr ) ), matsiz, zero, qstore( 1_${ik}$, submat ),ldqs )
                              
                 end if
                 iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$
                 curr = curr + 1_${ik}$
              end if
              k = 1_${ik}$
              do j = submat, iwork( i+1 )
                 iwork( indxq+j ) = k
                 k = k + 1_${ik}$
              end do
           end do loop_70
           ! successively merge eigensystems of adjacent submatrices
           ! into eigensystem for the corresponding larger matrix.
           ! while ( subpbs > 1 )
           curlvl = 1_${ik}$
           80 continue
           if( subpbs>1_${ik}$ ) then
              spm2 = subpbs - 2_${ik}$
              loop_90: do i = 0, spm2, 2
                 if( i==0_${ik}$ ) then
                    submat = 1_${ik}$
                    matsiz = iwork( 2_${ik}$ )
                    msd2 = iwork( 1_${ik}$ )
                    curprb = 0_${ik}$
                 else
                    submat = iwork( i ) + 1_${ik}$
                    matsiz = iwork( i+2 ) - iwork( i )
                    msd2 = matsiz / 2_${ik}$
                    curprb = curprb + 1_${ik}$
                 end if
           ! merge lower order eigensystems (of size msd2 and matsiz - msd2)
           ! into an eigensystem of size matsiz.
           ! stdlib${ii}$_dlaed1 is used only for the full eigensystem of a tridiagonal
           ! matrix.
           ! stdlib${ii}$_dlaed7 handles the cases in which eigenvalues only or eigenvalues
           ! and eigenvectors of a full symmetric matrix (which was reduced to
           ! tridiagonal form) are desired.
                 if( icompq==2_${ik}$ ) then
                    call stdlib${ii}$_dlaed1( matsiz, d( submat ), q( submat, submat ),ldq, iwork( &
                    indxq+submat ),e( submat+msd2-1 ), msd2, work,iwork( subpbs+1 ), info )
                              
                 else
                    call stdlib${ii}$_dlaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,d( submat ), &
                    qstore( 1_${ik}$, submat ), ldqs,iwork( indxq+submat ), e( submat+msd2-1 ),msd2, &
                    work( iq ), iwork( iqptr ),iwork( iprmpt ), iwork( iperm ),iwork( igivpt ), &
                    iwork( igivcl ),work( igivnm ), work( iwrem ),iwork( subpbs+1 ), info )
                              
                 end if
                 if( info/=0 )go to 130
                 iwork( i / 2_${ik}$+1 ) = iwork( i+2 )
              end do loop_90
              subpbs = subpbs / 2_${ik}$
              curlvl = curlvl + 1_${ik}$
              go to 80
           end if
           ! end while
           ! re-merge the eigenvalues/vectors which were deflated at the final
           ! merge step.
           if( icompq==1_${ik}$ ) then
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
                 call stdlib${ii}$_dcopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ )
              end do
              call stdlib${ii}$_dcopy( n, work, 1_${ik}$, d, 1_${ik}$ )
           else if( icompq==2_${ik}$ ) then
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
                 call stdlib${ii}$_dcopy( n, q( 1_${ik}$, j ), 1_${ik}$, work( n*i+1 ), 1_${ik}$ )
              end do
              call stdlib${ii}$_dcopy( n, work, 1_${ik}$, d, 1_${ik}$ )
              call stdlib${ii}$_dlacpy( 'A', n, n, work( n+1 ), n, q, ldq )
           else
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
              end do
              call stdlib${ii}$_dcopy( n, work, 1_${ik}$, d, 1_${ik}$ )
           end if
           go to 140
           130 continue
           info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
           140 continue
           return
     end subroutine stdlib${ii}$_dlaed0

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$laed0( icompq, qsiz, n, d, e, q, ldq, qstore, ldqs,work, iwork, info &
     !! DLAED0: computes all eigenvalues and corresponding eigenvectors of a
     !! symmetric tridiagonal matrix using the divide and conquer method.
               )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: icompq, ldq, ldqs, n, qsiz
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(inout) :: d(*), e(*), q(ldq,*)
           real(${rk}$), intent(out) :: qstore(ldqs,*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, &
           iq, iqptr, iwrem, j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, subpbs, &
                     tlvls
           real(${rk}$) :: temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<0_${ik}$ .or. icompq>2_${ik}$ ) then
              info = -1_${ik}$
           else if( ( icompq==1_${ik}$ ) .and. ( qsiz<max( 0_${ik}$, n ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldqs<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DLAED0', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'DLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ )
           ! determine the size and placement of the submatrices, and save in
           ! the leading elements of iwork.
           iwork( 1_${ik}$ ) = n
           subpbs = 1_${ik}$
           tlvls = 0_${ik}$
           10 continue
           if( iwork( subpbs )>smlsiz ) then
              do j = subpbs, 1, -1
                 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$
                 iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$
              end do
              tlvls = tlvls + 1_${ik}$
              subpbs = 2_${ik}$*subpbs
              go to 10
           end if
           do j = 2, subpbs
              iwork( j ) = iwork( j ) + iwork( j-1 )
           end do
           ! divide the matrix into subpbs submatrices of size at most smlsiz+1
           ! using rank-1 modifications (cuts).
           spm1 = subpbs - 1_${ik}$
           do i = 1, spm1
              submat = iwork( i ) + 1_${ik}$
              smm1 = submat - 1_${ik}$
              d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
              d( submat ) = d( submat ) - abs( e( smm1 ) )
           end do
           indxq = 4_${ik}$*n + 3_${ik}$
           if( icompq/=2_${ik}$ ) then
              ! set up workspaces for eigenvalues only/accumulate new vectors
              ! routine
              temp = log( real( n,KIND=${rk}$) ) / log( two )
              lgn = int( temp,KIND=${ik}$)
              if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
              if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
              iprmpt = indxq + n + 1_${ik}$
              iperm = iprmpt + n*lgn
              iqptr = iperm + n*lgn
              igivpt = iqptr + n + 2_${ik}$
              igivcl = igivpt + n*lgn
              igivnm = 1_${ik}$
              iq = igivnm + 2_${ik}$*n*lgn
              iwrem = iq + n**2_${ik}$ + 1_${ik}$
              ! initialize pointers
              do i = 0, subpbs
                 iwork( iprmpt+i ) = 1_${ik}$
                 iwork( igivpt+i ) = 1_${ik}$
              end do
              iwork( iqptr ) = 1_${ik}$
           end if
           ! solve each submatrix eigenproblem at the bottom of the divide and
           ! conquer tree.
           curr = 0_${ik}$
           loop_70: do i = 0, spm1
              if( i==0_${ik}$ ) then
                 submat = 1_${ik}$
                 matsiz = iwork( 1_${ik}$ )
              else
                 submat = iwork( i ) + 1_${ik}$
                 matsiz = iwork( i+1 ) - iwork( i )
              end if
              if( icompq==2_${ik}$ ) then
                 call stdlib${ii}$_${ri}$steqr( 'I', matsiz, d( submat ), e( submat ),q( submat, submat ), &
                           ldq, work, info )
                 if( info/=0 )go to 130
              else
                 call stdlib${ii}$_${ri}$steqr( 'I', matsiz, d( submat ), e( submat ),work( iq-1+iwork( &
                           iqptr+curr ) ), matsiz, work,info )
                 if( info/=0 )go to 130
                 if( icompq==1_${ik}$ ) then
                    call stdlib${ii}$_${ri}$gemm( 'N', 'N', qsiz, matsiz, matsiz, one,q( 1_${ik}$, submat ), ldq, &
                    work( iq-1+iwork( iqptr+curr ) ), matsiz, zero, qstore( 1_${ik}$, submat ),ldqs )
                              
                 end if
                 iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$
                 curr = curr + 1_${ik}$
              end if
              k = 1_${ik}$
              do j = submat, iwork( i+1 )
                 iwork( indxq+j ) = k
                 k = k + 1_${ik}$
              end do
           end do loop_70
           ! successively merge eigensystems of adjacent submatrices
           ! into eigensystem for the corresponding larger matrix.
           ! while ( subpbs > 1 )
           curlvl = 1_${ik}$
           80 continue
           if( subpbs>1_${ik}$ ) then
              spm2 = subpbs - 2_${ik}$
              loop_90: do i = 0, spm2, 2
                 if( i==0_${ik}$ ) then
                    submat = 1_${ik}$
                    matsiz = iwork( 2_${ik}$ )
                    msd2 = iwork( 1_${ik}$ )
                    curprb = 0_${ik}$
                 else
                    submat = iwork( i ) + 1_${ik}$
                    matsiz = iwork( i+2 ) - iwork( i )
                    msd2 = matsiz / 2_${ik}$
                    curprb = curprb + 1_${ik}$
                 end if
           ! merge lower order eigensystems (of size msd2 and matsiz - msd2)
           ! into an eigensystem of size matsiz.
           ! stdlib${ii}$_${ri}$laed1 is used only for the full eigensystem of a tridiagonal
           ! matrix.
           ! stdlib${ii}$_${ri}$laed7 handles the cases in which eigenvalues only or eigenvalues
           ! and eigenvectors of a full symmetric matrix (which was reduced to
           ! tridiagonal form) are desired.
                 if( icompq==2_${ik}$ ) then
                    call stdlib${ii}$_${ri}$laed1( matsiz, d( submat ), q( submat, submat ),ldq, iwork( &
                    indxq+submat ),e( submat+msd2-1 ), msd2, work,iwork( subpbs+1 ), info )
                              
                 else
                    call stdlib${ii}$_${ri}$laed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,d( submat ), &
                    qstore( 1_${ik}$, submat ), ldqs,iwork( indxq+submat ), e( submat+msd2-1 ),msd2, &
                    work( iq ), iwork( iqptr ),iwork( iprmpt ), iwork( iperm ),iwork( igivpt ), &
                    iwork( igivcl ),work( igivnm ), work( iwrem ),iwork( subpbs+1 ), info )
                              
                 end if
                 if( info/=0 )go to 130
                 iwork( i / 2_${ik}$+1 ) = iwork( i+2 )
              end do loop_90
              subpbs = subpbs / 2_${ik}$
              curlvl = curlvl + 1_${ik}$
              go to 80
           end if
           ! end while
           ! re-merge the eigenvalues/vectors which were deflated at the final
           ! merge step.
           if( icompq==1_${ik}$ ) then
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
                 call stdlib${ii}$_${ri}$copy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ )
              end do
              call stdlib${ii}$_${ri}$copy( n, work, 1_${ik}$, d, 1_${ik}$ )
           else if( icompq==2_${ik}$ ) then
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
                 call stdlib${ii}$_${ri}$copy( n, q( 1_${ik}$, j ), 1_${ik}$, work( n*i+1 ), 1_${ik}$ )
              end do
              call stdlib${ii}$_${ri}$copy( n, work, 1_${ik}$, d, 1_${ik}$ )
              call stdlib${ii}$_${ri}$lacpy( 'A', n, n, work( n+1 ), n, q, ldq )
           else
              do i = 1, n
                 j = iwork( indxq+i )
                 work( i ) = d( j )
              end do
              call stdlib${ii}$_${ri}$copy( n, work, 1_${ik}$, d, 1_${ik}$ )
           end if
           go to 140
           130 continue
           info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
           140 continue
           return
     end subroutine stdlib${ii}$_${ri}$laed0

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_claed0( qsiz, n, d, e, q, ldq, qstore, ldqs, rwork,iwork, info )
     !! Using the divide and conquer method, CLAED0: computes all eigenvalues
     !! of a symmetric tridiagonal matrix which is one diagonal block of
     !! those from reducing a dense or band Hermitian matrix and
     !! corresponding eigenvectors of the dense or band matrix.
               
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldq, ldqs, n, qsiz
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: d(*), e(*)
           real(sp), intent(out) :: rwork(*)
           complex(sp), intent(inout) :: q(ldq,*)
           complex(sp), intent(out) :: qstore(ldqs,*)
        ! =====================================================================
        ! warning:      n could be as big as qsiz!
           
           ! Local Scalars 
           integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, &
           iq, iqptr, iwrem, j, k, lgn, ll, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, &
                     subpbs, tlvls
           real(sp) :: temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           ! if( icompq < 0 .or. icompq > 2 ) then
              ! info = -1
           ! else if( ( icompq == 1 ) .and. ( qsiz < max( 0, n ) ) )
          ! $        then
           if( qsiz<max( 0_${ik}$, n ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldqs<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CLAED0', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'CLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ )
           ! determine the size and placement of the submatrices, and save in
           ! the leading elements of iwork.
           iwork( 1_${ik}$ ) = n
           subpbs = 1_${ik}$
           tlvls = 0_${ik}$
           10 continue
           if( iwork( subpbs )>smlsiz ) then
              do j = subpbs, 1, -1
                 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$
                 iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$
              end do
              tlvls = tlvls + 1_${ik}$
              subpbs = 2_${ik}$*subpbs
              go to 10
           end if
           do j = 2, subpbs
              iwork( j ) = iwork( j ) + iwork( j-1 )
           end do
           ! divide the matrix into subpbs submatrices of size at most smlsiz+1
           ! using rank-1 modifications (cuts).
           spm1 = subpbs - 1_${ik}$
           do i = 1, spm1
              submat = iwork( i ) + 1_${ik}$
              smm1 = submat - 1_${ik}$
              d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
              d( submat ) = d( submat ) - abs( e( smm1 ) )
           end do
           indxq = 4_${ik}$*n + 3_${ik}$
           ! set up workspaces for eigenvalues only/accumulate new vectors
           ! routine
           temp = log( real( n,KIND=sp) ) / log( two )
           lgn = int( temp,KIND=${ik}$)
           if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
           if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
           iprmpt = indxq + n + 1_${ik}$
           iperm = iprmpt + n*lgn
           iqptr = iperm + n*lgn
           igivpt = iqptr + n + 2_${ik}$
           igivcl = igivpt + n*lgn
           igivnm = 1_${ik}$
           iq = igivnm + 2_${ik}$*n*lgn
           iwrem = iq + n**2_${ik}$ + 1_${ik}$
           ! initialize pointers
           do i = 0, subpbs
              iwork( iprmpt+i ) = 1_${ik}$
              iwork( igivpt+i ) = 1_${ik}$
           end do
           iwork( iqptr ) = 1_${ik}$
           ! solve each submatrix eigenproblem at the bottom of the divide and
           ! conquer tree.
           curr = 0_${ik}$
           do i = 0, spm1
              if( i==0_${ik}$ ) then
                 submat = 1_${ik}$
                 matsiz = iwork( 1_${ik}$ )
              else
                 submat = iwork( i ) + 1_${ik}$
                 matsiz = iwork( i+1 ) - iwork( i )
              end if
              ll = iq - 1_${ik}$ + iwork( iqptr+curr )
              call stdlib${ii}$_ssteqr( 'I', matsiz, d( submat ), e( submat ),rwork( ll ), matsiz, &
                        rwork, info )
              call stdlib${ii}$_clacrm( qsiz, matsiz, q( 1_${ik}$, submat ), ldq, rwork( ll ),matsiz, qstore( &
                        1_${ik}$, submat ), ldqs,rwork( iwrem ) )
              iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$
              curr = curr + 1_${ik}$
              if( info>0_${ik}$ ) then
                 info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
                 return
              end if
              k = 1_${ik}$
              do j = submat, iwork( i+1 )
                 iwork( indxq+j ) = k
                 k = k + 1_${ik}$
              end do
           end do
           ! successively merge eigensystems of adjacent submatrices
           ! into eigensystem for the corresponding larger matrix.
           ! while ( subpbs > 1 )
           curlvl = 1_${ik}$
           80 continue
           if( subpbs>1_${ik}$ ) then
              spm2 = subpbs - 2_${ik}$
              do i = 0, spm2, 2
                 if( i==0_${ik}$ ) then
                    submat = 1_${ik}$
                    matsiz = iwork( 2_${ik}$ )
                    msd2 = iwork( 1_${ik}$ )
                    curprb = 0_${ik}$
                 else
                    submat = iwork( i ) + 1_${ik}$
                    matsiz = iwork( i+2 ) - iwork( i )
                    msd2 = matsiz / 2_${ik}$
                    curprb = curprb + 1_${ik}$
                 end if
           ! merge lower order eigensystems (of size msd2 and matsiz - msd2)
           ! into an eigensystem of size matsiz.  stdlib${ii}$_claed7 handles the case
           ! when the eigenvectors of a full or band hermitian matrix (which
           ! was reduced to tridiagonal form) are desired.
           ! i am free to use q as a valuable working space until loop 150.
                 call stdlib${ii}$_claed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,d( submat ), &
                 qstore( 1_${ik}$, submat ), ldqs,e( submat+msd2-1 ), iwork( indxq+submat ),rwork( iq ), &
                 iwork( iqptr ), iwork( iprmpt ),iwork( iperm ), iwork( igivpt ),iwork( igivcl ), &
                           rwork( igivnm ),q( 1_${ik}$, submat ), rwork( iwrem ),iwork( subpbs+1 ), info )
                 if( info>0_${ik}$ ) then
                    info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
                    return
                 end if
                 iwork( i / 2_${ik}$+1 ) = iwork( i+2 )
              end do
              subpbs = subpbs / 2_${ik}$
              curlvl = curlvl + 1_${ik}$
              go to 80
           end if
           ! end while
           ! re-merge the eigenvalues/vectors which were deflated at the final
           ! merge step.
           do i = 1, n
              j = iwork( indxq+i )
              rwork( i ) = d( j )
              call stdlib${ii}$_ccopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ )
           end do
           call stdlib${ii}$_scopy( n, rwork, 1_${ik}$, d, 1_${ik}$ )
           return
     end subroutine stdlib${ii}$_claed0

     pure module subroutine stdlib${ii}$_zlaed0( qsiz, n, d, e, q, ldq, qstore, ldqs, rwork,iwork, info )
     !! Using the divide and conquer method, ZLAED0: computes all eigenvalues
     !! of a symmetric tridiagonal matrix which is one diagonal block of
     !! those from reducing a dense or band Hermitian matrix and
     !! corresponding eigenvectors of the dense or band matrix.
               
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldq, ldqs, n, qsiz
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: d(*), e(*)
           real(dp), intent(out) :: rwork(*)
           complex(dp), intent(inout) :: q(ldq,*)
           complex(dp), intent(out) :: qstore(ldqs,*)
        ! =====================================================================
        ! warning:      n could be as big as qsiz!
           
           ! Local Scalars 
           integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, &
           iq, iqptr, iwrem, j, k, lgn, ll, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, &
                     subpbs, tlvls
           real(dp) :: temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           ! if( icompq < 0 .or. icompq > 2 ) then
              ! info = -1
           ! else if( ( icompq == 1 ) .and. ( qsiz < max( 0, n ) ) )
          ! $        then
           if( qsiz<max( 0_${ik}$, n ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldqs<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZLAED0', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'ZLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ )
           ! determine the size and placement of the submatrices, and save in
           ! the leading elements of iwork.
           iwork( 1_${ik}$ ) = n
           subpbs = 1_${ik}$
           tlvls = 0_${ik}$
           10 continue
           if( iwork( subpbs )>smlsiz ) then
              do j = subpbs, 1, -1
                 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$
                 iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$
              end do
              tlvls = tlvls + 1_${ik}$
              subpbs = 2_${ik}$*subpbs
              go to 10
           end if
           do j = 2, subpbs
              iwork( j ) = iwork( j ) + iwork( j-1 )
           end do
           ! divide the matrix into subpbs submatrices of size at most smlsiz+1
           ! using rank-1 modifications (cuts).
           spm1 = subpbs - 1_${ik}$
           do i = 1, spm1
              submat = iwork( i ) + 1_${ik}$
              smm1 = submat - 1_${ik}$
              d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
              d( submat ) = d( submat ) - abs( e( smm1 ) )
           end do
           indxq = 4_${ik}$*n + 3_${ik}$
           ! set up workspaces for eigenvalues only/accumulate new vectors
           ! routine
           temp = log( real( n,KIND=dp) ) / log( two )
           lgn = int( temp,KIND=${ik}$)
           if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
           if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
           iprmpt = indxq + n + 1_${ik}$
           iperm = iprmpt + n*lgn
           iqptr = iperm + n*lgn
           igivpt = iqptr + n + 2_${ik}$
           igivcl = igivpt + n*lgn
           igivnm = 1_${ik}$
           iq = igivnm + 2_${ik}$*n*lgn
           iwrem = iq + n**2_${ik}$ + 1_${ik}$
           ! initialize pointers
           do i = 0, subpbs
              iwork( iprmpt+i ) = 1_${ik}$
              iwork( igivpt+i ) = 1_${ik}$
           end do
           iwork( iqptr ) = 1_${ik}$
           ! solve each submatrix eigenproblem at the bottom of the divide and
           ! conquer tree.
           curr = 0_${ik}$
           do i = 0, spm1
              if( i==0_${ik}$ ) then
                 submat = 1_${ik}$
                 matsiz = iwork( 1_${ik}$ )
              else
                 submat = iwork( i ) + 1_${ik}$
                 matsiz = iwork( i+1 ) - iwork( i )
              end if
              ll = iq - 1_${ik}$ + iwork( iqptr+curr )
              call stdlib${ii}$_dsteqr( 'I', matsiz, d( submat ), e( submat ),rwork( ll ), matsiz, &
                        rwork, info )
              call stdlib${ii}$_zlacrm( qsiz, matsiz, q( 1_${ik}$, submat ), ldq, rwork( ll ),matsiz, qstore( &
                        1_${ik}$, submat ), ldqs,rwork( iwrem ) )
              iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$
              curr = curr + 1_${ik}$
              if( info>0_${ik}$ ) then
                 info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
                 return
              end if
              k = 1_${ik}$
              do j = submat, iwork( i+1 )
                 iwork( indxq+j ) = k
                 k = k + 1_${ik}$
              end do
           end do
           ! successively merge eigensystems of adjacent submatrices
           ! into eigensystem for the corresponding larger matrix.
           ! while ( subpbs > 1 )
           curlvl = 1_${ik}$
           80 continue
           if( subpbs>1_${ik}$ ) then
              spm2 = subpbs - 2_${ik}$
              do i = 0, spm2, 2
                 if( i==0_${ik}$ ) then
                    submat = 1_${ik}$
                    matsiz = iwork( 2_${ik}$ )
                    msd2 = iwork( 1_${ik}$ )
                    curprb = 0_${ik}$
                 else
                    submat = iwork( i ) + 1_${ik}$
                    matsiz = iwork( i+2 ) - iwork( i )
                    msd2 = matsiz / 2_${ik}$
                    curprb = curprb + 1_${ik}$
                 end if
           ! merge lower order eigensystems (of size msd2 and matsiz - msd2)
           ! into an eigensystem of size matsiz.  stdlib${ii}$_zlaed7 handles the case
           ! when the eigenvectors of a full or band hermitian matrix (which
           ! was reduced to tridiagonal form) are desired.
           ! i am free to use q as a valuable working space until loop 150.
                 call stdlib${ii}$_zlaed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,d( submat ), &
                 qstore( 1_${ik}$, submat ), ldqs,e( submat+msd2-1 ), iwork( indxq+submat ),rwork( iq ), &
                 iwork( iqptr ), iwork( iprmpt ),iwork( iperm ), iwork( igivpt ),iwork( igivcl ), &
                           rwork( igivnm ),q( 1_${ik}$, submat ), rwork( iwrem ),iwork( subpbs+1 ), info )
                 if( info>0_${ik}$ ) then
                    info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
                    return
                 end if
                 iwork( i / 2_${ik}$+1 ) = iwork( i+2 )
              end do
              subpbs = subpbs / 2_${ik}$
              curlvl = curlvl + 1_${ik}$
              go to 80
           end if
           ! end while
           ! re-merge the eigenvalues/vectors which were deflated at the final
           ! merge step.
           do i = 1, n
              j = iwork( indxq+i )
              rwork( i ) = d( j )
              call stdlib${ii}$_zcopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ )
           end do
           call stdlib${ii}$_dcopy( n, rwork, 1_${ik}$, d, 1_${ik}$ )
           return
     end subroutine stdlib${ii}$_zlaed0

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$laed0( qsiz, n, d, e, q, ldq, qstore, ldqs, rwork,iwork, info )
     !! Using the divide and conquer method, ZLAED0: computes all eigenvalues
     !! of a symmetric tridiagonal matrix which is one diagonal block of
     !! those from reducing a dense or band Hermitian matrix and
     !! corresponding eigenvectors of the dense or band matrix.
               
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldq, ldqs, n, qsiz
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(${ck}$), intent(inout) :: d(*), e(*)
           real(${ck}$), intent(out) :: rwork(*)
           complex(${ck}$), intent(inout) :: q(ldq,*)
           complex(${ck}$), intent(out) :: qstore(ldqs,*)
        ! =====================================================================
        ! warning:      n could be as big as qsiz!
           
           ! Local Scalars 
           integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, &
           iq, iqptr, iwrem, j, k, lgn, ll, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, &
                     subpbs, tlvls
           real(${ck}$) :: temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           ! if( icompq < 0 .or. icompq > 2 ) then
              ! info = -1
           ! else if( ( icompq == 1 ) .and. ( qsiz < max( 0, n ) ) )
          ! $        then
           if( qsiz<max( 0_${ik}$, n ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldqs<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZLAED0', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'ZLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ )
           ! determine the size and placement of the submatrices, and save in
           ! the leading elements of iwork.
           iwork( 1_${ik}$ ) = n
           subpbs = 1_${ik}$
           tlvls = 0_${ik}$
           10 continue
           if( iwork( subpbs )>smlsiz ) then
              do j = subpbs, 1, -1
                 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$
                 iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$
              end do
              tlvls = tlvls + 1_${ik}$
              subpbs = 2_${ik}$*subpbs
              go to 10
           end if
           do j = 2, subpbs
              iwork( j ) = iwork( j ) + iwork( j-1 )
           end do
           ! divide the matrix into subpbs submatrices of size at most smlsiz+1
           ! using rank-1 modifications (cuts).
           spm1 = subpbs - 1_${ik}$
           do i = 1, spm1
              submat = iwork( i ) + 1_${ik}$
              smm1 = submat - 1_${ik}$
              d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
              d( submat ) = d( submat ) - abs( e( smm1 ) )
           end do
           indxq = 4_${ik}$*n + 3_${ik}$
           ! set up workspaces for eigenvalues only/accumulate new vectors
           ! routine
           temp = log( real( n,KIND=${ck}$) ) / log( two )
           lgn = int( temp,KIND=${ik}$)
           if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
           if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$
           iprmpt = indxq + n + 1_${ik}$
           iperm = iprmpt + n*lgn
           iqptr = iperm + n*lgn
           igivpt = iqptr + n + 2_${ik}$
           igivcl = igivpt + n*lgn
           igivnm = 1_${ik}$
           iq = igivnm + 2_${ik}$*n*lgn
           iwrem = iq + n**2_${ik}$ + 1_${ik}$
           ! initialize pointers
           do i = 0, subpbs
              iwork( iprmpt+i ) = 1_${ik}$
              iwork( igivpt+i ) = 1_${ik}$
           end do
           iwork( iqptr ) = 1_${ik}$
           ! solve each submatrix eigenproblem at the bottom of the divide and
           ! conquer tree.
           curr = 0_${ik}$
           do i = 0, spm1
              if( i==0_${ik}$ ) then
                 submat = 1_${ik}$
                 matsiz = iwork( 1_${ik}$ )
              else
                 submat = iwork( i ) + 1_${ik}$
                 matsiz = iwork( i+1 ) - iwork( i )
              end if
              ll = iq - 1_${ik}$ + iwork( iqptr+curr )
              call stdlib${ii}$_${c2ri(ci)}$steqr( 'I', matsiz, d( submat ), e( submat ),rwork( ll ), matsiz, &
                        rwork, info )
              call stdlib${ii}$_${ci}$lacrm( qsiz, matsiz, q( 1_${ik}$, submat ), ldq, rwork( ll ),matsiz, qstore( &
                        1_${ik}$, submat ), ldqs,rwork( iwrem ) )
              iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$
              curr = curr + 1_${ik}$
              if( info>0_${ik}$ ) then
                 info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
                 return
              end if
              k = 1_${ik}$
              do j = submat, iwork( i+1 )
                 iwork( indxq+j ) = k
                 k = k + 1_${ik}$
              end do
           end do
           ! successively merge eigensystems of adjacent submatrices
           ! into eigensystem for the corresponding larger matrix.
           ! while ( subpbs > 1 )
           curlvl = 1_${ik}$
           80 continue
           if( subpbs>1_${ik}$ ) then
              spm2 = subpbs - 2_${ik}$
              do i = 0, spm2, 2
                 if( i==0_${ik}$ ) then
                    submat = 1_${ik}$
                    matsiz = iwork( 2_${ik}$ )
                    msd2 = iwork( 1_${ik}$ )
                    curprb = 0_${ik}$
                 else
                    submat = iwork( i ) + 1_${ik}$
                    matsiz = iwork( i+2 ) - iwork( i )
                    msd2 = matsiz / 2_${ik}$
                    curprb = curprb + 1_${ik}$
                 end if
           ! merge lower order eigensystems (of size msd2 and matsiz - msd2)
           ! into an eigensystem of size matsiz.  stdlib${ii}$_${ci}$laed7 handles the case
           ! when the eigenvectors of a full or band hermitian matrix (which
           ! was reduced to tridiagonal form) are desired.
           ! i am free to use q as a valuable working space until loop 150.
                 call stdlib${ii}$_${ci}$laed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,d( submat ), &
                 qstore( 1_${ik}$, submat ), ldqs,e( submat+msd2-1 ), iwork( indxq+submat ),rwork( iq ), &
                 iwork( iqptr ), iwork( iprmpt ),iwork( iperm ), iwork( igivpt ),iwork( igivcl ), &
                           rwork( igivnm ),q( 1_${ik}$, submat ), rwork( iwrem ),iwork( subpbs+1 ), info )
                 if( info>0_${ik}$ ) then
                    info = submat*( n+1 ) + submat + matsiz - 1_${ik}$
                    return
                 end if
                 iwork( i / 2_${ik}$+1 ) = iwork( i+2 )
              end do
              subpbs = subpbs / 2_${ik}$
              curlvl = curlvl + 1_${ik}$
              go to 80
           end if
           ! end while
           ! re-merge the eigenvalues/vectors which were deflated at the final
           ! merge step.
           do i = 1, n
              j = iwork( indxq+i )
              rwork( i ) = d( j )
              call stdlib${ii}$_${ci}$copy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ )
           end do
           call stdlib${ii}$_${c2ri(ci)}$copy( n, rwork, 1_${ik}$, d, 1_${ik}$ )
           return
     end subroutine stdlib${ii}$_${ci}$laed0

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slaed1( n, d, q, ldq, indxq, rho, cutpnt, work, iwork,info )
     !! SLAED1 computes the updated eigensystem of a diagonal
     !! matrix after modification by a rank-one symmetric matrix.  This
     !! routine is used only for the eigenproblem which requires all
     !! eigenvalues and eigenvectors of a tridiagonal matrix.  SLAED7 handles
     !! the case in which eigenvalues only or eigenvalues and eigenvectors
     !! of a full symmetric matrix (which was reduced to tridiagonal form)
     !! are desired.
     !! T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
     !! where Z = Q**T*u, u is a vector of length N with ones in the
     !! CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
     !! The eigenvectors of the original matrix are stored in Q, and the
     !! eigenvalues are in D.  The algorithm consists of three stages:
     !! The first stage consists of deflating the size of the problem
     !! when there are multiple eigenvalues or if there is a zero in
     !! the Z vector.  For each such occurrence the dimension of the
     !! secular equation problem is reduced by one.  This stage is
     !! performed by the routine SLAED2.
     !! The second stage consists of calculating the updated
     !! eigenvalues. This is done by finding the roots of the secular
     !! equation via the routine SLAED4 (as called by SLAED3).
     !! This routine also calculates the eigenvectors of the current
     !! problem.
     !! The final stage consists of computing the updated eigenvectors
     !! directly using the updated eigenvalues.  The eigenvectors for
     !! the current problem are multiplied with the eigenvectors from
     !! the overall problem.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: cutpnt, ldq, n
           integer(${ik}$), intent(out) :: info
           real(sp), intent(inout) :: rho
           ! Array Arguments 
           integer(${ik}$), intent(inout) :: indxq(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: d(*), q(ldq,*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           integer(${ik}$) :: coltyp, cpp1, i, idlmda, indx, indxc, indxp, iq2, is, iw, iz, k, n1, &
                     n2
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( n<0_${ik}$ ) then
              info = -1_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( min( 1_${ik}$, n / 2_${ik}$ )>cutpnt .or. ( n / 2_${ik}$ )<cutpnt ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SLAED1', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! the following values are integer pointers which indicate
           ! the portion of the workspace
           ! used by a particular array in stdlib${ii}$_slaed2 and stdlib${ii}$_slaed3.
           iz = 1_${ik}$
           idlmda = iz + n
           iw = idlmda + n
           iq2 = iw + n
           indx = 1_${ik}$
           indxc = indx + n
           coltyp = indxc + n
           indxp = coltyp + n
           ! form the z-vector which consists of the last row of q_1 and the
           ! first row of q_2.
           call stdlib${ii}$_scopy( cutpnt, q( cutpnt, 1_${ik}$ ), ldq, work( iz ), 1_${ik}$ )
           cpp1 = cutpnt + 1_${ik}$
           call stdlib${ii}$_scopy( n-cutpnt, q( cpp1, cpp1 ), ldq, work( iz+cutpnt ), 1_${ik}$ )
           ! deflate eigenvalues.
           call stdlib${ii}$_slaed2( k, n, cutpnt, d, q, ldq, indxq, rho, work( iz ),work( idlmda ), &
           work( iw ), work( iq2 ),iwork( indx ), iwork( indxc ), iwork( indxp ),iwork( coltyp ), &
                     info )
           if( info/=0 )go to 20
           ! solve secular equation.
           if( k/=0_${ik}$ ) then
              is = ( iwork( coltyp )+iwork( coltyp+1 ) )*cutpnt +( iwork( coltyp+1 )+iwork( &
                        coltyp+2 ) )*( n-cutpnt ) + iq2
              call stdlib${ii}$_slaed3( k, n, cutpnt, d, q, ldq, rho, work( idlmda ),work( iq2 ), iwork(&
                         indxc ), iwork( coltyp ),work( iw ), work( is ), info )
              if( info/=0 )go to 20
           ! prepare the indxq sorting permutation.
              n1 = k
              n2 = n - k
              call stdlib${ii}$_slamrg( n1, n2, d, 1_${ik}$, -1_${ik}$, indxq )
           else
              do i = 1, n
                 indxq( i ) = i
              end do
           end if
           20 continue
           return
     end subroutine stdlib${ii}$_slaed1

     pure module subroutine stdlib${ii}$_dlaed1( n, d, q, ldq, indxq, rho, cutpnt, work, iwork,info )
     !! DLAED1 computes the updated eigensystem of a diagonal
     !! matrix after modification by a rank-one symmetric matrix.  This
     !! routine is used only for the eigenproblem which requires all
     !! eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles
     !! the case in which eigenvalues only or eigenvalues and eigenvectors
     !! of a full symmetric matrix (which was reduced to tridiagonal form)
     !! are desired.
     !! T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
     !! where Z = Q**T*u, u is a vector of length N with ones in the
     !! CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
     !! The eigenvectors of the original matrix are stored in Q, and the
     !! eigenvalues are in D.  The algorithm consists of three stages:
     !! The first stage consists of deflating the size of the problem
     !! when there are multiple eigenvalues or if there is a zero in
     !! the Z vector.  For each such occurrence the dimension of the
     !! secular equation problem is reduced by one.  This stage is
     !! performed by the routine DLAED2.
     !! The second stage consists of calculating the updated
     !! eigenvalues. This is done by finding the roots of the secular
     !! equation via the routine DLAED4 (as called by DLAED3).
     !! This routine also calculates the eigenvectors of the current
     !! problem.
     !! The final stage consists of computing the updated eigenvectors
     !! directly using the updated eigenvalues.  The eigenvectors for
     !! the current problem are multiplied with the eigenvectors from
     !! the overall problem.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: cutpnt, ldq, n
           integer(${ik}$), intent(out) :: info
           real(dp), intent(inout) :: rho
           ! Array Arguments 
           integer(${ik}$), intent(inout) :: indxq(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: d(*), q(ldq,*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           integer(${ik}$) :: coltyp, i, idlmda, indx, indxc, indxp, iq2, is, iw, iz, k, n1, n2, &
                     zpp1
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( n<0_${ik}$ ) then
              info = -1_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( min( 1_${ik}$, n / 2_${ik}$ )>cutpnt .or. ( n / 2_${ik}$ )<cutpnt ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DLAED1', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! the following values are integer pointers which indicate
           ! the portion of the workspace
           ! used by a particular array in stdlib${ii}$_dlaed2 and stdlib${ii}$_dlaed3.
           iz = 1_${ik}$
           idlmda = iz + n
           iw = idlmda + n
           iq2 = iw + n
           indx = 1_${ik}$
           indxc = indx + n
           coltyp = indxc + n
           indxp = coltyp + n
           ! form the z-vector which consists of the last row of q_1 and the
           ! first row of q_2.
           call stdlib${ii}$_dcopy( cutpnt, q( cutpnt, 1_${ik}$ ), ldq, work( iz ), 1_${ik}$ )
           zpp1 = cutpnt + 1_${ik}$
           call stdlib${ii}$_dcopy( n-cutpnt, q( zpp1, zpp1 ), ldq, work( iz+cutpnt ), 1_${ik}$ )
           ! deflate eigenvalues.
           call stdlib${ii}$_dlaed2( k, n, cutpnt, d, q, ldq, indxq, rho, work( iz ),work( idlmda ), &
           work( iw ), work( iq2 ),iwork( indx ), iwork( indxc ), iwork( indxp ),iwork( coltyp ), &
                     info )
           if( info/=0 )go to 20
           ! solve secular equation.
           if( k/=0_${ik}$ ) then
              is = ( iwork( coltyp )+iwork( coltyp+1 ) )*cutpnt +( iwork( coltyp+1 )+iwork( &
                        coltyp+2 ) )*( n-cutpnt ) + iq2
              call stdlib${ii}$_dlaed3( k, n, cutpnt, d, q, ldq, rho, work( idlmda ),work( iq2 ), iwork(&
                         indxc ), iwork( coltyp ),work( iw ), work( is ), info )
              if( info/=0 )go to 20
           ! prepare the indxq sorting permutation.
              n1 = k
              n2 = n - k
              call stdlib${ii}$_dlamrg( n1, n2, d, 1_${ik}$, -1_${ik}$, indxq )
           else
              do i = 1, n
                 indxq( i ) = i
              end do
           end if
           20 continue
           return
     end subroutine stdlib${ii}$_dlaed1

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$laed1( n, d, q, ldq, indxq, rho, cutpnt, work, iwork,info )
     !! DLAED1: computes the updated eigensystem of a diagonal
     !! matrix after modification by a rank-one symmetric matrix.  This
     !! routine is used only for the eigenproblem which requires all
     !! eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles
     !! the case in which eigenvalues only or eigenvalues and eigenvectors
     !! of a full symmetric matrix (which was reduced to tridiagonal form)
     !! are desired.
     !! T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
     !! where Z = Q**T*u, u is a vector of length N with ones in the
     !! CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
     !! The eigenvectors of the original matrix are stored in Q, and the
     !! eigenvalues are in D.  The algorithm consists of three stages:
     !! The first stage consists of deflating the size of the problem
     !! when there are multiple eigenvalues or if there is a zero in
     !! the Z vector.  For each such occurrence the dimension of the
     !! secular equation problem is reduced by one.  This stage is
     !! performed by the routine DLAED2.
     !! The second stage consists of calculating the updated
     !! eigenvalues. This is done by finding the roots of the secular
     !! equation via the routine DLAED4 (as called by DLAED3).
     !! This routine also calculates the eigenvectors of the current
     !! problem.
     !! The final stage consists of computing the updated eigenvectors
     !! directly using the updated eigenvalues.  The eigenvectors for
     !! the current problem are multiplied with the eigenvectors from
     !! the overall problem.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(in) :: cutpnt, ldq, n
           integer(${ik}$), intent(out) :: info
           real(${rk}$), intent(inout) :: rho
           ! Array Arguments 
           integer(${ik}$), intent(inout) :: indxq(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(inout) :: d(*), q(ldq,*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           integer(${ik}$) :: coltyp, i, idlmda, indx, indxc, indxp, iq2, is, iw, iz, k, n1, n2, &
                     zpp1
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( n<0_${ik}$ ) then
              info = -1_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( min( 1_${ik}$, n / 2_${ik}$ )>cutpnt .or. ( n / 2_${ik}$ )<cutpnt ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DLAED1', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! the following values are integer pointers which indicate
           ! the portion of the workspace
           ! used by a particular array in stdlib${ii}$_${ri}$laed2 and stdlib${ii}$_${ri}$laed3.
           iz = 1_${ik}$
           idlmda = iz + n
           iw = idlmda + n
           iq2 = iw + n
           indx = 1_${ik}$
           indxc = indx + n
           coltyp = indxc + n
           indxp = coltyp + n
           ! form the z-vector which consists of the last row of q_1 and the
           ! first row of q_2.
           call stdlib${ii}$_${ri}$copy( cutpnt, q( cutpnt, 1_${ik}$ ), ldq, work( iz ), 1_${ik}$ )
           zpp1 = cutpnt + 1_${ik}$
           call stdlib${ii}$_${ri}$copy( n-cutpnt, q( zpp1, zpp1 ), ldq, work( iz+cutpnt ), 1_${ik}$ )
           ! deflate eigenvalues.
           call stdlib${ii}$_${ri}$laed2( k, n, cutpnt, d, q, ldq, indxq, rho, work( iz ),work( idlmda ), &
           work( iw ), work( iq2 ),iwork( indx ), iwork( indxc ), iwork( indxp ),iwork( coltyp ), &
                     info )
           if( info/=0 )go to 20
           ! solve secular equation.
           if( k/=0_${ik}$ ) then
              is = ( iwork( coltyp )+iwork( coltyp+1 ) )*cutpnt +( iwork( coltyp+1 )+iwork( &
                        coltyp+2 ) )*( n-cutpnt ) + iq2
              call stdlib${ii}$_${ri}$laed3( k, n, cutpnt, d, q, ldq, rho, work( idlmda ),work( iq2 ), iwork(&
                         indxc ), iwork( coltyp ),work( iw ), work( is ), info )
              if( info/=0 )go to 20
           ! prepare the indxq sorting permutation.
              n1 = k
              n2 = n - k
              call stdlib${ii}$_${ri}$lamrg( n1, n2, d, 1_${ik}$, -1_${ik}$, indxq )
           else
              do i = 1, n
                 indxq( i ) = i
              end do
           end if
           20 continue
           return
     end subroutine stdlib${ii}$_${ri}$laed1

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slaed2( k, n, n1, d, q, ldq, indxq, rho, z, dlamda, w,q2, indx, indxc,&
     !! SLAED2 merges the two sets of eigenvalues together into a single
     !! sorted set.  Then it tries to deflate the size of the problem.
     !! There are two ways in which deflation can occur:  when two or more
     !! eigenvalues are close together or if there is a tiny entry in the
     !! Z vector.  For each such occurrence the order of the related secular
     !! equation problem is reduced by one.
                indxp, coltyp, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info, k
           integer(${ik}$), intent(in) :: ldq, n, n1
           real(sp), intent(inout) :: rho
           ! Array Arguments 
           integer(${ik}$), intent(out) :: coltyp(*), indx(*), indxc(*), indxp(*)
           integer(${ik}$), intent(inout) :: indxq(*)
           real(sp), intent(inout) :: d(*), q(ldq,*), z(*)
           real(sp), intent(out) :: dlamda(*), q2(*), w(*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: mone = -1.0_sp
           
           ! Local Arrays 
           integer(${ik}$) :: ctot(4_${ik}$), psm(4_${ik}$)
           ! Local Scalars 
           integer(${ik}$) :: ct, i, imax, iq1, iq2, j, jmax, js, k2, n1p1, n2, nj, pj
           real(sp) :: c, eps, s, t, tau, tol
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ldq<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( min( 1_${ik}$, ( n / 2_${ik}$ ) )>n1 .or. ( n / 2_${ik}$ )<n1 ) then
              info = -3_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SLAED2', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           n2 = n - n1
           n1p1 = n1 + 1_${ik}$
           if( rho<zero ) then
              call stdlib${ii}$_sscal( n2, mone, z( n1p1 ), 1_${ik}$ )
           end if
           ! normalize z so that norm(z) = 1.  since z is the concatenation of
           ! two normalized vectors, norm2(z) = sqrt(2).
           t = one / sqrt( two )
           call stdlib${ii}$_sscal( n, t, z, 1_${ik}$ )
           ! rho = abs( norm(z)**2 * rho )
           rho = abs( two*rho )
           ! sort the eigenvalues into increasing order
           do i = n1p1, n
              indxq( i ) = indxq( i ) + n1
           end do
           ! re-integrate the deflated parts from the last pass
           do i = 1, n
              dlamda( i ) = d( indxq( i ) )
           end do
           call stdlib${ii}$_slamrg( n1, n2, dlamda, 1_${ik}$, 1_${ik}$, indxc )
           do i = 1, n
              indx( i ) = indxq( indxc( i ) )
           end do
           ! calculate the allowable deflation tolerance
           imax = stdlib${ii}$_isamax( n, z, 1_${ik}$ )
           jmax = stdlib${ii}$_isamax( n, d, 1_${ik}$ )
           eps = stdlib${ii}$_slamch( 'EPSILON' )
           tol = eight*eps*max( abs( d( jmax ) ), abs( z( imax ) ) )
           ! if the rank-1 modifier is small enough, no more needs to be done
           ! except to reorganize q so that its columns correspond with the
           ! elements in d.
           if( rho*abs( z( imax ) )<=tol ) then
              k = 0_${ik}$
              iq2 = 1_${ik}$
              do j = 1, n
                 i = indx( j )
                 call stdlib${ii}$_scopy( n, q( 1_${ik}$, i ), 1_${ik}$, q2( iq2 ), 1_${ik}$ )
                 dlamda( j ) = d( i )
                 iq2 = iq2 + n
              end do
              call stdlib${ii}$_slacpy( 'A', n, n, q2, n, q, ldq )
              call stdlib${ii}$_scopy( n, dlamda, 1_${ik}$, d, 1_${ik}$ )
              go to 190
           end if
           ! if there are multiple eigenvalues then the problem deflates.  here
           ! the number of equal eigenvalues are found.  as each equal
           ! eigenvalue is found, an elementary reflector is computed to rotate
           ! the corresponding eigensubspace so that the corresponding
           ! components of z are zero in this new basis.
           do i = 1, n1
              coltyp( i ) = 1_${ik}$
           end do
           do i = n1p1, n
              coltyp( i ) = 3_${ik}$
           end do
           k = 0_${ik}$
           k2 = n + 1_${ik}$
           do j = 1, n
              nj = indx( j )
              if( rho*abs( z( nj ) )<=tol ) then
                 ! deflate due to small z component.
                 k2 = k2 - 1_${ik}$
                 coltyp( nj ) = 4_${ik}$
                 indxp( k2 ) = nj
                 if( j==n )go to 100
              else
                 pj = nj
                 go to 80
              end if
           end do
           80 continue
           j = j + 1_${ik}$
           nj = indx( j )
           if( j>n )go to 100
           if( rho*abs( z( nj ) )<=tol ) then
              ! deflate due to small z component.
              k2 = k2 - 1_${ik}$
              coltyp( nj ) = 4_${ik}$
              indxp( k2 ) = nj
           else
              ! check if eigenvalues are close enough to allow deflation.
              s = z( pj )
              c = z( nj )
              ! find sqrt(a**2+b**2) without overflow or
              ! destructive underflow.
              tau = stdlib${ii}$_slapy2( c, s )
              t = d( nj ) - d( pj )
              c = c / tau
              s = -s / tau
              if( abs( t*c*s )<=tol ) then
                 ! deflation is possible.
                 z( nj ) = tau
                 z( pj ) = zero
                 if( coltyp( nj )/=coltyp( pj ) )coltyp( nj ) = 2_${ik}$
                 coltyp( pj ) = 4_${ik}$
                 call stdlib${ii}$_srot( n, q( 1_${ik}$, pj ), 1_${ik}$, q( 1_${ik}$, nj ), 1_${ik}$, c, s )
                 t = d( pj )*c**2_${ik}$ + d( nj )*s**2_${ik}$
                 d( nj ) = d( pj )*s**2_${ik}$ + d( nj )*c**2_${ik}$
                 d( pj ) = t
                 k2 = k2 - 1_${ik}$
                 i = 1_${ik}$
                 90 continue
                 if( k2+i<=n ) then
                    if( d( pj )<d( indxp( k2+i ) ) ) then
                       indxp( k2+i-1 ) = indxp( k2+i )
                       indxp( k2+i ) = pj
                       i = i + 1_${ik}$
                       go to 90
                    else
                       indxp( k2+i-1 ) = pj
                    end if
                 else
                    indxp( k2+i-1 ) = pj
                 end if
                 pj = nj
              else
                 k = k + 1_${ik}$
                 dlamda( k ) = d( pj )
                 w( k ) = z( pj )
                 indxp( k ) = pj
                 pj = nj
              end if
           end if
           go to 80
           100 continue
           ! record the last eigenvalue.
           k = k + 1_${ik}$
           dlamda( k ) = d( pj )
           w( k ) = z( pj )
           indxp( k ) = pj
           ! count up the total number of the various types of columns, then
           ! form a permutation which positions the four column types into
           ! four uniform groups (although one or more of these groups may be
           ! empty).
           do j = 1, 4
              ctot( j ) = 0_${ik}$
           end do
           do j = 1, n
              ct = coltyp( j )
              ctot( ct ) = ctot( ct ) + 1_${ik}$
           end do
           ! psm(*) = position in submatrix (of types 1 through 4)
           psm( 1_${ik}$ ) = 1_${ik}$
           psm( 2_${ik}$ ) = 1_${ik}$ + ctot( 1_${ik}$ )
           psm( 3_${ik}$ ) = psm( 2_${ik}$ ) + ctot( 2_${ik}$ )
           psm( 4_${ik}$ ) = psm( 3_${ik}$ ) + ctot( 3_${ik}$ )
           k = n - ctot( 4_${ik}$ )
           ! fill out the indxc array so that the permutation which it induces
           ! will place all type-1 columns first, all type-2 columns next,
           ! then all type-3's, and finally all type-4's.
           do j = 1, n
              js = indxp( j )
              ct = coltyp( js )
              indx( psm( ct ) ) = js
              indxc( psm( ct ) ) = j
              psm( ct ) = psm( ct ) + 1_${ik}$
           end do
           ! sort the eigenvalues and corresponding eigenvectors into dlamda
           ! and q2 respectively.  the eigenvalues/vectors which were not
           ! deflated go into the first k slots of dlamda and q2 respectively,
           ! while those which were deflated go into the last n - k slots.
           i = 1_${ik}$
           iq1 = 1_${ik}$
           iq2 = 1_${ik}$ + ( ctot( 1_${ik}$ )+ctot( 2_${ik}$ ) )*n1
           do j = 1, ctot( 1 )
              js = indx( i )
              call stdlib${ii}$_scopy( n1, q( 1_${ik}$, js ), 1_${ik}$, q2( iq1 ), 1_${ik}$ )
              z( i ) = d( js )
              i = i + 1_${ik}$
              iq1 = iq1 + n1
           end do
           do j = 1, ctot( 2 )
              js = indx( i )
              call stdlib${ii}$_scopy( n1, q( 1_${ik}$, js ), 1_${ik}$, q2( iq1 ), 1_${ik}$ )
              call stdlib${ii}$_scopy( n2, q( n1+1, js ), 1_${ik}$, q2( iq2 ), 1_${ik}$ )
              z( i ) = d( js )
              i = i + 1_${ik}$
              iq1 = iq1 + n1
              iq2 = iq2 + n2
           end do
           do j = 1, ctot( 3 )
              js = indx( i )
              call stdlib${ii}$_scopy( n2, q( n1+1, js ), 1_${ik}$, q2( iq2 ), 1_${ik}$ )
              z( i ) = d( js )
              i = i + 1_${ik}$
              iq2 = iq2 + n2
           end do
           iq1 = iq2
           do j = 1, ctot( 4 )
              js = indx( i )
              call stdlib${ii}$_scopy( n, q( 1_${ik}$, js ), 1_${ik}$, q2( iq2 ), 1_${ik}$ )
              iq2 = iq2 + n
              z( i ) = d( js )
              i = i + 1_${ik}$
           end do
           ! the deflated eigenvalues and their corresponding vectors go back
           ! into the last n - k slots of d and q respectively.
           if( k<n ) then
              call stdlib${ii}$_slacpy( 'A', n, ctot( 4_${ik}$ ), q2( iq1 ), n,q( 1_${ik}$, k+1 ), ldq )
              call stdlib${ii}$_scopy( n-k, z( k+1 ), 1_${ik}$, d( k+1 ), 1_${ik}$ )
           end if
           ! copy ctot into coltyp for referencing in stdlib${ii}$_slaed3.
           do j = 1, 4
              coltyp( j ) = ctot( j )
           end do
           190 continue
           return
     end subroutine stdlib${ii}$_slaed2

     pure module subroutine stdlib${ii}$_dlaed2( k, n, n1, d, q, ldq, indxq, rho, z, dlamda, w,q2, indx, indxc,&
     !! DLAED2 merges the two sets of eigenvalues together into a single
     !! sorted set.  Then it tries to deflate the size of the problem.
     !! There are two ways in which deflation can occur:  when two or more
     !! eigenvalues are close together or if there is a tiny entry in the
     !! Z vector.  For each such occurrence the order of the related secular
     !! equation problem is reduced by one.
                indxp, coltyp, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info, k
           integer(${ik}$), intent(in) :: ldq, n, n1
           real(dp), intent(inout) :: rho
           ! Array Arguments 
           integer(${ik}$), intent(out) :: coltyp(*), indx(*), indxc(*), indxp(*)
           integer(${ik}$), intent(inout) :: indxq(*)
           real(dp), intent(inout) :: d(*), q(ldq,*), z(*)
           real(dp), intent(out) :: dlamda(*), q2(*), w(*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: mone = -1.0_dp
           
           ! Local Arrays 
           integer(${ik}$) :: ctot(4_${ik}$