stdlib_lapack_eigv_sym.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     module subroutine stdlib${ii}$_ssygv( itype, jobz, uplo, n, a, lda, b, ldb, w, work,lwork, info )
     !! SSYGV computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
     !! Here A and B are assumed to be symmetric and B is also
     !! positive definite.
        ! -- 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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, lwork, n
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*), b(ldb,*)
           real(sp), intent(out) :: w(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: lwkmin, lwkopt, nb, neig
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info==0_${ik}$ ) then
              lwkmin = max( 1_${ik}$, 3_${ik}$*n - 1_${ik}$ )
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSYTRD', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = max( lwkmin, ( nb + 2_${ik}$ )*n )
              work( 1_${ik}$ ) = lwkopt
              if( lwork<lwkmin .and. .not.lquery ) then
                 info = -11_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSYGV ', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_spotrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_ssygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_ssyev( jobz, uplo, n, a, lda, w, work, lwork, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_strsm( 'LEFT', uplo, trans, 'NON-UNIT', n, neig, one,b, ldb, a, lda )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_strmm( 'LEFT', uplo, trans, 'NON-UNIT', n, neig, one,b, ldb, a, lda )
                           
              end if
           end if
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_ssygv

     module subroutine stdlib${ii}$_dsygv( itype, jobz, uplo, n, a, lda, b, ldb, w, work,lwork, info )
     !! DSYGV computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
     !! Here A and B are assumed to be symmetric and B is also
     !! positive definite.
        ! -- 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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, lwork, n
           ! Array Arguments 
           real(dp), intent(inout) :: a(lda,*), b(ldb,*)
           real(dp), intent(out) :: w(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: lwkmin, lwkopt, nb, neig
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info==0_${ik}$ ) then
              lwkmin = max( 1_${ik}$, 3_${ik}$*n - 1_${ik}$ )
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYTRD', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = max( lwkmin, ( nb + 2_${ik}$ )*n )
              work( 1_${ik}$ ) = lwkopt
              if( lwork<lwkmin .and. .not.lquery ) then
                 info = -11_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGV ', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_dpotrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_dsygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_dsyev( jobz, uplo, n, a, lda, w, work, lwork, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_dtrsm( 'LEFT', uplo, trans, 'NON-UNIT', n, neig, one,b, ldb, a, lda )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_dtrmm( 'LEFT', uplo, trans, 'NON-UNIT', n, neig, one,b, ldb, a, lda )
                           
              end if
           end if
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_dsygv

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$sygv( itype, jobz, uplo, n, a, lda, b, ldb, w, work,lwork, info )
     !! DSYGV: computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
     !! Here A and B are assumed to be symmetric and B is also
     !! positive definite.
        ! -- 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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, lwork, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*)
           real(${rk}$), intent(out) :: w(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: lwkmin, lwkopt, nb, neig
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info==0_${ik}$ ) then
              lwkmin = max( 1_${ik}$, 3_${ik}$*n - 1_${ik}$ )
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYTRD', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = max( lwkmin, ( nb + 2_${ik}$ )*n )
              work( 1_${ik}$ ) = lwkopt
              if( lwork<lwkmin .and. .not.lquery ) then
                 info = -11_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGV ', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_${ri}$potrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_${ri}$sygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_${ri}$syev( jobz, uplo, n, a, lda, w, work, lwork, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, trans, 'NON-UNIT', n, neig, one,b, ldb, a, lda )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_${ri}$trmm( 'LEFT', uplo, trans, 'NON-UNIT', n, neig, one,b, ldb, a, lda )
                           
              end if
           end if
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_${ri}$sygv

