#: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