#:include "common.fypp"
submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_tridiag2
implicit none
contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
pure module subroutine stdlib${ii}$_slamrg( n1, n2, a, strd1, strd2, index )
!! SLAMRG will create a permutation list which will merge the elements
!! of A (which is composed of two independently sorted sets) into a
!! single set which is sorted in ascending order.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: n1, n2, strd1, strd2
! Array Arguments
integer(${ik}$), intent(out) :: index(*)
real(sp), intent(in) :: a(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, ind1, ind2, n1sv, n2sv
! Executable Statements
n1sv = n1
n2sv = n2
if( strd1>0_${ik}$ ) then
ind1 = 1_${ik}$
else
ind1 = n1
end if
if( strd2>0_${ik}$ ) then
ind2 = 1_${ik}$ + n1
else
ind2 = n1 + n2
end if
i = 1_${ik}$
! while ( (n1sv > 0)
10 continue
if( n1sv>0_${ik}$ .and. n2sv>0_${ik}$ ) then
if( a( ind1 )<=a( ind2 ) ) then
index( i ) = ind1
i = i + 1_${ik}$
ind1 = ind1 + strd1
n1sv = n1sv - 1_${ik}$
else
index( i ) = ind2
i = i + 1_${ik}$
ind2 = ind2 + strd2
n2sv = n2sv - 1_${ik}$
end if
go to 10
end if
! end while
if( n1sv==0_${ik}$ ) then
do n1sv = 1, n2sv
index( i ) = ind2
i = i + 1_${ik}$
ind2 = ind2 + strd2
end do
else
! n2sv == 0
do n2sv = 1, n1sv
index( i ) = ind1
i = i + 1_${ik}$
ind1 = ind1 + strd1
end do
end if
return
end subroutine stdlib${ii}$_slamrg
pure module subroutine stdlib${ii}$_dlamrg( n1, n2, a, dtrd1, dtrd2, index )
!! DLAMRG will create a permutation list which will merge the elements
!! of A (which is composed of two independently sorted sets) into a
!! single set which is sorted in ascending order.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: dtrd1, dtrd2, n1, n2
! Array Arguments
integer(${ik}$), intent(out) :: index(*)
real(dp), intent(in) :: a(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, ind1, ind2, n1sv, n2sv
! Executable Statements
n1sv = n1
n2sv = n2
if( dtrd1>0_${ik}$ ) then
ind1 = 1_${ik}$
else
ind1 = n1
end if
if( dtrd2>0_${ik}$ ) then
ind2 = 1_${ik}$ + n1
else
ind2 = n1 + n2
end if
i = 1_${ik}$
! while ( (n1sv > 0)
10 continue
if( n1sv>0_${ik}$ .and. n2sv>0_${ik}$ ) then
if( a( ind1 )<=a( ind2 ) ) then
index( i ) = ind1
i = i + 1_${ik}$
ind1 = ind1 + dtrd1
n1sv = n1sv - 1_${ik}$
else
index( i ) = ind2
i = i + 1_${ik}$
ind2 = ind2 + dtrd2
n2sv = n2sv - 1_${ik}$
end if
go to 10
end if
! end while
if( n1sv==0_${ik}$ ) then
do n1sv = 1, n2sv
index( i ) = ind2
i = i + 1_${ik}$
ind2 = ind2 + dtrd2
end do
else
! n2sv == 0
do n2sv = 1, n1sv
index( i ) = ind1
i = i + 1_${ik}$
ind1 = ind1 + dtrd1
end do
end if
return
end subroutine stdlib${ii}$_dlamrg
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$lamrg( n1, n2, a, dtrd1, dtrd2, index )
!! DLAMRG: will create a permutation list which will merge the elements
!! of A (which is composed of two independently sorted sets) into a
!! single set which is sorted in ascending order.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: dtrd1, dtrd2, n1, n2
! Array Arguments
integer(${ik}$), intent(out) :: index(*)
real(${rk}$), intent(in) :: a(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, ind1, ind2, n1sv, n2sv
! Executable Statements
n1sv = n1
n2sv = n2
if( dtrd1>0_${ik}$ ) then
ind1 = 1_${ik}$
else
ind1 = n1
end if
if( dtrd2>0_${ik}$ ) then
ind2 = 1_${ik}$ + n1
else
ind2 = n1 + n2
end if
i = 1_${ik}$
! while ( (n1sv > 0)
10 continue
if( n1sv>0_${ik}$ .and. n2sv>0_${ik}$ ) then
if( a( ind1 )<=a( ind2 ) ) then
index( i ) = ind1
i = i + 1_${ik}$
ind1 = ind1 + dtrd1
n1sv = n1sv - 1_${ik}$
else
index( i ) = ind2
i = i + 1_${ik}$
ind2 = ind2 + dtrd2
n2sv = n2sv - 1_${ik}$
end if
go to 10
end if
! end while
if( n1sv==0_${ik}$ ) then
do n1sv = 1, n2sv
index( i ) = ind2
i = i + 1_${ik}$
ind2 = ind2 + dtrd2
end do
else
! n2sv == 0
do n2sv = 1, n1sv
index( i ) = ind1
i = i + 1_${ik}$
ind1 = ind1 + dtrd1
end do
end if
return
end subroutine stdlib${ii}$_${ri}$lamrg
#:endif
#:endfor
pure module subroutine stdlib${ii}$_slaeda( n, tlvls, curlvl, curpbm, prmptr, perm, givptr,givcol, givnum,&
!! SLAEDA computes the Z vector corresponding to the merge step in the
!! CURLVLth step of the merge process with TLVLS steps for the CURPBMth
!! problem.
q, qptr, z, ztemp, info )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: curlvl, curpbm, n, tlvls
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: givcol(2_${ik}$,*), givptr(*), perm(*), prmptr(*), qptr(*)
real(sp), intent(in) :: givnum(2_${ik}$,*), q(*)
real(sp), intent(out) :: z(*), ztemp(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: bsiz1, bsiz2, curr, i, k, mid, psiz1, psiz2, ptr, zptr1
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SLAEDA', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine location of first number in second half.
mid = n / 2_${ik}$ + 1_${ik}$
! gather last/first rows of appropriate eigenblocks into center of z
ptr = 1_${ik}$
! determine location of lowest level subproblem in the full storage
! scheme
curr = ptr + curpbm*2_${ik}$**curlvl + 2_${ik}$**( curlvl-1 ) - 1_${ik}$
! determine size of these matrices. we add half to the value of
! the sqrt in case the machine underestimates one of these square
! roots.
bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=sp) ),KIND=${ik}$)
bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=sp) ),KIND=${ik}$)
do k = 1, mid - bsiz1 - 1
z( k ) = zero
end do
call stdlib${ii}$_scopy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,z( mid-bsiz1 ), 1_${ik}$ )
call stdlib${ii}$_scopy( bsiz2, q( qptr( curr+1 ) ), bsiz2, z( mid ), 1_${ik}$ )
do k = mid + bsiz2, n
z( k ) = zero
end do
! loop through remaining levels 1 -> curlvl applying the givens
! rotations and permutation and then multiplying the center matrices
! against the current z.
ptr = 2_${ik}$**tlvls + 1_${ik}$
loop_70: do k = 1, curlvl - 1
curr = ptr + curpbm*2_${ik}$**( curlvl-k ) + 2_${ik}$**( curlvl-k-1 ) - 1_${ik}$
psiz1 = prmptr( curr+1 ) - prmptr( curr )
psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
zptr1 = mid - psiz1
! apply givens at curr and curr+1
do i = givptr( curr ), givptr( curr+1 ) - 1
call stdlib${ii}$_srot( 1_${ik}$, z( zptr1+givcol( 1_${ik}$, i )-1_${ik}$ ), 1_${ik}$,z( zptr1+givcol( 2_${ik}$, i )-1_${ik}$ ), &
1_${ik}$, givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
end do
do i = givptr( curr+1 ), givptr( curr+2 ) - 1
call stdlib${ii}$_srot( 1_${ik}$, z( mid-1+givcol( 1_${ik}$, i ) ), 1_${ik}$,z( mid-1+givcol( 2_${ik}$, i ) ), 1_${ik}$, &
givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
end do
psiz1 = prmptr( curr+1 ) - prmptr( curr )
psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
do i = 0, psiz1 - 1
ztemp( i+1 ) = z( zptr1+perm( prmptr( curr )+i )-1_${ik}$ )
end do
do i = 0, psiz2 - 1
ztemp( psiz1+i+1 ) = z( mid+perm( prmptr( curr+1 )+i )-1_${ik}$ )
end do
! multiply blocks at curr and curr+1
! determine size of these matrices. we add half to the value of
! the sqrt in case the machine underestimates one of these
! square roots.
bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=sp) ),KIND=${ik}$)
bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=sp) ),KIND=${ik}$)
if( bsiz1>0_${ik}$ ) then
call stdlib${ii}$_sgemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),bsiz1, ztemp( 1_${ik}$ ), &
1_${ik}$, zero, z( zptr1 ), 1_${ik}$ )
end if
call stdlib${ii}$_scopy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1_${ik}$, z( zptr1+bsiz1 ),1_${ik}$ )
if( bsiz2>0_${ik}$ ) then
call stdlib${ii}$_sgemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),bsiz2, ztemp( &
psiz1+1 ), 1_${ik}$, zero, z( mid ), 1_${ik}$ )
end if
call stdlib${ii}$_scopy( psiz2-bsiz2, ztemp( psiz1+bsiz2+1 ), 1_${ik}$,z( mid+bsiz2 ), 1_${ik}$ )
ptr = ptr + 2_${ik}$**( tlvls-k )
end do loop_70
return
end subroutine stdlib${ii}$_slaeda
pure module subroutine stdlib${ii}$_dlaeda( n, tlvls, curlvl, curpbm, prmptr, perm, givptr,givcol, givnum,&
!! DLAEDA computes the Z vector corresponding to the merge step in the
!! CURLVLth step of the merge process with TLVLS steps for the CURPBMth
!! problem.
q, qptr, z, ztemp, info )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: curlvl, curpbm, n, tlvls
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: givcol(2_${ik}$,*), givptr(*), perm(*), prmptr(*), qptr(*)
real(dp), intent(in) :: givnum(2_${ik}$,*), q(*)
real(dp), intent(out) :: z(*), ztemp(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: bsiz1, bsiz2, curr, i, k, mid, psiz1, psiz2, ptr, zptr1
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DLAEDA', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine location of first number in second half.
mid = n / 2_${ik}$ + 1_${ik}$
! gather last/first rows of appropriate eigenblocks into center of z
ptr = 1_${ik}$
! determine location of lowest level subproblem in the full storage
! scheme
curr = ptr + curpbm*2_${ik}$**curlvl + 2_${ik}$**( curlvl-1 ) - 1_${ik}$
! determine size of these matrices. we add half to the value of
! the sqrt in case the machine underestimates one of these square
! roots.
bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=dp) ),KIND=${ik}$)
bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=dp) ),KIND=${ik}$)
do k = 1, mid - bsiz1 - 1
z( k ) = zero
end do
call stdlib${ii}$_dcopy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,z( mid-bsiz1 ), 1_${ik}$ )
call stdlib${ii}$_dcopy( bsiz2, q( qptr( curr+1 ) ), bsiz2, z( mid ), 1_${ik}$ )
do k = mid + bsiz2, n
z( k ) = zero
end do
! loop through remaining levels 1 -> curlvl applying the givens
! rotations and permutation and then multiplying the center matrices
! against the current z.
ptr = 2_${ik}$**tlvls + 1_${ik}$
loop_70: do k = 1, curlvl - 1
curr = ptr + curpbm*2_${ik}$**( curlvl-k ) + 2_${ik}$**( curlvl-k-1 ) - 1_${ik}$
psiz1 = prmptr( curr+1 ) - prmptr( curr )
psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
zptr1 = mid - psiz1
! apply givens at curr and curr+1
do i = givptr( curr ), givptr( curr+1 ) - 1
call stdlib${ii}$_drot( 1_${ik}$, z( zptr1+givcol( 1_${ik}$, i )-1_${ik}$ ), 1_${ik}$,z( zptr1+givcol( 2_${ik}$, i )-1_${ik}$ ), &
1_${ik}$, givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
end do
do i = givptr( curr+1 ), givptr( curr+2 ) - 1
call stdlib${ii}$_drot( 1_${ik}$, z( mid-1+givcol( 1_${ik}$, i ) ), 1_${ik}$,z( mid-1+givcol( 2_${ik}$, i ) ), 1_${ik}$, &
givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
end do
psiz1 = prmptr( curr+1 ) - prmptr( curr )
psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
do i = 0, psiz1 - 1
ztemp( i+1 ) = z( zptr1+perm( prmptr( curr )+i )-1_${ik}$ )
end do
do i = 0, psiz2 - 1
ztemp( psiz1+i+1 ) = z( mid+perm( prmptr( curr+1 )+i )-1_${ik}$ )
end do
! multiply blocks at curr and curr+1
! determine size of these matrices. we add half to the value of
! the sqrt in case the machine underestimates one of these
! square roots.
bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=dp) ),KIND=${ik}$)
bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=dp) ),KIND=${ik}$)
if( bsiz1>0_${ik}$ ) then
call stdlib${ii}$_dgemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),bsiz1, ztemp( 1_${ik}$ ), &
1_${ik}$, zero, z( zptr1 ), 1_${ik}$ )
end if
call stdlib${ii}$_dcopy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1_${ik}$, z( zptr1+bsiz1 ),1_${ik}$ )
if( bsiz2>0_${ik}$ ) then
call stdlib${ii}$_dgemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),bsiz2, ztemp( &
psiz1+1 ), 1_${ik}$, zero, z( mid ), 1_${ik}$ )
end if
call stdlib${ii}$_dcopy( psiz2-bsiz2, ztemp( psiz1+bsiz2+1 ), 1_${ik}$,z( mid+bsiz2 ), 1_${ik}$ )
ptr = ptr + 2_${ik}$**( tlvls-k )
end do loop_70
return
end subroutine stdlib${ii}$_dlaeda
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$laeda( n, tlvls, curlvl, curpbm, prmptr, perm, givptr,givcol, givnum,&
!! DLAEDA: computes the Z vector corresponding to the merge step in the
!! CURLVLth step of the merge process with TLVLS steps for the CURPBMth
!! problem.
q, qptr, z, ztemp, info )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: curlvl, curpbm, n, tlvls
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: givcol(2_${ik}$,*), givptr(*), perm(*), prmptr(*), qptr(*)
real(${rk}$), intent(in) :: givnum(2_${ik}$,*), q(*)
real(${rk}$), intent(out) :: z(*), ztemp(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: bsiz1, bsiz2, curr, i, k, mid, psiz1, psiz2, ptr, zptr1
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DLAEDA', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine location of first number in second half.
mid = n / 2_${ik}$ + 1_${ik}$
! gather last/first rows of appropriate eigenblocks into center of z
ptr = 1_${ik}$
! determine location of lowest level subproblem in the full storage
! scheme
curr = ptr + curpbm*2_${ik}$**curlvl + 2_${ik}$**( curlvl-1 ) - 1_${ik}$
! determine size of these matrices. we add half to the value of
! the sqrt in case the machine underestimates one of these square
! roots.
bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=${rk}$) ),KIND=${ik}$)
bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=${rk}$) ),KIND=${ik}$)
do k = 1, mid - bsiz1 - 1
z( k ) = zero
end do
call stdlib${ii}$_${ri}$copy( bsiz1, q( qptr( curr )+bsiz1-1 ), bsiz1,z( mid-bsiz1 ), 1_${ik}$ )
call stdlib${ii}$_${ri}$copy( bsiz2, q( qptr( curr+1 ) ), bsiz2, z( mid ), 1_${ik}$ )
do k = mid + bsiz2, n
z( k ) = zero
end do
! loop through remaining levels 1 -> curlvl applying the givens
! rotations and permutation and then multiplying the center matrices
! against the current z.
ptr = 2_${ik}$**tlvls + 1_${ik}$
loop_70: do k = 1, curlvl - 1
curr = ptr + curpbm*2_${ik}$**( curlvl-k ) + 2_${ik}$**( curlvl-k-1 ) - 1_${ik}$
psiz1 = prmptr( curr+1 ) - prmptr( curr )
psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
zptr1 = mid - psiz1
! apply givens at curr and curr+1
do i = givptr( curr ), givptr( curr+1 ) - 1
call stdlib${ii}$_${ri}$rot( 1_${ik}$, z( zptr1+givcol( 1_${ik}$, i )-1_${ik}$ ), 1_${ik}$,z( zptr1+givcol( 2_${ik}$, i )-1_${ik}$ ), &
1_${ik}$, givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
end do
do i = givptr( curr+1 ), givptr( curr+2 ) - 1
call stdlib${ii}$_${ri}$rot( 1_${ik}$, z( mid-1+givcol( 1_${ik}$, i ) ), 1_${ik}$,z( mid-1+givcol( 2_${ik}$, i ) ), 1_${ik}$, &
givnum( 1_${ik}$, i ),givnum( 2_${ik}$, i ) )
end do
psiz1 = prmptr( curr+1 ) - prmptr( curr )
psiz2 = prmptr( curr+2 ) - prmptr( curr+1 )
do i = 0, psiz1 - 1
ztemp( i+1 ) = z( zptr1+perm( prmptr( curr )+i )-1_${ik}$ )
end do
do i = 0, psiz2 - 1
ztemp( psiz1+i+1 ) = z( mid+perm( prmptr( curr+1 )+i )-1_${ik}$ )
end do
! multiply blocks at curr and curr+1
! determine size of these matrices. we add half to the value of
! the sqrt in case the machine underestimates one of these
! square roots.
bsiz1 = int( half+sqrt( real( qptr( curr+1 )-qptr( curr ),KIND=${rk}$) ),KIND=${ik}$)
bsiz2 = int( half+sqrt( real( qptr( curr+2 )-qptr( curr+1 ),KIND=${rk}$) ),KIND=${ik}$)
if( bsiz1>0_${ik}$ ) then
call stdlib${ii}$_${ri}$gemv( 'T', bsiz1, bsiz1, one, q( qptr( curr ) ),bsiz1, ztemp( 1_${ik}$ ), &
1_${ik}$, zero, z( zptr1 ), 1_${ik}$ )
end if
call stdlib${ii}$_${ri}$copy( psiz1-bsiz1, ztemp( bsiz1+1 ), 1_${ik}$, z( zptr1+bsiz1 ),1_${ik}$ )
if( bsiz2>0_${ik}$ ) then
call stdlib${ii}$_${ri}$gemv( 'T', bsiz2, bsiz2, one, q( qptr( curr+1 ) ),bsiz2, ztemp( &
psiz1+1 ), 1_${ik}$, zero, z( mid ), 1_${ik}$ )
end if
call stdlib${ii}$_${ri}$copy( psiz2-bsiz2, ztemp( psiz1+bsiz2+1 ), 1_${ik}$,z( mid+bsiz2 ), 1_${ik}$ )
ptr = ptr + 2_${ik}$**( tlvls-k )
end do loop_70
return
end subroutine stdlib${ii}$_${ri}$laeda
#:endif
#:endfor
pure module subroutine stdlib${ii}$_slarra( n, d, e, e2, spltol, tnrm,nsplit, isplit, info )
!! Compute the splitting points with threshold SPLTOL.
!! SLARRA sets any "small" off-diagonal elements to zero.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(out) :: info, nsplit
integer(${ik}$), intent(in) :: n
real(sp), intent(in) :: spltol, tnrm
! Array Arguments
integer(${ik}$), intent(out) :: isplit(*)
real(sp), intent(in) :: d(*)
real(sp), intent(inout) :: e(*), e2(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i
real(sp) :: eabs, tmp1
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
! compute splitting points
nsplit = 1_${ik}$
if(spltol<zero) then
! criterion based on absolute off-diagonal value
tmp1 = abs(spltol)* tnrm
do i = 1, n-1
eabs = abs( e(i) )
if( eabs <= tmp1) then
e(i) = zero
e2(i) = zero
isplit( nsplit ) = i
nsplit = nsplit + 1_${ik}$
end if
end do
else
! criterion that guarantees relative accuracy
do i = 1, n-1
eabs = abs( e(i) )
if( eabs <= spltol * sqrt(abs(d(i)))*sqrt(abs(d(i+1))) )then
e(i) = zero
e2(i) = zero
isplit( nsplit ) = i
nsplit = nsplit + 1_${ik}$
end if
end do
endif
isplit( nsplit ) = n
return
end subroutine stdlib${ii}$_slarra
pure module subroutine stdlib${ii}$_dlarra( n, d, e, e2, spltol, tnrm,nsplit, isplit, info )
!! Compute the splitting points with threshold SPLTOL.
!! DLARRA sets any "small" off-diagonal elements to zero.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(out) :: info, nsplit
integer(${ik}$), intent(in) :: n
real(dp), intent(in) :: spltol, tnrm
! Array Arguments
integer(${ik}$), intent(out) :: isplit(*)
real(dp), intent(in) :: d(*)
real(dp), intent(inout) :: e(*), e2(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i
real(dp) :: eabs, tmp1
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
! compute splitting points
nsplit = 1_${ik}$
if(spltol<zero) then
! criterion based on absolute off-diagonal value
tmp1 = abs(spltol)* tnrm
do i = 1, n-1
eabs = abs( e(i) )
if( eabs <= tmp1) then
e(i) = zero
e2(i) = zero
isplit( nsplit ) = i
nsplit = nsplit + 1_${ik}$
end if
end do
else
! criterion that guarantees relative accuracy
do i = 1, n-1
eabs = abs( e(i) )
if( eabs <= spltol * sqrt(abs(d(i)))*sqrt(abs(d(i+1))) )then
e(i) = zero
e2(i) = zero
isplit( nsplit ) = i
nsplit = nsplit + 1_${ik}$
end if
end do
endif
isplit( nsplit ) = n
return
end subroutine stdlib${ii}$_dlarra
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$larra( n, d, e, e2, spltol, tnrm,nsplit, isplit, info )
!! Compute the splitting points with threshold SPLTOL.
!! DLARRA: sets any "small" off-diagonal elements to zero.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(out) :: info, nsplit
integer(${ik}$), intent(in) :: n
real(${rk}$), intent(in) :: spltol, tnrm
! Array Arguments
integer(${ik}$), intent(out) :: isplit(*)
real(${rk}$), intent(in) :: d(*)
real(${rk}$), intent(inout) :: e(*), e2(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i
real(${rk}$) :: eabs, tmp1
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
! compute splitting points
nsplit = 1_${ik}$
if(spltol<zero) then
! criterion based on absolute off-diagonal value
tmp1 = abs(spltol)* tnrm
do i = 1, n-1
eabs = abs( e(i) )
if( eabs <= tmp1) then
e(i) = zero
e2(i) = zero
isplit( nsplit ) = i
nsplit = nsplit + 1_${ik}$
end if
end do
else
! criterion that guarantees relative accuracy
do i = 1, n-1
eabs = abs( e(i) )
if( eabs <= spltol * sqrt(abs(d(i)))*sqrt(abs(d(i+1))) )then
e(i) = zero
e2(i) = zero
isplit( nsplit ) = i
nsplit = nsplit + 1_${ik}$
end if
end do
endif
isplit( nsplit ) = n
return
end subroutine stdlib${ii}$_${ri}$larra
#:endif
#:endfor
pure module subroutine stdlib${ii}$_slarrb( n, d, lld, ifirst, ilast, rtol1,rtol2, offset, w, wgap, werr, &
!! Given the relatively robust representation(RRR) L D L^T, SLARRB:
!! does "limited" bisection to refine the eigenvalues of L D L^T,
!! W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
!! guesses for these eigenvalues are input in W, the corresponding estimate
!! of the error in these guesses and their gaps are input in WERR
!! and WGAP, respectively. During bisection, intervals
!! [left, right] are maintained by storing their mid-points and
!! semi-widths in the arrays W and WERR respectively.
work, iwork,pivmin, spdiam, twist, info )
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: ifirst, ilast, n, offset, twist
integer(${ik}$), intent(out) :: info
real(sp), intent(in) :: pivmin, rtol1, rtol2, spdiam
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(in) :: d(*), lld(*)
real(sp), intent(inout) :: w(*), werr(*), wgap(*)
real(sp), intent(out) :: work(*)
! =====================================================================
integer(${ik}$) :: maxitr
! Local Scalars
integer(${ik}$) :: i, i1, ii, ip, iter, k, negcnt, next, nint, olnint, prev, r
real(sp) :: back, cvrgd, gap, left, lgap, mid, mnwdth, rgap, right, tmp, width
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
mnwdth = two * pivmin
r = twist
if((r<1_${ik}$).or.(r>n)) r = n
! initialize unconverged intervals in [ work(2*i-1), work(2*i) ].
! the sturm count, count( work(2*i-1) ) is arranged to be i-1, while
! count( work(2*i) ) is stored in iwork( 2*i ). the integer iwork( 2*i-1 )
! for an unconverged interval is set to the index of the next unconverged
! interval, and is -1 or 0 for a converged interval. thus a linked
! list of unconverged intervals is set up.
i1 = ifirst
! the number of unconverged intervals
nint = 0_${ik}$
! the last unconverged interval found
prev = 0_${ik}$
rgap = wgap( i1-offset )
loop_75: do i = i1, ilast
k = 2_${ik}$*i
ii = i - offset
left = w( ii ) - werr( ii )
right = w( ii ) + werr( ii )
lgap = rgap
rgap = wgap( ii )
gap = min( lgap, rgap )
! make sure that [left,right] contains the desired eigenvalue
! compute negcount from dstqds facto l+d+l+^t = l d l^t - left
! do while( negcnt(left)>i-1 )
back = werr( ii )
20 continue
negcnt = stdlib${ii}$_slaneg( n, d, lld, left, pivmin, r )
if( negcnt>i-1 ) then
left = left - back
back = two*back
go to 20
end if
! do while( negcnt(right)<i )
! compute negcount from dstqds facto l+d+l+^t = l d l^t - right
back = werr( ii )
50 continue
negcnt = stdlib${ii}$_slaneg( n, d, lld, right, pivmin, r )
if( negcnt<i ) then
right = right + back
back = two*back
go to 50
end if
width = half*abs( left - right )
tmp = max( abs( left ), abs( right ) )
cvrgd = max(rtol1*gap,rtol2*tmp)
if( width<=cvrgd .or. width<=mnwdth ) then
! this interval has already converged and does not need refinement.
! (note that the gaps might change through refining the
! eigenvalues, however, they can only get bigger.)
! remove it from the list.
iwork( k-1 ) = -1_${ik}$
! make sure that i1 always points to the first unconverged interval
if((i==i1).and.(i<ilast)) i1 = i + 1_${ik}$
if((prev>=i1).and.(i<=ilast)) iwork( 2_${ik}$*prev-1 ) = i + 1_${ik}$
else
! unconverged interval found
prev = i
nint = nint + 1_${ik}$
iwork( k-1 ) = i + 1_${ik}$
iwork( k ) = negcnt
end if
work( k-1 ) = left
work( k ) = right
end do loop_75
! do while( nint>0 ), i.e. there are still unconverged intervals
! and while (iter<maxitr)
iter = 0_${ik}$
80 continue
prev = i1 - 1_${ik}$
i = i1
olnint = nint
loop_100: do ip = 1, olnint
k = 2_${ik}$*i
ii = i - offset
rgap = wgap( ii )
lgap = rgap
if(ii>1_${ik}$) lgap = wgap( ii-1 )
gap = min( lgap, rgap )
next = iwork( k-1 )
left = work( k-1 )
right = work( k )
mid = half*( left + right )
! semiwidth of interval
width = right - mid
tmp = max( abs( left ), abs( right ) )
cvrgd = max(rtol1*gap,rtol2*tmp)
if( ( width<=cvrgd ) .or. ( width<=mnwdth ).or.( iter==maxitr ) )then
! reduce number of unconverged intervals
nint = nint - 1_${ik}$
! mark interval as converged.
iwork( k-1 ) = 0_${ik}$
if( i1==i ) then
i1 = next
else
! prev holds the last unconverged interval previously examined
if(prev>=i1) iwork( 2_${ik}$*prev-1 ) = next
end if
i = next
cycle loop_100
end if
prev = i
! perform one bisection step
negcnt = stdlib${ii}$_slaneg( n, d, lld, mid, pivmin, r )
if( negcnt<=i-1 ) then
work( k-1 ) = mid
else
work( k ) = mid
end if
i = next
end do loop_100
iter = iter + 1_${ik}$
! do another loop if there are still unconverged intervals
! however, in the last iteration, all intervals are accepted
! since this is the best we can do.
if( ( nint>0 ).and.(iter<=maxitr) ) go to 80
! at this point, all the intervals have converged
do i = ifirst, ilast
k = 2_${ik}$*i
ii = i - offset
! all intervals marked by '0' have been refined.
if( iwork( k-1 )==0_${ik}$ ) then
w( ii ) = half*( work( k-1 )+work( k ) )
werr( ii ) = work( k ) - w( ii )
end if
end do
do i = ifirst+1, ilast
k = 2_${ik}$*i
ii = i - offset
wgap( ii-1 ) = max( zero,w(ii) - werr (ii) - w( ii-1 ) - werr( ii-1 ))
end do
return
end subroutine stdlib${ii}$_slarrb
pure module subroutine stdlib${ii}$_dlarrb( n, d, lld, ifirst, ilast, rtol1,rtol2, offset, w, wgap, werr, &
!! Given the relatively robust representation(RRR) L D L^T, DLARRB:
!! does "limited" bisection to refine the eigenvalues of L D L^T,
!! W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
!! guesses for these eigenvalues are input in W, the corresponding estimate
!! of the error in these guesses and their gaps are input in WERR
!! and WGAP, respectively. During bisection, intervals
!! [left, right] are maintained by storing their mid-points and
!! semi-widths in the arrays W and WERR respectively.
work, iwork,pivmin, spdiam, twist, info )
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: ifirst, ilast, n, offset, twist
integer(${ik}$), intent(out) :: info
real(dp), intent(in) :: pivmin, rtol1, rtol2, spdiam
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(in) :: d(*), lld(*)
real(dp), intent(inout) :: w(*), werr(*), wgap(*)
real(dp), intent(out) :: work(*)
! =====================================================================
integer(${ik}$) :: maxitr
! Local Scalars
integer(${ik}$) :: i, i1, ii, ip, iter, k, negcnt, next, nint, olnint, prev, r
real(dp) :: back, cvrgd, gap, left, lgap, mid, mnwdth, rgap, right, tmp, width
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
mnwdth = two * pivmin
r = twist
if((r<1_${ik}$).or.(r>n)) r = n
! initialize unconverged intervals in [ work(2*i-1), work(2*i) ].
! the sturm count, count( work(2*i-1) ) is arranged to be i-1, while
! count( work(2*i) ) is stored in iwork( 2*i ). the integer iwork( 2*i-1 )
! for an unconverged interval is set to the index of the next unconverged
! interval, and is -1 or 0 for a converged interval. thus a linked
! list of unconverged intervals is set up.
i1 = ifirst
! the number of unconverged intervals
nint = 0_${ik}$
! the last unconverged interval found
prev = 0_${ik}$
rgap = wgap( i1-offset )
loop_75: do i = i1, ilast
k = 2_${ik}$*i
ii = i - offset
left = w( ii ) - werr( ii )
right = w( ii ) + werr( ii )
lgap = rgap
rgap = wgap( ii )
gap = min( lgap, rgap )
! make sure that [left,right] contains the desired eigenvalue
! compute negcount from dstqds facto l+d+l+^t = l d l^t - left
! do while( negcnt(left)>i-1 )
back = werr( ii )
20 continue
negcnt = stdlib${ii}$_dlaneg( n, d, lld, left, pivmin, r )
if( negcnt>i-1 ) then
left = left - back
back = two*back
go to 20
end if
! do while( negcnt(right)<i )
! compute negcount from dstqds facto l+d+l+^t = l d l^t - right
back = werr( ii )
50 continue
negcnt = stdlib${ii}$_dlaneg( n, d, lld, right, pivmin, r )
if( negcnt<i ) then
right = right + back
back = two*back
go to 50
end if
width = half*abs( left - right )
tmp = max( abs( left ), abs( right ) )
cvrgd = max(rtol1*gap,rtol2*tmp)
if( width<=cvrgd .or. width<=mnwdth ) then
! this interval has already converged and does not need refinement.
! (note that the gaps might change through refining the
! eigenvalues, however, they can only get bigger.)
! remove it from the list.
iwork( k-1 ) = -1_${ik}$
! make sure that i1 always points to the first unconverged interval
if((i==i1).and.(i<ilast)) i1 = i + 1_${ik}$
if((prev>=i1).and.(i<=ilast)) iwork( 2_${ik}$*prev-1 ) = i + 1_${ik}$
else
! unconverged interval found
prev = i
nint = nint + 1_${ik}$
iwork( k-1 ) = i + 1_${ik}$
iwork( k ) = negcnt
end if
work( k-1 ) = left
work( k ) = right
end do loop_75
! do while( nint>0 ), i.e. there are still unconverged intervals
! and while (iter<maxitr)
iter = 0_${ik}$
80 continue
prev = i1 - 1_${ik}$
i = i1
olnint = nint
loop_100: do ip = 1, olnint
k = 2_${ik}$*i
ii = i - offset
rgap = wgap( ii )
lgap = rgap
if(ii>1_${ik}$) lgap = wgap( ii-1 )
gap = min( lgap, rgap )
next = iwork( k-1 )
left = work( k-1 )
right = work( k )
mid = half*( left + right )
! semiwidth of interval
width = right - mid
tmp = max( abs( left ), abs( right ) )
cvrgd = max(rtol1*gap,rtol2*tmp)
if( ( width<=cvrgd ) .or. ( width<=mnwdth ).or.( iter==maxitr ) )then
! reduce number of unconverged intervals
nint = nint - 1_${ik}$
! mark interval as converged.
iwork( k-1 ) = 0_${ik}$
if( i1==i ) then
i1 = next
else
! prev holds the last unconverged interval previously examined
if(prev>=i1) iwork( 2_${ik}$*prev-1 ) = next
end if
i = next
cycle loop_100
end if
prev = i
! perform one bisection step
negcnt = stdlib${ii}$_dlaneg( n, d, lld, mid, pivmin, r )
if( negcnt<=i-1 ) then
work( k-1 ) = mid
else
work( k ) = mid
end if
i = next
end do loop_100
iter = iter + 1_${ik}$
! do another loop if there are still unconverged intervals
! however, in the last iteration, all intervals are accepted
! since this is the best we can do.
if( ( nint>0 ).and.(iter<=maxitr) ) go to 80
! at this point, all the intervals have converged
do i = ifirst, ilast
k = 2_${ik}$*i
ii = i - offset
! all intervals marked by '0' have been refined.
if( iwork( k-1 )==0_${ik}$ ) then
w( ii ) = half*( work( k-1 )+work( k ) )
werr( ii ) = work( k ) - w( ii )
end if
end do
do i = ifirst+1, ilast
k = 2_${ik}$*i
ii = i - offset
wgap( ii-1 ) = max( zero,w(ii) - werr (ii) - w( ii-1 ) - werr( ii-1 ))
end do
return
end subroutine stdlib${ii}$_dlarrb
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$larrb( n, d, lld, ifirst, ilast, rtol1,rtol2, offset, w, wgap, werr, &
!! Given the relatively robust representation(RRR) L D L^T, DLARRB:
!! does "limited" bisection to refine the eigenvalues of L D L^T,
!! W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
!! guesses for these eigenvalues are input in W, the corresponding estimate
!! of the error in these guesses and their gaps are input in WERR
!! and WGAP, respectively. During bisection, intervals
!! [left, right] are maintained by storing their mid-points and
!! semi-widths in the arrays W and WERR respectively.
work, iwork,pivmin, spdiam, twist, info )
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: ifirst, ilast, n, offset, twist
integer(${ik}$), intent(out) :: info
real(${rk}$), intent(in) :: pivmin, rtol1, rtol2, spdiam
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(in) :: d(*), lld(*)
real(${rk}$), intent(inout) :: w(*), werr(*), wgap(*)
real(${rk}$), intent(out) :: work(*)
! =====================================================================
integer(${ik}$) :: maxitr
! Local Scalars
integer(${ik}$) :: i, i1, ii, ip, iter, k, negcnt, next, nint, olnint, prev, r
real(${rk}$) :: back, cvrgd, gap, left, lgap, mid, mnwdth, rgap, right, tmp, width
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
mnwdth = two * pivmin
r = twist
if((r<1_${ik}$).or.(r>n)) r = n
! initialize unconverged intervals in [ work(2*i-1), work(2*i) ].
! the sturm count, count( work(2*i-1) ) is arranged to be i-1, while
! count( work(2*i) ) is stored in iwork( 2*i ). the integer iwork( 2*i-1 )
! for an unconverged interval is set to the index of the next unconverged
! interval, and is -1 or 0 for a converged interval. thus a linked
! list of unconverged intervals is set up.
i1 = ifirst
! the number of unconverged intervals
nint = 0_${ik}$
! the last unconverged interval found
prev = 0_${ik}$
rgap = wgap( i1-offset )
loop_75: do i = i1, ilast
k = 2_${ik}$*i
ii = i - offset
left = w( ii ) - werr( ii )
right = w( ii ) + werr( ii )
lgap = rgap
rgap = wgap( ii )
gap = min( lgap, rgap )
! make sure that [left,right] contains the desired eigenvalue
! compute negcount from dstqds facto l+d+l+^t = l d l^t - left
! do while( negcnt(left)>i-1 )
back = werr( ii )
20 continue
negcnt = stdlib${ii}$_${ri}$laneg( n, d, lld, left, pivmin, r )
if( negcnt>i-1 ) then
left = left - back
back = two*back
go to 20
end if
! do while( negcnt(right)<i )
! compute negcount from dstqds facto l+d+l+^t = l d l^t - right
back = werr( ii )
50 continue
negcnt = stdlib${ii}$_${ri}$laneg( n, d, lld, right, pivmin, r )
if( negcnt<i ) then
right = right + back
back = two*back
go to 50
end if
width = half*abs( left - right )
tmp = max( abs( left ), abs( right ) )
cvrgd = max(rtol1*gap,rtol2*tmp)
if( width<=cvrgd .or. width<=mnwdth ) then
! this interval has already converged and does not need refinement.
! (note that the gaps might change through refining the
! eigenvalues, however, they can only get bigger.)
! remove it from the list.
iwork( k-1 ) = -1_${ik}$
! make sure that i1 always points to the first unconverged interval
if((i==i1).and.(i<ilast)) i1 = i + 1_${ik}$
if((prev>=i1).and.(i<=ilast)) iwork( 2_${ik}$*prev-1 ) = i + 1_${ik}$
else
! unconverged interval found
prev = i
nint = nint + 1_${ik}$
iwork( k-1 ) = i + 1_${ik}$
iwork( k ) = negcnt
end if
work( k-1 ) = left
work( k ) = right
end do loop_75
! do while( nint>0 ), i.e. there are still unconverged intervals
! and while (iter<maxitr)
iter = 0_${ik}$
80 continue
prev = i1 - 1_${ik}$
i = i1
olnint = nint
loop_100: do ip = 1, olnint
k = 2_${ik}$*i
ii = i - offset
rgap = wgap( ii )
lgap = rgap
if(ii>1_${ik}$) lgap = wgap( ii-1 )
gap = min( lgap, rgap )
next = iwork( k-1 )
left = work( k-1 )
right = work( k )
mid = half*( left + right )
! semiwidth of interval
width = right - mid
tmp = max( abs( left ), abs( right ) )
cvrgd = max(rtol1*gap,rtol2*tmp)
if( ( width<=cvrgd ) .or. ( width<=mnwdth ).or.( iter==maxitr ) )then
! reduce number of unconverged intervals
nint = nint - 1_${ik}$
! mark interval as converged.
iwork( k-1 ) = 0_${ik}$
if( i1==i ) then
i1 = next
else
! prev holds the last unconverged interval previously examined
if(prev>=i1) iwork( 2_${ik}$*prev-1 ) = next
end if
i = next
cycle loop_100
end if
prev = i
! perform one bisection step
negcnt = stdlib${ii}$_${ri}$laneg( n, d, lld, mid, pivmin, r )
if( negcnt<=i-1 ) then
work( k-1 ) = mid
else
work( k ) = mid
end if
i = next
end do loop_100
iter = iter + 1_${ik}$
! do another loop if there are still unconverged intervals
! however, in the last iteration, all intervals are accepted
! since this is the best we can do.
if( ( nint>0 ).and.(iter<=maxitr) ) go to 80
! at this point, all the intervals have converged
do i = ifirst, ilast
k = 2_${ik}$*i
ii = i - offset
! all intervals marked by '0' have been refined.
if( iwork( k-1 )==0_${ik}$ ) then
w( ii ) = half*( work( k-1 )+work( k ) )
werr( ii ) = work( k ) - w( ii )
end if
end do
do i = ifirst+1, ilast
k = 2_${ik}$*i
ii = i - offset
wgap( ii-1 ) = max( zero,w(ii) - werr (ii) - w( ii-1 ) - werr( ii-1 ))
end do
return
end subroutine stdlib${ii}$_${ri}$larrb
#:endif
#:endfor
pure module subroutine stdlib${ii}$_slarrc( jobt, n, vl, vu, d, e, pivmin,eigcnt, lcnt, rcnt, info )
!! Find the number of eigenvalues of the symmetric tridiagonal matrix T
!! that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T
!! if JOBT = 'L'.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobt
integer(${ik}$), intent(out) :: eigcnt, info, lcnt, rcnt
integer(${ik}$), intent(in) :: n
real(sp), intent(in) :: pivmin, vl, vu
! Array Arguments
real(sp), intent(in) :: d(*), e(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i
logical(lk) :: matt
real(sp) :: lpivot, rpivot, sl, su, tmp, tmp2
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
lcnt = 0_${ik}$
rcnt = 0_${ik}$
eigcnt = 0_${ik}$
matt = stdlib_lsame( jobt, 'T' )
if (matt) then
! sturm sequence count on t
lpivot = d( 1_${ik}$ ) - vl
rpivot = d( 1_${ik}$ ) - vu
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
do i = 1, n-1
tmp = e(i)**2_${ik}$
lpivot = ( d( i+1 )-vl ) - tmp/lpivot
rpivot = ( d( i+1 )-vu ) - tmp/rpivot
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
end do
else
! sturm sequence count on l d l^t
sl = -vl
su = -vu
do i = 1, n - 1
lpivot = d( i ) + sl
rpivot = d( i ) + su
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
tmp = e(i) * d(i) * e(i)
tmp2 = tmp / lpivot
if( tmp2==zero ) then
sl = tmp - vl
else
sl = sl*tmp2 - vl
end if
tmp2 = tmp / rpivot
if( tmp2==zero ) then
su = tmp - vu
else
su = su*tmp2 - vu
end if
end do
lpivot = d( n ) + sl
rpivot = d( n ) + su
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
endif
eigcnt = rcnt - lcnt
return
end subroutine stdlib${ii}$_slarrc
pure module subroutine stdlib${ii}$_dlarrc( jobt, n, vl, vu, d, e, pivmin,eigcnt, lcnt, rcnt, info )
!! Find the number of eigenvalues of the symmetric tridiagonal matrix T
!! that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T
!! if JOBT = 'L'.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobt
integer(${ik}$), intent(out) :: eigcnt, info, lcnt, rcnt
integer(${ik}$), intent(in) :: n
real(dp), intent(in) :: pivmin, vl, vu
! Array Arguments
real(dp), intent(in) :: d(*), e(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i
logical(lk) :: matt
real(dp) :: lpivot, rpivot, sl, su, tmp, tmp2
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
lcnt = 0_${ik}$
rcnt = 0_${ik}$
eigcnt = 0_${ik}$
matt = stdlib_lsame( jobt, 'T' )
if (matt) then
! sturm sequence count on t
lpivot = d( 1_${ik}$ ) - vl
rpivot = d( 1_${ik}$ ) - vu
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
do i = 1, n-1
tmp = e(i)**2_${ik}$
lpivot = ( d( i+1 )-vl ) - tmp/lpivot
rpivot = ( d( i+1 )-vu ) - tmp/rpivot
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
end do
else
! sturm sequence count on l d l^t
sl = -vl
su = -vu
do i = 1, n - 1
lpivot = d( i ) + sl
rpivot = d( i ) + su
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
tmp = e(i) * d(i) * e(i)
tmp2 = tmp / lpivot
if( tmp2==zero ) then
sl = tmp - vl
else
sl = sl*tmp2 - vl
end if
tmp2 = tmp / rpivot
if( tmp2==zero ) then
su = tmp - vu
else
su = su*tmp2 - vu
end if
end do
lpivot = d( n ) + sl
rpivot = d( n ) + su
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
endif
eigcnt = rcnt - lcnt
return
end subroutine stdlib${ii}$_dlarrc
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$larrc( jobt, n, vl, vu, d, e, pivmin,eigcnt, lcnt, rcnt, info )
!! Find the number of eigenvalues of the symmetric tridiagonal matrix T
!! that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T
!! if JOBT = 'L'.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobt
integer(${ik}$), intent(out) :: eigcnt, info, lcnt, rcnt
integer(${ik}$), intent(in) :: n
real(${rk}$), intent(in) :: pivmin, vl, vu
! Array Arguments
real(${rk}$), intent(in) :: d(*), e(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i
logical(lk) :: matt
real(${rk}$) :: lpivot, rpivot, sl, su, tmp, tmp2
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
lcnt = 0_${ik}$
rcnt = 0_${ik}$
eigcnt = 0_${ik}$
matt = stdlib_lsame( jobt, 'T' )
if (matt) then
! sturm sequence count on t
lpivot = d( 1_${ik}$ ) - vl
rpivot = d( 1_${ik}$ ) - vu
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
do i = 1, n-1
tmp = e(i)**2_${ik}$
lpivot = ( d( i+1 )-vl ) - tmp/lpivot
rpivot = ( d( i+1 )-vu ) - tmp/rpivot
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
end do
else
! sturm sequence count on l d l^t
sl = -vl
su = -vu
do i = 1, n - 1
lpivot = d( i ) + sl
rpivot = d( i ) + su
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
tmp = e(i) * d(i) * e(i)
tmp2 = tmp / lpivot
if( tmp2==zero ) then
sl = tmp - vl
else
sl = sl*tmp2 - vl
end if
tmp2 = tmp / rpivot
if( tmp2==zero ) then
su = tmp - vu
else
su = su*tmp2 - vu
end if
end do
lpivot = d( n ) + sl
rpivot = d( n ) + su
if( lpivot<=zero ) then
lcnt = lcnt + 1_${ik}$
endif
if( rpivot<=zero ) then
rcnt = rcnt + 1_${ik}$
endif
endif
eigcnt = rcnt - lcnt
return
end subroutine stdlib${ii}$_${ri}$larrc
#:endif
#:endfor
pure module subroutine stdlib${ii}$_slarrd( range, order, n, vl, vu, il, iu, gers,reltol, d, e, e2, &
!! SLARRD computes the eigenvalues of a symmetric tridiagonal
!! matrix T to suitable accuracy. This is an auxiliary code to be
!! called from SSTEMR.
!! The user may ask for all eigenvalues, all eigenvalues
!! in the half-open interval (VL, VU], or the IL-th through IU-th
!! eigenvalues.
!! To avoid overflow, the matrix must be scaled so that its
!! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
!! accuracy, it should not be much smaller than that.
!! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
!! Matrix", Report CS41, Computer Science Dept., Stanford
!! University, July 21, 1966.
pivmin, nsplit, isplit,m, w, werr, wl, wu, iblock, indexw,work, iwork, info )
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: order, range
integer(${ik}$), intent(in) :: il, iu, n, nsplit
integer(${ik}$), intent(out) :: info, m
real(sp), intent(in) :: pivmin, reltol, vl, vu
real(sp), intent(out) :: wl, wu
! Array Arguments
integer(${ik}$), intent(out) :: iblock(*), indexw(*), iwork(*)
integer(${ik}$), intent(in) :: isplit(*)
real(sp), intent(in) :: d(*), e(*), e2(*), gers(*)
real(sp), intent(out) :: w(*), werr(*), work(*)
! =====================================================================
! Parameters
real(sp), parameter :: fudge = two
integer(${ik}$), parameter :: allrng = 1_${ik}$
integer(${ik}$), parameter :: valrng = 2_${ik}$
integer(${ik}$), parameter :: indrng = 3_${ik}$
! Local Scalars
logical(lk) :: ncnvrg, toofew
integer(${ik}$) :: i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iout, &
irange, itmax, itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb, nwl, nwu
real(sp) :: atoli, eps, gl, gu, rtoli, tmp1, tmp2, tnorm, uflow, wkill, wlu, &
wul
! Local Arrays
integer(${ik}$) :: idumma(1_${ik}$)
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
! decode range
if( stdlib_lsame( range, 'A' ) ) then
irange = allrng
else if( stdlib_lsame( range, 'V' ) ) then
irange = valrng
else if( stdlib_lsame( range, 'I' ) ) then
irange = indrng
else
irange = 0_${ik}$
end if
! check for errors
if( irange<=0_${ik}$ ) then
info = -1_${ik}$
else if( .not.(stdlib_lsame(order,'B').or.stdlib_lsame(order,'E')) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( irange==valrng ) then
if( vl>=vu )info = -5_${ik}$
else if( irange==indrng .and.( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) ) then
info = -6_${ik}$
else if( irange==indrng .and.( iu<min( n, il ) .or. iu>n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
return
end if
! initialize error flags
info = 0_${ik}$
ncnvrg = .false.
toofew = .false.
! quick return if possible
m = 0_${ik}$
if( n==0 ) return
! simplification:
if( irange==indrng .and. il==1_${ik}$ .and. iu==n ) irange = 1_${ik}$
! get machine constants
eps = stdlib${ii}$_slamch( 'P' )
uflow = stdlib${ii}$_slamch( 'U' )
! special case when n=1
! treat case of 1x1 matrix for quick return
if( n==1_${ik}$ ) then
if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
m = 1_${ik}$
w(1_${ik}$) = d(1_${ik}$)
! the computation error of the eigenvalue is zero
werr(1_${ik}$) = zero
iblock( 1_${ik}$ ) = 1_${ik}$
indexw( 1_${ik}$ ) = 1_${ik}$
endif
return
end if
! nb is the minimum vector length for vector bisection, or 0
! if only scalar is to be done.
nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
if( nb<=1_${ik}$ ) nb = 0_${ik}$
! find global spectral radius
gl = d(1_${ik}$)
gu = d(1_${ik}$)
do i = 1,n
gl = min( gl, gers( 2_${ik}$*i - 1_${ik}$))
gu = max( gu, gers(2_${ik}$*i) )
end do
! compute global gerschgorin bounds and spectral diameter
tnorm = max( abs( gl ), abs( gu ) )
gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
! [jan/28/2009] remove the line below since spdiam variable not use
! spdiam = gu - gl
! input arguments for stdlib${ii}$_slaebz:
! the relative tolerance. an interval (a,b] lies within
! "relative tolerance" if b-a < reltol*max(|a|,|b|),
rtoli = reltol
! set the absolute tolerance for interval convergence to zero to force
! interval convergence based on relative size of the interval.
! this is dangerous because intervals might not converge when reltol is
! small. but at least a very small number should be selected so that for
! strongly graded matrices, the code can get relatively accurate
! eigenvalues.
atoli = fudge*two*uflow + fudge*two*pivmin
if( irange==indrng ) then
! range='i': compute an interval containing eigenvalues
! il through iu. the initial interval [gl,gu] from the global
! gerschgorin bounds gl and gu is refined by stdlib${ii}$_slaebz.
itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
work( n+1 ) = gl
work( n+2 ) = gl
work( n+3 ) = gu
work( n+4 ) = gu
work( n+5 ) = gl
work( n+6 ) = gu
iwork( 1_${ik}$ ) = -1_${ik}$
iwork( 2_${ik}$ ) = -1_${ik}$
iwork( 3_${ik}$ ) = n + 1_${ik}$
iwork( 4_${ik}$ ) = n + 1_${ik}$
iwork( 5_${ik}$ ) = il - 1_${ik}$
iwork( 6_${ik}$ ) = iu
call stdlib${ii}$_slaebz( 3_${ik}$, itmax, n, 2_${ik}$, 2_${ik}$, nb, atoli, rtoli, pivmin,d, e, e2, iwork( 5_${ik}$ )&
, work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
! on exit, output intervals may not be ordered by ascending negcount
if( iwork( 6_${ik}$ )==iu ) then
wl = work( n+1 )
wlu = work( n+3 )
nwl = iwork( 1_${ik}$ )
wu = work( n+4 )
wul = work( n+2 )
nwu = iwork( 4_${ik}$ )
else
wl = work( n+2 )
wlu = work( n+4 )
nwl = iwork( 2_${ik}$ )
wu = work( n+3 )
wul = work( n+1 )
nwu = iwork( 3_${ik}$ )
end if
! on exit, the interval [wl, wlu] contains a value with negcount nwl,
! and [wul, wu] contains a value with negcount nwu.
if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then
info = 4_${ik}$
return
end if
elseif( irange==valrng ) then
wl = vl
wu = vu
elseif( irange==allrng ) then
wl = gl
wu = gu
endif
! find eigenvalues -- loop over blocks and recompute nwl and nwu.
! nwl accumulates the number of eigenvalues .le. wl,
! nwu accumulates the number of eigenvalues .le. wu
m = 0_${ik}$
iend = 0_${ik}$
info = 0_${ik}$
nwl = 0_${ik}$
nwu = 0_${ik}$
loop_70: do jblk = 1, nsplit
ioff = iend
ibegin = ioff + 1_${ik}$
iend = isplit( jblk )
in = iend - ioff
if( in==1_${ik}$ ) then
! 1x1 block
if( wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$
if( wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$
if( irange==allrng .or.( wl<d( ibegin )-pivmin.and. wu>= d( ibegin )-pivmin ) ) &
then
m = m + 1_${ik}$
w( m ) = d( ibegin )
werr(m) = zero
! the gap for a single block doesn't matter for the later
! algorithm and is assigned an arbitrary large value
iblock( m ) = jblk
indexw( m ) = 1_${ik}$
end if
! disabled 2x2 case because of a failure on the following matrix
! range = 'i', il = iu = 4
! original tridiagonal, d = [
! -0.150102010615740e+00_sp
! -0.849897989384260e+00_sp
! -0.128208148052635e-15_sp
! 0.128257718286320e-15_sp
! ];
! e = [
! -0.357171383266986e+00_sp
! -0.180411241501588e-15_sp
! -0.175152352710251e-15_sp
! ];
! else if( in==2 ) then
! * 2x2 block
! disc = sqrt( (half*(d(ibegin)-d(iend)))**2 + e(ibegin)**2 )
! tmp1 = half*(d(ibegin)+d(iend))
! l1 = tmp1 - disc
! if( wl>= l1-pivmin )
! $ nwl = nwl + 1
! if( wu>= l1-pivmin )
! $ nwu = nwu + 1
! if( irange==allrng .or. ( wl<l1-pivmin .and. wu>=
! $ l1-pivmin ) ) then
! m = m + 1
! w( m ) = l1
! * the uncertainty of eigenvalues of a 2x2 matrix is very small
! werr( m ) = eps * abs( w( m ) ) * two
! iblock( m ) = jblk
! indexw( m ) = 1
! endif
! l2 = tmp1 + disc
! if( wl>= l2-pivmin )
! $ nwl = nwl + 1
! if( wu>= l2-pivmin )
! $ nwu = nwu + 1
! if( irange==allrng .or. ( wl<l2-pivmin .and. wu>=
! $ l2-pivmin ) ) then
! m = m + 1
! w( m ) = l2
! * the uncertainty of eigenvalues of a 2x2 matrix is very small
! werr( m ) = eps * abs( w( m ) ) * two
! iblock( m ) = jblk
! indexw( m ) = 2
! endif
else
! general case - block of size in >= 2
! compute local gerschgorin interval and use it as the initial
! interval for stdlib${ii}$_slaebz
gu = d( ibegin )
gl = d( ibegin )
tmp1 = zero
do j = ibegin, iend
gl = min( gl, gers( 2_${ik}$*j - 1_${ik}$))
gu = max( gu, gers(2_${ik}$*j) )
end do
! [jan/28/2009]
! change spdiam by tnorm in lines 2 and 3 thereafter
! line 1: remove computation of spdiam (not useful anymore)
! spdiam = gu - gl
! gl = gl - fudge*spdiam*eps*in - fudge*pivmin
! gu = gu + fudge*spdiam*eps*in + fudge*pivmin
gl = gl - fudge*tnorm*eps*in - fudge*pivmin
gu = gu + fudge*tnorm*eps*in + fudge*pivmin
if( irange>1_${ik}$ ) then
if( gu<wl ) then
! the local block contains none of the wanted eigenvalues
nwl = nwl + in
nwu = nwu + in
cycle loop_70
end if
! refine search interval if possible, only range (wl,wu] matters
gl = max( gl, wl )
gu = min( gu, wu )
if( gl>=gu )cycle loop_70
end if
! find negcount of initial interval boundaries gl and gu
work( n+1 ) = gl
work( n+in+1 ) = gu
call stdlib${ii}$_slaebz( 1_${ik}$, 0_${ik}$, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e( &
ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 ),&
iblock( m+1 ), iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
nwl = nwl + iwork( 1_${ik}$ )
nwu = nwu + iwork( in+1 )
iwoff = m - iwork( 1_${ik}$ )
! compute eigenvalues
itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + &
2_${ik}$
call stdlib${ii}$_slaebz( 2_${ik}$, itmax, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e(&
ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( m+&
1_${ik}$ ), iblock( m+1 ), iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
! copy eigenvalues into w and iblock
! use -jblk for block number for unconverged eigenvalues.
! loop over the number of output intervals from stdlib${ii}$_slaebz
do j = 1, iout
! eigenvalue approximation is middle point of interval
tmp1 = half*( work( j+n )+work( j+in+n ) )
! semi length of error interval
tmp2 = half*abs( work( j+n )-work( j+in+n ) )
if( j>iout-iinfo ) then
! flag non-convergence.
ncnvrg = .true.
ib = -jblk
else
ib = jblk
end if
do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff
w( je ) = tmp1
werr( je ) = tmp2
indexw( je ) = je - iwoff
iblock( je ) = ib
end do
end do
m = m + im
end if
end do loop_70
! if range='i', then (wl,wu) contains eigenvalues nwl+1,...,nwu
! if nwl+1 < il or nwu > iu, discard extra eigenvalues.
if( irange==indrng ) then
idiscl = il - 1_${ik}$ - nwl
idiscu = nwu - iu
if( idiscl>0_${ik}$ ) then
im = 0_${ik}$
do je = 1, m
! remove some of the smallest eigenvalues from the left so that
! at the end idiscl =0. move all eigenvalues up to the left.
if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then
idiscl = idiscl - 1_${ik}$
else
im = im + 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
m = im
end if
if( idiscu>0_${ik}$ ) then
! remove some of the largest eigenvalues from the right so that
! at the end idiscu =0. move all eigenvalues up to the left.
im=m+1
do je = m, 1, -1
if( w( je )>=wul .and. idiscu>0_${ik}$ ) then
idiscu = idiscu - 1_${ik}$
else
im = im - 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
jee = 0_${ik}$
do je = im, m
jee = jee + 1_${ik}$
w( jee ) = w( je )
werr( jee ) = werr( je )
indexw( jee ) = indexw( je )
iblock( jee ) = iblock( je )
end do
m = m-im+1
end if
if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
! code to deal with effects of bad arithmetic. (if n(w) is
! monotone non-decreasing, this should never happen.)
! some low eigenvalues to be discarded are not in (wl,wlu],
! or high eigenvalues to be discarded are not in (wul,wu]
! so just kill off the smallest idiscl/largest idiscu
! eigenvalues, by marking the corresponding iblock = 0
if( idiscl>0_${ik}$ ) then
wkill = wu
do jdisc = 1, idiscl
iw = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ .and.( w( je )<wkill .or. iw==0_${ik}$ ) ) then
iw = je
wkill = w( je )
end if
end do
iblock( iw ) = 0_${ik}$
end do
end if
if( idiscu>0_${ik}$ ) then
wkill = wl
do jdisc = 1, idiscu
iw = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ .and.( w( je )>=wkill .or. iw==0_${ik}$ ) ) then
iw = je
wkill = w( je )
end if
end do
iblock( iw ) = 0_${ik}$
end do
end if
! now erase all eigenvalues with iblock set to zero
im = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ ) then
im = im + 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
m = im
end if
if( idiscl<0_${ik}$ .or. idiscu<0_${ik}$ ) then
toofew = .true.
end if
end if
if(( irange==allrng .and. m/=n ).or.( irange==indrng .and. m/=iu-il+1 ) ) then
toofew = .true.
end if
! if order='b', do nothing the eigenvalues are already sorted by
! block.
! if order='e', sort the eigenvalues from smallest to largest
if( stdlib_lsame(order,'E') .and. nsplit>1_${ik}$ ) then
do je = 1, m - 1
ie = 0_${ik}$
tmp1 = w( je )
do j = je + 1, m
if( w( j )<tmp1 ) then
ie = j
tmp1 = w( j )
end if
end do
if( ie/=0_${ik}$ ) then
tmp2 = werr( ie )
itmp1 = iblock( ie )
itmp2 = indexw( ie )
w( ie ) = w( je )
werr( ie ) = werr( je )
iblock( ie ) = iblock( je )
indexw( ie ) = indexw( je )
w( je ) = tmp1
werr( je ) = tmp2
iblock( je ) = itmp1
indexw( je ) = itmp2
end if
end do
end if
info = 0_${ik}$
if( ncnvrg )info = info + 1_${ik}$
if( toofew )info = info + 2_${ik}$
return
end subroutine stdlib${ii}$_slarrd
pure module subroutine stdlib${ii}$_dlarrd( range, order, n, vl, vu, il, iu, gers,reltol, d, e, e2, &
!! DLARRD computes the eigenvalues of a symmetric tridiagonal
!! matrix T to suitable accuracy. This is an auxiliary code to be
!! called from DSTEMR.
!! The user may ask for all eigenvalues, all eigenvalues
!! in the half-open interval (VL, VU], or the IL-th through IU-th
!! eigenvalues.
!! To avoid overflow, the matrix must be scaled so that its
!! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
!! accuracy, it should not be much smaller than that.
!! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
!! Matrix", Report CS41, Computer Science Dept., Stanford
!! University, July 21, 1966.
pivmin, nsplit, isplit,m, w, werr, wl, wu, iblock, indexw,work, iwork, info )
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: order, range
integer(${ik}$), intent(in) :: il, iu, n, nsplit
integer(${ik}$), intent(out) :: info, m
real(dp), intent(in) :: pivmin, reltol, vl, vu
real(dp), intent(out) :: wl, wu
! Array Arguments
integer(${ik}$), intent(out) :: iblock(*), indexw(*), iwork(*)
integer(${ik}$), intent(in) :: isplit(*)
real(dp), intent(in) :: d(*), e(*), e2(*), gers(*)
real(dp), intent(out) :: w(*), werr(*), work(*)
! =====================================================================
! Parameters
real(dp), parameter :: fudge = two
integer(${ik}$), parameter :: allrng = 1_${ik}$
integer(${ik}$), parameter :: valrng = 2_${ik}$
integer(${ik}$), parameter :: indrng = 3_${ik}$
! Local Scalars
logical(lk) :: ncnvrg, toofew
integer(${ik}$) :: i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iout, &
irange, itmax, itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb, nwl, nwu
real(dp) :: atoli, eps, gl, gu, rtoli, tmp1, tmp2, tnorm, uflow, wkill, wlu, &
wul
! Local Arrays
integer(${ik}$) :: idumma(1_${ik}$)
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
! decode range
if( stdlib_lsame( range, 'A' ) ) then
irange = allrng
else if( stdlib_lsame( range, 'V' ) ) then
irange = valrng
else if( stdlib_lsame( range, 'I' ) ) then
irange = indrng
else
irange = 0_${ik}$
end if
! check for errors
if( irange<=0_${ik}$ ) then
info = -1_${ik}$
else if( .not.(stdlib_lsame(order,'B').or.stdlib_lsame(order,'E')) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( irange==valrng ) then
if( vl>=vu )info = -5_${ik}$
else if( irange==indrng .and.( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) ) then
info = -6_${ik}$
else if( irange==indrng .and.( iu<min( n, il ) .or. iu>n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
return
end if
! initialize error flags
info = 0_${ik}$
ncnvrg = .false.
toofew = .false.
! quick return if possible
m = 0_${ik}$
if( n==0 ) return
! simplification:
if( irange==indrng .and. il==1_${ik}$ .and. iu==n ) irange = 1_${ik}$
! get machine constants
eps = stdlib${ii}$_dlamch( 'P' )
uflow = stdlib${ii}$_dlamch( 'U' )
! special case when n=1
! treat case of 1x1 matrix for quick return
if( n==1_${ik}$ ) then
if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
m = 1_${ik}$
w(1_${ik}$) = d(1_${ik}$)
! the computation error of the eigenvalue is zero
werr(1_${ik}$) = zero
iblock( 1_${ik}$ ) = 1_${ik}$
indexw( 1_${ik}$ ) = 1_${ik}$
endif
return
end if
! nb is the minimum vector length for vector bisection, or 0
! if only scalar is to be done.
nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
if( nb<=1_${ik}$ ) nb = 0_${ik}$
! find global spectral radius
gl = d(1_${ik}$)
gu = d(1_${ik}$)
do i = 1,n
gl = min( gl, gers( 2_${ik}$*i - 1_${ik}$))
gu = max( gu, gers(2_${ik}$*i) )
end do
! compute global gerschgorin bounds and spectral diameter
tnorm = max( abs( gl ), abs( gu ) )
gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
! [jan/28/2009] remove the line below since spdiam variable not use
! spdiam = gu - gl
! input arguments for stdlib${ii}$_dlaebz:
! the relative tolerance. an interval (a,b] lies within
! "relative tolerance" if b-a < reltol*max(|a|,|b|),
rtoli = reltol
! set the absolute tolerance for interval convergence to zero to force
! interval convergence based on relative size of the interval.
! this is dangerous because intervals might not converge when reltol is
! small. but at least a very small number should be selected so that for
! strongly graded matrices, the code can get relatively accurate
! eigenvalues.
atoli = fudge*two*uflow + fudge*two*pivmin
if( irange==indrng ) then
! range='i': compute an interval containing eigenvalues
! il through iu. the initial interval [gl,gu] from the global
! gerschgorin bounds gl and gu is refined by stdlib${ii}$_dlaebz.
itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
work( n+1 ) = gl
work( n+2 ) = gl
work( n+3 ) = gu
work( n+4 ) = gu
work( n+5 ) = gl
work( n+6 ) = gu
iwork( 1_${ik}$ ) = -1_${ik}$
iwork( 2_${ik}$ ) = -1_${ik}$
iwork( 3_${ik}$ ) = n + 1_${ik}$
iwork( 4_${ik}$ ) = n + 1_${ik}$
iwork( 5_${ik}$ ) = il - 1_${ik}$
iwork( 6_${ik}$ ) = iu
call stdlib${ii}$_dlaebz( 3_${ik}$, itmax, n, 2_${ik}$, 2_${ik}$, nb, atoli, rtoli, pivmin,d, e, e2, iwork( 5_${ik}$ )&
, work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
! on exit, output intervals may not be ordered by ascending negcount
if( iwork( 6_${ik}$ )==iu ) then
wl = work( n+1 )
wlu = work( n+3 )
nwl = iwork( 1_${ik}$ )
wu = work( n+4 )
wul = work( n+2 )
nwu = iwork( 4_${ik}$ )
else
wl = work( n+2 )
wlu = work( n+4 )
nwl = iwork( 2_${ik}$ )
wu = work( n+3 )
wul = work( n+1 )
nwu = iwork( 3_${ik}$ )
end if
! on exit, the interval [wl, wlu] contains a value with negcount nwl,
! and [wul, wu] contains a value with negcount nwu.
if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then
info = 4_${ik}$
return
end if
elseif( irange==valrng ) then
wl = vl
wu = vu
elseif( irange==allrng ) then
wl = gl
wu = gu
endif
! find eigenvalues -- loop over blocks and recompute nwl and nwu.
! nwl accumulates the number of eigenvalues .le. wl,
! nwu accumulates the number of eigenvalues .le. wu
m = 0_${ik}$
iend = 0_${ik}$
info = 0_${ik}$
nwl = 0_${ik}$
nwu = 0_${ik}$
loop_70: do jblk = 1, nsplit
ioff = iend
ibegin = ioff + 1_${ik}$
iend = isplit( jblk )
in = iend - ioff
if( in==1_${ik}$ ) then
! 1x1 block
if( wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$
if( wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$
if( irange==allrng .or.( wl<d( ibegin )-pivmin.and. wu>= d( ibegin )-pivmin ) ) &
then
m = m + 1_${ik}$
w( m ) = d( ibegin )
werr(m) = zero
! the gap for a single block doesn't matter for the later
! algorithm and is assigned an arbitrary large value
iblock( m ) = jblk
indexw( m ) = 1_${ik}$
end if
! disabled 2x2 case because of a failure on the following matrix
! range = 'i', il = iu = 4
! original tridiagonal, d = [
! -0.150102010615740e+00_dp
! -0.849897989384260e+00_dp
! -0.128208148052635e-15_dp
! 0.128257718286320e-15_dp
! ];
! e = [
! -0.357171383266986e+00_dp
! -0.180411241501588e-15_dp
! -0.175152352710251e-15_dp
! ];
! else if( in==2 ) then
! * 2x2 block
! disc = sqrt( (half*(d(ibegin)-d(iend)))**2 + e(ibegin)**2 )
! tmp1 = half*(d(ibegin)+d(iend))
! l1 = tmp1 - disc
! if( wl>= l1-pivmin )
! $ nwl = nwl + 1
! if( wu>= l1-pivmin )
! $ nwu = nwu + 1
! if( irange==allrng .or. ( wl<l1-pivmin .and. wu>=
! $ l1-pivmin ) ) then
! m = m + 1
! w( m ) = l1
! * the uncertainty of eigenvalues of a 2x2 matrix is very small
! werr( m ) = eps * abs( w( m ) ) * two
! iblock( m ) = jblk
! indexw( m ) = 1
! endif
! l2 = tmp1 + disc
! if( wl>= l2-pivmin )
! $ nwl = nwl + 1
! if( wu>= l2-pivmin )
! $ nwu = nwu + 1
! if( irange==allrng .or. ( wl<l2-pivmin .and. wu>=
! $ l2-pivmin ) ) then
! m = m + 1
! w( m ) = l2
! * the uncertainty of eigenvalues of a 2x2 matrix is very small
! werr( m ) = eps * abs( w( m ) ) * two
! iblock( m ) = jblk
! indexw( m ) = 2
! endif
else
! general case - block of size in >= 2
! compute local gerschgorin interval and use it as the initial
! interval for stdlib${ii}$_dlaebz
gu = d( ibegin )
gl = d( ibegin )
tmp1 = zero
do j = ibegin, iend
gl = min( gl, gers( 2_${ik}$*j - 1_${ik}$))
gu = max( gu, gers(2_${ik}$*j) )
end do
! [jan/28/2009]
! change spdiam by tnorm in lines 2 and 3 thereafter
! line 1: remove computation of spdiam (not useful anymore)
! spdiam = gu - gl
! gl = gl - fudge*spdiam*eps*in - fudge*pivmin
! gu = gu + fudge*spdiam*eps*in + fudge*pivmin
gl = gl - fudge*tnorm*eps*in - fudge*pivmin
gu = gu + fudge*tnorm*eps*in + fudge*pivmin
if( irange>1_${ik}$ ) then
if( gu<wl ) then
! the local block contains none of the wanted eigenvalues
nwl = nwl + in
nwu = nwu + in
cycle loop_70
end if
! refine search interval if possible, only range (wl,wu] matters
gl = max( gl, wl )
gu = min( gu, wu )
if( gl>=gu )cycle loop_70
end if
! find negcount of initial interval boundaries gl and gu
work( n+1 ) = gl
work( n+in+1 ) = gu
call stdlib${ii}$_dlaebz( 1_${ik}$, 0_${ik}$, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e( &
ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 ),&
iblock( m+1 ), iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
nwl = nwl + iwork( 1_${ik}$ )
nwu = nwu + iwork( in+1 )
iwoff = m - iwork( 1_${ik}$ )
! compute eigenvalues
itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + &
2_${ik}$
call stdlib${ii}$_dlaebz( 2_${ik}$, itmax, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e(&
ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( m+&
1_${ik}$ ), iblock( m+1 ), iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
! copy eigenvalues into w and iblock
! use -jblk for block number for unconverged eigenvalues.
! loop over the number of output intervals from stdlib${ii}$_dlaebz
do j = 1, iout
! eigenvalue approximation is middle point of interval
tmp1 = half*( work( j+n )+work( j+in+n ) )
! semi length of error interval
tmp2 = half*abs( work( j+n )-work( j+in+n ) )
if( j>iout-iinfo ) then
! flag non-convergence.
ncnvrg = .true.
ib = -jblk
else
ib = jblk
end if
do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff
w( je ) = tmp1
werr( je ) = tmp2
indexw( je ) = je - iwoff
iblock( je ) = ib
end do
end do
m = m + im
end if
end do loop_70
! if range='i', then (wl,wu) contains eigenvalues nwl+1,...,nwu
! if nwl+1 < il or nwu > iu, discard extra eigenvalues.
if( irange==indrng ) then
idiscl = il - 1_${ik}$ - nwl
idiscu = nwu - iu
if( idiscl>0_${ik}$ ) then
im = 0_${ik}$
do je = 1, m
! remove some of the smallest eigenvalues from the left so that
! at the end idiscl =0. move all eigenvalues up to the left.
if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then
idiscl = idiscl - 1_${ik}$
else
im = im + 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
m = im
end if
if( idiscu>0_${ik}$ ) then
! remove some of the largest eigenvalues from the right so that
! at the end idiscu =0. move all eigenvalues up to the left.
im=m+1
do je = m, 1, -1
if( w( je )>=wul .and. idiscu>0_${ik}$ ) then
idiscu = idiscu - 1_${ik}$
else
im = im - 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
jee = 0_${ik}$
do je = im, m
jee = jee + 1_${ik}$
w( jee ) = w( je )
werr( jee ) = werr( je )
indexw( jee ) = indexw( je )
iblock( jee ) = iblock( je )
end do
m = m-im+1
end if
if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
! code to deal with effects of bad arithmetic. (if n(w) is
! monotone non-decreasing, this should never happen.)
! some low eigenvalues to be discarded are not in (wl,wlu],
! or high eigenvalues to be discarded are not in (wul,wu]
! so just kill off the smallest idiscl/largest idiscu
! eigenvalues, by marking the corresponding iblock = 0
if( idiscl>0_${ik}$ ) then
wkill = wu
do jdisc = 1, idiscl
iw = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ .and.( w( je )<wkill .or. iw==0_${ik}$ ) ) then
iw = je
wkill = w( je )
end if
end do
iblock( iw ) = 0_${ik}$
end do
end if
if( idiscu>0_${ik}$ ) then
wkill = wl
do jdisc = 1, idiscu
iw = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ .and.( w( je )>=wkill .or. iw==0_${ik}$ ) ) then
iw = je
wkill = w( je )
end if
end do
iblock( iw ) = 0_${ik}$
end do
end if
! now erase all eigenvalues with iblock set to zero
im = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ ) then
im = im + 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
m = im
end if
if( idiscl<0_${ik}$ .or. idiscu<0_${ik}$ ) then
toofew = .true.
end if
end if
if(( irange==allrng .and. m/=n ).or.( irange==indrng .and. m/=iu-il+1 ) ) then
toofew = .true.
end if
! if order='b', do nothing the eigenvalues are already sorted by
! block.
! if order='e', sort the eigenvalues from smallest to largest
if( stdlib_lsame(order,'E') .and. nsplit>1_${ik}$ ) then
do je = 1, m - 1
ie = 0_${ik}$
tmp1 = w( je )
do j = je + 1, m
if( w( j )<tmp1 ) then
ie = j
tmp1 = w( j )
end if
end do
if( ie/=0_${ik}$ ) then
tmp2 = werr( ie )
itmp1 = iblock( ie )
itmp2 = indexw( ie )
w( ie ) = w( je )
werr( ie ) = werr( je )
iblock( ie ) = iblock( je )
indexw( ie ) = indexw( je )
w( je ) = tmp1
werr( je ) = tmp2
iblock( je ) = itmp1
indexw( je ) = itmp2
end if
end do
end if
info = 0_${ik}$
if( ncnvrg )info = info + 1_${ik}$
if( toofew )info = info + 2_${ik}$
return
end subroutine stdlib${ii}$_dlarrd
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$larrd( range, order, n, vl, vu, il, iu, gers,reltol, d, e, e2, &
!! DLARRD: computes the eigenvalues of a symmetric tridiagonal
!! matrix T to suitable accuracy. This is an auxiliary code to be
!! called from DSTEMR.
!! The user may ask for all eigenvalues, all eigenvalues
!! in the half-open interval (VL, VU], or the IL-th through IU-th
!! eigenvalues.
!! To avoid overflow, the matrix must be scaled so that its
!! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
!! accuracy, it should not be much smaller than that.
!! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
!! Matrix", Report CS41, Computer Science Dept., Stanford
!! University, July 21, 1966.
pivmin, nsplit, isplit,m, w, werr, wl, wu, iblock, indexw,work, iwork, info )
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: order, range
integer(${ik}$), intent(in) :: il, iu, n, nsplit
integer(${ik}$), intent(out) :: info, m
real(${rk}$), intent(in) :: pivmin, reltol, vl, vu
real(${rk}$), intent(out) :: wl, wu
! Array Arguments
integer(${ik}$), intent(out) :: iblock(*), indexw(*), iwork(*)
integer(${ik}$), intent(in) :: isplit(*)
real(${rk}$), intent(in) :: d(*), e(*), e2(*), gers(*)
real(${rk}$), intent(out) :: w(*), werr(*), work(*)
! =====================================================================
! Parameters
real(${rk}$), parameter :: fudge = two
integer(${ik}$), parameter :: allrng = 1_${ik}$
integer(${ik}$), parameter :: valrng = 2_${ik}$
integer(${ik}$), parameter :: indrng = 3_${ik}$
! Local Scalars
logical(lk) :: ncnvrg, toofew
integer(${ik}$) :: i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iout, &
irange, itmax, itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb, nwl, nwu
real(${rk}$) :: atoli, eps, gl, gu, rtoli, tmp1, tmp2, tnorm, uflow, wkill, wlu, &
wul
! Local Arrays
integer(${ik}$) :: idumma(1_${ik}$)
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
! decode range
if( stdlib_lsame( range, 'A' ) ) then
irange = allrng
else if( stdlib_lsame( range, 'V' ) ) then
irange = valrng
else if( stdlib_lsame( range, 'I' ) ) then
irange = indrng
else
irange = 0_${ik}$
end if
! check for errors
if( irange<=0_${ik}$ ) then
info = -1_${ik}$
else if( .not.(stdlib_lsame(order,'B').or.stdlib_lsame(order,'E')) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( irange==valrng ) then
if( vl>=vu )info = -5_${ik}$
else if( irange==indrng .and.( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) ) then
info = -6_${ik}$
else if( irange==indrng .and.( iu<min( n, il ) .or. iu>n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
return
end if
! initialize error flags
info = 0_${ik}$
ncnvrg = .false.
toofew = .false.
! quick return if possible
m = 0_${ik}$
if( n==0 ) return
! simplification:
if( irange==indrng .and. il==1_${ik}$ .and. iu==n ) irange = 1_${ik}$
! get machine constants
eps = stdlib${ii}$_${ri}$lamch( 'P' )
uflow = stdlib${ii}$_${ri}$lamch( 'U' )
! special case when n=1
! treat case of 1x1 matrix for quick return
if( n==1_${ik}$ ) then
if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
m = 1_${ik}$
w(1_${ik}$) = d(1_${ik}$)
! the computation error of the eigenvalue is zero
werr(1_${ik}$) = zero
iblock( 1_${ik}$ ) = 1_${ik}$
indexw( 1_${ik}$ ) = 1_${ik}$
endif
return
end if
! nb is the minimum vector length for vector bisection, or 0
! if only scalar is to be done.
nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
if( nb<=1_${ik}$ ) nb = 0_${ik}$
! find global spectral radius
gl = d(1_${ik}$)
gu = d(1_${ik}$)
do i = 1,n
gl = min( gl, gers( 2_${ik}$*i - 1_${ik}$))
gu = max( gu, gers(2_${ik}$*i) )
end do
! compute global gerschgorin bounds and spectral diameter
tnorm = max( abs( gl ), abs( gu ) )
gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
! [jan/28/2009] remove the line below since spdiam variable not use
! spdiam = gu - gl
! input arguments for stdlib${ii}$_${ri}$laebz:
! the relative tolerance. an interval (a,b] lies within
! "relative tolerance" if b-a < reltol*max(|a|,|b|),
rtoli = reltol
! set the absolute tolerance for interval convergence to zero to force
! interval convergence based on relative size of the interval.
! this is dangerous because intervals might not converge when reltol is
! small. but at least a very small number should be selected so that for
! strongly graded matrices, the code can get relatively accurate
! eigenvalues.
atoli = fudge*two*uflow + fudge*two*pivmin
if( irange==indrng ) then
! range='i': compute an interval containing eigenvalues
! il through iu. the initial interval [gl,gu] from the global
! gerschgorin bounds gl and gu is refined by stdlib${ii}$_${ri}$laebz.
itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$
work( n+1 ) = gl
work( n+2 ) = gl
work( n+3 ) = gu
work( n+4 ) = gu
work( n+5 ) = gl
work( n+6 ) = gu
iwork( 1_${ik}$ ) = -1_${ik}$
iwork( 2_${ik}$ ) = -1_${ik}$
iwork( 3_${ik}$ ) = n + 1_${ik}$
iwork( 4_${ik}$ ) = n + 1_${ik}$
iwork( 5_${ik}$ ) = il - 1_${ik}$
iwork( 6_${ik}$ ) = iu
call stdlib${ii}$_${ri}$laebz( 3_${ik}$, itmax, n, 2_${ik}$, 2_${ik}$, nb, atoli, rtoli, pivmin,d, e, e2, iwork( 5_${ik}$ )&
, work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
! on exit, output intervals may not be ordered by ascending negcount
if( iwork( 6_${ik}$ )==iu ) then
wl = work( n+1 )
wlu = work( n+3 )
nwl = iwork( 1_${ik}$ )
wu = work( n+4 )
wul = work( n+2 )
nwu = iwork( 4_${ik}$ )
else
wl = work( n+2 )
wlu = work( n+4 )
nwl = iwork( 2_${ik}$ )
wu = work( n+3 )
wul = work( n+1 )
nwu = iwork( 3_${ik}$ )
end if
! on exit, the interval [wl, wlu] contains a value with negcount nwl,
! and [wul, wu] contains a value with negcount nwu.
if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then
info = 4_${ik}$
return
end if
elseif( irange==valrng ) then
wl = vl
wu = vu
elseif( irange==allrng ) then
wl = gl
wu = gu
endif
! find eigenvalues -- loop over blocks and recompute nwl and nwu.
! nwl accumulates the number of eigenvalues .le. wl,
! nwu accumulates the number of eigenvalues .le. wu
m = 0_${ik}$
iend = 0_${ik}$
info = 0_${ik}$
nwl = 0_${ik}$
nwu = 0_${ik}$
loop_70: do jblk = 1, nsplit
ioff = iend
ibegin = ioff + 1_${ik}$
iend = isplit( jblk )
in = iend - ioff
if( in==1_${ik}$ ) then
! 1x1 block
if( wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$
if( wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$
if( irange==allrng .or.( wl<d( ibegin )-pivmin.and. wu>= d( ibegin )-pivmin ) ) &
then
m = m + 1_${ik}$
w( m ) = d( ibegin )
werr(m) = zero
! the gap for a single block doesn't matter for the later
! algorithm and is assigned an arbitrary large value
iblock( m ) = jblk
indexw( m ) = 1_${ik}$
end if
! disabled 2x2 case because of a failure on the following matrix
! range = 'i', il = iu = 4
! original tridiagonal, d = [
! -0.150102010615740e+00_${rk}$
! -0.849897989384260e+00_${rk}$
! -0.128208148052635e-15_${rk}$
! 0.128257718286320e-15_${rk}$
! ];
! e = [
! -0.357171383266986e+00_${rk}$
! -0.180411241501588e-15_${rk}$
! -0.175152352710251e-15_${rk}$
! ];
! else if( in==2 ) then
! * 2x2 block
! disc = sqrt( (half*(d(ibegin)-d(iend)))**2 + e(ibegin)**2 )
! tmp1 = half*(d(ibegin)+d(iend))
! l1 = tmp1 - disc
! if( wl>= l1-pivmin )
! $ nwl = nwl + 1
! if( wu>= l1-pivmin )
! $ nwu = nwu + 1
! if( irange==allrng .or. ( wl<l1-pivmin .and. wu>=
! $ l1-pivmin ) ) then
! m = m + 1
! w( m ) = l1
! * the uncertainty of eigenvalues of a 2x2 matrix is very small
! werr( m ) = eps * abs( w( m ) ) * two
! iblock( m ) = jblk
! indexw( m ) = 1
! endif
! l2 = tmp1 + disc
! if( wl>= l2-pivmin )
! $ nwl = nwl + 1
! if( wu>= l2-pivmin )
! $ nwu = nwu + 1
! if( irange==allrng .or. ( wl<l2-pivmin .and. wu>=
! $ l2-pivmin ) ) then
! m = m + 1
! w( m ) = l2
! * the uncertainty of eigenvalues of a 2x2 matrix is very small
! werr( m ) = eps * abs( w( m ) ) * two
! iblock( m ) = jblk
! indexw( m ) = 2
! endif
else
! general case - block of size in >= 2
! compute local gerschgorin interval and use it as the initial
! interval for stdlib${ii}$_${ri}$laebz
gu = d( ibegin )
gl = d( ibegin )
tmp1 = zero
do j = ibegin, iend
gl = min( gl, gers( 2_${ik}$*j - 1_${ik}$))
gu = max( gu, gers(2_${ik}$*j) )
end do
! [jan/28/2009]
! change spdiam by tnorm in lines 2 and 3 thereafter
! line 1: remove computation of spdiam (not useful anymore)
! spdiam = gu - gl
! gl = gl - fudge*spdiam*eps*in - fudge*pivmin
! gu = gu + fudge*spdiam*eps*in + fudge*pivmin
gl = gl - fudge*tnorm*eps*in - fudge*pivmin
gu = gu + fudge*tnorm*eps*in + fudge*pivmin
if( irange>1_${ik}$ ) then
if( gu<wl ) then
! the local block contains none of the wanted eigenvalues
nwl = nwl + in
nwu = nwu + in
cycle loop_70
end if
! refine search interval if possible, only range (wl,wu] matters
gl = max( gl, wl )
gu = min( gu, wu )
if( gl>=gu )cycle loop_70
end if
! find negcount of initial interval boundaries gl and gu
work( n+1 ) = gl
work( n+in+1 ) = gu
call stdlib${ii}$_${ri}$laebz( 1_${ik}$, 0_${ik}$, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e( &
ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 ),&
iblock( m+1 ), iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
nwl = nwl + iwork( 1_${ik}$ )
nwu = nwu + iwork( in+1 )
iwoff = m - iwork( 1_${ik}$ )
! compute eigenvalues
itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + &
2_${ik}$
call stdlib${ii}$_${ri}$laebz( 2_${ik}$, itmax, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e(&
ibegin ), e2( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( m+&
1_${ik}$ ), iblock( m+1 ), iinfo )
if( iinfo /= 0_${ik}$ ) then
info = iinfo
return
end if
! copy eigenvalues into w and iblock
! use -jblk for block number for unconverged eigenvalues.
! loop over the number of output intervals from stdlib${ii}$_${ri}$laebz
do j = 1, iout
! eigenvalue approximation is middle point of interval
tmp1 = half*( work( j+n )+work( j+in+n ) )
! semi length of error interval
tmp2 = half*abs( work( j+n )-work( j+in+n ) )
if( j>iout-iinfo ) then
! flag non-convergence.
ncnvrg = .true.
ib = -jblk
else
ib = jblk
end if
do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff
w( je ) = tmp1
werr( je ) = tmp2
indexw( je ) = je - iwoff
iblock( je ) = ib
end do
end do
m = m + im
end if
end do loop_70
! if range='i', then (wl,wu) contains eigenvalues nwl+1,...,nwu
! if nwl+1 < il or nwu > iu, discard extra eigenvalues.
if( irange==indrng ) then
idiscl = il - 1_${ik}$ - nwl
idiscu = nwu - iu
if( idiscl>0_${ik}$ ) then
im = 0_${ik}$
do je = 1, m
! remove some of the smallest eigenvalues from the left so that
! at the end idiscl =0. move all eigenvalues up to the left.
if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then
idiscl = idiscl - 1_${ik}$
else
im = im + 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
m = im
end if
if( idiscu>0_${ik}$ ) then
! remove some of the largest eigenvalues from the right so that
! at the end idiscu =0. move all eigenvalues up to the left.
im=m+1
do je = m, 1, -1
if( w( je )>=wul .and. idiscu>0_${ik}$ ) then
idiscu = idiscu - 1_${ik}$
else
im = im - 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
jee = 0_${ik}$
do je = im, m
jee = jee + 1_${ik}$
w( jee ) = w( je )
werr( jee ) = werr( je )
indexw( jee ) = indexw( je )
iblock( jee ) = iblock( je )
end do
m = m-im+1
end if
if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then
! code to deal with effects of bad arithmetic. (if n(w) is
! monotone non-decreasing, this should never happen.)
! some low eigenvalues to be discarded are not in (wl,wlu],
! or high eigenvalues to be discarded are not in (wul,wu]
! so just kill off the smallest idiscl/largest idiscu
! eigenvalues, by marking the corresponding iblock = 0
if( idiscl>0_${ik}$ ) then
wkill = wu
do jdisc = 1, idiscl
iw = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ .and.( w( je )<wkill .or. iw==0_${ik}$ ) ) then
iw = je
wkill = w( je )
end if
end do
iblock( iw ) = 0_${ik}$
end do
end if
if( idiscu>0_${ik}$ ) then
wkill = wl
do jdisc = 1, idiscu
iw = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ .and.( w( je )>=wkill .or. iw==0_${ik}$ ) ) then
iw = je
wkill = w( je )
end if
end do
iblock( iw ) = 0_${ik}$
end do
end if
! now erase all eigenvalues with iblock set to zero
im = 0_${ik}$
do je = 1, m
if( iblock( je )/=0_${ik}$ ) then
im = im + 1_${ik}$
w( im ) = w( je )
werr( im ) = werr( je )
indexw( im ) = indexw( je )
iblock( im ) = iblock( je )
end if
end do
m = im
end if
if( idiscl<0_${ik}$ .or. idiscu<0_${ik}$ ) then
toofew = .true.
end if
end if
if(( irange==allrng .and. m/=n ).or.( irange==indrng .and. m/=iu-il+1 ) ) then
toofew = .true.
end if
! if order='b', do nothing the eigenvalues are already sorted by
! block.
! if order='e', sort the eigenvalues from smallest to largest
if( stdlib_lsame(order,'E') .and. nsplit>1_${ik}$ ) then
do je = 1, m - 1
ie = 0_${ik}$
tmp1 = w( je )
do j = je + 1, m
if( w( j )<tmp1 ) then
ie = j
tmp1 = w( j )
end if
end do
if( ie/=0_${ik}$ ) then
tmp2 = werr( ie )
itmp1 = iblock( ie )
itmp2 = indexw( ie )
w( ie ) = w( je )
werr( ie ) = werr( je )
iblock( ie ) = iblock( je )
indexw( ie ) = indexw( je )
w( je ) = tmp1
werr( je ) = tmp2
iblock( je ) = itmp1
indexw( je ) = itmp2
end if
end do
end if
info = 0_${ik}$
if( ncnvrg )info = info + 1_${ik}$
if( toofew )info = info + 2_${ik}$
return
end subroutine stdlib${ii}$_${ri}$larrd
#:endif
#:endfor
pure module subroutine stdlib${ii}$_slarre( range, n, vl, vu, il, iu, d, e, e2,rtol1, rtol2, spltol, &
!! To find the desired eigenvalues of a given real symmetric
!! tridiagonal matrix T, SLARRE: sets any "small" off-diagonal
!! elements to zero, and for each unreduced block T_i, it finds
!! (a) a suitable shift at one end of the block's spectrum,
!! (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
!! (c) eigenvalues of each L_i D_i L_i^T.
!! The representations and eigenvalues found are then used by
!! SSTEMR to compute the eigenvectors of T.
!! The accuracy varies depending on whether bisection is used to
!! find a few eigenvalues or the dqds algorithm (subroutine SLASQ2) to
!! conpute all and then discard any unwanted one.
!! As an added benefit, SLARRE also outputs the n
!! Gerschgorin intervals for the matrices L_i D_i L_i^T.
nsplit, isplit, m,w, werr, wgap, iblock, indexw, gers, pivmin,work, iwork, info )
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: range
integer(${ik}$), intent(in) :: il, iu, n
integer(${ik}$), intent(out) :: info, m, nsplit
real(sp), intent(out) :: pivmin
real(sp), intent(in) :: rtol1, rtol2, spltol
real(sp), intent(inout) :: vl, vu
! Array Arguments
integer(${ik}$), intent(out) :: iblock(*), isplit(*), iwork(*), indexw(*)
real(sp), intent(inout) :: d(*), e(*), e2(*)
real(sp), intent(out) :: gers(*), w(*), werr(*), wgap(*), work(*)
! =====================================================================
! Parameters
real(sp), parameter :: hndrd = 100.0_sp
real(sp), parameter :: pert = 4.0_sp
real(sp), parameter :: fourth = one/four
real(sp), parameter :: fac = half
real(sp), parameter :: maxgrowth = 64.0_sp
real(sp), parameter :: fudge = two
integer(${ik}$), parameter :: maxtry = 6_${ik}$
integer(${ik}$), parameter :: allrng = 1_${ik}$
integer(${ik}$), parameter :: indrng = 2_${ik}$
integer(${ik}$), parameter :: valrng = 3_${ik}$
! Local Scalars
logical(lk) :: forceb, norep, usedqd
integer(${ik}$) :: cnt, cnt1, cnt2, i, ibegin, idum, iend, iinfo, in, indl, indu, irange, &
j, jblk, mb, mm, wbegin, wend
real(sp) :: avgap, bsrtol, clwdth, dmax, dpivot, eabs, emax, eold, eps, gl, gu, isleft,&
isrght, rtl, rtol, s1, s2, safmin, sgndef, sigma, spdiam, tau, tmp, tmp1
! Local Arrays
integer(${ik}$) :: iseed(4_${ik}$)
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
! quick return if possible
if( n<=0_${ik}$ ) then
return
end if
! decode range
if( stdlib_lsame( range, 'A' ) ) then
irange = allrng
else if( stdlib_lsame( range, 'V' ) ) then
irange = valrng
else if( stdlib_lsame( range, 'I' ) ) then
irange = indrng
end if
m = 0_${ik}$
! get machine constants
safmin = stdlib${ii}$_slamch( 'S' )
eps = stdlib${ii}$_slamch( 'P' )
! set parameters
rtl = hndrd*eps
! if one were ever to ask for less initial precision in bsrtol,
! one should keep in mind that for the subset case, the extremal
! eigenvalues must be at least as accurate as the current setting
! (eigenvalues in the middle need not as much accuracy)
bsrtol = sqrt(eps)*(0.5e-3_sp)
! treat case of 1x1 matrix for quick return
if( n==1_${ik}$ ) then
if( (irange==allrng).or.((irange==valrng).and.(d(1_${ik}$)>vl).and.(d(1_${ik}$)<=vu)).or.((&
irange==indrng).and.(il==1_${ik}$).and.(iu==1_${ik}$)) ) then
m = 1_${ik}$
w(1_${ik}$) = d(1_${ik}$)
! the computation error of the eigenvalue is zero
werr(1_${ik}$) = zero
wgap(1_${ik}$) = zero
iblock( 1_${ik}$ ) = 1_${ik}$
indexw( 1_${ik}$ ) = 1_${ik}$
gers(1_${ik}$) = d( 1_${ik}$ )
gers(2_${ik}$) = d( 1_${ik}$ )
endif
! store the shift for the initial rrr, which is zero in this case
e(1_${ik}$) = zero
return
end if
! general case: tridiagonal matrix of order > 1
! init werr, wgap. compute gerschgorin intervals and spectral diameter.
! compute maximum off-diagonal entry and pivmin.
gl = d(1_${ik}$)
gu = d(1_${ik}$)
eold = zero
emax = zero
e(n) = zero
do i = 1,n
werr(i) = zero
wgap(i) = zero
eabs = abs( e(i) )
if( eabs >= emax ) then
emax = eabs
end if
tmp1 = eabs + eold
gers( 2_${ik}$*i-1) = d(i) - tmp1
gl = min( gl, gers( 2_${ik}$*i - 1_${ik}$))
gers( 2_${ik}$*i ) = d(i) + tmp1
gu = max( gu, gers(2_${ik}$*i) )
eold = eabs
end do
! the minimum pivot allowed in the sturm sequence for t
pivmin = safmin * max( one, emax**2_${ik}$ )
! compute spectral diameter. the gerschgorin bounds give an
! estimate that is wrong by at most a factor of sqrt(2)
spdiam = gu - gl
! compute splitting points
call stdlib${ii}$_slarra( n, d, e, e2, spltol, spdiam,nsplit, isplit, iinfo )
! can force use of bisection instead of faster dqds.
! option left in the code for future multisection work.
forceb = .false.
! initialize usedqd, dqds should be used for allrng unless someone
! explicitly wants bisection.
usedqd = (( irange==allrng ) .and. (.not.forceb))
if( (irange==allrng) .and. (.not. forceb) ) then
! set interval [vl,vu] that contains all eigenvalues
vl = gl
vu = gu
else
! we call stdlib${ii}$_slarrd to find crude approximations to the eigenvalues
! in the desired range. in case irange = indrng, we also obtain the
! interval (vl,vu] that contains all the wanted eigenvalues.
! an interval [left,right] has converged if
! right-left<rtol*max(abs(left),abs(right))
! stdlib${ii}$_slarrd needs a work of size 4*n, iwork of size 3*n
call stdlib${ii}$_slarrd( range, 'B', n, vl, vu, il, iu, gers,bsrtol, d, e, e2, pivmin, &
nsplit, isplit,mm, w, werr, vl, vu, iblock, indexw,work, iwork, iinfo )
if( iinfo/=0_${ik}$ ) then
info = -1_${ik}$
return
endif
! make sure that the entries m+1 to n in w, werr, iblock, indexw are 0
do i = mm+1,n
w( i ) = zero
werr( i ) = zero
iblock( i ) = 0_${ik}$
indexw( i ) = 0_${ik}$
end do
end if
! **
! loop over unreduced blocks
ibegin = 1_${ik}$
wbegin = 1_${ik}$
loop_170: do jblk = 1, nsplit
iend = isplit( jblk )
in = iend - ibegin + 1