#:endif
#:endfor



     module subroutine stdlib${ii}$_ssygvd( itype, jobz, uplo, n, a, lda, b, ldb, w, work,lwork, iwork, liwork,&
     !! SSYGVD computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
     !! B are assumed to be symmetric and B is also positive definite.
     !! 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.
                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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: a(lda,*), b(ldb,*)
           real(sp), intent(out) :: w(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: liopt, liwmin, lopt, lwmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( n<=1_${ik}$ ) then
              liwmin = 1_${ik}$
              lwmin = 1_${ik}$
           else if( wantz ) then
              liwmin = 3_${ik}$ + 5_${ik}$*n
              lwmin = 1_${ik}$ + 6_${ik}$*n + 2_${ik}$*n**2_${ik}$
           else
              liwmin = 1_${ik}$
              lwmin = 2_${ik}$*n + 1_${ik}$
           end if
           lopt = lwmin
           liopt = liwmin
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lopt
              iwork( 1_${ik}$ ) = liopt
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -11_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -13_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSYGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_spotrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_ssygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_ssyevd( jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork,info )
           lopt = max( real( lopt,KIND=sp), real( work( 1_${ik}$ ),KIND=sp) )
           liopt = max( real( liopt,KIND=sp), real( iwork( 1_${ik}$ ),KIND=sp) )
           if( wantz .and. info==0_${ik}$ ) then
              ! backtransform eigenvectors to the original problem.
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_strsm( 'LEFT', uplo, trans, 'NON-UNIT', n, n, one,b, ldb, a, lda )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_strmm( 'LEFT', uplo, trans, 'NON-UNIT', n, n, one,b, ldb, a, lda )
                           
              end if
           end if
           work( 1_${ik}$ ) = lopt
           iwork( 1_${ik}$ ) = liopt
           return
     end subroutine stdlib${ii}$_ssygvd

     module subroutine stdlib${ii}$_dsygvd( itype, jobz, uplo, n, a, lda, b, ldb, w, work,lwork, iwork, liwork,&
     !! DSYGVD computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
     !! B are assumed to be symmetric and B is also positive definite.
     !! 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.
                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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: a(lda,*), b(ldb,*)
           real(dp), intent(out) :: w(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: liopt, liwmin, lopt, lwmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( n<=1_${ik}$ ) then
              liwmin = 1_${ik}$
              lwmin = 1_${ik}$
           else if( wantz ) then
              liwmin = 3_${ik}$ + 5_${ik}$*n
              lwmin = 1_${ik}$ + 6_${ik}$*n + 2_${ik}$*n**2_${ik}$
           else
              liwmin = 1_${ik}$
              lwmin = 2_${ik}$*n + 1_${ik}$
           end if
           lopt = lwmin
           liopt = liwmin
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lopt
              iwork( 1_${ik}$ ) = liopt
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -11_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -13_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_dpotrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_dsygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_dsyevd( jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork,info )
           lopt = max( real( lopt,KIND=dp), real( work( 1_${ik}$ ),KIND=dp) )
           liopt = max( real( liopt,KIND=dp), real( iwork( 1_${ik}$ ),KIND=dp) )
           if( wantz .and. info==0_${ik}$ ) then
              ! backtransform eigenvectors to the original problem.
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_dtrsm( 'LEFT', uplo, trans, 'NON-UNIT', n, n, one,b, ldb, a, lda )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_dtrmm( 'LEFT', uplo, trans, 'NON-UNIT', n, n, one,b, ldb, a, lda )
                           
              end if
           end if
           work( 1_${ik}$ ) = lopt
           iwork( 1_${ik}$ ) = liopt
           return
     end subroutine stdlib${ii}$_dsygvd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$sygvd( itype, jobz, uplo, n, a, lda, b, ldb, w, work,lwork, iwork, liwork,&
     !! DSYGVD: computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
     !! B are assumed to be symmetric and B is also positive definite.
     !! 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.
                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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*)
           real(${rk}$), intent(out) :: w(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: liopt, liwmin, lopt, lwmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( n<=1_${ik}$ ) then
              liwmin = 1_${ik}$
              lwmin = 1_${ik}$
           else if( wantz ) then
              liwmin = 3_${ik}$ + 5_${ik}$*n
              lwmin = 1_${ik}$ + 6_${ik}$*n + 2_${ik}$*n**2_${ik}$
           else
              liwmin = 1_${ik}$
              lwmin = 2_${ik}$*n + 1_${ik}$
           end if
           lopt = lwmin
           liopt = liwmin
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lopt
              iwork( 1_${ik}$ ) = liopt
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -11_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -13_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_${ri}$potrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_${ri}$sygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_${ri}$syevd( jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork,info )
           lopt = max( real( lopt,KIND=${rk}$), real( work( 1_${ik}$ ),KIND=${rk}$) )
           liopt = max( real( liopt,KIND=${rk}$), real( iwork( 1_${ik}$ ),KIND=${rk}$) )
           if( wantz .and. info==0_${ik}$ ) then
              ! backtransform eigenvectors to the original problem.
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, trans, 'NON-UNIT', n, n, one,b, ldb, a, lda )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_${ri}$trmm( 'LEFT', uplo, trans, 'NON-UNIT', n, n, one,b, ldb, a, lda )
                           
              end if
           end if
           work( 1_${ik}$ ) = lopt
           iwork( 1_${ik}$ ) = liopt
           return
     end subroutine stdlib${ii}$_${ri}$sygvd

#:endif
#:endfor



     module subroutine stdlib${ii}$_ssygvx( itype, jobz, range, uplo, n, a, lda, b, ldb,vl, vu, il, iu, abstol,&
     !! SSYGVX computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
     !! and B are assumed to be symmetric and B is also positive definite.
     !! Eigenvalues and eigenvectors can be selected by specifying either a
     !! range of values or a range of indices for the desired eigenvalues.
                m, w, z, ldz, work,lwork, 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, uplo
           integer(${ik}$), intent(in) :: il, itype, iu, lda, ldb, ldz, lwork, 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) :: a(lda,*), b(ldb,*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
       ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, lquery, upper, valeig, wantz
           character :: trans
           integer(${ik}$) :: lwkmin, lwkopt, nb
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           upper = stdlib_lsame( uplo, 'U' )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           lquery = ( lwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -3_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -4_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -11_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -12_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -13_${ik}$
                 end if
              end if
           end if
           if (info==0_${ik}$) then
              if (ldz<1_${ik}$ .or. (wantz .and. ldz<n)) then
                 info = -18_${ik}$
              end if
           end if
           if( info==0_${ik}$ ) then
              lwkmin = max( 1_${ik}$, 8_${ik}$*n )
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSYTRD', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = max( lwkmin, ( nb + 3_${ik}$ )*n )
              work( 1_${ik}$ ) = lwkopt
              if( lwork<lwkmin .and. .not.lquery ) then
                 info = -20_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSYGVX', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0_${ik}$ ) then
              return
           end if
           ! form a cholesky factorization of b.
           call stdlib${ii}$_spotrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_ssygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_ssyevx( jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol,m, w, z, ldz, &
                     work, lwork, iwork, ifail, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              if( info>0_${ik}$ )m = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_strsm( 'LEFT', uplo, trans, 'NON-UNIT', n, m, one, b,ldb, z, ldz )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_strmm( 'LEFT', uplo, trans, 'NON-UNIT', n, m, one, b,ldb, z, ldz )
                           
              end if
           end if
           ! set work(1) to optimal workspace size.
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_ssygvx

     module subroutine stdlib${ii}$_dsygvx( itype, jobz, range, uplo, n, a, lda, b, ldb,vl, vu, il, iu, abstol,&
     !! DSYGVX computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
     !! and B are assumed to be symmetric and B is also positive definite.
     !! Eigenvalues and eigenvectors can be selected by specifying either a
     !! range of values or a range of indices for the desired eigenvalues.
                m, w, z, ldz, work,lwork, 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, uplo
           integer(${ik}$), intent(in) :: il, itype, iu, lda, ldb, ldz, lwork, 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) :: a(lda,*), b(ldb,*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
       ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, lquery, upper, valeig, wantz
           character :: trans
           integer(${ik}$) :: lwkmin, lwkopt, nb
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           upper = stdlib_lsame( uplo, 'U' )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           lquery = ( lwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -3_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -4_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -11_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -12_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -13_${ik}$
                 end if
              end if
           end if
           if (info==0_${ik}$) then
              if (ldz<1_${ik}$ .or. (wantz .and. ldz<n)) then
                 info = -18_${ik}$
              end if
           end if
           if( info==0_${ik}$ ) then
              lwkmin = max( 1_${ik}$, 8_${ik}$*n )
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYTRD', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = max( lwkmin, ( nb + 3_${ik}$ )*n )
              work( 1_${ik}$ ) = lwkopt
              if( lwork<lwkmin .and. .not.lquery ) then
                 info = -20_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGVX', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0_${ik}$ ) then
              return
           end if
           ! form a cholesky factorization of b.
           call stdlib${ii}$_dpotrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_dsygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_dsyevx( jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol,m, w, z, ldz, &
                     work, lwork, iwork, ifail, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              if( info>0_${ik}$ )m = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_dtrsm( 'LEFT', uplo, trans, 'NON-UNIT', n, m, one, b,ldb, z, ldz )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_dtrmm( 'LEFT', uplo, trans, 'NON-UNIT', n, m, one, b,ldb, z, ldz )
                           
              end if
           end if
           ! set work(1) to optimal workspace size.
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_dsygvx

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$sygvx( itype, jobz, range, uplo, n, a, lda, b, ldb,vl, vu, il, iu, abstol,&
     !! DSYGVX: computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
     !! and B are assumed to be symmetric and B is also positive definite.
     !! Eigenvalues and eigenvectors can be selected by specifying either a
     !! range of values or a range of indices for the desired eigenvalues.
                m, w, z, ldz, work,lwork, 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, uplo
           integer(${ik}$), intent(in) :: il, itype, iu, lda, ldb, ldz, lwork, 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) :: a(lda,*), b(ldb,*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
       ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, lquery, upper, valeig, wantz
           character :: trans
           integer(${ik}$) :: lwkmin, lwkopt, nb
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           upper = stdlib_lsame( uplo, 'U' )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           lquery = ( lwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -3_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -4_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -11_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -12_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -13_${ik}$
                 end if
              end if
           end if
           if (info==0_${ik}$) then
              if (ldz<1_${ik}$ .or. (wantz .and. ldz<n)) then
                 info = -18_${ik}$
              end if
           end if
           if( info==0_${ik}$ ) then
              lwkmin = max( 1_${ik}$, 8_${ik}$*n )
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYTRD', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = max( lwkmin, ( nb + 3_${ik}$ )*n )
              work( 1_${ik}$ ) = lwkopt
              if( lwork<lwkmin .and. .not.lquery ) then
                 info = -20_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGVX', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0_${ik}$ ) then
              return
           end if
           ! form a cholesky factorization of b.
           call stdlib${ii}$_${ri}$potrf( uplo, n, b, ldb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_${ri}$sygst( itype, uplo, n, a, lda, b, ldb, info )
           call stdlib${ii}$_${ri}$syevx( jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol,m, w, z, ldz, &
                     work, lwork, iwork, ifail, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              if( info>0_${ik}$ )m = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, trans, 'NON-UNIT', n, m, one, b,ldb, z, ldz )
                           
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 call stdlib${ii}$_${ri}$trmm( 'LEFT', uplo, trans, 'NON-UNIT', n, m, one, b,ldb, z, ldz )
                           
              end if
           end if
           ! set work(1) to optimal workspace size.
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_${ri}$sygvx

#:endif
#:endfor



     module subroutine stdlib${ii}$_sspgv( itype, jobz, uplo, n, ap, bp, w, z, ldz, work,info )
     !! SSPGV computes all the eigenvalues and, optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
     !! Here A and B are assumed to be symmetric, stored in packed format,
     !! and B is also positive definite.
        ! -- 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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, ldz, n
           ! Array Arguments 
           real(sp), intent(inout) :: ap(*), bp(*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: upper, wantz
           character :: trans
           integer(${ik}$) :: j, neig
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -9_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSPGV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_spptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_sspgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_sspev( jobz, uplo, n, ap, w, z, ldz, work, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_stpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_stpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_sspgv

     module subroutine stdlib${ii}$_dspgv( itype, jobz, uplo, n, ap, bp, w, z, ldz, work,info )
     !! DSPGV computes all the eigenvalues and, optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
     !! Here A and B are assumed to be symmetric, stored in packed format,
     !! and B is also positive definite.
        ! -- 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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, ldz, n
           ! Array Arguments 
           real(dp), intent(inout) :: ap(*), bp(*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: upper, wantz
           character :: trans
           integer(${ik}$) :: j, neig
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -9_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_dpptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_dspgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_dspev( jobz, uplo, n, ap, w, z, ldz, work, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_dtpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_dtpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_dspgv

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$spgv( itype, jobz, uplo, n, ap, bp, w, z, ldz, work,info )
     !! DSPGV: computes all the eigenvalues and, optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
     !! Here A and B are assumed to be symmetric, stored in packed format,
     !! and B is also positive definite.
        ! -- 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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, ldz, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: ap(*), bp(*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: upper, wantz
           character :: trans
           integer(${ik}$) :: j, neig
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -9_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_${ri}$pptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_${ri}$spgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_${ri}$spev( jobz, uplo, n, ap, w, z, ldz, work, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_${ri}$tpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_${ri}$tpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_${ri}$spgv

#:endif
#:endfor



     module subroutine stdlib${ii}$_sspgvd( itype, jobz, uplo, n, ap, bp, w, z, ldz, work,lwork, iwork, liwork,&
     !! SSPGVD computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
     !! B are assumed to be symmetric, stored in packed format, and B is also
     !! positive definite.
     !! 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.
                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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: ap(*), bp(*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: j, liwmin, lwmin, neig
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -9_${ik}$
           end if
           if( info==0_${ik}$ ) then
              if( n<=1_${ik}$ ) then
                 liwmin = 1_${ik}$
                 lwmin = 1_${ik}$
              else
                 if( wantz ) then
                    liwmin = 3_${ik}$ + 5_${ik}$*n
                    lwmin = 1_${ik}$ + 6_${ik}$*n + 2_${ik}$*n**2_${ik}$
                 else
                    liwmin = 1_${ik}$
                    lwmin = 2_${ik}$*n
                 end if
              end if
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -11_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -13_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSPGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of bp.
           call stdlib${ii}$_spptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_sspgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_sspevd( jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork,liwork, info )
                     
           lwmin = max( real( lwmin,KIND=sp), real( work( 1_${ik}$ ),KIND=sp) )
           liwmin = max( real( liwmin,KIND=sp), real( iwork( 1_${ik}$ ),KIND=sp) )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t *y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_stpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t *y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_stpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_sspgvd

     module subroutine stdlib${ii}$_dspgvd( itype, jobz, uplo, n, ap, bp, w, z, ldz, work,lwork, iwork, liwork,&
     !! DSPGVD computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
     !! B are assumed to be symmetric, stored in packed format, and B is also
     !! positive definite.
     !! 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.
                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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: ap(*), bp(*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: j, liwmin, lwmin, neig
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -9_${ik}$
           end if
           if( info==0_${ik}$ ) then
              if( n<=1_${ik}$ ) then
                 liwmin = 1_${ik}$
                 lwmin = 1_${ik}$
              else
                 if( wantz ) then
                    liwmin = 3_${ik}$ + 5_${ik}$*n
                    lwmin = 1_${ik}$ + 6_${ik}$*n + 2_${ik}$*n**2_${ik}$
                 else
                    liwmin = 1_${ik}$
                    lwmin = 2_${ik}$*n
                 end if
              end if
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -11_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -13_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of bp.
           call stdlib${ii}$_dpptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_dspgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_dspevd( jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork,liwork, info )
                     
           lwmin = max( real( lwmin,KIND=dp), real( work( 1_${ik}$ ),KIND=dp) )
           liwmin = max( real( liwmin,KIND=dp), real( iwork( 1_${ik}$ ),KIND=dp) )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t *y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_dtpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t *y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_dtpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_dspgvd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$spgvd( itype, jobz, uplo, n, ap, bp, w, z, ldz, work,lwork, iwork, liwork,&
     !! DSPGVD: computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
     !! B are assumed to be symmetric, stored in packed format, and B is also
     !! positive definite.
     !! 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.
                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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(inout) :: ap(*), bp(*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: trans
           integer(${ik}$) :: j, liwmin, lwmin, neig
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -9_${ik}$
           end if
           if( info==0_${ik}$ ) then
              if( n<=1_${ik}$ ) then
                 liwmin = 1_${ik}$
                 lwmin = 1_${ik}$
              else
                 if( wantz ) then
                    liwmin = 3_${ik}$ + 5_${ik}$*n
                    lwmin = 1_${ik}$ + 6_${ik}$*n + 2_${ik}$*n**2_${ik}$
                 else
                    liwmin = 1_${ik}$
                    lwmin = 2_${ik}$*n
                 end if
              end if
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -11_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -13_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a cholesky factorization of bp.
           call stdlib${ii}$_${ri}$pptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_${ri}$spgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_${ri}$spevd( jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork,liwork, info )
                     
           lwmin = max( real( lwmin,KIND=${rk}$), real( work( 1_${ik}$ ),KIND=${rk}$) )
           liwmin = max( real( liwmin,KIND=${rk}$), real( iwork( 1_${ik}$ ),KIND=${rk}$) )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              neig = n
              if( info>0_${ik}$ )neig = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t *y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_${ri}$tpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t *y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, neig
                    call stdlib${ii}$_${ri}$tpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_${ri}$spgvd

#:endif
#:endfor



     module subroutine stdlib${ii}$_sspgvx( itype, jobz, range, uplo, n, ap, bp, vl, vu,il, iu, abstol, m, w, &
     !! SSPGVX computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
     !! and B are assumed to be symmetric, stored in packed storage, and B
     !! is also positive definite.  Eigenvalues and eigenvectors can be
     !! selected by specifying either a range of values or a range of indices
     !! for the desired eigenvalues.
               z, ldz, 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, uplo
           integer(${ik}$), intent(in) :: il, itype, 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) :: ap(*), bp(*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
       ! =====================================================================
           ! Local Scalars 
           logical(lk) :: alleig, indeig, upper, valeig, wantz
           character :: trans
           integer(${ik}$) :: j
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           upper = stdlib_lsame( uplo, 'U' )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -3_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -4_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl ) then
                    info = -9_${ik}$
                 end if
              else if( indeig ) then
                 if( il<1_${ik}$ ) then
                    info = -10_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -11_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -16_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSPGVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_spptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_sspgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_sspevx( jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m,w, z, ldz, &
                     work, iwork, ifail, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              if( info>0_${ik}$ )m = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, m
                    call stdlib${ii}$_stpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, m
                    call stdlib${ii}$_stpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_sspgvx

     module subroutine stdlib${ii}$_dspgvx( itype, jobz, range, uplo, n, ap, bp, vl, vu,il, iu, abstol, m, w, &
     !! DSPGVX computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
     !! and B are assumed to be symmetric, stored in packed storage, and B
     !! is also positive definite.  Eigenvalues and eigenvectors can be
     !! selected by specifying either a range of values or a range of indices
     !! for the desired eigenvalues.
               z, ldz, 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, uplo
           integer(${ik}$), intent(in) :: il, itype, 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) :: ap(*), bp(*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
       ! =====================================================================
           ! Local Scalars 
           logical(lk) :: alleig, indeig, upper, valeig, wantz
           character :: trans
           integer(${ik}$) :: j
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           upper = stdlib_lsame( uplo, 'U' )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -3_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -4_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl ) then
                    info = -9_${ik}$
                 end if
              else if( indeig ) then
                 if( il<1_${ik}$ ) then
                    info = -10_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -11_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -16_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_dpptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_dspgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_dspevx( jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m,w, z, ldz, &
                     work, iwork, ifail, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              if( info>0_${ik}$ )m = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, m
                    call stdlib${ii}$_dtpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, m
                    call stdlib${ii}$_dtpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_dspgvx

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     module subroutine stdlib${ii}$_${ri}$spgvx( itype, jobz, range, uplo, n, ap, bp, vl, vu,il, iu, abstol, m, w, &
     !! DSPGVX: computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite eigenproblem, of the form
     !! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
     !! and B are assumed to be symmetric, stored in packed storage, and B
     !! is also positive definite.  Eigenvalues and eigenvectors can be
     !! selected by specifying either a range of values or a range of indices
     !! for the desired eigenvalues.
               z, ldz, 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, uplo
           integer(${ik}$), intent(in) :: il, itype, 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) :: ap(*), bp(*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
       ! =====================================================================
           ! Local Scalars 
           logical(lk) :: alleig, indeig, upper, valeig, wantz
           character :: trans
           integer(${ik}$) :: j
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           upper = stdlib_lsame( uplo, 'U' )
           wantz = stdlib_lsame( jobz, 'V' )
           alleig = stdlib_lsame( range, 'A' )
           valeig = stdlib_lsame( range, 'V' )
           indeig = stdlib_lsame( range, 'I' )
           info = 0_${ik}$
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -2_${ik}$
           else if( .not.( alleig .or. valeig .or. indeig ) ) then
              info = -3_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -4_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl ) then
                    info = -9_${ik}$
                 end if
              else if( indeig ) then
                 if( il<1_${ik}$ ) then
                    info = -10_${ik}$
                 else if( iu<min( n, il ) .or. iu>n ) then
                    info = -11_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$ ) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -16_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           ! form a cholesky factorization of b.
           call stdlib${ii}$_${ri}$pptrf( uplo, n, bp, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem and solve.
           call stdlib${ii}$_${ri}$spgst( itype, uplo, n, ap, bp, info )
           call stdlib${ii}$_${ri}$spevx( jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m,w, z, ldz, &
                     work, iwork, ifail, info )
           if( wantz ) then
              ! backtransform eigenvectors to the original problem.
              if( info>0_${ik}$ )m = info - 1_${ik}$
              if( itype==1_${ik}$ .or. itype==2_${ik}$ ) then
                 ! for a*x=(lambda)*b*x and a*b*x=(lambda)*x;
                 ! backtransform eigenvectors: x = inv(l)**t*y or inv(u)*y
                 if( upper ) then
                    trans = 'N'
                 else
                    trans = 'T'
                 end if
                 do j = 1, m
                    call stdlib${ii}$_${ri}$tpsv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              else if( itype==3_${ik}$ ) then
                 ! for b*a*x=(lambda)*x;
                 ! backtransform eigenvectors: x = l*y or u**t*y
                 if( upper ) then
                    trans = 'T'
                 else
                    trans = 'N'
                 end if
                 do j = 1, m
                    call stdlib${ii}$_${ri}$tpmv( uplo, trans, 'NON-UNIT', n, bp, z( 1_${ik}$, j ),1_${ik}$ )
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_${ri}$spgvx

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_ssbgv( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z,ldz, work, &
     !! SSBGV computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of
     !! the form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric
     !! and banded, and B is also positive definite.
               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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldz, n
           ! Array Arguments 
           real(sp), intent(inout) :: ab(ldab,*), bb(ldbb,*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: upper, wantz
           character :: vect
           integer(${ik}$) :: iinfo, inde, indwrk
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -5_${ik}$
           else if( ldab<ka+1 ) then
              info = -7_${ik}$
           else if( ldbb<kb+1 ) then
              info = -9_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -12_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSBGV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_spbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           inde = 1_${ik}$
           indwrk = inde + n
           call stdlib${ii}$_ssbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,work( indwrk ), &
                     iinfo )
           ! reduce to tridiagonal form.
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_ssbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,work( indwrk )&
                     , iinfo )
           ! for eigenvalues only, call stdlib${ii}$_ssterf.  for eigenvectors, call stdlib${ii}$_ssteqr.
           if( .not.wantz ) then
              call stdlib${ii}$_ssterf( n, w, work( inde ), info )
           else
              call stdlib${ii}$_ssteqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),info )
           end if
           return
     end subroutine stdlib${ii}$_ssbgv

     pure module subroutine stdlib${ii}$_dsbgv( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z,ldz, work, &
     !! DSBGV computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of
     !! the form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric
     !! and banded, and B is also positive definite.
               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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldz, n
           ! Array Arguments 
           real(dp), intent(inout) :: ab(ldab,*), bb(ldbb,*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: upper, wantz
           character :: vect
           integer(${ik}$) :: iinfo, inde, indwrk
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -5_${ik}$
           else if( ldab<ka+1 ) then
              info = -7_${ik}$
           else if( ldbb<kb+1 ) then
              info = -9_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -12_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSBGV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_dpbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           inde = 1_${ik}$
           indwrk = inde + n
           call stdlib${ii}$_dsbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,work( indwrk ), &
                     iinfo )
           ! reduce to tridiagonal form.
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_dsbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,work( indwrk )&
                     , iinfo )
           ! for eigenvalues only, call stdlib${ii}$_dsterf.  for eigenvectors, call stdlib${ii}$_ssteqr.
           if( .not.wantz ) then
              call stdlib${ii}$_dsterf( n, w, work( inde ), info )
           else
              call stdlib${ii}$_dsteqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),info )
           end if
           return
     end subroutine stdlib${ii}$_dsbgv

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sbgv( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z,ldz, work, &
     !! DSBGV: computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of
     !! the form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric
     !! and banded, and B is also positive definite.
               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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldz, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: ab(ldab,*), bb(ldbb,*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: upper, wantz
           character :: vect
           integer(${ik}$) :: iinfo, inde, indwrk
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           info = 0_${ik}$
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -5_${ik}$
           else if( ldab<ka+1 ) then
              info = -7_${ik}$
           else if( ldbb<kb+1 ) then
              info = -9_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -12_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSBGV ', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_${ri}$pbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           inde = 1_${ik}$
           indwrk = inde + n
           call stdlib${ii}$_${ri}$sbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,work( indwrk ), &
                     iinfo )
           ! reduce to tridiagonal form.
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_${ri}$sbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,work( indwrk )&
                     , iinfo )
           ! for eigenvalues only, call stdlib${ii}$_${ri}$sterf.  for eigenvectors, call stdlib${ii}$_dsteqr.
           if( .not.wantz ) then
              call stdlib${ii}$_${ri}$sterf( n, w, work( inde ), info )
           else
              call stdlib${ii}$_${ri}$steqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),info )
           end if
           return
     end subroutine stdlib${ii}$_${ri}$sbgv

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_ssbgvd( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w,z, ldz, work, &
     !! SSBGVD computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of the
     !! form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric and
     !! banded, and B is also positive definite.  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.
               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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: ab(ldab,*), bb(ldbb,*)
           real(sp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: vect
           integer(${ik}$) :: iinfo, inde, indwk2, indwrk, liwmin, llwrk2, lwmin
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( n<=1_${ik}$ ) then
              liwmin = 1_${ik}$
              lwmin = 1_${ik}$
           else if( wantz ) then
              liwmin = 3_${ik}$ + 5_${ik}$*n
              lwmin = 1_${ik}$ + 5_${ik}$*n + 2_${ik}$*n**2_${ik}$
           else
              liwmin = 1_${ik}$
              lwmin = 2_${ik}$*n
           end if
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -5_${ik}$
           else if( ldab<ka+1 ) then
              info = -7_${ik}$
           else if( ldbb<kb+1 ) then
              info = -9_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -12_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -14_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -16_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSBGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_spbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           inde = 1_${ik}$
           indwrk = inde + n
           indwk2 = indwrk + n*n
           llwrk2 = lwork - indwk2 + 1_${ik}$
           call stdlib${ii}$_ssbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,work, iinfo )
                     
           ! reduce to tridiagonal form.
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_ssbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,work( indwrk )&
                     , iinfo )
           ! for eigenvalues only, call stdlib${ii}$_ssterf. for eigenvectors, call stdlib${ii}$_sstedc.
           if( .not.wantz ) then
              call stdlib${ii}$_ssterf( n, w, work( inde ), info )
           else
              call stdlib${ii}$_sstedc( 'I', n, w, work( inde ), work( indwrk ), n,work( indwk2 ), &
                        llwrk2, iwork, liwork, info )
              call stdlib${ii}$_sgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,zero, work( &
                        indwk2 ), n )
              call stdlib${ii}$_slacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_ssbgvd

     pure module subroutine stdlib${ii}$_dsbgvd( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w,z, ldz, work, &
     !! DSBGVD computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of the
     !! form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric and
     !! banded, and B is also positive definite.  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.
               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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: ab(ldab,*), bb(ldbb,*)
           real(dp), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: vect
           integer(${ik}$) :: iinfo, inde, indwk2, indwrk, liwmin, llwrk2, lwmin
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( n<=1_${ik}$ ) then
              liwmin = 1_${ik}$
              lwmin = 1_${ik}$
           else if( wantz ) then
              liwmin = 3_${ik}$ + 5_${ik}$*n
              lwmin = 1_${ik}$ + 5_${ik}$*n + 2_${ik}$*n**2_${ik}$
           else
              liwmin = 1_${ik}$
              lwmin = 2_${ik}$*n
           end if
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -5_${ik}$
           else if( ldab<ka+1 ) then
              info = -7_${ik}$
           else if( ldbb<kb+1 ) then
              info = -9_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -12_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -14_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -16_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSBGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_dpbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           inde = 1_${ik}$
           indwrk = inde + n
           indwk2 = indwrk + n*n
           llwrk2 = lwork - indwk2 + 1_${ik}$
           call stdlib${ii}$_dsbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,work, iinfo )
                     
           ! reduce to tridiagonal form.
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_dsbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,work( indwrk )&
                     , iinfo )
           ! for eigenvalues only, call stdlib${ii}$_dsterf. for eigenvectors, call stdlib${ii}$_sstedc.
           if( .not.wantz ) then
              call stdlib${ii}$_dsterf( n, w, work( inde ), info )
           else
              call stdlib${ii}$_dstedc( 'I', n, w, work( inde ), work( indwrk ), n,work( indwk2 ), &
                        llwrk2, iwork, liwork, info )
              call stdlib${ii}$_dgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,zero, work( &
                        indwk2 ), n )
              call stdlib${ii}$_dlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_dsbgvd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sbgvd( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w,z, ldz, work, &
     !! DSBGVD: computes all the eigenvalues, and optionally, the eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of the
     !! form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric and
     !! banded, and B is also positive definite.  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.
               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, uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldz, liwork, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(inout) :: ab(ldab,*), bb(ldbb,*)
           real(${rk}$), intent(out) :: w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery, upper, wantz
           character :: vect
           integer(${ik}$) :: iinfo, inde, indwk2, indwrk, liwmin, llwrk2, lwmin
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           info = 0_${ik}$
           if( n<=1_${ik}$ ) then
              liwmin = 1_${ik}$
              lwmin = 1_${ik}$
           else if( wantz ) then
              liwmin = 3_${ik}$ + 5_${ik}$*n
              lwmin = 1_${ik}$ + 5_${ik}$*n + 2_${ik}$*n**2_${ik}$
           else
              liwmin = 1_${ik}$
              lwmin = 2_${ik}$*n
           end if
           if( .not.( wantz .or. stdlib_lsame( jobz, 'N' ) ) ) then
              info = -1_${ik}$
           else if( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -5_${ik}$
           else if( ldab<ka+1 ) then
              info = -7_${ik}$
           else if( ldbb<kb+1 ) then
              info = -9_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -12_${ik}$
           end if
           if( info==0_${ik}$ ) then
              work( 1_${ik}$ ) = lwmin
              iwork( 1_${ik}$ ) = liwmin
              if( lwork<lwmin .and. .not.lquery ) then
                 info = -14_${ik}$
              else if( liwork<liwmin .and. .not.lquery ) then
                 info = -16_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSBGVD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_${ri}$pbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           inde = 1_${ik}$
           indwrk = inde + n
           indwk2 = indwrk + n*n
           llwrk2 = lwork - indwk2 + 1_${ik}$
           call stdlib${ii}$_${ri}$sbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,work, iinfo )
                     
           ! reduce to tridiagonal form.
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_${ri}$sbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,work( indwrk )&
                     , iinfo )
           ! for eigenvalues only, call stdlib${ii}$_${ri}$sterf. for eigenvectors, call stdlib${ii}$_dstedc.
           if( .not.wantz ) then
              call stdlib${ii}$_${ri}$sterf( n, w, work( inde ), info )
           else
              call stdlib${ii}$_${ri}$stedc( 'I', n, w, work( inde ), work( indwrk ), n,work( indwk2 ), &
                        llwrk2, iwork, liwork, info )
              call stdlib${ii}$_${ri}$gemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,zero, work( &
                        indwk2 ), n )
              call stdlib${ii}$_${ri}$lacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_${ri}$sbgvd

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_ssbgvx( jobz, range, uplo, n, ka, kb, ab, ldab, bb,ldbb, q, ldq, vl, &
     !! SSBGVX computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of
     !! the form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric
     !! and banded, and B is also positive definite.  Eigenvalues and
     !! eigenvectors can be selected by specifying either all eigenvalues,
     !! a range of values or a range of indices for the desired eigenvalues.
               vu, il, iu, abstol, m, w, z,ldz, 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, uplo
           integer(${ik}$), intent(in) :: il, iu, ka, kb, ldab, ldbb, ldq, 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) :: ab(ldab,*), bb(ldbb,*)
           real(sp), intent(out) :: q(ldq,*), w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, upper, valeig, wantz
           character :: order, vect
           integer(${ik}$) :: i, iinfo, indd, inde, indee, indibl, indisp, indiwo, indwrk, itmp1, j, &
                     jj, nsplit
           real(sp) :: tmp1
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           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( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -5_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -6_${ik}$
           else if( ldab<ka+1 ) then
              info = -8_${ik}$
           else if( ldbb<kb+1 ) then
              info = -10_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantz .and. ldq<n ) ) then
              info = -12_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -14_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -15_${ik}$
                 else if ( iu<min( n, il ) .or. iu>n ) then
                    info = -16_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -21_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSBGVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_spbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           call stdlib${ii}$_ssbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, ldq,work, iinfo )
                     
           ! reduce symmetric band matrix to tridiagonal form.
           indd = 1_${ik}$
           inde = indd + n
           indwrk = inde + n
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_ssbtrd( vect, uplo, n, ka, ab, ldab, work( indd ),work( inde ), q, ldq, &
                     work( indwrk ), iinfo )
           ! if all eigenvalues are desired and abstol is less than or equal
           ! to 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, work( indd ), 1_${ik}$, w, 1_${ik}$ )
              indee = indwrk + 2_${ik}$*n
              call stdlib${ii}$_scopy( n-1, work( inde ), 1_${ik}$, work( indee ), 1_${ik}$ )
              if( .not.wantz ) then
                 call stdlib${ii}$_ssterf( n, w, work( indee ), info )
              else
                 call stdlib${ii}$_slacpy( 'A', n, n, q, ldq, z, ldz )
                 call stdlib${ii}$_ssteqr( jobz, n, w, work( indee ), 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 30
              end if
              info = 0_${ik}$
           end if
           ! otherwise, call stdlib${ii}$_sstebz and, if eigenvectors are desired,
           ! call stdlib${ii}$_sstein.
           if( wantz ) then
              order = 'B'
           else
              order = 'E'
           end if
           indibl = 1_${ik}$
           indisp = indibl + n
           indiwo = indisp + n
           call stdlib${ii}$_sstebz( range, order, n, vl, vu, il, iu, abstol,work( indd ), work( inde ),&
            m, nsplit, w,iwork( indibl ), iwork( indisp ), work( indwrk ),iwork( indiwo ), info )
                      
           if( wantz ) then
              call stdlib${ii}$_sstein( n, work( indd ), work( inde ), m, w,iwork( indibl ), iwork( &
                        indisp ), z, ldz,work( indwrk ), iwork( indiwo ), ifail, info )
              ! apply transformation matrix used in reduction to tridiagonal
              ! form to eigenvectors returned by stdlib${ii}$_sstein.
              do j = 1, m
                 call stdlib${ii}$_scopy( n, z( 1_${ik}$, j ), 1_${ik}$, work( 1_${ik}$ ), 1_${ik}$ )
                 call stdlib${ii}$_sgemv( 'N', n, n, one, q, ldq, work, 1_${ik}$, zero,z( 1_${ik}$, j ), 1_${ik}$ )
              end do
           end if
           30 continue
           ! 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}$_ssbgvx

     pure module subroutine stdlib${ii}$_dsbgvx( jobz, range, uplo, n, ka, kb, ab, ldab, bb,ldbb, q, ldq, vl, &
     !! DSBGVX computes selected eigenvalues, and optionally, eigenvectors
     !! of a real generalized symmetric-definite banded eigenproblem, of
     !! the form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric
     !! and banded, and B is also positive definite.  Eigenvalues and
     !! eigenvectors can be selected by specifying either all eigenvalues,
     !! a range of values or a range of indices for the desired eigenvalues.
               vu, il, iu, abstol, m, w, z,ldz, 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, uplo
           integer(${ik}$), intent(in) :: il, iu, ka, kb, ldab, ldbb, ldq, 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) :: ab(ldab,*), bb(ldbb,*)
           real(dp), intent(out) :: q(ldq,*), w(*), work(*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: alleig, indeig, test, upper, valeig, wantz
           character :: order, vect
           integer(${ik}$) :: i, iinfo, indd, inde, indee, indibl, indisp, indiwo, indwrk, itmp1, j, &
                     jj, nsplit
           real(dp) :: tmp1
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           wantz = stdlib_lsame( jobz, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           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( .not.( upper .or. stdlib_lsame( uplo, 'L' ) ) ) then
              info = -3_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ka<0_${ik}$ ) then
              info = -5_${ik}$
           else if( kb<0_${ik}$ .or. kb>ka ) then
              info = -6_${ik}$
           else if( ldab<ka+1 ) then
              info = -8_${ik}$
           else if( ldbb<kb+1 ) then
              info = -10_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantz .and. ldq<n ) ) then
              info = -12_${ik}$
           else
              if( valeig ) then
                 if( n>0_${ik}$ .and. vu<=vl )info = -14_${ik}$
              else if( indeig ) then
                 if( il<1_${ik}$ .or. il>max( 1_${ik}$, n ) ) then
                    info = -15_${ik}$
                 else if ( iu<min( n, il ) .or. iu>n ) then
                    info = -16_${ik}$
                 end if
              end if
           end if
           if( info==0_${ik}$) then
              if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
                 info = -21_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSBGVX', -info )
              return
           end if
           ! quick return if possible
           m = 0_${ik}$
           if( n==0 )return
           ! form a split cholesky factorization of b.
           call stdlib${ii}$_dpbstf( uplo, n, kb, bb, ldbb, info )
           if( info/=0_${ik}$ ) then
              info = n + info
              return
           end if
           ! transform problem to standard eigenvalue problem.
           call stdlib${ii}$_dsbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, ldq,work, iinfo )
                     
           ! reduce symmetric band matrix to tridiagonal form.
           indd = 1_${ik}$
           inde = indd + n
           indwrk = inde + n
           if( wantz ) then
              vect = 'U'
           else
              vect = 'N'
           end if
           call stdlib${ii}$_dsbtrd( vect, uplo, n, ka, ab, ldab, work( indd ),work( inde ), q, ldq, &
                     work( indwrk ), iinfo )
           ! if all eigenvalues are desired and abstol is less than or equal
           ! to 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, work( indd ), 1_${ik}$, w, 1_${ik}$ )
              indee = indwrk + 2_${ik}$*n
              call stdlib${ii}$_dcopy( n-1, work( inde ), 1_${ik}$, work( indee ),