stdlib_lapack_eigv_sym_comp.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_ssygst( itype, uplo, n, a, lda, b, ldb, info )
     !! SSYGST reduces a real symmetric-definite generalized eigenproblem
     !! to standard form.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
     !! B must have been previously factorized as U**T*U or L*L**T by SPOTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, n
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*)
           real(sp), intent(in) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: k, kb, nb
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -5_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSYGST', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! determine the block size for this environment.
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSYGST', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ .or. nb>=n ) then
              ! use unblocked code
              call stdlib${ii}$_ssygs2( itype, uplo, n, a, lda, b, ldb, info )
           else
              ! use blocked code
              if( itype==1_${ik}$ ) then
                 if( upper ) then
                    ! compute inv(u**t)*a*inv(u)
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the upper triangle of a(k:n,k:n)
                       call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                       if( k+kb<=n ) then
                          call stdlib${ii}$_strsm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT',kb, n-k-kb+1, &
                                    one, b( k, k ), ldb,a( k, k+kb ), lda )
                          call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( &
                                    k, k+kb ), ldb, one,a( k, k+kb ), lda )
                          call stdlib${ii}$_ssyr2k( uplo, 'TRANSPOSE', n-k-kb+1, kb, -one,a( k, k+kb ), &
                                    lda, b( k, k+kb ), ldb,one, a( k+kb, k+kb ), lda )
                          call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( &
                                    k, k+kb ), ldb, one,a( k, k+kb ), lda )
                          call stdlib${ii}$_strsm( 'RIGHT', uplo, 'NO TRANSPOSE','NON-UNIT', kb, n-k-kb+&
                                    1_${ik}$, one,b( k+kb, k+kb ), ldb, a( k, k+kb ),lda )
                       end if
                    end do
                 else
                    ! compute inv(l)*a*inv(l**t)
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the lower triangle of a(k:n,k:n)
                       call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                       if( k+kb<=n ) then
                          call stdlib${ii}$_strsm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',n-k-kb+1, kb, &
                                    one, b( k, k ), ldb,a( k+kb, k ), lda )
                          call stdlib${ii}$_ssymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(&
                                     k+kb, k ), ldb, one,a( k+kb, k ), lda )
                          call stdlib${ii}$_ssyr2k( uplo, 'NO TRANSPOSE', n-k-kb+1, kb,-one, a( k+kb, k &
                                    ), lda, b( k+kb, k ),ldb, one, a( k+kb, k+kb ), lda )
                          call stdlib${ii}$_ssymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(&
                                     k+kb, k ), ldb, one,a( k+kb, k ), lda )
                          call stdlib${ii}$_strsm( 'LEFT', uplo, 'NO TRANSPOSE','NON-UNIT', n-k-kb+1, &
                                    kb, one,b( k+kb, k+kb ), ldb, a( k+kb, k ),lda )
                       end if
                    end do
                 end if
              else
                 if( upper ) then
                    ! compute u*a*u**t
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the upper triangle of a(1:k+kb-1,1:k+kb-1)
                       call stdlib${ii}$_strmm( 'LEFT', uplo, 'NO TRANSPOSE', 'NON-UNIT',k-1, kb, one, &
                                 b, ldb, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_ssymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), &
                                 ldb, one, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_ssyr2k( uplo, 'NO TRANSPOSE', k-1, kb, one,a( 1_${ik}$, k ), lda, b( &
                                 1_${ik}$, k ), ldb, one, a,lda )
                       call stdlib${ii}$_ssymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), &
                                 ldb, one, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_strmm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',k-1, kb, one, b( &
                                 k, k ), ldb, a( 1_${ik}$, k ),lda )
                       call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                    end do
                 else
                    ! compute l**t*a*l
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the lower triangle of a(1:k+kb-1,1:k+kb-1)
                       call stdlib${ii}$_strmm( 'RIGHT', uplo, 'NO TRANSPOSE', 'NON-UNIT',kb, k-1, one, &
                                 b, ldb, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), &
                                 ldb, one, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_ssyr2k( uplo, 'TRANSPOSE', k-1, kb, one,a( k, 1_${ik}$ ), lda, b( k, &
                                 1_${ik}$ ), ldb, one, a,lda )
                       call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), &
                                 ldb, one, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_strmm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT', kb,k-1, one, b( &
                                 k, k ), ldb, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                    end do
                 end if
              end if
           end if
           return
     end subroutine stdlib${ii}$_ssygst

     pure module subroutine stdlib${ii}$_dsygst( itype, uplo, n, a, lda, b, ldb, info )
     !! DSYGST reduces a real symmetric-definite generalized eigenproblem
     !! to standard form.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
     !! B must have been previously factorized as U**T*U or L*L**T by DPOTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, n
           ! Array Arguments 
           real(dp), intent(inout) :: a(lda,*)
           real(dp), intent(in) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: k, kb, nb
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -5_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGST', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! determine the block size for this environment.
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYGST', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ .or. nb>=n ) then
              ! use unblocked code
              call stdlib${ii}$_dsygs2( itype, uplo, n, a, lda, b, ldb, info )
           else
              ! use blocked code
              if( itype==1_${ik}$ ) then
                 if( upper ) then
                    ! compute inv(u**t)*a*inv(u)
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the upper triangle of a(k:n,k:n)
                       call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                       if( k+kb<=n ) then
                          call stdlib${ii}$_dtrsm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT',kb, n-k-kb+1, &
                                    one, b( k, k ), ldb,a( k, k+kb ), lda )
                          call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( &
                                    k, k+kb ), ldb, one,a( k, k+kb ), lda )
                          call stdlib${ii}$_dsyr2k( uplo, 'TRANSPOSE', n-k-kb+1, kb, -one,a( k, k+kb ), &
                                    lda, b( k, k+kb ), ldb,one, a( k+kb, k+kb ), lda )
                          call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( &
                                    k, k+kb ), ldb, one,a( k, k+kb ), lda )
                          call stdlib${ii}$_dtrsm( 'RIGHT', uplo, 'NO TRANSPOSE','NON-UNIT', kb, n-k-kb+&
                                    1_${ik}$, one,b( k+kb, k+kb ), ldb, a( k, k+kb ),lda )
                       end if
                    end do
                 else
                    ! compute inv(l)*a*inv(l**t)
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the lower triangle of a(k:n,k:n)
                       call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                       if( k+kb<=n ) then
                          call stdlib${ii}$_dtrsm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',n-k-kb+1, kb, &
                                    one, b( k, k ), ldb,a( k+kb, k ), lda )
                          call stdlib${ii}$_dsymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(&
                                     k+kb, k ), ldb, one,a( k+kb, k ), lda )
                          call stdlib${ii}$_dsyr2k( uplo, 'NO TRANSPOSE', n-k-kb+1, kb,-one, a( k+kb, k &
                                    ), lda, b( k+kb, k ),ldb, one, a( k+kb, k+kb ), lda )
                          call stdlib${ii}$_dsymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(&
                                     k+kb, k ), ldb, one,a( k+kb, k ), lda )
                          call stdlib${ii}$_dtrsm( 'LEFT', uplo, 'NO TRANSPOSE','NON-UNIT', n-k-kb+1, &
                                    kb, one,b( k+kb, k+kb ), ldb, a( k+kb, k ),lda )
                       end if
                    end do
                 end if
              else
                 if( upper ) then
                    ! compute u*a*u**t
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the upper triangle of a(1:k+kb-1,1:k+kb-1)
                       call stdlib${ii}$_dtrmm( 'LEFT', uplo, 'NO TRANSPOSE', 'NON-UNIT',k-1, kb, one, &
                                 b, ldb, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_dsymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), &
                                 ldb, one, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_dsyr2k( uplo, 'NO TRANSPOSE', k-1, kb, one,a( 1_${ik}$, k ), lda, b( &
                                 1_${ik}$, k ), ldb, one, a,lda )
                       call stdlib${ii}$_dsymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), &
                                 ldb, one, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_dtrmm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',k-1, kb, one, b( &
                                 k, k ), ldb, a( 1_${ik}$, k ),lda )
                       call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                    end do
                 else
                    ! compute l**t*a*l
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the lower triangle of a(1:k+kb-1,1:k+kb-1)
                       call stdlib${ii}$_dtrmm( 'RIGHT', uplo, 'NO TRANSPOSE', 'NON-UNIT',kb, k-1, one, &
                                 b, ldb, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), &
                                 ldb, one, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_dsyr2k( uplo, 'TRANSPOSE', k-1, kb, one,a( k, 1_${ik}$ ), lda, b( k, &
                                 1_${ik}$ ), ldb, one, a,lda )
                       call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), &
                                 ldb, one, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_dtrmm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT', kb,k-1, one, b( &
                                 k, k ), ldb, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                    end do
                 end if
              end if
           end if
           return
     end subroutine stdlib${ii}$_dsygst

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sygst( itype, uplo, n, a, lda, b, ldb, info )
     !! DSYGST: reduces a real symmetric-definite generalized eigenproblem
     !! to standard form.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
     !! B must have been previously factorized as U**T*U or L*L**T by DPOTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: a(lda,*)
           real(${rk}$), intent(in) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: k, kb, nb
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -5_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGST', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! determine the block size for this environment.
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYGST', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
           if( nb<=1_${ik}$ .or. nb>=n ) then
              ! use unblocked code
              call stdlib${ii}$_${ri}$sygs2( itype, uplo, n, a, lda, b, ldb, info )
           else
              ! use blocked code
              if( itype==1_${ik}$ ) then
                 if( upper ) then
                    ! compute inv(u**t)*a*inv(u)
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the upper triangle of a(k:n,k:n)
                       call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                       if( k+kb<=n ) then
                          call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT',kb, n-k-kb+1, &
                                    one, b( k, k ), ldb,a( k, k+kb ), lda )
                          call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( &
                                    k, k+kb ), ldb, one,a( k, k+kb ), lda )
                          call stdlib${ii}$_${ri}$syr2k( uplo, 'TRANSPOSE', n-k-kb+1, kb, -one,a( k, k+kb ), &
                                    lda, b( k, k+kb ), ldb,one, a( k+kb, k+kb ), lda )
                          call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( &
                                    k, k+kb ), ldb, one,a( k, k+kb ), lda )
                          call stdlib${ii}$_${ri}$trsm( 'RIGHT', uplo, 'NO TRANSPOSE','NON-UNIT', kb, n-k-kb+&
                                    1_${ik}$, one,b( k+kb, k+kb ), ldb, a( k, k+kb ),lda )
                       end if
                    end do
                 else
                    ! compute inv(l)*a*inv(l**t)
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the lower triangle of a(k:n,k:n)
                       call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                       if( k+kb<=n ) then
                          call stdlib${ii}$_${ri}$trsm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',n-k-kb+1, kb, &
                                    one, b( k, k ), ldb,a( k+kb, k ), lda )
                          call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(&
                                     k+kb, k ), ldb, one,a( k+kb, k ), lda )
                          call stdlib${ii}$_${ri}$syr2k( uplo, 'NO TRANSPOSE', n-k-kb+1, kb,-one, a( k+kb, k &
                                    ), lda, b( k+kb, k ),ldb, one, a( k+kb, k+kb ), lda )
                          call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(&
                                     k+kb, k ), ldb, one,a( k+kb, k ), lda )
                          call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, 'NO TRANSPOSE','NON-UNIT', n-k-kb+1, &
                                    kb, one,b( k+kb, k+kb ), ldb, a( k+kb, k ),lda )
                       end if
                    end do
                 end if
              else
                 if( upper ) then
                    ! compute u*a*u**t
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the upper triangle of a(1:k+kb-1,1:k+kb-1)
                       call stdlib${ii}$_${ri}$trmm( 'LEFT', uplo, 'NO TRANSPOSE', 'NON-UNIT',k-1, kb, one, &
                                 b, ldb, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), &
                                 ldb, one, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_${ri}$syr2k( uplo, 'NO TRANSPOSE', k-1, kb, one,a( 1_${ik}$, k ), lda, b( &
                                 1_${ik}$, k ), ldb, one, a,lda )
                       call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), &
                                 ldb, one, a( 1_${ik}$, k ), lda )
                       call stdlib${ii}$_${ri}$trmm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',k-1, kb, one, b( &
                                 k, k ), ldb, a( 1_${ik}$, k ),lda )
                       call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                    end do
                 else
                    ! compute l**t*a*l
                    do k = 1, n, nb
                       kb = min( n-k+1, nb )
                       ! update the lower triangle of a(1:k+kb-1,1:k+kb-1)
                       call stdlib${ii}$_${ri}$trmm( 'RIGHT', uplo, 'NO TRANSPOSE', 'NON-UNIT',kb, k-1, one, &
                                 b, ldb, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), &
                                 ldb, one, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_${ri}$syr2k( uplo, 'TRANSPOSE', k-1, kb, one,a( k, 1_${ik}$ ), lda, b( k, &
                                 1_${ik}$ ), ldb, one, a,lda )
                       call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), &
                                 ldb, one, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_${ri}$trmm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT', kb,k-1, one, b( &
                                 k, k ), ldb, a( k, 1_${ik}$ ), lda )
                       call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info )
                                 
                    end do
                 end if
              end if
           end if
           return
     end subroutine stdlib${ii}$_${ri}$sygst

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_ssygs2( itype, uplo, n, a, lda, b, ldb, info )
     !! SSYGS2 reduces a real symmetric-definite generalized eigenproblem
     !! to standard form.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T *A*L.
     !! B must have been previously factorized as U**T *U or L*L**T by SPOTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, n
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*)
           real(sp), intent(in) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: k
           real(sp) :: akk, bkk, ct
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -5_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSYGS2', -info )
              return
           end if
           if( itype==1_${ik}$ ) then
              if( upper ) then
                 ! compute inv(u**t)*a*inv(u)
                 do k = 1, n
                    ! update the upper triangle of a(k:n,k:n)
                    akk = a( k, k )
                    bkk = b( k, k )
                    akk = akk / bkk**2_${ik}$
                    a( k, k ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_sscal( n-k, one / bkk, a( k, k+1 ), lda )
                       ct = -half*akk
                       call stdlib${ii}$_saxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda )
                       call stdlib${ii}$_ssyr2( uplo, n-k, -one, a( k, k+1 ), lda,b( k, k+1 ), ldb, a( &
                                 k+1, k+1 ), lda )
                       call stdlib${ii}$_saxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda )
                       call stdlib${ii}$_strsv( uplo, 'TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), ldb, &
                                 a( k, k+1 ), lda )
                    end if
                 end do
              else
                 ! compute inv(l)*a*inv(l**t)
                 do k = 1, n
                    ! update the lower triangle of a(k:n,k:n)
                    akk = a( k, k )
                    bkk = b( k, k )
                    akk = akk / bkk**2_${ik}$
                    a( k, k ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_sscal( n-k, one / bkk, a( k+1, k ), 1_${ik}$ )
                       ct = -half*akk
                       call stdlib${ii}$_saxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ )
                       call stdlib${ii}$_ssyr2( uplo, n-k, -one, a( k+1, k ), 1_${ik}$,b( k+1, k ), 1_${ik}$, a( k+1, &
                                 k+1 ), lda )
                       call stdlib${ii}$_saxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ )
                       call stdlib${ii}$_strsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), &
                                 ldb, a( k+1, k ), 1_${ik}$ )
                    end if
                 end do
              end if
           else
              if( upper ) then
                 ! compute u*a*u**t
                 do k = 1, n
                    ! update the upper triangle of a(1:k,1:k)
                    akk = a( k, k )
                    bkk = b( k, k )
                    call stdlib${ii}$_strmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, b,ldb, a( 1_${ik}$, k ), 1_${ik}$ &
                              )
                    ct = half*akk
                    call stdlib${ii}$_saxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    call stdlib${ii}$_ssyr2( uplo, k-1, one, a( 1_${ik}$, k ), 1_${ik}$, b( 1_${ik}$, k ), 1_${ik}$,a, lda )
                    call stdlib${ii}$_saxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    call stdlib${ii}$_sscal( k-1, bkk, a( 1_${ik}$, k ), 1_${ik}$ )
                    a( k, k ) = akk*bkk**2_${ik}$
                 end do
              else
                 ! compute l**t *a*l
                 do k = 1, n
                    ! update the lower triangle of a(1:k,1:k)
                    akk = a( k, k )
                    bkk = b( k, k )
                    call stdlib${ii}$_strmv( uplo, 'TRANSPOSE', 'NON-UNIT', k-1, b, ldb,a( k, 1_${ik}$ ), lda )
                              
                    ct = half*akk
                    call stdlib${ii}$_saxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda )
                    call stdlib${ii}$_ssyr2( uplo, k-1, one, a( k, 1_${ik}$ ), lda, b( k, 1_${ik}$ ),ldb, a, lda )
                              
                    call stdlib${ii}$_saxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda )
                    call stdlib${ii}$_sscal( k-1, bkk, a( k, 1_${ik}$ ), lda )
                    a( k, k ) = akk*bkk**2_${ik}$
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_ssygs2

     pure module subroutine stdlib${ii}$_dsygs2( itype, uplo, n, a, lda, b, ldb, info )
     !! DSYGS2 reduces a real symmetric-definite generalized eigenproblem
     !! to standard form.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T *A*L.
     !! B must have been previously factorized as U**T *U or L*L**T by DPOTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, n
           ! Array Arguments 
           real(dp), intent(inout) :: a(lda,*)
           real(dp), intent(in) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: k
           real(dp) :: akk, bkk, ct
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -5_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGS2', -info )
              return
           end if
           if( itype==1_${ik}$ ) then
              if( upper ) then
                 ! compute inv(u**t)*a*inv(u)
                 do k = 1, n
                    ! update the upper triangle of a(k:n,k:n)
                    akk = a( k, k )
                    bkk = b( k, k )
                    akk = akk / bkk**2_${ik}$
                    a( k, k ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_dscal( n-k, one / bkk, a( k, k+1 ), lda )
                       ct = -half*akk
                       call stdlib${ii}$_daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda )
                       call stdlib${ii}$_dsyr2( uplo, n-k, -one, a( k, k+1 ), lda,b( k, k+1 ), ldb, a( &
                                 k+1, k+1 ), lda )
                       call stdlib${ii}$_daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda )
                       call stdlib${ii}$_dtrsv( uplo, 'TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), ldb, &
                                 a( k, k+1 ), lda )
                    end if
                 end do
              else
                 ! compute inv(l)*a*inv(l**t)
                 do k = 1, n
                    ! update the lower triangle of a(k:n,k:n)
                    akk = a( k, k )
                    bkk = b( k, k )
                    akk = akk / bkk**2_${ik}$
                    a( k, k ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_dscal( n-k, one / bkk, a( k+1, k ), 1_${ik}$ )
                       ct = -half*akk
                       call stdlib${ii}$_daxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ )
                       call stdlib${ii}$_dsyr2( uplo, n-k, -one, a( k+1, k ), 1_${ik}$,b( k+1, k ), 1_${ik}$, a( k+1, &
                                 k+1 ), lda )
                       call stdlib${ii}$_daxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ )
                       call stdlib${ii}$_dtrsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), &
                                 ldb, a( k+1, k ), 1_${ik}$ )
                    end if
                 end do
              end if
           else
              if( upper ) then
                 ! compute u*a*u**t
                 do k = 1, n
                    ! update the upper triangle of a(1:k,1:k)
                    akk = a( k, k )
                    bkk = b( k, k )
                    call stdlib${ii}$_dtrmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, b,ldb, a( 1_${ik}$, k ), 1_${ik}$ &
                              )
                    ct = half*akk
                    call stdlib${ii}$_daxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    call stdlib${ii}$_dsyr2( uplo, k-1, one, a( 1_${ik}$, k ), 1_${ik}$, b( 1_${ik}$, k ), 1_${ik}$,a, lda )
                    call stdlib${ii}$_daxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    call stdlib${ii}$_dscal( k-1, bkk, a( 1_${ik}$, k ), 1_${ik}$ )
                    a( k, k ) = akk*bkk**2_${ik}$
                 end do
              else
                 ! compute l**t *a*l
                 do k = 1, n
                    ! update the lower triangle of a(1:k,1:k)
                    akk = a( k, k )
                    bkk = b( k, k )
                    call stdlib${ii}$_dtrmv( uplo, 'TRANSPOSE', 'NON-UNIT', k-1, b, ldb,a( k, 1_${ik}$ ), lda )
                              
                    ct = half*akk
                    call stdlib${ii}$_daxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda )
                    call stdlib${ii}$_dsyr2( uplo, k-1, one, a( k, 1_${ik}$ ), lda, b( k, 1_${ik}$ ),ldb, a, lda )
                              
                    call stdlib${ii}$_daxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda )
                    call stdlib${ii}$_dscal( k-1, bkk, a( k, 1_${ik}$ ), lda )
                    a( k, k ) = akk*bkk**2_${ik}$
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_dsygs2

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sygs2( itype, uplo, n, a, lda, b, ldb, info )
     !! DSYGS2: reduces a real symmetric-definite generalized eigenproblem
     !! to standard form.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T *A*L.
     !! B must have been previously factorized as U**T *U or L*L**T by DPOTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, lda, ldb, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: a(lda,*)
           real(${rk}$), intent(in) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: k
           real(${rk}$) :: akk, bkk, ct
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -5_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYGS2', -info )
              return
           end if
           if( itype==1_${ik}$ ) then
              if( upper ) then
                 ! compute inv(u**t)*a*inv(u)
                 do k = 1, n
                    ! update the upper triangle of a(k:n,k:n)
                    akk = a( k, k )
                    bkk = b( k, k )
                    akk = akk / bkk**2_${ik}$
                    a( k, k ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_${ri}$scal( n-k, one / bkk, a( k, k+1 ), lda )
                       ct = -half*akk
                       call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda )
                       call stdlib${ii}$_${ri}$syr2( uplo, n-k, -one, a( k, k+1 ), lda,b( k, k+1 ), ldb, a( &
                                 k+1, k+1 ), lda )
                       call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda )
                       call stdlib${ii}$_${ri}$trsv( uplo, 'TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), ldb, &
                                 a( k, k+1 ), lda )
                    end if
                 end do
              else
                 ! compute inv(l)*a*inv(l**t)
                 do k = 1, n
                    ! update the lower triangle of a(k:n,k:n)
                    akk = a( k, k )
                    bkk = b( k, k )
                    akk = akk / bkk**2_${ik}$
                    a( k, k ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_${ri}$scal( n-k, one / bkk, a( k+1, k ), 1_${ik}$ )
                       ct = -half*akk
                       call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ )
                       call stdlib${ii}$_${ri}$syr2( uplo, n-k, -one, a( k+1, k ), 1_${ik}$,b( k+1, k ), 1_${ik}$, a( k+1, &
                                 k+1 ), lda )
                       call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ )
                       call stdlib${ii}$_${ri}$trsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), &
                                 ldb, a( k+1, k ), 1_${ik}$ )
                    end if
                 end do
              end if
           else
              if( upper ) then
                 ! compute u*a*u**t
                 do k = 1, n
                    ! update the upper triangle of a(1:k,1:k)
                    akk = a( k, k )
                    bkk = b( k, k )
                    call stdlib${ii}$_${ri}$trmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, b,ldb, a( 1_${ik}$, k ), 1_${ik}$ &
                              )
                    ct = half*akk
                    call stdlib${ii}$_${ri}$axpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$syr2( uplo, k-1, one, a( 1_${ik}$, k ), 1_${ik}$, b( 1_${ik}$, k ), 1_${ik}$,a, lda )
                    call stdlib${ii}$_${ri}$axpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$scal( k-1, bkk, a( 1_${ik}$, k ), 1_${ik}$ )
                    a( k, k ) = akk*bkk**2_${ik}$
                 end do
              else
                 ! compute l**t *a*l
                 do k = 1, n
                    ! update the lower triangle of a(1:k,1:k)
                    akk = a( k, k )
                    bkk = b( k, k )
                    call stdlib${ii}$_${ri}$trmv( uplo, 'TRANSPOSE', 'NON-UNIT', k-1, b, ldb,a( k, 1_${ik}$ ), lda )
                              
                    ct = half*akk
                    call stdlib${ii}$_${ri}$axpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda )
                    call stdlib${ii}$_${ri}$syr2( uplo, k-1, one, a( k, 1_${ik}$ ), lda, b( k, 1_${ik}$ ),ldb, a, lda )
                              
                    call stdlib${ii}$_${ri}$axpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda )
                    call stdlib${ii}$_${ri}$scal( k-1, bkk, a( k, 1_${ik}$ ), lda )
                    a( k, k ) = akk*bkk**2_${ik}$
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_${ri}$sygs2

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sspgst( itype, uplo, n, ap, bp, info )
     !! SSPGST reduces a real symmetric-definite generalized eigenproblem
     !! to standard form, using packed storage.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
     !! B must have been previously factorized as U**T*U or L*L**T by SPPTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, n
           ! Array Arguments 
           real(sp), intent(inout) :: ap(*)
           real(sp), intent(in) :: bp(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, j1, j1j1, jj, k, k1, k1k1, kk
           real(sp) :: ajj, akk, bjj, bkk, ct
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSPGST', -info )
              return
           end if
           if( itype==1_${ik}$ ) then
              if( upper ) then
                 ! compute inv(u**t)*a*inv(u)
                 ! j1 and jj are the indices of a(1,j) and a(j,j)
                 jj = 0_${ik}$
                 do j = 1, n
                    j1 = jj + 1_${ik}$
                    jj = jj + j
                    ! compute the j-th column of the upper triangle of a
                    bjj = bp( jj )
                    call stdlib${ii}$_stpsv( uplo, 'TRANSPOSE', 'NONUNIT', j, bp,ap( j1 ), 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, j-1, -one, ap, bp( j1 ), 1_${ik}$, one,ap( j1 ), 1_${ik}$ )
                    call stdlib${ii}$_sscal( j-1, one / bjj, ap( j1 ), 1_${ik}$ )
                    ap( jj ) = ( ap( jj )-stdlib${ii}$_sdot( j-1, ap( j1 ), 1_${ik}$, bp( j1 ),1_${ik}$ ) ) / &
                              bjj
                 end do
              else
                 ! compute inv(l)*a*inv(l**t)
                 ! kk and k1k1 are the indices of a(k,k) and a(k+1,k+1)
                 kk = 1_${ik}$
                 do k = 1, n
                    k1k1 = kk + n - k + 1_${ik}$
                    ! update the lower triangle of a(k:n,k:n)
                    akk = ap( kk )
                    bkk = bp( kk )
                    akk = akk / bkk**2_${ik}$
                    ap( kk ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_sscal( n-k, one / bkk, ap( kk+1 ), 1_${ik}$ )
                       ct = -half*akk
                       call stdlib${ii}$_saxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ )
                       call stdlib${ii}$_sspr2( uplo, n-k, -one, ap( kk+1 ), 1_${ik}$,bp( kk+1 ), 1_${ik}$, ap( k1k1 )&
                                  )
                       call stdlib${ii}$_saxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ )
                       call stdlib${ii}$_stpsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,bp( k1k1 ), ap( &
                                 kk+1 ), 1_${ik}$ )
                    end if
                    kk = k1k1
                 end do
              end if
           else
              if( upper ) then
                 ! compute u*a*u**t
                 ! k1 and kk are the indices of a(1,k) and a(k,k)
                 kk = 0_${ik}$
                 do k = 1, n
                    k1 = kk + 1_${ik}$
                    kk = kk + k
                    ! update the upper triangle of a(1:k,1:k)
                    akk = ap( kk )
                    bkk = bp( kk )
                    call stdlib${ii}$_stpmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, bp,ap( k1 ), 1_${ik}$ )
                              
                    ct = half*akk
                    call stdlib${ii}$_saxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ )
                    call stdlib${ii}$_sspr2( uplo, k-1, one, ap( k1 ), 1_${ik}$, bp( k1 ), 1_${ik}$,ap )
                    call stdlib${ii}$_saxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ )
                    call stdlib${ii}$_sscal( k-1, bkk, ap( k1 ), 1_${ik}$ )
                    ap( kk ) = akk*bkk**2_${ik}$
                 end do
              else
                 ! compute l**t *a*l
                 ! jj and j1j1 are the indices of a(j,j) and a(j+1,j+1)
                 jj = 1_${ik}$
                 do j = 1, n
                    j1j1 = jj + n - j + 1_${ik}$
                    ! compute the j-th column of the lower triangle of a
                    ajj = ap( jj )
                    bjj = bp( jj )
                    ap( jj ) = ajj*bjj + stdlib${ii}$_sdot( n-j, ap( jj+1 ), 1_${ik}$,bp( jj+1 ), 1_${ik}$ )
                    call stdlib${ii}$_sscal( n-j, bjj, ap( jj+1 ), 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, n-j, one, ap( j1j1 ), bp( jj+1 ), 1_${ik}$,one, ap( jj+1 ), &
                              1_${ik}$ )
                    call stdlib${ii}$_stpmv( uplo, 'TRANSPOSE', 'NON-UNIT', n-j+1,bp( jj ), ap( jj ), 1_${ik}$ &
                              )
                    jj = j1j1
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_sspgst

     pure module subroutine stdlib${ii}$_dspgst( itype, uplo, n, ap, bp, info )
     !! DSPGST reduces a real symmetric-definite generalized eigenproblem
     !! to standard form, using packed storage.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
     !! B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, n
           ! Array Arguments 
           real(dp), intent(inout) :: ap(*)
           real(dp), intent(in) :: bp(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, j1, j1j1, jj, k, k1, k1k1, kk
           real(dp) :: ajj, akk, bjj, bkk, ct
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGST', -info )
              return
           end if
           if( itype==1_${ik}$ ) then
              if( upper ) then
                 ! compute inv(u**t)*a*inv(u)
                 ! j1 and jj are the indices of a(1,j) and a(j,j)
                 jj = 0_${ik}$
                 do j = 1, n
                    j1 = jj + 1_${ik}$
                    jj = jj + j
                    ! compute the j-th column of the upper triangle of a
                    bjj = bp( jj )
                    call stdlib${ii}$_dtpsv( uplo, 'TRANSPOSE', 'NONUNIT', j, bp,ap( j1 ), 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, j-1, -one, ap, bp( j1 ), 1_${ik}$, one,ap( j1 ), 1_${ik}$ )
                    call stdlib${ii}$_dscal( j-1, one / bjj, ap( j1 ), 1_${ik}$ )
                    ap( jj ) = ( ap( jj )-stdlib${ii}$_ddot( j-1, ap( j1 ), 1_${ik}$, bp( j1 ),1_${ik}$ ) ) / &
                              bjj
                 end do
              else
                 ! compute inv(l)*a*inv(l**t)
                 ! kk and k1k1 are the indices of a(k,k) and a(k+1,k+1)
                 kk = 1_${ik}$
                 do k = 1, n
                    k1k1 = kk + n - k + 1_${ik}$
                    ! update the lower triangle of a(k:n,k:n)
                    akk = ap( kk )
                    bkk = bp( kk )
                    akk = akk / bkk**2_${ik}$
                    ap( kk ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_dscal( n-k, one / bkk, ap( kk+1 ), 1_${ik}$ )
                       ct = -half*akk
                       call stdlib${ii}$_daxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ )
                       call stdlib${ii}$_dspr2( uplo, n-k, -one, ap( kk+1 ), 1_${ik}$,bp( kk+1 ), 1_${ik}$, ap( k1k1 )&
                                  )
                       call stdlib${ii}$_daxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ )
                       call stdlib${ii}$_dtpsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,bp( k1k1 ), ap( &
                                 kk+1 ), 1_${ik}$ )
                    end if
                    kk = k1k1
                 end do
              end if
           else
              if( upper ) then
                 ! compute u*a*u**t
                 ! k1 and kk are the indices of a(1,k) and a(k,k)
                 kk = 0_${ik}$
                 do k = 1, n
                    k1 = kk + 1_${ik}$
                    kk = kk + k
                    ! update the upper triangle of a(1:k,1:k)
                    akk = ap( kk )
                    bkk = bp( kk )
                    call stdlib${ii}$_dtpmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, bp,ap( k1 ), 1_${ik}$ )
                              
                    ct = half*akk
                    call stdlib${ii}$_daxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ )
                    call stdlib${ii}$_dspr2( uplo, k-1, one, ap( k1 ), 1_${ik}$, bp( k1 ), 1_${ik}$,ap )
                    call stdlib${ii}$_daxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ )
                    call stdlib${ii}$_dscal( k-1, bkk, ap( k1 ), 1_${ik}$ )
                    ap( kk ) = akk*bkk**2_${ik}$
                 end do
              else
                 ! compute l**t *a*l
                 ! jj and j1j1 are the indices of a(j,j) and a(j+1,j+1)
                 jj = 1_${ik}$
                 do j = 1, n
                    j1j1 = jj + n - j + 1_${ik}$
                    ! compute the j-th column of the lower triangle of a
                    ajj = ap( jj )
                    bjj = bp( jj )
                    ap( jj ) = ajj*bjj + stdlib${ii}$_ddot( n-j, ap( jj+1 ), 1_${ik}$,bp( jj+1 ), 1_${ik}$ )
                    call stdlib${ii}$_dscal( n-j, bjj, ap( jj+1 ), 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, n-j, one, ap( j1j1 ), bp( jj+1 ), 1_${ik}$,one, ap( jj+1 ), &
                              1_${ik}$ )
                    call stdlib${ii}$_dtpmv( uplo, 'TRANSPOSE', 'NON-UNIT', n-j+1,bp( jj ), ap( jj ), 1_${ik}$ &
                              )
                    jj = j1j1
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_dspgst

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$spgst( itype, uplo, n, ap, bp, info )
     !! DSPGST: reduces a real symmetric-definite generalized eigenproblem
     !! to standard form, using packed storage.
     !! If ITYPE = 1, the problem is A*x = lambda*B*x,
     !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
     !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
     !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
     !! B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: itype, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: ap(*)
           real(${rk}$), intent(in) :: bp(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, j1, j1j1, jj, k, k1, k1k1, kk
           real(${rk}$) :: ajj, akk, bjj, bkk, ct
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPGST', -info )
              return
           end if
           if( itype==1_${ik}$ ) then
              if( upper ) then
                 ! compute inv(u**t)*a*inv(u)
                 ! j1 and jj are the indices of a(1,j) and a(j,j)
                 jj = 0_${ik}$
                 do j = 1, n
                    j1 = jj + 1_${ik}$
                    jj = jj + j
                    ! compute the j-th column of the upper triangle of a
                    bjj = bp( jj )
                    call stdlib${ii}$_${ri}$tpsv( uplo, 'TRANSPOSE', 'NONUNIT', j, bp,ap( j1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, j-1, -one, ap, bp( j1 ), 1_${ik}$, one,ap( j1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$scal( j-1, one / bjj, ap( j1 ), 1_${ik}$ )
                    ap( jj ) = ( ap( jj )-stdlib${ii}$_${ri}$dot( j-1, ap( j1 ), 1_${ik}$, bp( j1 ),1_${ik}$ ) ) / &
                              bjj
                 end do
              else
                 ! compute inv(l)*a*inv(l**t)
                 ! kk and k1k1 are the indices of a(k,k) and a(k+1,k+1)
                 kk = 1_${ik}$
                 do k = 1, n
                    k1k1 = kk + n - k + 1_${ik}$
                    ! update the lower triangle of a(k:n,k:n)
                    akk = ap( kk )
                    bkk = bp( kk )
                    akk = akk / bkk**2_${ik}$
                    ap( kk ) = akk
                    if( k<n ) then
                       call stdlib${ii}$_${ri}$scal( n-k, one / bkk, ap( kk+1 ), 1_${ik}$ )
                       ct = -half*akk
                       call stdlib${ii}$_${ri}$axpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ )
                       call stdlib${ii}$_${ri}$spr2( uplo, n-k, -one, ap( kk+1 ), 1_${ik}$,bp( kk+1 ), 1_${ik}$, ap( k1k1 )&
                                  )
                       call stdlib${ii}$_${ri}$axpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ )
                       call stdlib${ii}$_${ri}$tpsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,bp( k1k1 ), ap( &
                                 kk+1 ), 1_${ik}$ )
                    end if
                    kk = k1k1
                 end do
              end if
           else
              if( upper ) then
                 ! compute u*a*u**t
                 ! k1 and kk are the indices of a(1,k) and a(k,k)
                 kk = 0_${ik}$
                 do k = 1, n
                    k1 = kk + 1_${ik}$
                    kk = kk + k
                    ! update the upper triangle of a(1:k,1:k)
                    akk = ap( kk )
                    bkk = bp( kk )
                    call stdlib${ii}$_${ri}$tpmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, bp,ap( k1 ), 1_${ik}$ )
                              
                    ct = half*akk
                    call stdlib${ii}$_${ri}$axpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spr2( uplo, k-1, one, ap( k1 ), 1_${ik}$, bp( k1 ), 1_${ik}$,ap )
                    call stdlib${ii}$_${ri}$axpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$scal( k-1, bkk, ap( k1 ), 1_${ik}$ )
                    ap( kk ) = akk*bkk**2_${ik}$
                 end do
              else
                 ! compute l**t *a*l
                 ! jj and j1j1 are the indices of a(j,j) and a(j+1,j+1)
                 jj = 1_${ik}$
                 do j = 1, n
                    j1j1 = jj + n - j + 1_${ik}$
                    ! compute the j-th column of the lower triangle of a
                    ajj = ap( jj )
                    bjj = bp( jj )
                    ap( jj ) = ajj*bjj + stdlib${ii}$_${ri}$dot( n-j, ap( jj+1 ), 1_${ik}$,bp( jj+1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$scal( n-j, bjj, ap( jj+1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, n-j, one, ap( j1j1 ), bp( jj+1 ), 1_${ik}$,one, ap( jj+1 ), &
                              1_${ik}$ )
                    call stdlib${ii}$_${ri}$tpmv( uplo, 'TRANSPOSE', 'NON-UNIT', n-j+1,bp( jj ), ap( jj ), 1_${ik}$ &
                              )
                    jj = j1j1
                 end do
              end if
           end if
           return
     end subroutine stdlib${ii}$_${ri}$spgst

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_ssbgst( vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x,ldx, work, info )
     !! SSBGST reduces a real symmetric-definite banded generalized
     !! eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
     !! such that C has the same bandwidth as A.
     !! B must have been previously factorized as S**T*S by SPBSTF, using a
     !! split Cholesky factorization. A is overwritten by C = X**T*A*X, where
     !! X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
     !! bandwidth of A.
               
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo, vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldx, n
           ! Array Arguments 
           real(sp), intent(inout) :: ab(ldab,*)
           real(sp), intent(in) :: bb(ldbb,*)
           real(sp), intent(out) :: work(*), x(ldx,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: update, upper, wantx
           integer(${ik}$) :: i, i0, i1, i2, inca, j, j1, j1t, j2, j2t, k, ka1, kb1, kbt, l, m, nr, &
                     nrt, nx
           real(sp) :: bii, ra, ra1, t
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           wantx = stdlib_lsame( vect, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           ka1 = ka + 1_${ik}$
           kb1 = kb + 1_${ik}$
           info = 0_${ik}$
           if( .not.wantx .and. .not.stdlib_lsame( vect, 'N' ) ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.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( ldx<1_${ik}$ .or. wantx .and. ldx<max( 1_${ik}$, n ) ) then
              info = -11_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSBGST', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           inca = ldab*ka1
           ! initialize x to the unit matrix, if needed
           if( wantx )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, x, ldx )
           ! set m to the splitting point m. it must be the same value as is
           ! used in stdlib${ii}$_spbstf. the chosen value allows the arrays work and rwork
           ! to be of dimension (n).
           m = ( n+kb ) / 2_${ik}$
           ! the routine works in two phases, corresponding to the two halves
           ! of the split cholesky factorization of b as s**t*s where
           ! s = ( u    )
               ! ( m  l )
           ! with u upper triangular of order m, and l lower triangular of
           ! order n-m. s has the same bandwidth as b.
           ! s is treated as a product of elementary matrices:
           ! s = s(m)*s(m-1)*...*s(2)*s(1)*s(m+1)*s(m+2)*...*s(n-1)*s(n)
           ! where s(i) is determined by the i-th row of s.
           ! in phase 1, the index i takes the values n, n-1, ... , m+1;
           ! in phase 2, it takes the values 1, 2, ... , m.
           ! for each value of i, the current matrix a is updated by forming
           ! inv(s(i))**t*a*inv(s(i)). this creates a triangular bulge outside
           ! the band of a. the bulge is then pushed down toward the bottom of
           ! a in phase 1, and up toward the top of a in phase 2, by applying
           ! plane rotations.
           ! there are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
           ! of them are linearly independent, so annihilating a bulge requires
           ! only 2*kb-1 plane rotations. the rotations are divided into a 1st
           ! set of kb-1 rotations, and a 2nd set of kb rotations.
           ! wherever possible, rotations are generated and applied in vector
           ! operations of length nr between the indices j1 and j2 (sometimes
           ! replaced by modified values nrt, j1t or j2t).
           ! the cosines and sines of the rotations are stored in the array
           ! work. the cosines of the 1st set of rotations are stored in
           ! elements n+2:n+m-kb-1 and the sines of the 1st set in elements
           ! 2:m-kb-1; the cosines of the 2nd set are stored in elements
           ! n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
           ! the bulges are not formed explicitly; nonzero elements outside the
           ! band are created only when they are required for generating new
           ! rotations; they are stored in the array work, in positions where
           ! they are later overwritten by the sines of the rotations which
           ! annihilate them.
           ! **************************** phase 1 *****************************
           ! the logical structure of this phase is:
           ! update = .true.
           ! do i = n, m + 1, -1
              ! use s(i) to update a and create a new bulge
              ! apply rotations to push all bulges ka positions downward
           ! end do
           ! update = .false.
           ! do i = m + ka + 1, n - 1
              ! apply rotations to push all bulges ka positions downward
           ! end do
           ! to avoid duplicating code, the two loops are merged.
           update = .true.
           i = n + 1_${ik}$
           10 continue
           if( update ) then
              i = i - 1_${ik}$
              kbt = min( kb, i-1 )
              i0 = i - 1_${ik}$
              i1 = min( n, i+ka )
              i2 = i - kbt + ka1
              if( i<m+1 ) then
                 update = .false.
                 i = i + 1_${ik}$
                 i0 = m
                 if( ka==0 )go to 480
                 go to 10
              end if
           else
              i = i + ka
              if( i>n-1 )go to 480
           end if
           if( upper ) then
              ! transform a, working with the upper triangle
              if( update ) then
                 ! form  inv(s(i))**t * a * inv(s(i))
                 bii = bb( kb1, i )
                 do j = i, i1
                    ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
                 end do
                 do j = max( 1, i-ka ), i
                    ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
                 end do
                 do k = i - kbt, i - 1
                    do j = i - kbt, k
                       ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( j-i+kb1, i )*ab( k-i+ka1, i ) -bb(&
                        k-i+kb1, i )*ab( j-i+ka1, i ) +ab( ka1, i )*bb( j-i+kb1, i )*bb( k-i+kb1, &
                                  i )
                    end do
                    do j = max( 1, i-ka ), i - kbt - 1
                       ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( k-i+kb1, i )*ab( j-i+ka1, i )
                                 
                    end do
                 end do
                 do j = i, i1
                    do k = max( j-ka, i-kbt ), i - 1
                       ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( k-i+kb1, i )*ab( i-j+ka1, j )
                                 
                    end do
                 end do
                 if( wantx ) then
                    ! post-multiply x by inv(s(i))
                    call stdlib${ii}$_sscal( n-m, one / bii, x( m+1, i ), 1_${ik}$ )
                    if( kbt>0_${ik}$ )call stdlib${ii}$_sger( n-m, kbt, -one, x( m+1, i ), 1_${ik}$,bb( kb1-kbt, i ), &
                              1_${ik}$, x( m+1, i-kbt ), ldx )
                 end if
                 ! store a(i,i1) in ra1 for use in next loop over k
                 ra1 = ab( i-i1+ka1, i1 )
              end if
              ! generate and apply vectors of rotations to chase all the
              ! existing bulges ka positions down toward the bottom of the
              ! band
              loop_130: do k = 1, kb - 1
                 if( update ) then
                    ! determine the rotations which would annihilate the bulge
                    ! which has in theory just been created
                    if( i-k+ka<n .and. i-k>1_${ik}$ ) then
                       ! generate rotation to annihilate a(i,i-k+ka+1)
                       call stdlib${ii}$_slartg( ab( k+1, i-k+ka ), ra1,work( n+i-k+ka-m ), work( i-k+&
                                 ka-m ),ra )
                       ! create nonzero element a(i-k,i-k+ka+1) outside the
                       ! band and store it in work(i-k)
                       t = -bb( kb1-k, i )*ra1
                       work( i-k ) = work( n+i-k+ka-m )*t -work( i-k+ka-m )*ab( 1_${ik}$, i-k+ka )
                                 
                       ab( 1_${ik}$, i-k+ka ) = work( i-k+ka-m )*t +work( n+i-k+ka-m )*ab( 1_${ik}$, i-k+ka )
                                 
                       ra1 = ra
                    end if
                 end if
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 if( update ) then
                    j2t = max( j2, i+2*ka-k+1 )
                 else
                    j2t = j2
                 end if
                 nrt = ( n-j2t+ka ) / ka1
                 do j = j2t, j1, ka1
                    ! create nonzero element a(j-ka,j+1) outside the band
                    ! and store it in work(j-m)
                    work( j-m ) = work( j-m )*ab( 1_${ik}$, j+1 )
                    ab( 1_${ik}$, j+1 ) = work( n+j-m )*ab( 1_${ik}$, j+1 )
                 end do
                 ! generate rotations in 1st set to annihilate elements which
                 ! have been created outside the band
                 if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( 1_${ik}$, j2t ), inca, work( j2t-m ), ka1,work( &
                           n+j2t-m ), ka1 )
                 if( nr>0_${ik}$ ) then
                    ! apply rotations in 1st set from the right
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( ka1-l, j2 ), inca,ab( ka-l, j2+1 ), inca, work(&
                                  n+j2-m ),work( j2-m ), ka1 )
                    end do
                    ! apply rotations in 1st set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),ab( ka, j2+1 ), inca, &
                              work( n+j2-m ),work( j2-m ), ka1 )
                 end if
                 ! start applying rotations in 1st set from the left
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( n-j2+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l &
                              ), inca,work( n+j2-m ), work( j2-m ), ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 1st set
                    do j = j2, j1, ka1
                       call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j-m ), &
                                 work( j-m ) )
                    end do
                 end if
              end do loop_130
              if( update ) then
                 if( i2<=n .and. kbt>0_${ik}$ ) then
                    ! create nonzero element a(i-kbt,i-kbt+ka+1) outside the
                    ! band and store it in work(i-kbt)
                    work( i-kbt ) = -bb( kb1-kbt, i )*ra1
                 end if
              end if
              loop_170: do k = kb, 1, -1
                 if( update ) then
                    j2 = i - k - 1_${ik}$ + max( 2_${ik}$, k-i0+1 )*ka1
                 else
                    j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1
                 end if
                 ! finish applying rotations in 2nd set from the left
                 do l = kb - k, 1, -1
                    nrt = ( n-j2+ka+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2-l+1 ), inca,ab( l+1, j2-l+1 ), &
                              inca, work( n+j2-ka ),work( j2-ka ), ka1 )
                 end do
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 do j = j1, j2, -ka1
                    work( j ) = work( j-ka )
                    work( n+j ) = work( n+j-ka )
                 end do
                 do j = j2, j1, ka1
                    ! create nonzero element a(j-ka,j+1) outside the band
                    ! and store it in work(j)
                    work( j ) = work( j )*ab( 1_${ik}$, j+1 )
                    ab( 1_${ik}$, j+1 ) = work( n+j )*ab( 1_${ik}$, j+1 )
                 end do
                 if( update ) then
                    if( i-k<n-ka .and. k<=kbt )work( i-k+ka ) = work( i-k )
                 end if
              end do loop_170
              loop_210: do k = kb, 1, -1
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 if( nr>0_${ik}$ ) then
                    ! generate rotations in 2nd set to annihilate elements
                    ! which have been created outside the band
                    call stdlib${ii}$_slargv( nr, ab( 1_${ik}$, j2 ), inca, work( j2 ), ka1,work( n+j2 ), ka1 )
                              
                    ! apply rotations in 2nd set from the right
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( ka1-l, j2 ), inca,ab( ka-l, j2+1 ), inca, work(&
                                  n+j2 ),work( j2 ), ka1 )
                    end do
                    ! apply rotations in 2nd set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),ab( ka, j2+1 ), inca, &
                              work( n+j2 ),work( j2 ), ka1 )
                 end if
                 ! start applying rotations in 2nd set from the left
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( n-j2+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l &
                              ), inca, work( n+j2 ),work( j2 ), ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 2nd set
                    do j = j2, j1, ka1
                       call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j ), work( &
                                 j ) )
                    end do
                 end if
              end do loop_210
              do k = 1, kb - 1
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1
                 ! finish applying rotations in 1st set from the left
                 do l = kb - k, 1, -1
                    nrt = ( n-j2+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l &
                              ), inca,work( n+j2-m ), work( j2-m ), ka1 )
                 end do
              end do
              if( kb>1_${ik}$ ) then
                 do j = n - 1, i - kb + 2*ka + 1, -1
                    work( n+j-m ) = work( n+j-ka-m )
                    work( j-m ) = work( j-ka-m )
                 end do
              end if
           else
              ! transform a, working with the lower triangle
              if( update ) then
                 ! form  inv(s(i))**t * a * inv(s(i))
                 bii = bb( 1_${ik}$, i )
                 do j = i, i1
                    ab( j-i+1, i ) = ab( j-i+1, i ) / bii
                 end do
                 do j = max( 1, i-ka ), i
                    ab( i-j+1, j ) = ab( i-j+1, j ) / bii
                 end do
                 do k = i - kbt, i - 1
                    do j = i - kbt, k
                       ab( k-j+1, j ) = ab( k-j+1, j ) -bb( i-j+1, j )*ab( i-k+1, k ) -bb( i-k+1, &
                                 k )*ab( i-j+1, j ) +ab( 1_${ik}$, i )*bb( i-j+1, j )*bb( i-k+1, k )
                    end do
                    do j = max( 1, i-ka ), i - kbt - 1
                       ab( k-j+1, j ) = ab( k-j+1, j ) -bb( i-k+1, k )*ab( i-j+1, j )
                    end do
                 end do
                 do j = i, i1
                    do k = max( j-ka, i-kbt ), i - 1
                       ab( j-k+1, k ) = ab( j-k+1, k ) -bb( i-k+1, k )*ab( j-i+1, i )
                    end do
                 end do
                 if( wantx ) then
                    ! post-multiply x by inv(s(i))
                    call stdlib${ii}$_sscal( n-m, one / bii, x( m+1, i ), 1_${ik}$ )
                    if( kbt>0_${ik}$ )call stdlib${ii}$_sger( n-m, kbt, -one, x( m+1, i ), 1_${ik}$,bb( kbt+1, i-kbt )&
                              , ldbb-1,x( m+1, i-kbt ), ldx )
                 end if
                 ! store a(i1,i) in ra1 for use in next loop over k
                 ra1 = ab( i1-i+1, i )
              end if
              ! generate and apply vectors of rotations to chase all the
              ! existing bulges ka positions down toward the bottom of the
              ! band
              loop_360: do k = 1, kb - 1
                 if( update ) then
                    ! determine the rotations which would annihilate the bulge
                    ! which has in theory just been created
                    if( i-k+ka<n .and. i-k>1_${ik}$ ) then
                       ! generate rotation to annihilate a(i-k+ka+1,i)
                       call stdlib${ii}$_slartg( ab( ka1-k, i ), ra1, work( n+i-k+ka-m ),work( i-k+ka-m &
                                 ), ra )
                       ! create nonzero element a(i-k+ka+1,i-k) outside the
                       ! band and store it in work(i-k)
                       t = -bb( k+1, i-k )*ra1
                       work( i-k ) = work( n+i-k+ka-m )*t -work( i-k+ka-m )*ab( ka1, i-k )
                       ab( ka1, i-k ) = work( i-k+ka-m )*t +work( n+i-k+ka-m )*ab( ka1, i-k )
                                 
                       ra1 = ra
                    end if
                 end if
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 if( update ) then
                    j2t = max( j2, i+2*ka-k+1 )
                 else
                    j2t = j2
                 end if
                 nrt = ( n-j2t+ka ) / ka1
                 do j = j2t, j1, ka1
                    ! create nonzero element a(j+1,j-ka) outside the band
                    ! and store it in work(j-m)
                    work( j-m ) = work( j-m )*ab( ka1, j-ka+1 )
                    ab( ka1, j-ka+1 ) = work( n+j-m )*ab( ka1, j-ka+1 )
                 end do
                 ! generate rotations in 1st set to annihilate elements which
                 ! have been created outside the band
                 if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( ka1, j2t-ka ), inca, work( j2t-m ),ka1, &
                           work( n+j2t-m ), ka1 )
                 if( nr>0_${ik}$ ) then
                    ! apply rotations in 1st set from the left
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( l+1, j2-l ), inca,ab( l+2, j2-l ), inca, work( &
                                 n+j2-m ),work( j2-m ), ka1 )
                    end do
                    ! apply rotations in 1st set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j2 ), ab( 1_${ik}$, j2+1 ), ab( 2_${ik}$, j2 ),inca, work( n+&
                              j2-m ), work( j2-m ), ka1 )
                 end if
                 ! start applying rotations in 1st set from the right
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( n-j2+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2 ), inca,ab( ka1-l, j2+1 ),&
                               inca, work( n+j2-m ),work( j2-m ), ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 1st set
                    do j = j2, j1, ka1
                       call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j-m ), &
                                 work( j-m ) )
                    end do
                 end if
              end do loop_360
              if( update ) then
                 if( i2<=n .and. kbt>0_${ik}$ ) then
                    ! create nonzero element a(i-kbt+ka+1,i-kbt) outside the
                    ! band and store it in work(i-kbt)
                    work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1
                 end if
              end if
              loop_400: do k = kb, 1, -1
                 if( update ) then
                    j2 = i - k - 1_${ik}$ + max( 2_${ik}$, k-i0+1 )*ka1
                 else
                    j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1
                 end if
                 ! finish applying rotations in 2nd set from the right
                 do l = kb - k, 1, -1
                    nrt = ( n-j2+ka+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2-ka ), inca,ab( ka1-l, j2-&
                              ka+1 ), inca,work( n+j2-ka ), work( j2-ka ), ka1 )
                 end do
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 do j = j1, j2, -ka1
                    work( j ) = work( j-ka )
                    work( n+j ) = work( n+j-ka )
                 end do
                 do j = j2, j1, ka1
                    ! create nonzero element a(j+1,j-ka) outside the band
                    ! and store it in work(j)
                    work( j ) = work( j )*ab( ka1, j-ka+1 )
                    ab( ka1, j-ka+1 ) = work( n+j )*ab( ka1, j-ka+1 )
                 end do
                 if( update ) then
                    if( i-k<n-ka .and. k<=kbt )work( i-k+ka ) = work( i-k )
                 end if
              end do loop_400
              loop_440: do k = kb, 1, -1
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 if( nr>0_${ik}$ ) then
                    ! generate rotations in 2nd set to annihilate elements
                    ! which have been created outside the band
                    call stdlib${ii}$_slargv( nr, ab( ka1, j2-ka ), inca, work( j2 ), ka1,work( n+j2 ), &
                              ka1 )
                    ! apply rotations in 2nd set from the left
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( l+1, j2-l ), inca,ab( l+2, j2-l ), inca, work( &
                                 n+j2 ),work( j2 ), ka1 )
                    end do
                    ! apply rotations in 2nd set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j2 ), ab( 1_${ik}$, j2+1 ), ab( 2_${ik}$, j2 ),inca, work( n+&
                              j2 ), work( j2 ), ka1 )
                 end if
                 ! start applying rotations in 2nd set from the right
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( n-j2+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2 ), inca,ab( ka1-l, j2+1 ),&
                               inca, work( n+j2 ),work( j2 ), ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 2nd set
                    do j = j2, j1, ka1
                       call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j ), work( &
                                 j ) )
                    end do
                 end if
              end do loop_440
              do k = 1, kb - 1
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1
                 ! finish applying rotations in 1st set from the right
                 do l = kb - k, 1, -1
                    nrt = ( n-j2+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2 ), inca,ab( ka1-l, j2+1 ),&
                               inca, work( n+j2-m ),work( j2-m ), ka1 )
                 end do
              end do
              if( kb>1_${ik}$ ) then
                 do j = n - 1, i - kb + 2*ka + 1, -1
                    work( n+j-m ) = work( n+j-ka-m )
                    work( j-m ) = work( j-ka-m )
                 end do
              end if
           end if
           go to 10
           480 continue
           ! **************************** phase 2 *****************************
           ! the logical structure of this phase is:
           ! update = .true.
           ! do i = 1, m
              ! use s(i) to update a and create a new bulge
              ! apply rotations to push all bulges ka positions upward
           ! end do
           ! update = .false.
           ! do i = m - ka - 1, 2, -1
              ! apply rotations to push all bulges ka positions upward
           ! end do
           ! to avoid duplicating code, the two loops are merged.
           update = .true.
           i = 0_${ik}$
           490 continue
           if( update ) then
              i = i + 1_${ik}$
              kbt = min( kb, m-i )
              i0 = i + 1_${ik}$
              i1 = max( 1_${ik}$, i-ka )
              i2 = i + kbt - ka1
              if( i>m ) then
                 update = .false.
                 i = i - 1_${ik}$
                 i0 = m + 1_${ik}$
                 if( ka==0 )return
                 go to 490
              end if
           else
              i = i - ka
              if( i<2 )return
           end if
           if( i<m-kbt ) then
              nx = m
           else
              nx = n
           end if
           if( upper ) then
              ! transform a, working with the upper triangle
              if( update ) then
                 ! form  inv(s(i))**t * a * inv(s(i))
                 bii = bb( kb1, i )
                 do j = i1, i
                    ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
                 end do
                 do j = i, min( n, i+ka )
                    ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
                 end do
                 do k = i + 1, i + kbt
                    do j = k, i + kbt
                       ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( i-j+kb1, j )*ab( i-k+ka1, k ) -bb(&
                        i-k+kb1, k )*ab( i-j+ka1, j ) +ab( ka1, i )*bb( i-j+kb1, j )*bb( i-k+kb1, &
                                  k )
                    end do
                    do j = i + kbt + 1, min( n, i+ka )
                       ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( i-k+kb1, k )*ab( i-j+ka1, j )
                                 
                    end do
                 end do
                 do j = i1, i
                    do k = i + 1, min( j+ka, i+kbt )
                       ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( i-k+kb1, k )*ab( j-i+ka1, i )
                                 
                    end do
                 end do
                 if( wantx ) then
                    ! post-multiply x by inv(s(i))
                    call stdlib${ii}$_sscal( nx, one / bii, x( 1_${ik}$, i ), 1_${ik}$ )
                    if( kbt>0_${ik}$ )call stdlib${ii}$_sger( nx, kbt, -one, x( 1_${ik}$, i ), 1_${ik}$, bb( kb, i+1 ),ldbb-&
                              1_${ik}$, x( 1_${ik}$, i+1 ), ldx )
                 end if
                 ! store a(i1,i) in ra1 for use in next loop over k
                 ra1 = ab( i1-i+ka1, i )
              end if
              ! generate and apply vectors of rotations to chase all the
              ! existing bulges ka positions up toward the top of the band
              loop_610: do k = 1, kb - 1
                 if( update ) then
                    ! determine the rotations which would annihilate the bulge
                    ! which has in theory just been created
                    if( i+k-ka1>0_${ik}$ .and. i+k<m ) then
                       ! generate rotation to annihilate a(i+k-ka-1,i)
                       call stdlib${ii}$_slartg( ab( k+1, i ), ra1, work( n+i+k-ka ),work( i+k-ka ), ra &
                                 )
                       ! create nonzero element a(i+k-ka-1,i+k) outside the
                       ! band and store it in work(m-kb+i+k)
                       t = -bb( kb1-k, i+k )*ra1
                       work( m-kb+i+k ) = work( n+i+k-ka )*t -work( i+k-ka )*ab( 1_${ik}$, i+k )
                       ab( 1_${ik}$, i+k ) = work( i+k-ka )*t +work( n+i+k-ka )*ab( 1_${ik}$, i+k )
                       ra1 = ra
                    end if
                 end if
                 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1
                 nr = ( j2+ka-1 ) / ka1
                 j1 = j2 - ( nr-1 )*ka1
                 if( update ) then
                    j2t = min( j2, i-2*ka+k-1 )
                 else
                    j2t = j2
                 end if
                 nrt = ( j2t+ka-1 ) / ka1
                 do j = j1, j2t, ka1
                    ! create nonzero element a(j-1,j+ka) outside the band
                    ! and store it in work(j)
                    work( j ) = work( j )*ab( 1_${ik}$, j+ka-1 )
                    ab( 1_${ik}$, j+ka-1 ) = work( n+j )*ab( 1_${ik}$, j+ka-1 )
                 end do
                 ! generate rotations in 1st set to annihilate elements which
                 ! have been created outside the band
                 if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( 1_${ik}$, j1+ka ), inca, work( j1 ), ka1,work( &
                           n+j1 ), ka1 )
                 if( nr>0_${ik}$ ) then
                    ! apply rotations in 1st set from the left
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( ka1-l, j1+l ), inca,ab( ka-l, j1+l ), inca, &
                                 work( n+j1 ),work( j1 ), ka1 )
                    end do
                    ! apply rotations in 1st set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),ab( ka, j1 ), inca, &
                              work( n+j1 ),work( j1 ), ka1 )
                 end if
                 ! start applying rotations in 1st set from the right
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( j2+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t ), inca,ab( l+1, j1t-1 ), inca,&
                               work( n+j1t ),work( j1t ), ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 1st set
                    do j = j1, j2, ka1
                       call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+j ), work( j ) )
                                 
                    end do
                 end if
              end do loop_610
              if( update ) then
                 if( i2>0_${ik}$ .and. kbt>0_${ik}$ ) then
                    ! create nonzero element a(i+kbt-ka-1,i+kbt) outside the
                    ! band and store it in work(m-kb+i+kbt)
                    work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1
                 end if
              end if
              loop_650: do k = kb, 1, -1
                 if( update ) then
                    j2 = i + k + 1_${ik}$ - max( 2_${ik}$, k+i0-m )*ka1
                 else
                    j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1
                 end if
                 ! finish applying rotations in 2nd set from the right
                 do l = kb - k, 1, -1
                    nrt = ( j2+ka+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t+ka ), inca,ab( l+1, j1t+ka-1 ),&
                               inca,work( n+m-kb+j1t+ka ),work( m-kb+j1t+ka ), ka1 )
                 end do
                 nr = ( j2+ka-1 ) / ka1
                 j1 = j2 - ( nr-1 )*ka1
                 do j = j1, j2, ka1
                    work( m-kb+j ) = work( m-kb+j+ka )
                    work( n+m-kb+j ) = work( n+m-kb+j+ka )
                 end do
                 do j = j1, j2, ka1
                    ! create nonzero element a(j-1,j+ka) outside the band
                    ! and store it in work(m-kb+j)
                    work( m-kb+j ) = work( m-kb+j )*ab( 1_${ik}$, j+ka-1 )
                    ab( 1_${ik}$, j+ka-1 ) = work( n+m-kb+j )*ab( 1_${ik}$, j+ka-1 )
                 end do
                 if( update ) then
                    if( i+k>ka1 .and. k<=kbt )work( m-kb+i+k-ka ) = work( m-kb+i+k )
                 end if
              end do loop_650
              loop_690: do k = kb, 1, -1
                 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1
                 nr = ( j2+ka-1 ) / ka1
                 j1 = j2 - ( nr-1 )*ka1
                 if( nr>0_${ik}$ ) then
                    ! generate rotations in 2nd set to annihilate elements
                    ! which have been created outside the band
                    call stdlib${ii}$_slargv( nr, ab( 1_${ik}$, j1+ka ), inca, work( m-kb+j1 ),ka1, work( n+m-&
                              kb+j1 ), ka1 )
                    ! apply rotations in 2nd set from the left
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( ka1-l, j1+l ), inca,ab( ka-l, j1+l ), inca,&
                                 work( n+m-kb+j1 ), work( m-kb+j1 ), ka1 )
                    end do
                    ! apply rotations in 2nd set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),ab( ka, j1 ), inca, &
                              work( n+m-kb+j1 ),work( m-kb+j1 ), ka1 )
                 end if
                 ! start applying rotations in 2nd set from the right
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( j2+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t ), inca,ab( l+1, j1t-1 ), inca,&
                              work( n+m-kb+j1t ), work( m-kb+j1t ),ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 2nd set
                    do j = j1, j2, ka1
                       call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+m-kb+j ), work( &
                                 m-kb+j ) )
                    end do
                 end if
              end do loop_690
              do k = 1, kb - 1
                 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1
                 ! finish applying rotations in 1st set from the right
                 do l = kb - k, 1, -1
                    nrt = ( j2+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t ), inca,ab( l+1, j1t-1 ), inca,&
                               work( n+j1t ),work( j1t ), ka1 )
                 end do
              end do
              if( kb>1_${ik}$ ) then
                 do j = 2, min( i+kb, m ) - 2*ka - 1
                    work( n+j ) = work( n+j+ka )
                    work( j ) = work( j+ka )
                 end do
              end if
           else
              ! transform a, working with the lower triangle
              if( update ) then
                 ! form  inv(s(i))**t * a * inv(s(i))
                 bii = bb( 1_${ik}$, i )
                 do j = i1, i
                    ab( i-j+1, j ) = ab( i-j+1, j ) / bii
                 end do
                 do j = i, min( n, i+ka )
                    ab( j-i+1, i ) = ab( j-i+1, i ) / bii
                 end do
                 do k = i + 1, i + kbt
                    do j = k, i + kbt
                       ab( j-k+1, k ) = ab( j-k+1, k ) -bb( j-i+1, i )*ab( k-i+1, i ) -bb( k-i+1, &
                                 i )*ab( j-i+1, i ) +ab( 1_${ik}$, i )*bb( j-i+1, i )*bb( k-i+1, i )
                    end do
                    do j = i + kbt + 1, min( n, i+ka )
                       ab( j-k+1, k ) = ab( j-k+1, k ) -bb( k-i+1, i )*ab( j-i+1, i )
                    end do
                 end do
                 do j = i1, i
                    do k = i + 1, min( j+ka, i+kbt )
                       ab( k-j+1, j ) = ab( k-j+1, j ) -bb( k-i+1, i )*ab( i-j+1, j )
                    end do
                 end do
                 if( wantx ) then
                    ! post-multiply x by inv(s(i))
                    call stdlib${ii}$_sscal( nx, one / bii, x( 1_${ik}$, i ), 1_${ik}$ )
                    if( kbt>0_${ik}$ )call stdlib${ii}$_sger( nx, kbt, -one, x( 1_${ik}$, i ), 1_${ik}$, bb( 2_${ik}$, i ), 1_${ik}$,x( 1_${ik}$, &
                              i+1 ), ldx )
                 end if
                 ! store a(i,i1) in ra1 for use in next loop over k
                 ra1 = ab( i-i1+1, i1 )
              end if
              ! generate and apply vectors of rotations to chase all the
              ! existing bulges ka positions up toward the top of the band
              loop_840: do k = 1, kb - 1
                 if( update ) then
                    ! determine the rotations which would annihilate the bulge
                    ! which has in theory just been created
                    if( i+k-ka1>0_${ik}$ .and. i+k<m ) then
                       ! generate rotation to annihilate a(i,i+k-ka-1)
                       call stdlib${ii}$_slartg( ab( ka1-k, i+k-ka ), ra1,work( n+i+k-ka ), work( i+k-&
                                 ka ), ra )
                       ! create nonzero element a(i+k,i+k-ka-1) outside the
                       ! band and store it in work(m-kb+i+k)
                       t = -bb( k+1, i )*ra1
                       work( m-kb+i+k ) = work( n+i+k-ka )*t -work( i+k-ka )*ab( ka1, i+k-ka )
                                 
                       ab( ka1, i+k-ka ) = work( i+k-ka )*t +work( n+i+k-ka )*ab( ka1, i+k-ka )
                                 
                       ra1 = ra
                    end if
                 end if
                 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1
                 nr = ( j2+ka-1 ) / ka1
                 j1 = j2 - ( nr-1 )*ka1
                 if( update ) then
                    j2t = min( j2, i-2*ka+k-1 )
                 else
                    j2t = j2
                 end if
                 nrt = ( j2t+ka-1 ) / ka1
                 do j = j1, j2t, ka1
                    ! create nonzero element a(j+ka,j-1) outside the band
                    ! and store it in work(j)
                    work( j ) = work( j )*ab( ka1, j-1 )
                    ab( ka1, j-1 ) = work( n+j )*ab( ka1, j-1 )
                 end do
                 ! generate rotations in 1st set to annihilate elements which
                 ! have been created outside the band
                 if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( ka1, j1 ), inca, work( j1 ), ka1,work( n+&
                           j1 ), ka1 )
                 if( nr>0_${ik}$ ) then
                    ! apply rotations in 1st set from the right
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),inca, work( n+&
                                 j1 ), work( j1 ), ka1 )
                    end do
                    ! apply rotations in 1st set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j1 ), ab( 1_${ik}$, j1-1 ),ab( 2_${ik}$, j1-1 ), inca, work( &
                              n+j1 ),work( j1 ), ka1 )
                 end if
                 ! start applying rotations in 1st set from the left
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( j2+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,ab( ka1-l, &
                              j1t-ka1+l ), inca,work( n+j1t ), work( j1t ), ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 1st set
                    do j = j1, j2, ka1
                       call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+j ), work( j ) )
                                 
                    end do
                 end if
              end do loop_840
              if( update ) then
                 if( i2>0_${ik}$ .and. kbt>0_${ik}$ ) then
                    ! create nonzero element a(i+kbt,i+kbt-ka-1) outside the
                    ! band and store it in work(m-kb+i+kbt)
                    work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1
                 end if
              end if
              loop_880: do k = kb, 1, -1
                 if( update ) then
                    j2 = i + k + 1_${ik}$ - max( 2_${ik}$, k+i0-m )*ka1
                 else
                    j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1
                 end if
                 ! finish applying rotations in 2nd set from the left
                 do l = kb - k, 1, -1
                    nrt = ( j2+ka+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,ab( ka1-l, &
                              j1t+l-1 ), inca,work( n+m-kb+j1t+ka ),work( m-kb+j1t+ka ), ka1 )
                 end do
                 nr = ( j2+ka-1 ) / ka1
                 j1 = j2 - ( nr-1 )*ka1
                 do j = j1, j2, ka1
                    work( m-kb+j ) = work( m-kb+j+ka )
                    work( n+m-kb+j ) = work( n+m-kb+j+ka )
                 end do
                 do j = j1, j2, ka1
                    ! create nonzero element a(j+ka,j-1) outside the band
                    ! and store it in work(m-kb+j)
                    work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 )
                    ab( ka1, j-1 ) = work( n+m-kb+j )*ab( ka1, j-1 )
                 end do
                 if( update ) then
                    if( i+k>ka1 .and. k<=kbt )work( m-kb+i+k-ka ) = work( m-kb+i+k )
                 end if
              end do loop_880
              loop_920: do k = kb, 1, -1
                 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1
                 nr = ( j2+ka-1 ) / ka1
                 j1 = j2 - ( nr-1 )*ka1
                 if( nr>0_${ik}$ ) then
                    ! generate rotations in 2nd set to annihilate elements
                    ! which have been created outside the band
                    call stdlib${ii}$_slargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),ka1, work( n+m-&
                              kb+j1 ), ka1 )
                    ! apply rotations in 2nd set from the right
                    do l = 1, ka - 1
                       call stdlib${ii}$_slartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),inca, work( n+&
                                 m-kb+j1 ), work( m-kb+j1 ),ka1 )
                    end do
                    ! apply rotations in 2nd set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j1 ), ab( 1_${ik}$, j1-1 ),ab( 2_${ik}$, j1-1 ), inca, work( &
                              n+m-kb+j1 ),work( m-kb+j1 ), ka1 )
                 end if
                 ! start applying rotations in 2nd set from the left
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( j2+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,ab( ka1-l, &
                              j1t-ka1+l ), inca,work( n+m-kb+j1t ), work( m-kb+j1t ),ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 2nd set
                    do j = j1, j2, ka1
                       call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+m-kb+j ), work( &
                                 m-kb+j ) )
                    end do
                 end if
              end do loop_920
              do k = 1, kb - 1
                 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1
                 ! finish applying rotations in 1st set from the left
                 do l = kb - k, 1, -1
                    nrt = ( j2+l-1 ) / ka1
                    j1t = j2 - ( nrt-1 )*ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,ab( ka1-l, &
                              j1t-ka1+l ), inca,work( n+j1t ), work( j1t ), ka1 )
                 end do
              end do
              if( kb>1_${ik}$ ) then
                 do j = 2, min( i+kb, m ) - 2*ka - 1
                    work( n+j ) = work( n+j+ka )
                    work( j ) = work( j+ka )
                 end do
              end if
           end if
           go to 490
     end subroutine stdlib${ii}$_ssbgst

     pure module subroutine stdlib${ii}$_dsbgst( vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x,ldx, work, info )
     !! DSBGST reduces a real symmetric-definite banded generalized
     !! eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
     !! such that C has the same bandwidth as A.
     !! B must have been previously factorized as S**T*S by DPBSTF, using a
     !! split Cholesky factorization. A is overwritten by C = X**T*A*X, where
     !! X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
     !! bandwidth of A.
               
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: uplo, vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldx, n
           ! Array Arguments 
           real(dp), intent(inout) :: ab(ldab,*)
           real(dp), intent(in) :: bb(ldbb,*)
           real(dp), intent(out) :: work(*), x(ldx,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: update, upper, wantx
           integer(${ik}$) :: i, i0, i1, i2, inca, j, j1, j1t, j2, j2t, k, ka1, kb1, kbt, l, m, nr, &
                     nrt, nx
           real(dp) :: bii, ra, ra1, t
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           wantx = stdlib_lsame( vect, 'V' )
           upper = stdlib_lsame( uplo, 'U' )
           ka1 = ka + 1_${ik}$
           kb1 = kb + 1_${ik}$
           info = 0_${ik}$
           if( .not.wantx .and. .not.stdlib_lsame( vect, 'N' ) ) then
              info = -1_${ik}$
           else if( .not.upper .and. .not.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( ldx<1_${ik}$ .or. wantx .and. ldx<max( 1_${ik}$, n ) ) then
              info = -11_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSBGST', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           inca = ldab*ka1
           ! initialize x to the unit matrix, if needed
           if( wantx )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, x, ldx )
           ! set m to the splitting point m. it must be the same value as is
           ! used in stdlib${ii}$_dpbstf. the chosen value allows the arrays work and rwork
           ! to be of dimension (n).
           m = ( n+kb ) / 2_${ik}$
           ! the routine works in two phases, corresponding to the two halves
           ! of the split cholesky factorization of b as s**t*s where
           ! s = ( u    )
               ! ( m  l )
           ! with u upper triangular of order m, and l lower triangular of
           ! order n-m. s has the same bandwidth as b.
           ! s is treated as a product of elementary matrices:
           ! s = s(m)*s(m-1)*...*s(2)*s(1)*s(m+1)*s(m+2)*...*s(n-1)*s(n)
           ! where s(i) is determined by the i-th row of s.
           ! in phase 1, the index i takes the values n, n-1, ... , m+1;
           ! in phase 2, it takes the values 1, 2, ... , m.
           ! for each value of i, the current matrix a is updated by forming
           ! inv(s(i))**t*a*inv(s(i)). this creates a triangular bulge outside
           ! the band of a. the bulge is then pushed down toward the bottom of
           ! a in phase 1, and up toward the top of a in phase 2, by applying
           ! plane rotations.
           ! there are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
           ! of them are linearly independent, so annihilating a bulge requires
           ! only 2*kb-1 plane rotations. the rotations are divided into a 1st
           ! set of kb-1 rotations, and a 2nd set of kb rotations.
           ! wherever possible, rotations are generated and applied in vector
           ! operations of length nr between the indices j1 and j2 (sometimes
           ! replaced by modified values nrt, j1t or j2t).
           ! the cosines and sines of the rotations are stored in the array
           ! work. the cosines of the 1st set of rotations are stored in
           ! elements n+2:n+m-kb-1 and the sines of the 1st set in elements
           ! 2:m-kb-1; the cosines of the 2nd set are stored in elements
           ! n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
           ! the bulges are not formed explicitly; nonzero elements outside the
           ! band are created only when they are required for generating new
           ! rotations; they are stored in the array work, in positions where
           ! they are later overwritten by the sines of the rotations which
           ! annihilate them.
           ! **************************** phase 1 *****************************
           ! the logical structure of this phase is:
           ! update = .true.
           ! do i = n, m + 1, -1
              ! use s(i) to update a and create a new bulge
              ! apply rotations to push all bulges ka positions downward
           ! end do
           ! update = .false.
           ! do i = m + ka + 1, n - 1
              ! apply rotations to push all bulges ka positions downward
           ! end do
           ! to avoid duplicating code, the two loops are merged.
           update = .true.
           i = n + 1_${ik}$
           10 continue
           if( update ) then
              i = i - 1_${ik}$
              kbt = min( kb, i-1 )
              i0 = i - 1_${ik}$
              i1 = min( n, i+ka )
              i2 = i - kbt + ka1
              if( i<m+1 ) then
                 update = .false.
                 i = i + 1_${ik}$
                 i0 = m
                 if( ka==0 )go to 480
                 go to 10
              end if
           else
              i = i + ka
              if( i>n-1 )go to 480
           end if
           if( upper ) then
              ! transform a, working with the upper triangle
              if( update ) then
                 ! form  inv(s(i))**t * a * inv(s(i))
                 bii = bb( kb1, i )
                 do j = i, i1
                    ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
                 end do
                 do j = max( 1, i-ka ), i
                    ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
                 end do
                 do k = i - kbt, i - 1
                    do j = i - kbt, k
                       ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( j-i+kb1, i )*ab( k-i+ka1, i ) -bb(&
                        k-i+kb1, i )*ab( j-i+ka1, i ) +ab( ka1, i )*bb( j-i+kb1, i )*bb( k-i+kb1, &
                                  i )
                    end do
                    do j = max( 1, i-ka ), i - kbt - 1
                       ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( k-i+kb1, i )*ab( j-i+ka1, i )
                                 
                    end do
                 end do
                 do j = i, i1
                    do k = max( j-ka, i-kbt ), i - 1
                       ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( k-i+kb1, i )*ab( i-j+ka1, j )
                                 
                    end do
                 end do
                 if( wantx ) then
                    ! post-multiply x by inv(s(i))
                    call stdlib${ii}$_dscal( n-m, one / bii, x( m+1, i ), 1_${ik}$ )
                    if( kbt>0_${ik}$ )call stdlib${ii}$_dger( n-m, kbt, -one, x( m+1, i ), 1_${ik}$,bb( kb1-kbt, i ), &
                              1_${ik}$, x( m+1, i-kbt ), ldx )
                 end if
                 ! store a(i,i1) in ra1 for use in next loop over k
                 ra1 = ab( i-i1+ka1, i1 )
              end if
              ! generate and apply vectors of rotations to chase all the
              ! existing bulges ka positions down toward the bottom of the
              ! band
              loop_130: do k = 1, kb - 1
                 if( update ) then
                    ! determine the rotations which would annihilate the bulge
                    ! which has in theory just been created
                    if( i-k+ka<n .and. i-k>1_${ik}$ ) then
                       ! generate rotation to annihilate a(i,i-k+ka+1)
                       call stdlib${ii}$_dlartg( ab( k+1, i-k+ka ), ra1,work( n+i-k+ka-m ), work( i-k+&
                                 ka-m ),ra )
                       ! create nonzero element a(i-k,i-k+ka+1) outside the
                       ! band and store it in work(i-k)
                       t = -bb( kb1-k, i )*ra1
                       work( i-k ) = work( n+i-k+ka-m )*t -work( i-k+ka-m )*ab( 1_${ik}$, i-k+ka )
                                 
                       ab( 1_${ik}$, i-k+ka ) = work( i-k+ka-m )*t +work( n+i-k+ka-m )*ab( 1_${ik}$, i-k+ka )
                                 
                       ra1 = ra
                    end if
                 end if
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 if( update ) then
                    j2t = max( j2, i+2*ka-k+1 )
                 else
                    j2t = j2
                 end if
                 nrt = ( n-j2t+ka ) / ka1
                 do j = j2t, j1, ka1
                    ! create nonzero element a(j-ka,j+1) outside the band
                    ! and store it in work(j-m)
                    work( j-m ) = work( j-m )*ab( 1_${ik}$, j+1 )
                    ab( 1_${ik}$, j+1 ) = work( n+j-m )*ab( 1_${ik}$, j+1 )
                 end do
                 ! generate rotations in 1st set to annihilate elements which
                 ! have been created outside the band
                 if( nrt>0_${ik}$ )call stdlib${ii}$_dlargv( nrt, ab( 1_${ik}$, j2t ), inca, work( j2t-m ), ka1,work( &
                           n+j2t-m ), ka1 )
                 if( nr>0_${ik}$ ) then
                    ! apply rotations in 1st set from the right
                    do l = 1, ka - 1
                       call stdlib${ii}$_dlartv( nr, ab( ka1-l, j2 ), inca,ab( ka-l, j2+1 ), inca, work(&
                                  n+j2-m ),work( j2-m ), ka1 )
                    end do
                    ! apply rotations in 1st set from both sides to diagonal
                    ! blocks
                    call stdlib${ii}$_dlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),ab( ka, j2+1 ), inca, &
                              work( n+j2-m ),work( j2-m ), ka1 )
                 end if
                 ! start applying rotations in 1st set from the left
                 do l = ka - 1, kb - k + 1, -1
                    nrt = ( n-j2+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l &
                              ), inca,work( n+j2-m ), work( j2-m ), ka1 )
                 end do
                 if( wantx ) then
                    ! post-multiply x by product of rotations in 1st set
                    do j = j2, j1, ka1
                       call stdlib${ii}$_drot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j-m ), &
                                 work( j-m ) )
                    end do
                 end if
              end do loop_130
              if( update ) then
                 if( i2<=n .and. kbt>0_${ik}$ ) then
                    ! create nonzero element a(i-kbt,i-kbt+ka+1) outside the
                    ! band and store it in work(i-kbt)
                    work( i-kbt ) = -bb( kb1-kbt, i )*ra1
                 end if
              end if
              loop_170: do k = kb, 1, -1
                 if( update ) then
                    j2 = i - k - 1_${ik}$ + max( 2_${ik}$, k-i0+1 )*ka1
                 else
                    j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1
                 end if
                 ! finish applying rotations in 2nd set from the left
                 do l = kb - k, 1, -1
                    nrt = ( n-j2+ka+l ) / ka1
                    if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( l, j2-l+1 ), inca,ab( l+1, j2-l+1 ), &
                              inca, work( n+j2-ka ),work( j2-ka ), ka1 )
                 end do
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 do j = j1, j2, -ka1
                    work( j ) = work( j-ka )
                    work( n+j ) = work( n+j-ka )
                 end do
                 do j = j2, j1, ka1
                    ! create nonzero element a(j-ka,j+1) outside the band
                    ! and store it in work(j)
                    work( j ) = work( j )*ab( 1_${ik}$, j+1 )
                    ab( 1_${ik}$, j+1 ) = work( n+j )*ab( 1_${ik}$, j+1 )
                 end do
                 if( update ) then
                    if( i-k<n-ka .and. k<=kbt )work( i-k+ka ) = work( i-k )
                 end if
              end do loop_170
              loop_210: do k = kb, 1, -1
                 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1
                 nr = ( n-j2+ka ) / ka1
                 j1 = j2 + ( nr-1 )*ka1
                 if( nr>0_${ik}$ ) then
                    ! generate rotations in 2nd set to annihilate elements
                    ! which have been created outside the band
                    call stdlib${ii}$_dlargv( nr, ab( 1_${ik}$, j2