#:include "common.fypp" submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_tridiag3 implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_sstev( jobz, n, d, e, z, ldz, work, info ) !! SSTEV computes all eigenvalues and, optionally, eigenvectors of a !! real symmetric tridiagonal matrix A. ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(sp), intent(inout) :: d(*), e(*) real(sp), intent(out) :: work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: wantz integer(${ik}$) :: imax, iscale real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSTEV ', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = sqrt( bignum ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ tnrm = stdlib${ii}$_slanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) end if ! for eigenvalues only, call stdlib${ii}$_ssterf. for eigenvalues and ! eigenvectors, call stdlib${ii}$_ssteqr. if( .not.wantz ) then call stdlib${ii}$_ssterf( n, d, e, info ) else call stdlib${ii}$_ssteqr( 'I', n, d, e, z, ldz, work, info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = n else imax = info - 1_${ik}$ end if call stdlib${ii}$_sscal( imax, one / sigma, d, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_sstev pure module subroutine stdlib${ii}$_dstev( jobz, n, d, e, z, ldz, work, info ) !! DSTEV computes all eigenvalues and, optionally, eigenvectors of a !! real symmetric tridiagonal matrix A. ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(dp), intent(inout) :: d(*), e(*) real(dp), intent(out) :: work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: wantz integer(${ik}$) :: imax, iscale real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEV ', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_dlamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = sqrt( bignum ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ tnrm = stdlib${ii}$_dlanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) end if ! for eigenvalues only, call stdlib${ii}$_dsterf. for eigenvalues and ! eigenvectors, call stdlib${ii}$_dsteqr. if( .not.wantz ) then call stdlib${ii}$_dsterf( n, d, e, info ) else call stdlib${ii}$_dsteqr( 'I', n, d, e, z, ldz, work, info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = n else imax = info - 1_${ik}$ end if call stdlib${ii}$_dscal( imax, one / sigma, d, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_dstev #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$stev( jobz, n, d, e, z, ldz, work, info ) !! DSTEV: computes all eigenvalues and, optionally, eigenvectors of a !! real symmetric tridiagonal matrix A. ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(${rk}$), intent(inout) :: d(*), e(*) real(${rk}$), intent(out) :: work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: wantz integer(${ik}$) :: imax, iscale real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEV ', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = sqrt( bignum ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) end if ! for eigenvalues only, call stdlib${ii}$_${ri}$sterf. for eigenvalues and ! eigenvectors, call stdlib${ii}$_${ri}$steqr. if( .not.wantz ) then call stdlib${ii}$_${ri}$sterf( n, d, e, info ) else call stdlib${ii}$_${ri}$steqr( 'I', n, d, e, z, ldz, work, info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = n else imax = info - 1_${ik}$ end if call stdlib${ii}$_${ri}$scal( imax, one / sigma, d, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_${ri}$stev #:endif #:endfor pure module subroutine stdlib${ii}$_sstevd( jobz, n, d, e, z, ldz, work, lwork, iwork,liwork, info ) !! SSTEVD computes all eigenvalues and, optionally, eigenvectors of a !! real symmetric tridiagonal matrix. If eigenvectors are desired, it !! uses a divide and conquer algorithm. !! The divide and conquer algorithm makes very mild assumptions about !! floating point arithmetic. It will work on machines with a guard !! digit in add/subtract, or on those binary machines without guard !! digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or !! Cray-2. It could conceivably fail on hexadecimal or decimal machines !! without guard digits, but we know of none. ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, liwork, lwork, n ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(inout) :: d(*), e(*) real(sp), intent(out) :: work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery, wantz integer(${ik}$) :: iscale, liwmin, lwmin real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) info = 0_${ik}$ liwmin = 1_${ik}$ lwmin = 1_${ik}$ if( n>1_${ik}$ .and. wantz ) then lwmin = 1_${ik}$ + 4_${ik}$*n + n**2_${ik}$ liwmin = 3_${ik}$ + 5_${ik}$*n end if if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -6_${ik}$ end if if( info==0_${ik}$ ) then work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -8_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -10_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSTEVD', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = sqrt( bignum ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ tnrm = stdlib${ii}$_slanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) end if ! for eigenvalues only, call stdlib${ii}$_ssterf. for eigenvalues and ! eigenvectors, call stdlib${ii}$_sstedc. if( .not.wantz ) then call stdlib${ii}$_ssterf( n, d, e, info ) else call stdlib${ii}$_sstedc( 'I', n, d, e, z, ldz, work, lwork, iwork, liwork,info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. if( iscale==1_${ik}$ )call stdlib${ii}$_sscal( n, one / sigma, d, 1_${ik}$ ) work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_sstevd pure module subroutine stdlib${ii}$_dstevd( jobz, n, d, e, z, ldz, work, lwork, iwork,liwork, info ) !! DSTEVD computes all eigenvalues and, optionally, eigenvectors of a !! real symmetric tridiagonal matrix. If eigenvectors are desired, it !! uses a divide and conquer algorithm. !! The divide and conquer algorithm makes very mild assumptions about !! floating point arithmetic. It will work on machines with a guard !! digit in add/subtract, or on those binary machines without guard !! digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or !! Cray-2. It could conceivably fail on hexadecimal or decimal machines !! without guard digits, but we know of none. ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, liwork, lwork, n ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(inout) :: d(*), e(*) real(dp), intent(out) :: work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery, wantz integer(${ik}$) :: iscale, liwmin, lwmin real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) info = 0_${ik}$ liwmin = 1_${ik}$ lwmin = 1_${ik}$ if( n>1_${ik}$ .and. wantz ) then lwmin = 1_${ik}$ + 4_${ik}$*n + n**2_${ik}$ liwmin = 3_${ik}$ + 5_${ik}$*n end if if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -6_${ik}$ end if if( info==0_${ik}$ ) then work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -8_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -10_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEVD', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_dlamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = sqrt( bignum ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ tnrm = stdlib${ii}$_dlanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) end if ! for eigenvalues only, call stdlib${ii}$_dsterf. for eigenvalues and ! eigenvectors, call stdlib${ii}$_dstedc. if( .not.wantz ) then call stdlib${ii}$_dsterf( n, d, e, info ) else call stdlib${ii}$_dstedc( 'I', n, d, e, z, ldz, work, lwork, iwork, liwork,info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. if( iscale==1_${ik}$ )call stdlib${ii}$_dscal( n, one / sigma, d, 1_${ik}$ ) work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_dstevd #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$stevd( jobz, n, d, e, z, ldz, work, lwork, iwork,liwork, info ) !! DSTEVD: computes all eigenvalues and, optionally, eigenvectors of a !! real symmetric tridiagonal matrix. If eigenvectors are desired, it !! uses a divide and conquer algorithm. !! The divide and conquer algorithm makes very mild assumptions about !! floating point arithmetic. It will work on machines with a guard !! digit in add/subtract, or on those binary machines without guard !! digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or !! Cray-2. It could conceivably fail on hexadecimal or decimal machines !! without guard digits, but we know of none. ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, liwork, lwork, n ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(${rk}$), intent(inout) :: d(*), e(*) real(${rk}$), intent(out) :: work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery, wantz integer(${ik}$) :: iscale, liwmin, lwmin real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tnrm ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) info = 0_${ik}$ liwmin = 1_${ik}$ lwmin = 1_${ik}$ if( n>1_${ik}$ .and. wantz ) then lwmin = 1_${ik}$ + 4_${ik}$*n + n**2_${ik}$ liwmin = 3_${ik}$ + 5_${ik}$*n end if if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -6_${ik}$ end if if( info==0_${ik}$ ) then work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -8_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -10_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEVD', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = sqrt( bignum ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) end if ! for eigenvalues only, call stdlib${ii}$_${ri}$sterf. for eigenvalues and ! eigenvectors, call stdlib${ii}$_${ri}$stedc. if( .not.wantz ) then call stdlib${ii}$_${ri}$sterf( n, d, e, info ) else call stdlib${ii}$_${ri}$stedc( 'I', n, d, e, z, ldz, work, lwork, iwork, liwork,info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. if( iscale==1_${ik}$ )call stdlib${ii}$_${ri}$scal( n, one / sigma, d, 1_${ik}$ ) work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_${ri}$stevd #:endif #:endfor pure module subroutine stdlib${ii}$_sstevr( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, & !! SSTEVR computes selected eigenvalues and, optionally, eigenvectors !! of a real symmetric tridiagonal matrix T. Eigenvalues and !! eigenvectors can be selected by specifying either a range of values !! or a range of indices for the desired eigenvalues. !! Whenever possible, SSTEVR calls SSTEMR to compute the !! eigenspectrum using Relatively Robust Representations. SSTEMR !! computes eigenvalues by the dqds algorithm, while orthogonal !! eigenvectors are computed from various "good" L D L^T representations !! (also known as Relatively Robust Representations). Gram-Schmidt !! orthogonalization is avoided as far as possible. More specifically, !! the various steps of the algorithm are as follows. For the i-th !! unreduced block of T, !! (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T !! is a relatively robust representation, !! (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high !! relative accuracy by the dqds algorithm, !! (c) If there is a cluster of close eigenvalues, "choose" sigma_i !! close to the cluster, and go to step (a), !! (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, !! compute the corresponding eigenvector by forming a !! rank-revealing twisted factorization. !! The desired accuracy of the output can be specified by the input !! parameter ABSTOL. !! For more details, see "A new O(n^2) algorithm for the symmetric !! tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, !! Computer Science Division Technical Report No. UCB//CSD-97-971, !! UC Berkeley, May 1997. !! Note 1 : SSTEVR calls SSTEMR when the full spectrum is requested !! on machines which conform to the ieee-754 floating point standard. !! SSTEVR calls SSTEBZ and SSTEIN on non-ieee machines and !! when partial spectrum requests are made. !! Normal execution of SSTEMR may create NaNs and infinities and !! hence may abort due to a floating point exception in environments !! which do not handle NaNs and infinities in the ieee standard default !! manner. isuppz, work, lwork, iwork,liwork, info ) ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz, range integer(${ik}$), intent(in) :: il, iu, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(sp), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: isuppz(*), iwork(*) real(sp), intent(inout) :: d(*), e(*) real(sp), intent(out) :: w(*), work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: alleig, indeig, test, lquery, valeig, wantz, tryrac character :: order integer(${ik}$) :: i, ieeeok, imax, indibl, indifl, indisp, indiwo, iscale, j, jj, liwmin,& lwmin, nsplit real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, & vuu ! Intrinsic Functions ! Executable Statements ! test the input parameters. ieeeok = stdlib${ii}$_ilaenv( 10_${ik}$, 'SSTEVR', 'N', 1_${ik}$, 2_${ik}$, 3_${ik}$, 4_${ik}$ ) wantz = stdlib_lsame( jobz, 'V' ) alleig = stdlib_lsame( range, 'A' ) valeig = stdlib_lsame( range, 'V' ) indeig = stdlib_lsame( range, 'I' ) lquery = ( ( lwork==-1_${ik}$ ) .or. ( liwork==-1_${ik}$ ) ) lwmin = max( 1_${ik}$, 20_${ik}$*n ) liwmin = max(1_${ik}$, 10_${ik}$*n ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( .not.( alleig .or. valeig .or. indeig ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( valeig ) then if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$ else if( indeig ) then if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( iu<min( n, il ) .or. iu>n ) then info = -9_${ik}$ end if end if end if if( info==0_${ik}$ ) then if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -14_${ik}$ end if end if if( info==0_${ik}$ ) then work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -17_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -19_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSTEVR', -info ) return else if( lquery ) then return end if ! quick return if possible m = 0_${ik}$ if( n==0 )return if( n==1_${ik}$ ) then if( alleig .or. indeig ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) else if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) end if end if if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ if( valeig ) then vll = vl vuu = vu end if tnrm = stdlib${ii}$_slanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) if( valeig ) then vll = vl*sigma vuu = vu*sigma end if end if ! initialize indices into workspaces. note: these indices are used only ! if stdlib${ii}$_ssterf or stdlib${ii}$_sstemr fail. ! iwork(indibl:indibl+m-1) corresponds to iblock in stdlib${ii}$_sstebz and ! stores the block indices of each of the m<=n eigenvalues. indibl = 1_${ik}$ ! iwork(indisp:indisp+nsplit-1) corresponds to isplit in stdlib${ii}$_sstebz and ! stores the starting and finishing indices of each block. indisp = indibl + n ! iwork(indifl:indifl+n-1) stores the indices of eigenvectors ! that corresponding to eigenvectors that fail to converge in ! stdlib${ii}$_sstein. this information is discarded; if any fail, the driver ! returns info > 0. indifl = indisp + n ! indiwo is the offset of the remaining integer workspace. indiwo = indisp + n ! if all eigenvalues are desired, then ! call stdlib${ii}$_ssterf or stdlib${ii}$_sstemr. if this fails for some eigenvalue, then ! try stdlib${ii}$_sstebz. test = .false. if( indeig ) then if( il==1_${ik}$ .and. iu==n ) then test = .true. end if end if if( ( alleig .or. test ) .and. ieeeok==1_${ik}$ ) then call stdlib${ii}$_scopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ ) if( .not.wantz ) then call stdlib${ii}$_scopy( n, d, 1_${ik}$, w, 1_${ik}$ ) call stdlib${ii}$_ssterf( n, w, work, info ) else call stdlib${ii}$_scopy( n, d, 1_${ik}$, work( n+1 ), 1_${ik}$ ) if (abstol <= two*n*eps) then tryrac = .true. else tryrac = .false. end if call stdlib${ii}$_sstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,iu, m, w, z, ldz,& n, isuppz, tryrac,work( 2_${ik}$*n+1 ), lwork-2*n, iwork, liwork, info ) end if if( info==0_${ik}$ ) then m = n go to 10 end if info = 0_${ik}$ end if ! otherwise, call stdlib${ii}$_sstebz and, if eigenvectors are desired, stdlib${ii}$_sstein. if( wantz ) then order = 'B' else order = 'E' end if call stdlib${ii}$_sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, & iwork( indibl ), iwork( indisp ), work,iwork( indiwo ), info ) if( wantz ) then call stdlib${ii}$_sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work, & iwork( indiwo ), iwork( indifl ),info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. 10 continue if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = m else imax = info - 1_${ik}$ end if call stdlib${ii}$_sscal( imax, one / sigma, w, 1_${ik}$ ) end if ! if eigenvalues are not in order, then sort them, along with ! eigenvectors. if( wantz ) then do j = 1, m - 1 i = 0_${ik}$ tmp1 = w( j ) do jj = j + 1, m if( w( jj )<tmp1 ) then i = jj tmp1 = w( jj ) end if end do if( i/=0_${ik}$ ) then w( i ) = w( j ) w( j ) = tmp1 call stdlib${ii}$_sswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ ) end if end do end if ! causes problems with tests 19 ! if (wantz .and. indeig ) z( 1,1) = z(1,1) / 1.002_sp + .002 work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_sstevr pure module subroutine stdlib${ii}$_dstevr( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, & !! DSTEVR computes selected eigenvalues and, optionally, eigenvectors !! of a real symmetric tridiagonal matrix T. Eigenvalues and !! eigenvectors can be selected by specifying either a range of values !! or a range of indices for the desired eigenvalues. !! Whenever possible, DSTEVR calls DSTEMR to compute the !! eigenspectrum using Relatively Robust Representations. DSTEMR !! computes eigenvalues by the dqds algorithm, while orthogonal !! eigenvectors are computed from various "good" L D L^T representations !! (also known as Relatively Robust Representations). Gram-Schmidt !! orthogonalization is avoided as far as possible. More specifically, !! the various steps of the algorithm are as follows. For the i-th !! unreduced block of T, !! (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T !! is a relatively robust representation, !! (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high !! relative accuracy by the dqds algorithm, !! (c) If there is a cluster of close eigenvalues, "choose" sigma_i !! close to the cluster, and go to step (a), !! (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, !! compute the corresponding eigenvector by forming a !! rank-revealing twisted factorization. !! The desired accuracy of the output can be specified by the input !! parameter ABSTOL. !! For more details, see "A new O(n^2) algorithm for the symmetric !! tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, !! Computer Science Division Technical Report No. UCB//CSD-97-971, !! UC Berkeley, May 1997. !! Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested !! on machines which conform to the ieee-754 floating point standard. !! DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and !! when partial spectrum requests are made. !! Normal execution of DSTEMR may create NaNs and infinities and !! hence may abort due to a floating point exception in environments !! which do not handle NaNs and infinities in the ieee standard default !! manner. isuppz, work, lwork, iwork,liwork, info ) ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz, range integer(${ik}$), intent(in) :: il, iu, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(dp), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: isuppz(*), iwork(*) real(dp), intent(inout) :: d(*), e(*) real(dp), intent(out) :: w(*), work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: alleig, indeig, test, lquery, valeig, wantz, tryrac character :: order integer(${ik}$) :: i, ieeeok, imax, indibl, indifl, indisp, indiwo, iscale, itmp1, j, jj, & liwmin, lwmin, nsplit real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, & vuu ! Intrinsic Functions ! Executable Statements ! test the input parameters. ieeeok = stdlib${ii}$_ilaenv( 10_${ik}$, 'DSTEVR', 'N', 1_${ik}$, 2_${ik}$, 3_${ik}$, 4_${ik}$ ) wantz = stdlib_lsame( jobz, 'V' ) alleig = stdlib_lsame( range, 'A' ) valeig = stdlib_lsame( range, 'V' ) indeig = stdlib_lsame( range, 'I' ) lquery = ( ( lwork==-1_${ik}$ ) .or. ( liwork==-1_${ik}$ ) ) lwmin = max( 1_${ik}$, 20_${ik}$*n ) liwmin = max( 1_${ik}$, 10_${ik}$*n ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( .not.( alleig .or. valeig .or. indeig ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( valeig ) then if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$ else if( indeig ) then if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( iu<min( n, il ) .or. iu>n ) then info = -9_${ik}$ end if end if end if if( info==0_${ik}$ ) then if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -14_${ik}$ end if end if if( info==0_${ik}$ ) then work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -17_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -19_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEVR', -info ) return else if( lquery ) then return end if ! quick return if possible m = 0_${ik}$ if( n==0 )return if( n==1_${ik}$ ) then if( alleig .or. indeig ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) else if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) end if end if if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_dlamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ if( valeig ) then vll = vl vuu = vu end if tnrm = stdlib${ii}$_dlanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) if( valeig ) then vll = vl*sigma vuu = vu*sigma end if end if ! initialize indices into workspaces. note: these indices are used only ! if stdlib${ii}$_dsterf or stdlib${ii}$_dstemr fail. ! iwork(indibl:indibl+m-1) corresponds to iblock in stdlib${ii}$_dstebz and ! stores the block indices of each of the m<=n eigenvalues. indibl = 1_${ik}$ ! iwork(indisp:indisp+nsplit-1) corresponds to isplit in stdlib${ii}$_dstebz and ! stores the starting and finishing indices of each block. indisp = indibl + n ! iwork(indifl:indifl+n-1) stores the indices of eigenvectors ! that corresponding to eigenvectors that fail to converge in ! stdlib${ii}$_dstein. this information is discarded; if any fail, the driver ! returns info > 0. indifl = indisp + n ! indiwo is the offset of the remaining integer workspace. indiwo = indisp + n ! if all eigenvalues are desired, then ! call stdlib${ii}$_dsterf or stdlib${ii}$_dstemr. if this fails for some eigenvalue, then ! try stdlib${ii}$_dstebz. test = .false. if( indeig ) then if( il==1_${ik}$ .and. iu==n ) then test = .true. end if end if if( ( alleig .or. test ) .and. ieeeok==1_${ik}$ ) then call stdlib${ii}$_dcopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ ) if( .not.wantz ) then call stdlib${ii}$_dcopy( n, d, 1_${ik}$, w, 1_${ik}$ ) call stdlib${ii}$_dsterf( n, w, work, info ) else call stdlib${ii}$_dcopy( n, d, 1_${ik}$, work( n+1 ), 1_${ik}$ ) if (abstol <= two*n*eps) then tryrac = .true. else tryrac = .false. end if call stdlib${ii}$_dstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,iu, m, w, z, ldz,& n, isuppz, tryrac,work( 2_${ik}$*n+1 ), lwork-2*n, iwork, liwork, info ) end if if( info==0_${ik}$ ) then m = n go to 10 end if info = 0_${ik}$ end if ! otherwise, call stdlib${ii}$_dstebz and, if eigenvectors are desired, stdlib${ii}$_dstein. if( wantz ) then order = 'B' else order = 'E' end if call stdlib${ii}$_dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, & iwork( indibl ), iwork( indisp ), work,iwork( indiwo ), info ) if( wantz ) then call stdlib${ii}$_dstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work, & iwork( indiwo ), iwork( indifl ),info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. 10 continue if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = m else imax = info - 1_${ik}$ end if call stdlib${ii}$_dscal( imax, one / sigma, w, 1_${ik}$ ) end if ! if eigenvalues are not in order, then sort them, along with ! eigenvectors. if( wantz ) then do j = 1, m - 1 i = 0_${ik}$ tmp1 = w( j ) do jj = j + 1, m if( w( jj )<tmp1 ) then i = jj tmp1 = w( jj ) end if end do if( i/=0_${ik}$ ) then itmp1 = iwork( i ) w( i ) = w( j ) iwork( i ) = iwork( j ) w( j ) = tmp1 iwork( j ) = itmp1 call stdlib${ii}$_dswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ ) end if end do end if ! causes problems with tests 19 ! if (wantz .and. indeig ) z( 1,1) = z(1,1) / 1.002_dp + .002 work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_dstevr #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$stevr( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, & !! DSTEVR: computes selected eigenvalues and, optionally, eigenvectors !! of a real symmetric tridiagonal matrix T. Eigenvalues and !! eigenvectors can be selected by specifying either a range of values !! or a range of indices for the desired eigenvalues. !! Whenever possible, DSTEVR calls DSTEMR to compute the !! eigenspectrum using Relatively Robust Representations. DSTEMR !! computes eigenvalues by the dqds algorithm, while orthogonal !! eigenvectors are computed from various "good" L D L^T representations !! (also known as Relatively Robust Representations). Gram-Schmidt !! orthogonalization is avoided as far as possible. More specifically, !! the various steps of the algorithm are as follows. For the i-th !! unreduced block of T, !! (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T !! is a relatively robust representation, !! (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high !! relative accuracy by the dqds algorithm, !! (c) If there is a cluster of close eigenvalues, "choose" sigma_i !! close to the cluster, and go to step (a), !! (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, !! compute the corresponding eigenvector by forming a !! rank-revealing twisted factorization. !! The desired accuracy of the output can be specified by the input !! parameter ABSTOL. !! For more details, see "A new O(n^2) algorithm for the symmetric !! tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, !! Computer Science Division Technical Report No. UCB//CSD-97-971, !! UC Berkeley, May 1997. !! Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested !! on machines which conform to the ieee-754 floating point standard. !! DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and !! when partial spectrum requests are made. !! Normal execution of DSTEMR may create NaNs and infinities and !! hence may abort due to a floating point exception in environments !! which do not handle NaNs and infinities in the ieee standard default !! manner. isuppz, work, lwork, iwork,liwork, info ) ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz, range integer(${ik}$), intent(in) :: il, iu, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(${rk}$), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: isuppz(*), iwork(*) real(${rk}$), intent(inout) :: d(*), e(*) real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: alleig, indeig, test, lquery, valeig, wantz, tryrac character :: order integer(${ik}$) :: i, ieeeok, imax, indibl, indifl, indisp, indiwo, iscale, itmp1, j, jj, & liwmin, lwmin, nsplit real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, & vuu ! Intrinsic Functions ! Executable Statements ! test the input parameters. ieeeok = stdlib${ii}$_ilaenv( 10_${ik}$, 'DSTEVR', 'N', 1_${ik}$, 2_${ik}$, 3_${ik}$, 4_${ik}$ ) wantz = stdlib_lsame( jobz, 'V' ) alleig = stdlib_lsame( range, 'A' ) valeig = stdlib_lsame( range, 'V' ) indeig = stdlib_lsame( range, 'I' ) lquery = ( ( lwork==-1_${ik}$ ) .or. ( liwork==-1_${ik}$ ) ) lwmin = max( 1_${ik}$, 20_${ik}$*n ) liwmin = max( 1_${ik}$, 10_${ik}$*n ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( .not.( alleig .or. valeig .or. indeig ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( valeig ) then if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$ else if( indeig ) then if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( iu<min( n, il ) .or. iu>n ) then info = -9_${ik}$ end if end if end if if( info==0_${ik}$ ) then if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -14_${ik}$ end if end if if( info==0_${ik}$ ) then work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -17_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -19_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEVR', -info ) return else if( lquery ) then return end if ! quick return if possible m = 0_${ik}$ if( n==0 )return if( n==1_${ik}$ ) then if( alleig .or. indeig ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) else if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) end if end if if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ if( valeig ) then vll = vl vuu = vu end if tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) if( valeig ) then vll = vl*sigma vuu = vu*sigma end if end if ! initialize indices into workspaces. note: these indices are used only ! if stdlib${ii}$_${ri}$sterf or stdlib${ii}$_${ri}$stemr fail. ! iwork(indibl:indibl+m-1) corresponds to iblock in stdlib${ii}$_${ri}$stebz and ! stores the block indices of each of the m<=n eigenvalues. indibl = 1_${ik}$ ! iwork(indisp:indisp+nsplit-1) corresponds to isplit in stdlib${ii}$_${ri}$stebz and ! stores the starting and finishing indices of each block. indisp = indibl + n ! iwork(indifl:indifl+n-1) stores the indices of eigenvectors ! that corresponding to eigenvectors that fail to converge in ! stdlib${ii}$_${ri}$stein. this information is discarded; if any fail, the driver ! returns info > 0. indifl = indisp + n ! indiwo is the offset of the remaining integer workspace. indiwo = indisp + n ! if all eigenvalues are desired, then ! call stdlib${ii}$_${ri}$sterf or stdlib${ii}$_${ri}$stemr. if this fails for some eigenvalue, then ! try stdlib${ii}$_${ri}$stebz. test = .false. if( indeig ) then if( il==1_${ik}$ .and. iu==n ) then test = .true. end if end if if( ( alleig .or. test ) .and. ieeeok==1_${ik}$ ) then call stdlib${ii}$_${ri}$copy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ ) if( .not.wantz ) then call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, w, 1_${ik}$ ) call stdlib${ii}$_${ri}$sterf( n, w, work, info ) else call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, work( n+1 ), 1_${ik}$ ) if (abstol <= two*n*eps) then tryrac = .true. else tryrac = .false. end if call stdlib${ii}$_${ri}$stemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,iu, m, w, z, ldz,& n, isuppz, tryrac,work( 2_${ik}$*n+1 ), lwork-2*n, iwork, liwork, info ) end if if( info==0_${ik}$ ) then m = n go to 10 end if info = 0_${ik}$ end if ! otherwise, call stdlib${ii}$_${ri}$stebz and, if eigenvectors are desired, stdlib${ii}$_${ri}$stein. if( wantz ) then order = 'B' else order = 'E' end if call stdlib${ii}$_${ri}$stebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, & iwork( indibl ), iwork( indisp ), work,iwork( indiwo ), info ) if( wantz ) then call stdlib${ii}$_${ri}$stein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work, & iwork( indiwo ), iwork( indifl ),info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. 10 continue if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = m else imax = info - 1_${ik}$ end if call stdlib${ii}$_${ri}$scal( imax, one / sigma, w, 1_${ik}$ ) end if ! if eigenvalues are not in order, then sort them, along with ! eigenvectors. if( wantz ) then do j = 1, m - 1 i = 0_${ik}$ tmp1 = w( j ) do jj = j + 1, m if( w( jj )<tmp1 ) then i = jj tmp1 = w( jj ) end if end do if( i/=0_${ik}$ ) then itmp1 = iwork( i ) w( i ) = w( j ) iwork( i ) = iwork( j ) w( j ) = tmp1 iwork( j ) = itmp1 call stdlib${ii}$_${ri}$swap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ ) end if end do end if ! causes problems with tests 19 ! if (wantz .and. indeig ) z( 1,1) = z(1,1) / 1.002_${rk}$ + .002 work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_${ri}$stevr #:endif #:endfor pure module subroutine stdlib${ii}$_sstevx( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, & !! SSTEVX computes selected eigenvalues and, optionally, eigenvectors !! of a real symmetric tridiagonal matrix A. Eigenvalues and !! eigenvectors can be selected by specifying either a range of values !! or a range of indices for the desired eigenvalues. work, iwork, ifail, info ) ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz, range integer(${ik}$), intent(in) :: il, iu, ldz, n integer(${ik}$), intent(out) :: info, m real(sp), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: ifail(*), iwork(*) real(sp), intent(inout) :: d(*), e(*) real(sp), intent(out) :: w(*), work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: alleig, indeig, test, valeig, wantz character :: order integer(${ik}$) :: i, imax, indibl, indisp, indiwo, indwrk, iscale, itmp1, j, jj, & nsplit real(sp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, & vuu ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) alleig = stdlib_lsame( range, 'A' ) valeig = stdlib_lsame( range, 'V' ) indeig = stdlib_lsame( range, 'I' ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( .not.( alleig .or. valeig .or. indeig ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( valeig ) then if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$ else if( indeig ) then if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( iu<min( n, il ) .or. iu>n ) then info = -9_${ik}$ end if end if end if if( info==0_${ik}$ ) then if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) )info = -14_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSTEVX', -info ) return end if ! quick return if possible m = 0_${ik}$ if( n==0 )return if( n==1_${ik}$ ) then if( alleig .or. indeig ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) else if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) end if end if if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ if ( valeig ) then vll = vl vuu = vu else vll = zero vuu = zero endif tnrm = stdlib${ii}$_slanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_sscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_sscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) if( valeig ) then vll = vl*sigma vuu = vu*sigma end if end if ! if all eigenvalues are desired and abstol is less than zero, then ! call stdlib${ii}$_ssterf or stdlib${ii}$_ssteqr. if this fails for some eigenvalue, then ! try stdlib${ii}$_sstebz. test = .false. if( indeig ) then if( il==1_${ik}$ .and. iu==n ) then test = .true. end if end if if( ( alleig .or. test ) .and. ( abstol<=zero ) ) then call stdlib${ii}$_scopy( n, d, 1_${ik}$, w, 1_${ik}$ ) call stdlib${ii}$_scopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ ) indwrk = n + 1_${ik}$ if( .not.wantz ) then call stdlib${ii}$_ssterf( n, w, work, info ) else call stdlib${ii}$_ssteqr( 'I', n, w, work, z, ldz, work( indwrk ), info ) if( info==0_${ik}$ ) then do i = 1, n ifail( i ) = 0_${ik}$ end do end if end if if( info==0_${ik}$ ) then m = n go to 20 end if info = 0_${ik}$ end if ! otherwise, call stdlib${ii}$_sstebz and, if eigenvectors are desired, stdlib${ii}$_sstein. if( wantz ) then order = 'B' else order = 'E' end if indwrk = 1_${ik}$ indibl = 1_${ik}$ indisp = indibl + n indiwo = indisp + n call stdlib${ii}$_sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, & iwork( indibl ), iwork( indisp ),work( indwrk ), iwork( indiwo ), info ) if( wantz ) then call stdlib${ii}$_sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work( & indwrk ), iwork( indiwo ), ifail,info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. 20 continue if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = m else imax = info - 1_${ik}$ end if call stdlib${ii}$_sscal( imax, one / sigma, w, 1_${ik}$ ) end if ! if eigenvalues are not in order, then sort them, along with ! eigenvectors. if( wantz ) then do j = 1, m - 1 i = 0_${ik}$ tmp1 = w( j ) do jj = j + 1, m if( w( jj )<tmp1 ) then i = jj tmp1 = w( jj ) end if end do if( i/=0_${ik}$ ) then itmp1 = iwork( indibl+i-1 ) w( i ) = w( j ) iwork( indibl+i-1 ) = iwork( indibl+j-1 ) w( j ) = tmp1 iwork( indibl+j-1 ) = itmp1 call stdlib${ii}$_sswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ ) if( info/=0_${ik}$ ) then itmp1 = ifail( i ) ifail( i ) = ifail( j ) ifail( j ) = itmp1 end if end if end do end if return end subroutine stdlib${ii}$_sstevx pure module subroutine stdlib${ii}$_dstevx( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, & !! DSTEVX computes selected eigenvalues and, optionally, eigenvectors !! of a real symmetric tridiagonal matrix A. Eigenvalues and !! eigenvectors can be selected by specifying either a range of values !! or a range of indices for the desired eigenvalues. work, iwork, ifail, info ) ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz, range integer(${ik}$), intent(in) :: il, iu, ldz, n integer(${ik}$), intent(out) :: info, m real(dp), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: ifail(*), iwork(*) real(dp), intent(inout) :: d(*), e(*) real(dp), intent(out) :: w(*), work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: alleig, indeig, test, valeig, wantz character :: order integer(${ik}$) :: i, imax, indibl, indisp, indiwo, indwrk, iscale, itmp1, j, jj, & nsplit real(dp) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, & vuu ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) alleig = stdlib_lsame( range, 'A' ) valeig = stdlib_lsame( range, 'V' ) indeig = stdlib_lsame( range, 'I' ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( .not.( alleig .or. valeig .or. indeig ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( valeig ) then if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$ else if( indeig ) then if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( iu<min( n, il ) .or. iu>n ) then info = -9_${ik}$ end if end if end if if( info==0_${ik}$ ) then if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) )info = -14_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEVX', -info ) return end if ! quick return if possible m = 0_${ik}$ if( n==0 )return if( n==1_${ik}$ ) then if( alleig .or. indeig ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) else if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) end if end if if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_dlamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ if( valeig ) then vll = vl vuu = vu else vll = zero vuu = zero end if tnrm = stdlib${ii}$_dlanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_dscal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_dscal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) if( valeig ) then vll = vl*sigma vuu = vu*sigma end if end if ! if all eigenvalues are desired and abstol is less than zero, then ! call stdlib${ii}$_dsterf or stdlib${ii}$_ssteqr. if this fails for some eigenvalue, then ! try stdlib${ii}$_dstebz. test = .false. if( indeig ) then if( il==1_${ik}$ .and. iu==n ) then test = .true. end if end if if( ( alleig .or. test ) .and. ( abstol<=zero ) ) then call stdlib${ii}$_dcopy( n, d, 1_${ik}$, w, 1_${ik}$ ) call stdlib${ii}$_dcopy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ ) indwrk = n + 1_${ik}$ if( .not.wantz ) then call stdlib${ii}$_dsterf( n, w, work, info ) else call stdlib${ii}$_dsteqr( 'I', n, w, work, z, ldz, work( indwrk ), info ) if( info==0_${ik}$ ) then do i = 1, n ifail( i ) = 0_${ik}$ end do end if end if if( info==0_${ik}$ ) then m = n go to 20 end if info = 0_${ik}$ end if ! otherwise, call stdlib${ii}$_dstebz and, if eigenvectors are desired, stdlib${ii}$_sstein. if( wantz ) then order = 'B' else order = 'E' end if indwrk = 1_${ik}$ indibl = 1_${ik}$ indisp = indibl + n indiwo = indisp + n call stdlib${ii}$_dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, & iwork( indibl ), iwork( indisp ),work( indwrk ), iwork( indiwo ), info ) if( wantz ) then call stdlib${ii}$_dstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work( & indwrk ), iwork( indiwo ), ifail,info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. 20 continue if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = m else imax = info - 1_${ik}$ end if call stdlib${ii}$_dscal( imax, one / sigma, w, 1_${ik}$ ) end if ! if eigenvalues are not in order, then sort them, along with ! eigenvectors. if( wantz ) then do j = 1, m - 1 i = 0_${ik}$ tmp1 = w( j ) do jj = j + 1, m if( w( jj )<tmp1 ) then i = jj tmp1 = w( jj ) end if end do if( i/=0_${ik}$ ) then itmp1 = iwork( indibl+i-1 ) w( i ) = w( j ) iwork( indibl+i-1 ) = iwork( indibl+j-1 ) w( j ) = tmp1 iwork( indibl+j-1 ) = itmp1 call stdlib${ii}$_dswap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ ) if( info/=0_${ik}$ ) then itmp1 = ifail( i ) ifail( i ) = ifail( j ) ifail( j ) = itmp1 end if end if end do end if return end subroutine stdlib${ii}$_dstevx #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$stevx( jobz, range, n, d, e, vl, vu, il, iu, abstol,m, w, z, ldz, & !! DSTEVX: computes selected eigenvalues and, optionally, eigenvectors !! of a real symmetric tridiagonal matrix A. Eigenvalues and !! eigenvectors can be selected by specifying either a range of values !! or a range of indices for the desired eigenvalues. work, iwork, ifail, info ) ! -- lapack driver routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: jobz, range integer(${ik}$), intent(in) :: il, iu, ldz, n integer(${ik}$), intent(out) :: info, m real(${rk}$), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: ifail(*), iwork(*) real(${rk}$), intent(inout) :: d(*), e(*) real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*) ! ===================================================================== ! Local Scalars logical(lk) :: alleig, indeig, test, valeig, wantz character :: order integer(${ik}$) :: i, imax, indibl, indisp, indiwo, indwrk, iscale, itmp1, j, jj, & nsplit real(${rk}$) :: bignum, eps, rmax, rmin, safmin, sigma, smlnum, tmp1, tnrm, vll, & vuu ! Intrinsic Functions ! Executable Statements ! test the input parameters. wantz = stdlib_lsame( jobz, 'V' ) alleig = stdlib_lsame( range, 'A' ) valeig = stdlib_lsame( range, 'V' ) indeig = stdlib_lsame( range, 'I' ) info = 0_${ik}$ if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then info = -1_${ik}$ else if( .not.( alleig .or. valeig .or. indeig ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( valeig ) then if( n>0_${ik}$ .and. vu<=vl )info = -7_${ik}$ else if( indeig ) then if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( iu<min( n, il ) .or. iu>n ) then info = -9_${ik}$ end if end if end if if( info==0_${ik}$ ) then if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) )info = -14_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEVX', -info ) return end if ! quick return if possible m = 0_${ik}$ if( n==0 )return if( n==1_${ik}$ ) then if( alleig .or. indeig ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) else if( vl<d( 1_${ik}$ ) .and. vu>=d( 1_${ik}$ ) ) then m = 1_${ik}$ w( 1_${ik}$ ) = d( 1_${ik}$ ) end if end if if( wantz )z( 1_${ik}$, 1_${ik}$ ) = one return end if ! get machine constants. safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) eps = stdlib${ii}$_${ri}$lamch( 'PRECISION' ) smlnum = safmin / eps bignum = one / smlnum rmin = sqrt( smlnum ) rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) ) ! scale matrix to allowable range, if necessary. iscale = 0_${ik}$ if( valeig ) then vll = vl vuu = vu else vll = zero vuu = zero end if tnrm = stdlib${ii}$_${ri}$lanst( 'M', n, d, e ) if( tnrm>zero .and. tnrm<rmin ) then iscale = 1_${ik}$ sigma = rmin / tnrm else if( tnrm>rmax ) then iscale = 1_${ik}$ sigma = rmax / tnrm end if if( iscale==1_${ik}$ ) then call stdlib${ii}$_${ri}$scal( n, sigma, d, 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( n-1, sigma, e( 1_${ik}$ ), 1_${ik}$ ) if( valeig ) then vll = vl*sigma vuu = vu*sigma end if end if ! if all eigenvalues are desired and abstol is less than zero, then ! call stdlib${ii}$_${ri}$sterf or stdlib${ii}$_dsteqr. if this fails for some eigenvalue, then ! try stdlib${ii}$_${ri}$stebz. test = .false. if( indeig ) then if( il==1_${ik}$ .and. iu==n ) then test = .true. end if end if if( ( alleig .or. test ) .and. ( abstol<=zero ) ) then call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, w, 1_${ik}$ ) call stdlib${ii}$_${ri}$copy( n-1, e( 1_${ik}$ ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ ) indwrk = n + 1_${ik}$ if( .not.wantz ) then call stdlib${ii}$_${ri}$sterf( n, w, work, info ) else call stdlib${ii}$_${ri}$steqr( 'I', n, w, work, z, ldz, work( indwrk ), info ) if( info==0_${ik}$ ) then do i = 1, n ifail( i ) = 0_${ik}$ end do end if end if if( info==0_${ik}$ ) then m = n go to 20 end if info = 0_${ik}$ end if ! otherwise, call stdlib${ii}$_${ri}$stebz and, if eigenvectors are desired, stdlib${ii}$_dstein. if( wantz ) then order = 'B' else order = 'E' end if indwrk = 1_${ik}$ indibl = 1_${ik}$ indisp = indibl + n indiwo = indisp + n call stdlib${ii}$_${ri}$stebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,nsplit, w, & iwork( indibl ), iwork( indisp ),work( indwrk ), iwork( indiwo ), info ) if( wantz ) then call stdlib${ii}$_${ri}$stein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),z, ldz, work( & indwrk ), iwork( indiwo ), ifail,info ) end if ! if matrix was scaled, then rescale eigenvalues appropriately. 20 continue if( iscale==1_${ik}$ ) then if( info==0_${ik}$ ) then imax = m else imax = info - 1_${ik}$ end if call stdlib${ii}$_${ri}$scal( imax, one / sigma, w, 1_${ik}$ ) end if ! if eigenvalues are not in order, then sort them, along with ! eigenvectors. if( wantz ) then do j = 1, m - 1 i = 0_${ik}$ tmp1 = w( j ) do jj = j + 1, m if( w( jj )<tmp1 ) then i = jj tmp1 = w( jj ) end if end do if( i/=0_${ik}$ ) then itmp1 = iwork( indibl+i-1 ) w( i ) = w( j ) iwork( indibl+i-1 ) = iwork( indibl+j-1 ) w( j ) = tmp1 iwork( indibl+j-1 ) = itmp1 call stdlib${ii}$_${ri}$swap( n, z( 1_${ik}$, i ), 1_${ik}$, z( 1_${ik}$, j ), 1_${ik}$ ) if( info/=0_${ik}$ ) then itmp1 = ifail( i ) ifail( i ) = ifail( j ) ifail( j ) = itmp1 end if end if end do end if return end subroutine stdlib${ii}$_${ri}$stevx #:endif #:endfor pure module subroutine stdlib${ii}$_spteqr( compz, n, d, e, z, ldz, work, info ) !! SPTEQR computes all eigenvalues and, optionally, eigenvectors of a !! symmetric positive definite tridiagonal matrix by first factoring the !! matrix using SPTTRF, and then calling SBDSQR to compute the singular !! values of the bidiagonal factor. !! This routine computes the eigenvalues of the positive definite !! tridiagonal matrix to high relative accuracy. This means that if the !! eigenvalues range over many orders of magnitude in size, then the !! small eigenvalues and corresponding eigenvectors will be computed !! more accurately than, for example, with the standard QR method. !! The eigenvectors of a full or band symmetric positive definite matrix !! can also be found if SSYTRD, SSPTRD, or SSBTRD has been used to !! reduce this matrix to tridiagonal form. (The reduction to tridiagonal !! form, however, may preclude the possibility of obtaining high !! relative accuracy in the small eigenvalues of the original matrix, if !! these eigenvalues range over many orders of magnitude.) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: compz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(sp), intent(inout) :: d(*), e(*), z(ldz,*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Arrays real(sp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$) ! Local Scalars integer(${ik}$) :: i, icompz, nru ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( stdlib_lsame( compz, 'N' ) ) then icompz = 0_${ik}$ else if( stdlib_lsame( compz, 'V' ) ) then icompz = 1_${ik}$ else if( stdlib_lsame( compz, 'I' ) ) then icompz = 2_${ik}$ else icompz = -1_${ik}$ end if if( icompz<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SPTEQR', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = one return end if if( icompz==2_${ik}$ )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, z, ldz ) ! call stdlib${ii}$_spttrf to factor the matrix. call stdlib${ii}$_spttrf( n, d, e, info ) if( info/=0 )return do i = 1, n d( i ) = sqrt( d( i ) ) end do do i = 1, n - 1 e( i ) = e( i )*d( i ) end do ! call stdlib${ii}$_sbdsqr to compute the singular values/vectors of the ! bidiagonal factor. if( icompz>0_${ik}$ ) then nru = n else nru = 0_${ik}$ end if call stdlib${ii}$_sbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info ) ! square the singular values. if( info==0_${ik}$ ) then do i = 1, n d( i ) = d( i )*d( i ) end do else info = n + info end if return end subroutine stdlib${ii}$_spteqr pure module subroutine stdlib${ii}$_dpteqr( compz, n, d, e, z, ldz, work, info ) !! DPTEQR computes all eigenvalues and, optionally, eigenvectors of a !! symmetric positive definite tridiagonal matrix by first factoring the !! matrix using DPTTRF, and then calling DBDSQR to compute the singular !! values of the bidiagonal factor. !! This routine computes the eigenvalues of the positive definite !! tridiagonal matrix to high relative accuracy. This means that if the !! eigenvalues range over many orders of magnitude in size, then the !! small eigenvalues and corresponding eigenvectors will be computed !! more accurately than, for example, with the standard QR method. !! The eigenvectors of a full or band symmetric positive definite matrix !! can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to !! reduce this matrix to tridiagonal form. (The reduction to tridiagonal !! form, however, may preclude the possibility of obtaining high !! relative accuracy in the small eigenvalues of the original matrix, if !! these eigenvalues range over many orders of magnitude.) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: compz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(dp), intent(inout) :: d(*), e(*), z(ldz,*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Arrays real(dp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$) ! Local Scalars integer(${ik}$) :: i, icompz, nru ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( stdlib_lsame( compz, 'N' ) ) then icompz = 0_${ik}$ else if( stdlib_lsame( compz, 'V' ) ) then icompz = 1_${ik}$ else if( stdlib_lsame( compz, 'I' ) ) then icompz = 2_${ik}$ else icompz = -1_${ik}$ end if if( icompz<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPTEQR', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = one return end if if( icompz==2_${ik}$ )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, z, ldz ) ! call stdlib${ii}$_dpttrf to factor the matrix. call stdlib${ii}$_dpttrf( n, d, e, info ) if( info/=0 )return do i = 1, n d( i ) = sqrt( d( i ) ) end do do i = 1, n - 1 e( i ) = e( i )*d( i ) end do ! call stdlib${ii}$_dbdsqr to compute the singular values/vectors of the ! bidiagonal factor. if( icompz>0_${ik}$ ) then nru = n else nru = 0_${ik}$ end if call stdlib${ii}$_dbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info ) ! square the singular values. if( info==0_${ik}$ ) then do i = 1, n d( i ) = d( i )*d( i ) end do else info = n + info end if return end subroutine stdlib${ii}$_dpteqr #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$pteqr( compz, n, d, e, z, ldz, work, info ) !! DPTEQR: computes all eigenvalues and, optionally, eigenvectors of a !! symmetric positive definite tridiagonal matrix by first factoring the !! matrix using DPTTRF, and then calling DBDSQR to compute the singular !! values of the bidiagonal factor. !! This routine computes the eigenvalues of the positive definite !! tridiagonal matrix to high relative accuracy. This means that if the !! eigenvalues range over many orders of magnitude in size, then the !! small eigenvalues and corresponding eigenvectors will be computed !! more accurately than, for example, with the standard QR method. !! The eigenvectors of a full or band symmetric positive definite matrix !! can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to !! reduce this matrix to tridiagonal form. (The reduction to tridiagonal !! form, however, may preclude the possibility of obtaining high !! relative accuracy in the small eigenvalues of the original matrix, if !! these eigenvalues range over many orders of magnitude.) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: compz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(${rk}$), intent(inout) :: d(*), e(*), z(ldz,*) real(${rk}$), intent(out) :: work(*) ! ===================================================================== ! Local Arrays real(${rk}$) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$) ! Local Scalars integer(${ik}$) :: i, icompz, nru ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( stdlib_lsame( compz, 'N' ) ) then icompz = 0_${ik}$ else if( stdlib_lsame( compz, 'V' ) ) then icompz = 1_${ik}$ else if( stdlib_lsame( compz, 'I' ) ) then icompz = 2_${ik}$ else icompz = -1_${ik}$ end if if( icompz<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPTEQR', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = one return end if if( icompz==2_${ik}$ )call stdlib${ii}$_${ri}$laset( 'FULL', n, n, zero, one, z, ldz ) ! call stdlib${ii}$_${ri}$pttrf to factor the matrix. call stdlib${ii}$_${ri}$pttrf( n, d, e, info ) if( info/=0 )return do i = 1, n d( i ) = sqrt( d( i ) ) end do do i = 1, n - 1 e( i ) = e( i )*d( i ) end do ! call stdlib${ii}$_${ri}$bdsqr to compute the singular values/vectors of the ! bidiagonal factor. if( icompz>0_${ik}$ ) then nru = n else nru = 0_${ik}$ end if call stdlib${ii}$_${ri}$bdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info ) ! square the singular values. if( info==0_${ik}$ ) then do i = 1, n d( i ) = d( i )*d( i ) end do else info = n + info end if return end subroutine stdlib${ii}$_${ri}$pteqr #:endif #:endfor pure module subroutine stdlib${ii}$_cpteqr( compz, n, d, e, z, ldz, work, info ) !! CPTEQR computes all eigenvalues and, optionally, eigenvectors of a !! symmetric positive definite tridiagonal matrix by first factoring the !! matrix using SPTTRF and then calling CBDSQR to compute the singular !! values of the bidiagonal factor. !! This routine computes the eigenvalues of the positive definite !! tridiagonal matrix to high relative accuracy. This means that if the !! eigenvalues range over many orders of magnitude in size, then the !! small eigenvalues and corresponding eigenvectors will be computed !! more accurately than, for example, with the standard QR method. !! The eigenvectors of a full or band positive definite Hermitian matrix !! can also be found if CHETRD, CHPTRD, or CHBTRD has been used to !! reduce this matrix to tridiagonal form. (The reduction to !! tridiagonal form, however, may preclude the possibility of obtaining !! high relative accuracy in the small eigenvalues of the original !! matrix, if these eigenvalues range over many orders of magnitude.) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: compz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(sp), intent(inout) :: d(*), e(*) real(sp), intent(out) :: work(*) complex(sp), intent(inout) :: z(ldz,*) ! ==================================================================== ! Local Arrays complex(sp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$) ! Local Scalars integer(${ik}$) :: i, icompz, nru ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( stdlib_lsame( compz, 'N' ) ) then icompz = 0_${ik}$ else if( stdlib_lsame( compz, 'V' ) ) then icompz = 1_${ik}$ else if( stdlib_lsame( compz, 'I' ) ) then icompz = 2_${ik}$ else icompz = -1_${ik}$ end if if( icompz<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CPTEQR', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = cone return end if if( icompz==2_${ik}$ )call stdlib${ii}$_claset( 'FULL', n, n, czero, cone, z, ldz ) ! call stdlib${ii}$_spttrf to factor the matrix. call stdlib${ii}$_spttrf( n, d, e, info ) if( info/=0 )return do i = 1, n d( i ) = sqrt( d( i ) ) end do do i = 1, n - 1 e( i ) = e( i )*d( i ) end do ! call stdlib${ii}$_cbdsqr to compute the singular values/vectors of the ! bidiagonal factor. if( icompz>0_${ik}$ ) then nru = n else nru = 0_${ik}$ end if call stdlib${ii}$_cbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info ) ! square the singular values. if( info==0_${ik}$ ) then do i = 1, n d( i ) = d( i )*d( i ) end do else info = n + info end if return end subroutine stdlib${ii}$_cpteqr pure module subroutine stdlib${ii}$_zpteqr( compz, n, d, e, z, ldz, work, info ) !! ZPTEQR computes all eigenvalues and, optionally, eigenvectors of a !! symmetric positive definite tridiagonal matrix by first factoring the !! matrix using DPTTRF and then calling ZBDSQR to compute the singular !! values of the bidiagonal factor. !! This routine computes the eigenvalues of the positive definite !! tridiagonal matrix to high relative accuracy. This means that if the !! eigenvalues range over many orders of magnitude in size, then the !! small eigenvalues and corresponding eigenvectors will be computed !! more accurately than, for example, with the standard QR method. !! The eigenvectors of a full or band positive definite Hermitian matrix !! can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to !! reduce this matrix to tridiagonal form. (The reduction to !! tridiagonal form, however, may preclude the possibility of obtaining !! high relative accuracy in the small eigenvalues of the original !! matrix, if these eigenvalues range over many orders of magnitude.) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: compz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(dp), intent(inout) :: d(*), e(*) real(dp), intent(out) :: work(*) complex(dp), intent(inout) :: z(ldz,*) ! ==================================================================== ! Local Arrays complex(dp) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$) ! Local Scalars integer(${ik}$) :: i, icompz, nru ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( stdlib_lsame( compz, 'N' ) ) then icompz = 0_${ik}$ else if( stdlib_lsame( compz, 'V' ) ) then icompz = 1_${ik}$ else if( stdlib_lsame( compz, 'I' ) ) then icompz = 2_${ik}$ else icompz = -1_${ik}$ end if if( icompz<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPTEQR', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = cone return end if if( icompz==2_${ik}$ )call stdlib${ii}$_zlaset( 'FULL', n, n, czero, cone, z, ldz ) ! call stdlib${ii}$_dpttrf to factor the matrix. call stdlib${ii}$_dpttrf( n, d, e, info ) if( info/=0 )return do i = 1, n d( i ) = sqrt( d( i ) ) end do do i = 1, n - 1 e( i ) = e( i )*d( i ) end do ! call stdlib${ii}$_zbdsqr to compute the singular values/vectors of the ! bidiagonal factor. if( icompz>0_${ik}$ ) then nru = n else nru = 0_${ik}$ end if call stdlib${ii}$_zbdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info ) ! square the singular values. if( info==0_${ik}$ ) then do i = 1, n d( i ) = d( i )*d( i ) end do else info = n + info end if return end subroutine stdlib${ii}$_zpteqr #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$pteqr( compz, n, d, e, z, ldz, work, info ) !! ZPTEQR: computes all eigenvalues and, optionally, eigenvectors of a !! symmetric positive definite tridiagonal matrix by first factoring the !! matrix using DPTTRF and then calling ZBDSQR to compute the singular !! values of the bidiagonal factor. !! This routine computes the eigenvalues of the positive definite !! tridiagonal matrix to high relative accuracy. This means that if the !! eigenvalues range over many orders of magnitude in size, then the !! small eigenvalues and corresponding eigenvectors will be computed !! more accurately than, for example, with the standard QR method. !! The eigenvectors of a full or band positive definite Hermitian matrix !! can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to !! reduce this matrix to tridiagonal form. (The reduction to !! tridiagonal form, however, may preclude the possibility of obtaining !! high relative accuracy in the small eigenvalues of the original !! matrix, if these eigenvalues range over many orders of magnitude.) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: compz integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldz, n ! Array Arguments real(${ck}$), intent(inout) :: d(*), e(*) real(${ck}$), intent(out) :: work(*) complex(${ck}$), intent(inout) :: z(ldz,*) ! ==================================================================== ! Local Arrays complex(${ck}$) :: c(1_${ik}$,1_${ik}$), vt(1_${ik}$,1_${ik}$) ! Local Scalars integer(${ik}$) :: i, icompz, nru ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( stdlib_lsame( compz, 'N' ) ) then icompz = 0_${ik}$ else if( stdlib_lsame( compz, 'V' ) ) then icompz = 1_${ik}$ else if( stdlib_lsame( compz, 'I' ) ) then icompz = 2_${ik}$ else icompz = -1_${ik}$ end if if( icompz<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ( ldz<1_${ik}$ ) .or. ( icompz>0_${ik}$ .and. ldz<max( 1_${ik}$,n ) ) ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPTEQR', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( icompz>0_${ik}$ )z( 1_${ik}$, 1_${ik}$ ) = cone return end if if( icompz==2_${ik}$ )call stdlib${ii}$_${ci}$laset( 'FULL', n, n, czero, cone, z, ldz ) ! call stdlib${ii}$_${c2ri(ci)}$pttrf to factor the matrix. call stdlib${ii}$_${c2ri(ci)}$pttrf( n, d, e, info ) if( info/=0 )return do i = 1, n d( i ) = sqrt( d( i ) ) end do do i = 1, n - 1 e( i ) = e( i )*d( i ) end do ! call stdlib${ii}$_${ci}$bdsqr to compute the singular values/vectors of the ! bidiagonal factor. if( icompz>0_${ik}$ ) then nru = n else nru = 0_${ik}$ end if call stdlib${ii}$_${ci}$bdsqr( 'LOWER', n, 0_${ik}$, nru, 0_${ik}$, d, e, vt, 1_${ik}$, z, ldz, c, 1_${ik}$,work, info ) ! square the singular values. if( info==0_${ik}$ ) then do i = 1, n d( i ) = d( i )*d( i ) end do else info = n + info end if return end subroutine stdlib${ii}$_${ci}$pteqr #:endif #:endfor pure module subroutine stdlib${ii}$_sstebz( range, order, n, vl, vu, il, iu, abstol, d, e,m, nsplit, w, & !! SSTEBZ computes the eigenvalues of a symmetric tridiagonal !! matrix T. The user may ask for all eigenvalues, all eigenvalues !! in the half-open interval (VL, VU], or the IL-th through IU-th !! eigenvalues. !! To avoid overflow, the matrix must be scaled so that its !! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest !! accuracy, it should not be much smaller than that. !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal !! Matrix", Report CS41, Computer Science Dept., Stanford !! University, July 21, 1966. iblock, isplit, work, iwork,info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: order, range integer(${ik}$), intent(in) :: il, iu, n integer(${ik}$), intent(out) :: info, m, nsplit real(sp), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: iblock(*), isplit(*), iwork(*) real(sp), intent(in) :: d(*), e(*) real(sp), intent(out) :: w(*), work(*) ! ===================================================================== ! Parameters real(sp), parameter :: fudge = 2.1_sp real(sp), parameter :: relfac = two ! Local Scalars logical(lk) :: ncnvrg, toofew integer(${ik}$) :: ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iorder, & iout, irange, itmax, itmp1, iw, iwoff, j, jb, jdisc, je, nb, nwl, nwu real(sp) :: atoli, bnorm, gl, gu, pivmin, rtoli, safemn, tmp1, tmp2, tnorm, ulp, wkill,& wl, wlu, wu, wul ! Local Arrays integer(${ik}$) :: idumma(1_${ik}$) ! Intrinsic Functions ! Executable Statements info = 0_${ik}$ ! decode range if( stdlib_lsame( range, 'A' ) ) then irange = 1_${ik}$ else if( stdlib_lsame( range, 'V' ) ) then irange = 2_${ik}$ else if( stdlib_lsame( range, 'I' ) ) then irange = 3_${ik}$ else irange = 0_${ik}$ end if ! decode order if( stdlib_lsame( order, 'B' ) ) then iorder = 2_${ik}$ else if( stdlib_lsame( order, 'E' ) ) then iorder = 1_${ik}$ else iorder = 0_${ik}$ end if ! check for errors if( irange<=0_${ik}$ ) then info = -1_${ik}$ else if( iorder<=0_${ik}$ ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( irange==2_${ik}$ ) then if( vl>=vu ) info = -5_${ik}$ else if( irange==3_${ik}$ .and. ( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) )then info = -6_${ik}$ else if( irange==3_${ik}$ .and. ( iu<min( n, il ) .or. iu>n ) )then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSTEBZ', -info ) return end if ! initialize error flags info = 0_${ik}$ ncnvrg = .false. toofew = .false. ! quick return if possible m = 0_${ik}$ if( n==0 )return ! simplifications: if( irange==3_${ik}$ .and. il==1_${ik}$ .and. iu==n )irange = 1_${ik}$ ! get machine constants ! nb is the minimum vector length for vector bisection, or 0 ! if only scalar is to be done. safemn = stdlib${ii}$_slamch( 'S' ) ulp = stdlib${ii}$_slamch( 'P' ) rtoli = ulp*relfac nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ )nb = 0_${ik}$ ! special case when n=1 if( n==1_${ik}$ ) then nsplit = 1_${ik}$ isplit( 1_${ik}$ ) = 1_${ik}$ if( irange==2_${ik}$ .and. ( vl>=d( 1_${ik}$ ) .or. vu<d( 1_${ik}$ ) ) ) then m = 0_${ik}$ else w( 1_${ik}$ ) = d( 1_${ik}$ ) iblock( 1_${ik}$ ) = 1_${ik}$ m = 1_${ik}$ end if return end if ! compute splitting points nsplit = 1_${ik}$ work( n ) = zero pivmin = one do j = 2, n tmp1 = e( j-1 )**2_${ik}$ if( abs( d( j )*d( j-1 ) )*ulp**2_${ik}$+safemn>tmp1 ) then isplit( nsplit ) = j - 1_${ik}$ nsplit = nsplit + 1_${ik}$ work( j-1 ) = zero else work( j-1 ) = tmp1 pivmin = max( pivmin, tmp1 ) end if end do isplit( nsplit ) = n pivmin = pivmin*safemn ! compute interval and atoli if( irange==3_${ik}$ ) then ! range='i': compute the interval containing eigenvalues ! il through iu. ! compute gershgorin interval for entire (split) matrix ! and use it as the initial interval gu = d( 1_${ik}$ ) gl = d( 1_${ik}$ ) tmp1 = zero do j = 1, n - 1 tmp2 = sqrt( work( j ) ) gu = max( gu, d( j )+tmp1+tmp2 ) gl = min( gl, d( j )-tmp1-tmp2 ) tmp1 = tmp2 end do gu = max( gu, d( n )+tmp1 ) gl = min( gl, d( n )-tmp1 ) tnorm = max( abs( gl ), abs( gu ) ) gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin gu = gu + fudge*tnorm*ulp*n + fudge*pivmin ! compute iteration parameters itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$ if( abstol<=zero ) then atoli = ulp*tnorm else atoli = abstol end if work( n+1 ) = gl work( n+2 ) = gl work( n+3 ) = gu work( n+4 ) = gu work( n+5 ) = gl work( n+6 ) = gu iwork( 1_${ik}$ ) = -1_${ik}$ iwork( 2_${ik}$ ) = -1_${ik}$ iwork( 3_${ik}$ ) = n + 1_${ik}$ iwork( 4_${ik}$ ) = n + 1_${ik}$ iwork( 5_${ik}$ ) = il - 1_${ik}$ iwork( 6_${ik}$ ) = iu call stdlib${ii}$_slaebz( 3_${ik}$, itmax, n, 2_${ik}$, 2_${ik}$, nb, atoli, rtoli, pivmin, d, e,work, iwork( & 5_${ik}$ ), work( n+1 ), work( n+5 ), iout,iwork, w, iblock, iinfo ) if( iwork( 6_${ik}$ )==iu ) then wl = work( n+1 ) wlu = work( n+3 ) nwl = iwork( 1_${ik}$ ) wu = work( n+4 ) wul = work( n+2 ) nwu = iwork( 4_${ik}$ ) else wl = work( n+2 ) wlu = work( n+4 ) nwl = iwork( 2_${ik}$ ) wu = work( n+3 ) wul = work( n+1 ) nwu = iwork( 3_${ik}$ ) end if if( nwl<0_${ik}$ .or. nwl>=n .or. nwu<1_${ik}$ .or. nwu>n ) then info = 4_${ik}$ return end if else ! range='a' or 'v' -- set atoli tnorm = max( abs( d( 1_${ik}$ ) )+abs( e( 1_${ik}$ ) ),abs( d( n ) )+abs( e( n-1 ) ) ) do j = 2, n - 1 tnorm = max( tnorm, abs( d( j ) )+abs( e( j-1 ) )+abs( e( j ) ) ) end do if( abstol<=zero ) then atoli = ulp*tnorm else atoli = abstol end if if( irange==2_${ik}$ ) then wl = vl wu = vu else wl = zero wu = zero end if end if ! find eigenvalues -- loop over blocks and recompute nwl and nwu. ! nwl accumulates the number of eigenvalues .le. wl, ! nwu accumulates the number of eigenvalues .le. wu m = 0_${ik}$ iend = 0_${ik}$ info = 0_${ik}$ nwl = 0_${ik}$ nwu = 0_${ik}$ loop_70: do jb = 1, nsplit ioff = iend ibegin = ioff + 1_${ik}$ iend = isplit( jb ) in = iend - ioff if( in==1_${ik}$ ) then ! special case -- in=1 if( irange==1_${ik}$ .or. wl>=d( ibegin )-pivmin )nwl = nwl + 1_${ik}$ if( irange==1_${ik}$ .or. wu>=d( ibegin )-pivmin )nwu = nwu + 1_${ik}$ if( irange==1_${ik}$ .or. ( wl<d( ibegin )-pivmin .and. wu>=d( ibegin )-pivmin ) ) & then m = m + 1_${ik}$ w( m ) = d( ibegin ) iblock( m ) = jb end if else ! general case -- in > 1 ! compute gershgorin interval ! and use it as the initial interval gu = d( ibegin ) gl = d( ibegin ) tmp1 = zero do j = ibegin, iend - 1 tmp2 = abs( e( j ) ) gu = max( gu, d( j )+tmp1+tmp2 ) gl = min( gl, d( j )-tmp1-tmp2 ) tmp1 = tmp2 end do gu = max( gu, d( iend )+tmp1 ) gl = min( gl, d( iend )-tmp1 ) bnorm = max( abs( gl ), abs( gu ) ) gl = gl - fudge*bnorm*ulp*in - fudge*pivmin gu = gu + fudge*bnorm*ulp*in + fudge*pivmin ! compute atoli for the current submatrix if( abstol<=zero ) then atoli = ulp*max( abs( gl ), abs( gu ) ) else atoli = abstol end if if( irange>1_${ik}$ ) then if( gu<wl ) then nwl = nwl + in nwu = nwu + in cycle loop_70 end if gl = max( gl, wl ) gu = min( gu, wu ) if( gl>=gu )cycle loop_70 end if ! set up initial interval work( n+1 ) = gl work( n+in+1 ) = gu call stdlib${ii}$_slaebz( 1_${ik}$, 0_${ik}$, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e( & ibegin ), work( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), im,iwork, w( m+1 & ), iblock( m+1 ), iinfo ) nwl = nwl + iwork( 1_${ik}$ ) nwu = nwu + iwork( in+1 ) iwoff = m - iwork( 1_${ik}$ ) ! compute eigenvalues itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + & 2_${ik}$ call stdlib${ii}$_slaebz( 2_${ik}$, itmax, in, in, 1_${ik}$, nb, atoli, rtoli, pivmin,d( ibegin ), e(& ibegin ), work( ibegin ),idumma, work( n+1 ), work( n+2*in+1 ), iout,iwork, w( & m+1 ), iblock( m+1 ), iinfo ) ! copy eigenvalues into w and iblock ! use -jb for block number for unconverged eigenvalues. do j = 1, iout tmp1 = half*( work( j+n )+work( j+in+n ) ) ! flag non-convergence. if( j>iout-iinfo ) then ncnvrg = .true. ib = -jb else ib = jb end if do je = iwork( j ) + 1 + iwoff,iwork( j+in ) + iwoff w( je ) = tmp1 iblock( je ) = ib end do end do m = m + im end if end do loop_70 ! if range='i', then (wl,wu) contains eigenvalues nwl+1,...,nwu ! if nwl+1 < il or nwu > iu, discard extra eigenvalues. if( irange==3_${ik}$ ) then im = 0_${ik}$ idiscl = il - 1_${ik}$ - nwl idiscu = nwu - iu if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then do je = 1, m if( w( je )<=wlu .and. idiscl>0_${ik}$ ) then idiscl = idiscl - 1_${ik}$ else if( w( je )>=wul .and. idiscu>0_${ik}$ ) then idiscu = idiscu - 1_${ik}$ else im = im + 1_${ik}$ w( im ) = w( je ) iblock( im ) = iblock( je ) end if end do m = im end if if( idiscl>0_${ik}$ .or. idiscu>0_${ik}$ ) then ! code to deal with effects of bad arithmetic: ! some low eigenvalues to be discarded are not in (wl,wlu], ! or high eigenvalues to be discarded are not in (wul,wu] ! so just kill off the smallest idiscl/largest idiscu ! eigenvalues, by simply finding the smallest/largest ! eigenvalue(s). ! (if n(w) is monotone non-decreasing, this should never ! happen.) if( idiscl>0_${ik}$ ) then wkill = wu do jdisc = 1, idiscl iw = 0_${ik}$ do je = 1, m if( iblock( je )/=0_${ik}$ .and.( w( je )<wkill .or. iw==0_${ik}$ ) ) then iw = je wkill = w( je ) end if end do iblock( iw ) = 0_${ik}$ end do end if if( idiscu>0_${ik}$ ) then wkill = wl do jdisc = 1, idiscu iw = 0_${ik}$ do je = 1, m if( iblock( je )/=0_${ik}$ .and.( w( je )>wkill .or. iw==0_${ik}$ ) ) then iw = je wkill = w( je ) end if end do iblock( iw ) = 0_${ik}$ end do end if im = 0_${ik}$ do je = 1, m if( iblock( je )/=0_${ik}$ ) then im = im + 1_${ik}$ w( im ) = w( je ) iblock( im ) = iblock( je ) end if end do m = im end if if( idiscl<0_${ik}$ .or. idiscu<0_${ik}$ ) then toofew = .true. end if end if ! if order='b', do nothing -- the eigenvalues are already sorted ! by block. ! if order='e', sort the eigenvalues from smallest to largest if( iorder==1_${ik}$ .and. nsplit>1_${ik}$ ) then do je = 1, m - 1 ie = 0_${ik}$ tmp1 = w( je ) do j = je + 1, m if( w( j )<tmp1 ) then ie = j tmp1 = w( j ) end if end do if( ie/=0_${ik}$ ) then itmp1 = iblock( ie ) w( ie ) = w( je ) iblock( ie ) = iblock( je ) w( je ) = tmp1 iblock( je ) = itmp1 end if end do end if info = 0_${ik}$ if( ncnvrg )info = info + 1_${ik}$ if( toofew )info = info + 2_${ik}$ return end subroutine stdlib${ii}$_sstebz pure module subroutine stdlib${ii}$_dstebz( range, order, n, vl, vu, il, iu, abstol, d, e,m, nsplit, w, & !! DSTEBZ computes the eigenvalues of a symmetric tridiagonal !! matrix T. The user may ask for all eigenvalues, all eigenvalues !! in the half-open interval (VL, VU], or the IL-th through IU-th !! eigenvalues. !! To avoid overflow, the matrix must be scaled so that its !! largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest !! accuracy, it should not be much smaller than that. !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal !! Matrix", Report CS41, Computer Science Dept., Stanford !! University, July 21, 1966. iblock, isplit, work, iwork,info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: order, range integer(${ik}$), intent(in) :: il, iu, n integer(${ik}$), intent(out) :: info, m, nsplit real(dp), intent(in) :: abstol, vl, vu ! Array Arguments integer(${ik}$), intent(out) :: iblock(*), isplit(*), iwork(*) real(dp), intent(in) :: d(*), e(*) real(dp), intent(out) :: w(*), work(*) ! ===================================================================== ! Parameters real(dp), parameter :: fudge = 2.1_dp real(dp), parameter :: relfac = two ! Local Scalars logical(lk) :: ncnvrg, toofew integer(${ik}$) :: ib, ibegin, idiscl, idiscu, ie, iend, iinfo, im, in, ioff, iorder, & iout, irange, itmax, itmp1, iw, iwoff, j, jb, jdisc, je, nb, nwl, nwu real(dp) :: atoli, bnorm, gl, gu, pivmin, rtoli, safemn, tmp1, tmp2, tnorm, ulp, wkill,& wl, wlu, wu, wul ! Local Arrays integer(${ik}$) :: idumma(1_${ik}$) ! Intrinsic Functions ! Executable Statements info = 0_${ik}$ ! decode range if( stdlib_lsame( range, 'A' ) ) then irange = 1_${ik}$ else if( stdlib_lsame( range, 'V' ) ) then irange = 2_${ik}$ else if( stdlib_lsame( range, 'I' ) ) then irange = 3_${ik}$ else irange = 0_${ik}$ end if ! decode order if( stdlib_lsame( order, 'B' ) ) then iorder = 2_${ik}$ else if( stdlib_lsame( order, 'E' ) ) then iorder = 1_${ik}$ else iorder = 0_${ik}$ end if ! check for errors if( irange<=0_${ik}$ ) then info = -1_${ik}$ else if( iorder<=0_${ik}$ ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( irange==2_${ik}$ ) then if( vl>=vu )info = -5_${ik}$ else if( irange==3_${ik}$ .and. ( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) )then info = -6_${ik}$ else if( irange==3_${ik}$ .and. ( iu<min( n, il ) .or. iu>n ) )then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSTEBZ', -info ) return end if ! initialize error flags info = 0_${ik}$ ncnvrg = .false. toofew = .false. ! quick return if possible m = 0_${ik}$ if( n==0 )return ! simplifications: if( irange==3_${ik}$ .and. il==1_${ik}$ .and. iu==n )irange = 1_${ik}$ ! get machine constants ! nb is the minimum vector length for vector bisection, or 0 ! if only scalar is to be done. safemn = stdlib${ii}$_dlamch( 'S' ) ulp = stdlib${ii}$_dlamch( 'P' ) rtoli = ulp*relfac nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSTEBZ', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ )nb = 0_${ik}$ ! special case when n=1 if( n==1_${ik}$ ) then nsplit = 1_${ik}$ isplit( 1_${ik}$ ) = 1_${ik}$ if( irange==2_${ik}$ .and. ( vl>=d( 1_${ik}$ ) .or. vu<d( 1_${ik}$ ) ) ) then m = 0_${ik}$ else w( 1_${ik}$ ) = d( 1_${ik}$ ) iblock( 1_${ik}$ ) = 1_${ik}$ m = 1_${ik}$ end if return end if ! compute splitting points nsplit = 1_${ik}$ work( n ) = zero pivmin = one do j = 2, n tmp1 = e( j-1 )**2_${ik}$ if( abs( d( j )*d( j-1 ) )*ulp**2_${ik}$+safemn>tmp1 ) then isplit( nsplit ) = j - 1_${ik}$ nsplit = nsplit + 1_${ik}$ work( j-1 ) = zero else work( j-1 ) = tmp1 pivmin = max( pivmin, tmp1 ) end if end do isplit( nsplit ) = n pivmin = pivmin*safemn ! compute interval and atoli if( irange==3_${ik}$ ) then ! range='i': compute the interval containing eigenvalues ! il through iu. ! compute gershgorin interval for entire (split) matrix ! and use it as the initial interval gu = d( 1_${ik}$ ) gl = d( 1_${ik}$ ) tmp1 = zero do j = 1, n - 1 tmp2 = sqrt( work( j ) ) gu = max( gu, d( j )+tmp1+tmp2 ) gl = min( gl, d( j )-tmp1-tmp2 ) tmp1 = tmp2 end do gu = max( gu, d( n )+tmp1 ) gl = min( gl, d( n )-tmp1 ) tnorm = max( abs( gl ), abs( gu ) ) gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin gu = gu + fudge*tnorm*ulp*n + fudge*pivmin ! compute iteration parameters itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /log( two ),KIND=${ik}$) + 2_${ik}$ if( abstol<=zero ) then atoli = ulp*tnorm else atoli = abstol end if work( n+1 ) = gl work( n+2 ) = gl work( n+3 ) = gu work( n+