stdlib_lapack_svd_comp.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_sgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info )
     !! SGEBRD reduces a general real M-by-N matrix A to upper or lower
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, lwork, m, n
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*)
           real(sp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery
           integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'SGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
           lwkopt = ( m+n )*nb
           work( 1_${ik}$ ) = real( lwkopt,KIND=sp)
           lquery = ( lwork==-1_${ik}$ )
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then
              info = -10_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SGEBRD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           minmn = min( m, n )
           if( minmn==0_${ik}$ ) then
              work( 1_${ik}$ ) = 1_${ik}$
              return
           end if
           ws = max( m, n )
           ldwrkx = m
           ldwrky = n
           if( nb>1_${ik}$ .and. nb<minmn ) then
              ! set the crossover point nx.
              nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'SGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
              ! determine when to switch from blocked to unblocked code.
              if( nx<minmn ) then
                 ws = ( m+n )*nb
                 if( lwork<ws ) then
                    ! not enough work space for the optimal nb, consider using
                    ! a smaller block size.
                    nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'SGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ )
                    if( lwork>=( m+n )*nbmin ) then
                       nb = lwork / ( m+n )
                    else
                       nb = 1_${ik}$
                       nx = minmn
                    end if
                 end if
              end if
           else
              nx = minmn
           end if
           do i = 1, minmn - nx, nb
              ! reduce rows and columns i:i+nb-1 to bidiagonal form and return
              ! the matrices x and y which are needed to update the unreduced
              ! part of the matrix
              call stdlib${ii}$_slabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), &
                        taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky )
              ! update the trailing submatrix a(i+nb:m,i+nb:n), using an update
              ! of the form  a := a - v*y**t - x*u**t
              call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, a( i+&
                        nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, one,a( i+nb, i+nb ), lda )
              call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, &
                        work( nb+1 ), ldwrkx, a( i, i+nb ), lda,one, a( i+nb, i+nb ), lda )
              ! copy diagonal and off-diagonal elements of b back into a
              if( m>=n ) then
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j, j+1 ) = e( j )
                 end do
              else
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j+1, j ) = e( j )
                 end do
              end if
           end do
           ! use unblocked code to reduce the remainder of the matrix
           call stdlib${ii}$_sgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), &
                     work, iinfo )
           work( 1_${ik}$ ) = ws
           return
     end subroutine stdlib${ii}$_sgebrd

     pure module subroutine stdlib${ii}$_dgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info )
     !! DGEBRD reduces a general real M-by-N matrix A to upper or lower
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, lwork, m, n
           ! Array Arguments 
           real(dp), intent(inout) :: a(lda,*)
           real(dp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery
           integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
           lwkopt = ( m+n )*nb
           work( 1_${ik}$ ) = real( lwkopt,KIND=dp)
           lquery = ( lwork==-1_${ik}$ )
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then
              info = -10_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGEBRD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           minmn = min( m, n )
           if( minmn==0_${ik}$ ) then
              work( 1_${ik}$ ) = 1_${ik}$
              return
           end if
           ws = max( m, n )
           ldwrkx = m
           ldwrky = n
           if( nb>1_${ik}$ .and. nb<minmn ) then
              ! set the crossover point nx.
              nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
              ! determine when to switch from blocked to unblocked code.
              if( nx<minmn ) then
                 ws = ( m+n )*nb
                 if( lwork<ws ) then
                    ! not enough work space for the optimal nb, consider using
                    ! a smaller block size.
                    nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ )
                    if( lwork>=( m+n )*nbmin ) then
                       nb = lwork / ( m+n )
                    else
                       nb = 1_${ik}$
                       nx = minmn
                    end if
                 end if
              end if
           else
              nx = minmn
           end if
           do i = 1, minmn - nx, nb
              ! reduce rows and columns i:i+nb-1 to bidiagonal form and return
              ! the matrices x and y which are needed to update the unreduced
              ! part of the matrix
              call stdlib${ii}$_dlabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), &
                        taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky )
              ! update the trailing submatrix a(i+nb:m,i+nb:n), using an update
              ! of the form  a := a - v*y**t - x*u**t
              call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, a( i+&
                        nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, one,a( i+nb, i+nb ), lda )
              call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, &
                        work( nb+1 ), ldwrkx, a( i, i+nb ), lda,one, a( i+nb, i+nb ), lda )
              ! copy diagonal and off-diagonal elements of b back into a
              if( m>=n ) then
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j, j+1 ) = e( j )
                 end do
              else
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j+1, j ) = e( j )
                 end do
              end if
           end do
           ! use unblocked code to reduce the remainder of the matrix
           call stdlib${ii}$_dgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), &
                     work, iinfo )
           work( 1_${ik}$ ) = ws
           return
     end subroutine stdlib${ii}$_dgebrd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$gebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info )
     !! DGEBRD: reduces a general real M-by-N matrix A to upper or lower
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, lwork, m, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: a(lda,*)
           real(${rk}$), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery
           integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
           lwkopt = ( m+n )*nb
           work( 1_${ik}$ ) = real( lwkopt,KIND=${rk}$)
           lquery = ( lwork==-1_${ik}$ )
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then
              info = -10_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGEBRD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           minmn = min( m, n )
           if( minmn==0_${ik}$ ) then
              work( 1_${ik}$ ) = 1_${ik}$
              return
           end if
           ws = max( m, n )
           ldwrkx = m
           ldwrky = n
           if( nb>1_${ik}$ .and. nb<minmn ) then
              ! set the crossover point nx.
              nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
              ! determine when to switch from blocked to unblocked code.
              if( nx<minmn ) then
                 ws = ( m+n )*nb
                 if( lwork<ws ) then
                    ! not enough work space for the optimal nb, consider using
                    ! a smaller block size.
                    nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ )
                    if( lwork>=( m+n )*nbmin ) then
                       nb = lwork / ( m+n )
                    else
                       nb = 1_${ik}$
                       nx = minmn
                    end if
                 end if
              end if
           else
              nx = minmn
           end if
           do i = 1, minmn - nx, nb
              ! reduce rows and columns i:i+nb-1 to bidiagonal form and return
              ! the matrices x and y which are needed to update the unreduced
              ! part of the matrix
              call stdlib${ii}$_${ri}$labrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), &
                        taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky )
              ! update the trailing submatrix a(i+nb:m,i+nb:n), using an update
              ! of the form  a := a - v*y**t - x*u**t
              call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, a( i+&
                        nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, one,a( i+nb, i+nb ), lda )
              call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, &
                        work( nb+1 ), ldwrkx, a( i, i+nb ), lda,one, a( i+nb, i+nb ), lda )
              ! copy diagonal and off-diagonal elements of b back into a
              if( m>=n ) then
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j, j+1 ) = e( j )
                 end do
              else
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j+1, j ) = e( j )
                 end do
              end if
           end do
           ! use unblocked code to reduce the remainder of the matrix
           call stdlib${ii}$_${ri}$gebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), &
                     work, iinfo )
           work( 1_${ik}$ ) = ws
           return
     end subroutine stdlib${ii}$_${ri}$gebrd

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_cgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info )
     !! CGEBRD reduces a general complex M-by-N matrix A to upper or lower
     !! bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, lwork, m, n
           ! Array Arguments 
           real(sp), intent(out) :: d(*), e(*)
           complex(sp), intent(inout) :: a(lda,*)
           complex(sp), intent(out) :: taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery
           integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'CGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
           lwkopt = ( m+n )*nb
           work( 1_${ik}$ ) = real( lwkopt,KIND=sp)
           lquery = ( lwork==-1_${ik}$ )
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then
              info = -10_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CGEBRD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           minmn = min( m, n )
           if( minmn==0_${ik}$ ) then
              work( 1_${ik}$ ) = 1_${ik}$
              return
           end if
           ws = max( m, n )
           ldwrkx = m
           ldwrky = n
           if( nb>1_${ik}$ .and. nb<minmn ) then
              ! set the crossover point nx.
              nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'CGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
              ! determine when to switch from blocked to unblocked code.
              if( nx<minmn ) then
                 ws = ( m+n )*nb
                 if( lwork<ws ) then
                    ! not enough work space for the optimal nb, consider using
                    ! a smaller block size.
                    nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'CGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ )
                    if( lwork>=( m+n )*nbmin ) then
                       nb = lwork / ( m+n )
                    else
                       nb = 1_${ik}$
                       nx = minmn
                    end if
                 end if
              end if
           else
              nx = minmn
           end if
           do i = 1, minmn - nx, nb
              ! reduce rows and columns i:i+ib-1 to bidiagonal form and return
              ! the matrices x and y which are needed to update the unreduced
              ! part of the matrix
              call stdlib${ii}$_clabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), &
                        taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky )
              ! update the trailing submatrix a(i+ib:m,i+ib:n), using
              ! an update of the form  a := a - v*y**h - x*u**h
              call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE', m-i-nb+1,n-i-nb+1, nb, -&
              cone, a( i+nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, cone,a( i+nb, i+nb ), lda )
                        
              call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -cone, &
                        work( nb+1 ), ldwrkx, a( i, i+nb ), lda,cone, a( i+nb, i+nb ), lda )
              ! copy diagonal and off-diagonal elements of b back into a
              if( m>=n ) then
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j, j+1 ) = e( j )
                 end do
              else
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j+1, j ) = e( j )
                 end do
              end if
           end do
           ! use unblocked code to reduce the remainder of the matrix
           call stdlib${ii}$_cgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), &
                     work, iinfo )
           work( 1_${ik}$ ) = ws
           return
     end subroutine stdlib${ii}$_cgebrd

     pure module subroutine stdlib${ii}$_zgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info )
     !! ZGEBRD reduces a general complex M-by-N matrix A to upper or lower
     !! bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, lwork, m, n
           ! Array Arguments 
           real(dp), intent(out) :: d(*), e(*)
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(out) :: taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery
           integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
           lwkopt = ( m+n )*nb
           work( 1_${ik}$ ) = real( lwkopt,KIND=dp)
           lquery = ( lwork==-1_${ik}$ )
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then
              info = -10_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGEBRD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           minmn = min( m, n )
           if( minmn==0_${ik}$ ) then
              work( 1_${ik}$ ) = 1_${ik}$
              return
           end if
           ws = max( m, n )
           ldwrkx = m
           ldwrky = n
           if( nb>1_${ik}$ .and. nb<minmn ) then
              ! set the crossover point nx.
              nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
              ! determine when to switch from blocked to unblocked code.
              if( nx<minmn ) then
                 ws = ( m+n )*nb
                 if( lwork<ws ) then
                    ! not enough work space for the optimal nb, consider using
                    ! a smaller block size.
                    nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ )
                    if( lwork>=( m+n )*nbmin ) then
                       nb = lwork / ( m+n )
                    else
                       nb = 1_${ik}$
                       nx = minmn
                    end if
                 end if
              end if
           else
              nx = minmn
           end if
           do i = 1, minmn - nx, nb
              ! reduce rows and columns i:i+ib-1 to bidiagonal form and return
              ! the matrices x and y which are needed to update the unreduced
              ! part of the matrix
              call stdlib${ii}$_zlabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), &
                        taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky )
              ! update the trailing submatrix a(i+ib:m,i+ib:n), using
              ! an update of the form  a := a - v*y**h - x*u**h
              call stdlib${ii}$_zgemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE', m-i-nb+1,n-i-nb+1, nb, -&
              cone, a( i+nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, cone,a( i+nb, i+nb ), lda )
                        
              call stdlib${ii}$_zgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -cone, &
                        work( nb+1 ), ldwrkx, a( i, i+nb ), lda,cone, a( i+nb, i+nb ), lda )
              ! copy diagonal and off-diagonal elements of b back into a
              if( m>=n ) then
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j, j+1 ) = e( j )
                 end do
              else
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j+1, j ) = e( j )
                 end do
              end if
           end do
           ! use unblocked code to reduce the remainder of the matrix
           call stdlib${ii}$_zgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), &
                     work, iinfo )
           work( 1_${ik}$ ) = ws
           return
     end subroutine stdlib${ii}$_zgebrd

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$gebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info )
     !! ZGEBRD: reduces a general complex M-by-N matrix A to upper or lower
     !! bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, lwork, m, n
           ! Array Arguments 
           real(${ck}$), intent(out) :: d(*), e(*)
           complex(${ck}$), intent(inout) :: a(lda,*)
           complex(${ck}$), intent(out) :: taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: lquery
           integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
           lwkopt = ( m+n )*nb
           work( 1_${ik}$ ) = real( lwkopt,KIND=${ck}$)
           lquery = ( lwork==-1_${ik}$ )
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then
              info = -10_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGEBRD', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible
           minmn = min( m, n )
           if( minmn==0_${ik}$ ) then
              work( 1_${ik}$ ) = 1_${ik}$
              return
           end if
           ws = max( m, n )
           ldwrkx = m
           ldwrky = n
           if( nb>1_${ik}$ .and. nb<minmn ) then
              ! set the crossover point nx.
              nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) )
              ! determine when to switch from blocked to unblocked code.
              if( nx<minmn ) then
                 ws = ( m+n )*nb
                 if( lwork<ws ) then
                    ! not enough work space for the optimal nb, consider using
                    ! a smaller block size.
                    nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ )
                    if( lwork>=( m+n )*nbmin ) then
                       nb = lwork / ( m+n )
                    else
                       nb = 1_${ik}$
                       nx = minmn
                    end if
                 end if
              end if
           else
              nx = minmn
           end if
           do i = 1, minmn - nx, nb
              ! reduce rows and columns i:i+ib-1 to bidiagonal form and return
              ! the matrices x and y which are needed to update the unreduced
              ! part of the matrix
              call stdlib${ii}$_${ci}$labrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), &
                        taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky )
              ! update the trailing submatrix a(i+ib:m,i+ib:n), using
              ! an update of the form  a := a - v*y**h - x*u**h
              call stdlib${ii}$_${ci}$gemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE', m-i-nb+1,n-i-nb+1, nb, -&
              cone, a( i+nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, cone,a( i+nb, i+nb ), lda )
                        
              call stdlib${ii}$_${ci}$gemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -cone, &
                        work( nb+1 ), ldwrkx, a( i, i+nb ), lda,cone, a( i+nb, i+nb ), lda )
              ! copy diagonal and off-diagonal elements of b back into a
              if( m>=n ) then
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j, j+1 ) = e( j )
                 end do
              else
                 do j = i, i + nb - 1
                    a( j, j ) = d( j )
                    a( j+1, j ) = e( j )
                 end do
              end if
           end do
           ! use unblocked code to reduce the remainder of the matrix
           call stdlib${ii}$_${ci}$gebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), &
                     work, iinfo )
           work( 1_${ik}$ ) = ws
           return
     end subroutine stdlib${ii}$_${ci}$gebrd

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sgebd2( m, n, a, lda, d, e, tauq, taup, work, info )
     !! SGEBD2 reduces a real general m by n matrix A to upper or lower
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, m, n
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*)
           real(sp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SGEBD2', -info )
              return
           end if
           if( m>=n ) then
              ! reduce to upper bidiagonal form
              do i = 1, n
                 ! generate elementary reflector h(i) to annihilate a(i+1:m,i)
                 call stdlib${ii}$_slarfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) )
                           
                 d( i ) = a( i, i )
                 a( i, i ) = one
                 ! apply h(i) to a(i:m,i+1:n) from the left
                 if( i<n )call stdlib${ii}$_slarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$, tauq( i ),a( i, i+&
                           1_${ik}$ ), lda, work )
                 a( i, i ) = d( i )
                 if( i<n ) then
                    ! generate elementary reflector g(i) to annihilate
                    ! a(i,i+2:n)
                    call stdlib${ii}$_slarfg( n-i, a( i, i+1 ), a( i, min( i+2, n ) ),lda, taup( i ) )
                              
                    e( i ) = a( i, i+1 )
                    a( i, i+1 ) = one
                    ! apply g(i) to a(i+1:m,i+1:n) from the right
                    call stdlib${ii}$_slarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 &
                              ), lda, work )
                    a( i, i+1 ) = e( i )
                 else
                    taup( i ) = zero
                 end if
              end do
           else
              ! reduce to lower bidiagonal form
              do i = 1, m
                 ! generate elementary reflector g(i) to annihilate a(i,i+1:n)
                 call stdlib${ii}$_slarfg( n-i+1, a( i, i ), a( i, min( i+1, n ) ), lda,taup( i ) )
                           
                 d( i ) = a( i, i )
                 a( i, i ) = one
                 ! apply g(i) to a(i+1:m,i:n) from the right
                 if( i<m )call stdlib${ii}$_slarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+&
                           1_${ik}$, i ), lda, work )
                 a( i, i ) = d( i )
                 if( i<m ) then
                    ! generate elementary reflector h(i) to annihilate
                    ! a(i+2:m,i)
                    call stdlib${ii}$_slarfg( m-i, a( i+1, i ), a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) )
                              
                    e( i ) = a( i+1, i )
                    a( i+1, i ) = one
                    ! apply h(i) to a(i+1:m,i+1:n) from the left
                    call stdlib${ii}$_slarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$, tauq( i ),a( i+1, i+1 ), &
                              lda, work )
                    a( i+1, i ) = e( i )
                 else
                    tauq( i ) = zero
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_sgebd2

     pure module subroutine stdlib${ii}$_dgebd2( m, n, a, lda, d, e, tauq, taup, work, info )
     !! DGEBD2 reduces a real general m by n matrix A to upper or lower
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, m, n
           ! Array Arguments 
           real(dp), intent(inout) :: a(lda,*)
           real(dp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGEBD2', -info )
              return
           end if
           if( m>=n ) then
              ! reduce to upper bidiagonal form
              do i = 1, n
                 ! generate elementary reflector h(i) to annihilate a(i+1:m,i)
                 call stdlib${ii}$_dlarfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) )
                           
                 d( i ) = a( i, i )
                 a( i, i ) = one
                 ! apply h(i) to a(i:m,i+1:n) from the left
                 if( i<n )call stdlib${ii}$_dlarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$, tauq( i ),a( i, i+&
                           1_${ik}$ ), lda, work )
                 a( i, i ) = d( i )
                 if( i<n ) then
                    ! generate elementary reflector g(i) to annihilate
                    ! a(i,i+2:n)
                    call stdlib${ii}$_dlarfg( n-i, a( i, i+1 ), a( i, min( i+2, n ) ),lda, taup( i ) )
                              
                    e( i ) = a( i, i+1 )
                    a( i, i+1 ) = one
                    ! apply g(i) to a(i+1:m,i+1:n) from the right
                    call stdlib${ii}$_dlarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 &
                              ), lda, work )
                    a( i, i+1 ) = e( i )
                 else
                    taup( i ) = zero
                 end if
              end do
           else
              ! reduce to lower bidiagonal form
              do i = 1, m
                 ! generate elementary reflector g(i) to annihilate a(i,i+1:n)
                 call stdlib${ii}$_dlarfg( n-i+1, a( i, i ), a( i, min( i+1, n ) ), lda,taup( i ) )
                           
                 d( i ) = a( i, i )
                 a( i, i ) = one
                 ! apply g(i) to a(i+1:m,i:n) from the right
                 if( i<m )call stdlib${ii}$_dlarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+&
                           1_${ik}$, i ), lda, work )
                 a( i, i ) = d( i )
                 if( i<m ) then
                    ! generate elementary reflector h(i) to annihilate
                    ! a(i+2:m,i)
                    call stdlib${ii}$_dlarfg( m-i, a( i+1, i ), a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) )
                              
                    e( i ) = a( i+1, i )
                    a( i+1, i ) = one
                    ! apply h(i) to a(i+1:m,i+1:n) from the left
                    call stdlib${ii}$_dlarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$, tauq( i ),a( i+1, i+1 ), &
                              lda, work )
                    a( i+1, i ) = e( i )
                 else
                    tauq( i ) = zero
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_dgebd2

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$gebd2( m, n, a, lda, d, e, tauq, taup, work, info )
     !! DGEBD2: reduces a real general m by n matrix A to upper or lower
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, m, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: a(lda,*)
           real(${rk}$), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGEBD2', -info )
              return
           end if
           if( m>=n ) then
              ! reduce to upper bidiagonal form
              do i = 1, n
                 ! generate elementary reflector h(i) to annihilate a(i+1:m,i)
                 call stdlib${ii}$_${ri}$larfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) )
                           
                 d( i ) = a( i, i )
                 a( i, i ) = one
                 ! apply h(i) to a(i:m,i+1:n) from the left
                 if( i<n )call stdlib${ii}$_${ri}$larf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$, tauq( i ),a( i, i+&
                           1_${ik}$ ), lda, work )
                 a( i, i ) = d( i )
                 if( i<n ) then
                    ! generate elementary reflector g(i) to annihilate
                    ! a(i,i+2:n)
                    call stdlib${ii}$_${ri}$larfg( n-i, a( i, i+1 ), a( i, min( i+2, n ) ),lda, taup( i ) )
                              
                    e( i ) = a( i, i+1 )
                    a( i, i+1 ) = one
                    ! apply g(i) to a(i+1:m,i+1:n) from the right
                    call stdlib${ii}$_${ri}$larf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 &
                              ), lda, work )
                    a( i, i+1 ) = e( i )
                 else
                    taup( i ) = zero
                 end if
              end do
           else
              ! reduce to lower bidiagonal form
              do i = 1, m
                 ! generate elementary reflector g(i) to annihilate a(i,i+1:n)
                 call stdlib${ii}$_${ri}$larfg( n-i+1, a( i, i ), a( i, min( i+1, n ) ), lda,taup( i ) )
                           
                 d( i ) = a( i, i )
                 a( i, i ) = one
                 ! apply g(i) to a(i+1:m,i:n) from the right
                 if( i<m )call stdlib${ii}$_${ri}$larf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+&
                           1_${ik}$, i ), lda, work )
                 a( i, i ) = d( i )
                 if( i<m ) then
                    ! generate elementary reflector h(i) to annihilate
                    ! a(i+2:m,i)
                    call stdlib${ii}$_${ri}$larfg( m-i, a( i+1, i ), a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) )
                              
                    e( i ) = a( i+1, i )
                    a( i+1, i ) = one
                    ! apply h(i) to a(i+1:m,i+1:n) from the left
                    call stdlib${ii}$_${ri}$larf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$, tauq( i ),a( i+1, i+1 ), &
                              lda, work )
                    a( i+1, i ) = e( i )
                 else
                    tauq( i ) = zero
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_${ri}$gebd2

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_cgebd2( m, n, a, lda, d, e, tauq, taup, work, info )
     !! CGEBD2 reduces a complex general m by n matrix A to upper or lower
     !! real bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, m, n
           ! Array Arguments 
           real(sp), intent(out) :: d(*), e(*)
           complex(sp), intent(inout) :: a(lda,*)
           complex(sp), intent(out) :: taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           complex(sp) :: alpha
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CGEBD2', -info )
              return
           end if
           if( m>=n ) then
              ! reduce to upper bidiagonal form
              do i = 1, n
                 ! generate elementary reflector h(i) to annihilate a(i+1:m,i)
                 alpha = a( i, i )
                 call stdlib${ii}$_clarfg( m-i+1, alpha, a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) )
                 d( i ) = real( alpha,KIND=sp)
                 a( i, i ) = cone
                 ! apply h(i)**h to a(i:m,i+1:n) from the left
                 if( i<n )call stdlib${ii}$_clarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$,conjg( tauq( i ) ), &
                           a( i, i+1 ), lda, work )
                 a( i, i ) = d( i )
                 if( i<n ) then
                    ! generate elementary reflector g(i) to annihilate
                    ! a(i,i+2:n)
                    call stdlib${ii}$_clacgv( n-i, a( i, i+1 ), lda )
                    alpha = a( i, i+1 )
                    call stdlib${ii}$_clarfg( n-i, alpha, a( i, min( i+2, n ) ),lda, taup( i ) )
                    e( i ) = real( alpha,KIND=sp)
                    a( i, i+1 ) = cone
                    ! apply g(i) to a(i+1:m,i+1:n) from the right
                    call stdlib${ii}$_clarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 &
                              ), lda, work )
                    call stdlib${ii}$_clacgv( n-i, a( i, i+1 ), lda )
                    a( i, i+1 ) = e( i )
                 else
                    taup( i ) = czero
                 end if
              end do
           else
              ! reduce to lower bidiagonal form
              do i = 1, m
                 ! generate elementary reflector g(i) to annihilate a(i,i+1:n)
                 call stdlib${ii}$_clacgv( n-i+1, a( i, i ), lda )
                 alpha = a( i, i )
                 call stdlib${ii}$_clarfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,taup( i ) )
                 d( i ) = real( alpha,KIND=sp)
                 a( i, i ) = cone
                 ! apply g(i) to a(i+1:m,i:n) from the right
                 if( i<m )call stdlib${ii}$_clarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+&
                           1_${ik}$, i ), lda, work )
                 call stdlib${ii}$_clacgv( n-i+1, a( i, i ), lda )
                 a( i, i ) = d( i )
                 if( i<m ) then
                    ! generate elementary reflector h(i) to annihilate
                    ! a(i+2:m,i)
                    alpha = a( i+1, i )
                    call stdlib${ii}$_clarfg( m-i, alpha, a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) )
                    e( i ) = real( alpha,KIND=sp)
                    a( i+1, i ) = cone
                    ! apply h(i)**h to a(i+1:m,i+1:n) from the left
                    call stdlib${ii}$_clarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$,conjg( tauq( i ) ), a( i+&
                              1_${ik}$, i+1 ), lda,work )
                    a( i+1, i ) = e( i )
                 else
                    tauq( i ) = czero
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_cgebd2

     pure module subroutine stdlib${ii}$_zgebd2( m, n, a, lda, d, e, tauq, taup, work, info )
     !! ZGEBD2 reduces a complex general m by n matrix A to upper or lower
     !! real bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, m, n
           ! Array Arguments 
           real(dp), intent(out) :: d(*), e(*)
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(out) :: taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           complex(dp) :: alpha
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGEBD2', -info )
              return
           end if
           if( m>=n ) then
              ! reduce to upper bidiagonal form
              do i = 1, n
                 ! generate elementary reflector h(i) to annihilate a(i+1:m,i)
                 alpha = a( i, i )
                 call stdlib${ii}$_zlarfg( m-i+1, alpha, a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) )
                 d( i ) = real( alpha,KIND=dp)
                 a( i, i ) = cone
                 ! apply h(i)**h to a(i:m,i+1:n) from the left
                 if( i<n )call stdlib${ii}$_zlarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$,conjg( tauq( i ) ), &
                           a( i, i+1 ), lda, work )
                 a( i, i ) = d( i )
                 if( i<n ) then
                    ! generate elementary reflector g(i) to annihilate
                    ! a(i,i+2:n)
                    call stdlib${ii}$_zlacgv( n-i, a( i, i+1 ), lda )
                    alpha = a( i, i+1 )
                    call stdlib${ii}$_zlarfg( n-i, alpha, a( i, min( i+2, n ) ), lda,taup( i ) )
                    e( i ) = real( alpha,KIND=dp)
                    a( i, i+1 ) = cone
                    ! apply g(i) to a(i+1:m,i+1:n) from the right
                    call stdlib${ii}$_zlarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 &
                              ), lda, work )
                    call stdlib${ii}$_zlacgv( n-i, a( i, i+1 ), lda )
                    a( i, i+1 ) = e( i )
                 else
                    taup( i ) = czero
                 end if
              end do
           else
              ! reduce to lower bidiagonal form
              do i = 1, m
                 ! generate elementary reflector g(i) to annihilate a(i,i+1:n)
                 call stdlib${ii}$_zlacgv( n-i+1, a( i, i ), lda )
                 alpha = a( i, i )
                 call stdlib${ii}$_zlarfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,taup( i ) )
                 d( i ) = real( alpha,KIND=dp)
                 a( i, i ) = cone
                 ! apply g(i) to a(i+1:m,i:n) from the right
                 if( i<m )call stdlib${ii}$_zlarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+&
                           1_${ik}$, i ), lda, work )
                 call stdlib${ii}$_zlacgv( n-i+1, a( i, i ), lda )
                 a( i, i ) = d( i )
                 if( i<m ) then
                    ! generate elementary reflector h(i) to annihilate
                    ! a(i+2:m,i)
                    alpha = a( i+1, i )
                    call stdlib${ii}$_zlarfg( m-i, alpha, a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) )
                    e( i ) = real( alpha,KIND=dp)
                    a( i+1, i ) = cone
                    ! apply h(i)**h to a(i+1:m,i+1:n) from the left
                    call stdlib${ii}$_zlarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$,conjg( tauq( i ) ), a( i+&
                              1_${ik}$, i+1 ), lda,work )
                    a( i+1, i ) = e( i )
                 else
                    tauq( i ) = czero
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_zgebd2

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$gebd2( m, n, a, lda, d, e, tauq, taup, work, info )
     !! ZGEBD2: reduces a complex general m by n matrix A to upper or lower
     !! real bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: lda, m, n
           ! Array Arguments 
           real(${ck}$), intent(out) :: d(*), e(*)
           complex(${ck}$), intent(inout) :: a(lda,*)
           complex(${ck}$), intent(out) :: taup(*), tauq(*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(${ik}$) :: i
           complex(${ck}$) :: alpha
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( m<0_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, m ) ) then
              info = -4_${ik}$
           end if
           if( info<0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGEBD2', -info )
              return
           end if
           if( m>=n ) then
              ! reduce to upper bidiagonal form
              do i = 1, n
                 ! generate elementary reflector h(i) to annihilate a(i+1:m,i)
                 alpha = a( i, i )
                 call stdlib${ii}$_${ci}$larfg( m-i+1, alpha, a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) )
                 d( i ) = real( alpha,KIND=${ck}$)
                 a( i, i ) = cone
                 ! apply h(i)**h to a(i:m,i+1:n) from the left
                 if( i<n )call stdlib${ii}$_${ci}$larf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$,conjg( tauq( i ) ), &
                           a( i, i+1 ), lda, work )
                 a( i, i ) = d( i )
                 if( i<n ) then
                    ! generate elementary reflector g(i) to annihilate
                    ! a(i,i+2:n)
                    call stdlib${ii}$_${ci}$lacgv( n-i, a( i, i+1 ), lda )
                    alpha = a( i, i+1 )
                    call stdlib${ii}$_${ci}$larfg( n-i, alpha, a( i, min( i+2, n ) ), lda,taup( i ) )
                    e( i ) = real( alpha,KIND=${ck}$)
                    a( i, i+1 ) = cone
                    ! apply g(i) to a(i+1:m,i+1:n) from the right
                    call stdlib${ii}$_${ci}$larf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 &
                              ), lda, work )
                    call stdlib${ii}$_${ci}$lacgv( n-i, a( i, i+1 ), lda )
                    a( i, i+1 ) = e( i )
                 else
                    taup( i ) = czero
                 end if
              end do
           else
              ! reduce to lower bidiagonal form
              do i = 1, m
                 ! generate elementary reflector g(i) to annihilate a(i,i+1:n)
                 call stdlib${ii}$_${ci}$lacgv( n-i+1, a( i, i ), lda )
                 alpha = a( i, i )
                 call stdlib${ii}$_${ci}$larfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,taup( i ) )
                 d( i ) = real( alpha,KIND=${ck}$)
                 a( i, i ) = cone
                 ! apply g(i) to a(i+1:m,i:n) from the right
                 if( i<m )call stdlib${ii}$_${ci}$larf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+&
                           1_${ik}$, i ), lda, work )
                 call stdlib${ii}$_${ci}$lacgv( n-i+1, a( i, i ), lda )
                 a( i, i ) = d( i )
                 if( i<m ) then
                    ! generate elementary reflector h(i) to annihilate
                    ! a(i+2:m,i)
                    alpha = a( i+1, i )
                    call stdlib${ii}$_${ci}$larfg( m-i, alpha, a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) )
                    e( i ) = real( alpha,KIND=${ck}$)
                    a( i+1, i ) = cone
                    ! apply h(i)**h to a(i+1:m,i+1:n) from the left
                    call stdlib${ii}$_${ci}$larf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$,conjg( tauq( i ) ), a( i+&
                              1_${ik}$, i+1 ), lda,work )
                    a( i+1, i ) = e( i )
                 else
                    tauq( i ) = czero
                 end if
              end do
           end if
           return
     end subroutine stdlib${ii}$_${ci}$gebd2

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, &
     !! SGBBRD reduces a real general m-by-n band matrix A to upper
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! The routine computes B, and optionally forms Q or P**T, or computes
     !! Q**T*C for a given matrix C.
               ldc, work, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc
           ! Array Arguments 
           real(sp), intent(inout) :: ab(ldab,*), c(ldc,*)
           real(sp), intent(out) :: d(*), e(*), pt(ldpt,*), q(ldq,*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: wantb, wantc, wantpt, wantq
           integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mn,&
                      mu, mu0, nr, nrt
           real(sp) :: ra, rb, rc, rs
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           wantb = stdlib_lsame( vect, 'B' )
           wantq = stdlib_lsame( vect, 'Q' ) .or. wantb
           wantpt = stdlib_lsame( vect, 'P' ) .or. wantb
           wantc = ncc>0_${ik}$
           klu1 = kl + ku + 1_${ik}$
           info = 0_${ik}$
           if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then
              info = -1_${ik}$
           else if( m<0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kl<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ku<0_${ik}$ ) then
              info = -6_${ik}$
           else if( ldab<klu1 ) then
              info = -8_${ik}$
           else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then
              info = -12_${ik}$
           else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then
              info = -14_${ik}$
           else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SGBBRD', -info )
              return
           end if
           ! initialize q and p**t to the unit matrix, if needed
           if( wantq )call stdlib${ii}$_slaset( 'FULL', m, m, zero, one, q, ldq )
           if( wantpt )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, pt, ldpt )
           ! quick return if possible.
           if( m==0 .or. n==0 )return
           minmn = min( m, n )
           if( kl+ku>1_${ik}$ ) then
              ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce
              ! first to lower bidiagonal form and then transform to upper
              ! bidiagonal
              if( ku>0_${ik}$ ) then
                 ml0 = 1_${ik}$
                 mu0 = 2_${ik}$
              else
                 ml0 = 2_${ik}$
                 mu0 = 1_${ik}$
              end if
              ! wherever possible, plane rotations are generated and applied in
              ! vector operations of length nr over the index set j1:j2:klu1.
              ! the sines of the plane rotations are stored in work(1:max(m,n))
              ! and the cosines in work(max(m,n)+1:2*max(m,n)).
              mn = max( m, n )
              klm = min( m-1, kl )
              kun = min( n-1, ku )
              kb = klm + kun
              kb1 = kb + 1_${ik}$
              inca = kb1*ldab
              nr = 0_${ik}$
              j1 = klm + 2_${ik}$
              j2 = 1_${ik}$ - kun
              loop_90: do i = 1, minmn
                 ! reduce i-th column and i-th row of matrix to bidiagonal form
                 ml = klm + 1_${ik}$
                 mu = kun + 1_${ik}$
                 loop_80: do kk = 1, kb
                    j1 = j1 + kb
                    j2 = j2 + kb
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been created below the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_slargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, &
                              work( mn+j1 ), kb1 )
                    ! apply plane rotations from the left
                    do l = 1, kb
                       if( j2-klm+l-1>n ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( &
                                 klu1-l+1, j1-klm+l-1 ), inca,work( mn+j1 ), work( j1 ), kb1 )
                    end do
                    if( ml>ml0 ) then
                       if( ml<=m-i+1 ) then
                          ! generate plane rotation to annihilate a(i+ml-1,i)
                          ! within the band, and apply rotation from the left
                          call stdlib${ii}$_slartg( ab( ku+ml-1, i ), ab( ku+ml, i ),work( mn+i+ml-1 ), &
                                    work( i+ml-1 ),ra )
                          ab( ku+ml-1, i ) = ra
                          if( i<n )call stdlib${ii}$_srot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-&
                                    1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,work( mn+i+ml-1 ), work( i+ml-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantq ) then
                       ! accumulate product of plane rotations in q
                       do j = j1, j2, kb1
                          call stdlib${ii}$_srot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,work( mn+j ), work( j &
                                    ) )
                       end do
                    end if
                    if( wantc ) then
                       ! apply plane rotations to c
                       do j = j1, j2, kb1
                          call stdlib${ii}$_srot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,work( mn+j ), &
                                    work( j ) )
                       end do
                    end if
                    if( j2+kun>n ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j-1,j+ku) above the band
                       ! and store it in work(n+1:2*n)
                       work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun )
                       ab( 1_${ik}$, j+kun ) = work( mn+j )*ab( 1_${ik}$, j+kun )
                    end do
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been generated above the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_slargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,&
                               work( mn+j1+kun ),kb1 )
                    ! apply plane rotations from the right
                    do l = 1, kb
                       if( j2+l-1>m ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+&
                                 kun ), inca,work( mn+j1+kun ), work( j1+kun ),kb1 )
                    end do
                    if( ml==ml0 .and. mu>mu0 ) then
                       if( mu<=n-i+1 ) then
                          ! generate plane rotation to annihilate a(i,i+mu-1)
                          ! within the band, and apply rotation from the right
                          call stdlib${ii}$_slartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),work( &
                                    mn+i+mu-1 ), work( i+mu-1 ),ra )
                          ab( ku-mu+3, i+mu-2 ) = ra
                          call stdlib${ii}$_srot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-&
                                    mu+3, i+mu-1 ), 1_${ik}$,work( mn+i+mu-1 ), work( i+mu-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantpt ) then
                       ! accumulate product of plane rotations in p**t
                       do j = j1, j2, kb1
                          call stdlib${ii}$_srot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, work( &
                                    mn+j+kun ),work( j+kun ) )
                       end do
                    end if
                    if( j2+kb>m ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j+kl+ku,j+ku-1) below the
                       ! band and store it in work(1:n)
                       work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
                       ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun )
                    end do
                    if( ml>ml0 ) then
                       ml = ml - 1_${ik}$
                    else
                       mu = mu - 1_${ik}$
                    end if
                 end do loop_80
              end do loop_90
           end if
           if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then
              ! a has been reduced to lower bidiagonal form
              ! transform lower bidiagonal form to upper bidiagonal by applying
              ! plane rotations from the left, storing diagonal elements in d
              ! and off-diagonal elements in e
              do i = 1, min( m-1, n )
                 call stdlib${ii}$_slartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra )
                 d( i ) = ra
                 if( i<n ) then
                    e( i ) = rs*ab( 1_${ik}$, i+1 )
                    ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 )
                 end if
                 if( wantq )call stdlib${ii}$_srot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc, rs )
                 if( wantc )call stdlib${ii}$_srot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs )
                           
              end do
              if( m<=n )d( m ) = ab( 1_${ik}$, m )
           else if( ku>0_${ik}$ ) then
              ! a has been reduced to upper bidiagonal form
              if( m<n ) then
                 ! annihilate a(m,m+1) by applying plane rotations from the
                 ! right, storing diagonal elements in d and off-diagonal
                 ! elements in e
                 rb = ab( ku, m+1 )
                 do i = m, 1, -1
                    call stdlib${ii}$_slartg( ab( ku+1, i ), rb, rc, rs, ra )
                    d( i ) = ra
                    if( i>1_${ik}$ ) then
                       rb = -rs*ab( ku, i )
                       e( i-1 ) = rc*ab( ku, i )
                    end if
                    if( wantpt )call stdlib${ii}$_srot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, rs )
                              
                 end do
              else
                 ! copy off-diagonal elements to e and diagonal elements to d
                 do i = 1, minmn - 1
                    e( i ) = ab( ku, i+1 )
                 end do
                 do i = 1, minmn
                    d( i ) = ab( ku+1, i )
                 end do
              end if
           else
              ! a is diagonal. set elements of e to zero and copy diagonal
              ! elements to d.
              do i = 1, minmn - 1
                 e( i ) = zero
              end do
              do i = 1, minmn
                 d( i ) = ab( 1_${ik}$, i )
              end do
           end if
           return
     end subroutine stdlib${ii}$_sgbbrd

     pure module subroutine stdlib${ii}$_dgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, &
     !! DGBBRD reduces a real general m-by-n band matrix A to upper
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! The routine computes B, and optionally forms Q or P**T, or computes
     !! Q**T*C for a given matrix C.
               ldc, work, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc
           ! Array Arguments 
           real(dp), intent(inout) :: ab(ldab,*), c(ldc,*)
           real(dp), intent(out) :: d(*), e(*), pt(ldpt,*), q(ldq,*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: wantb, wantc, wantpt, wantq
           integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mn,&
                      mu, mu0, nr, nrt
           real(dp) :: ra, rb, rc, rs
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           wantb = stdlib_lsame( vect, 'B' )
           wantq = stdlib_lsame( vect, 'Q' ) .or. wantb
           wantpt = stdlib_lsame( vect, 'P' ) .or. wantb
           wantc = ncc>0_${ik}$
           klu1 = kl + ku + 1_${ik}$
           info = 0_${ik}$
           if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then
              info = -1_${ik}$
           else if( m<0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kl<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ku<0_${ik}$ ) then
              info = -6_${ik}$
           else if( ldab<klu1 ) then
              info = -8_${ik}$
           else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then
              info = -12_${ik}$
           else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then
              info = -14_${ik}$
           else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGBBRD', -info )
              return
           end if
           ! initialize q and p**t to the unit matrix, if needed
           if( wantq )call stdlib${ii}$_dlaset( 'FULL', m, m, zero, one, q, ldq )
           if( wantpt )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, pt, ldpt )
           ! quick return if possible.
           if( m==0 .or. n==0 )return
           minmn = min( m, n )
           if( kl+ku>1_${ik}$ ) then
              ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce
              ! first to lower bidiagonal form and then transform to upper
              ! bidiagonal
              if( ku>0_${ik}$ ) then
                 ml0 = 1_${ik}$
                 mu0 = 2_${ik}$
              else
                 ml0 = 2_${ik}$
                 mu0 = 1_${ik}$
              end if
              ! wherever possible, plane rotations are generated and applied in
              ! vector operations of length nr over the index set j1:j2:klu1.
              ! the sines of the plane rotations are stored in work(1:max(m,n))
              ! and the cosines in work(max(m,n)+1:2*max(m,n)).
              mn = max( m, n )
              klm = min( m-1, kl )
              kun = min( n-1, ku )
              kb = klm + kun
              kb1 = kb + 1_${ik}$
              inca = kb1*ldab
              nr = 0_${ik}$
              j1 = klm + 2_${ik}$
              j2 = 1_${ik}$ - kun
              loop_90: do i = 1, minmn
                 ! reduce i-th column and i-th row of matrix to bidiagonal form
                 ml = klm + 1_${ik}$
                 mu = kun + 1_${ik}$
                 loop_80: do kk = 1, kb
                    j1 = j1 + kb
                    j2 = j2 + kb
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been created below the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_dlargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, &
                              work( mn+j1 ), kb1 )
                    ! apply plane rotations from the left
                    do l = 1, kb
                       if( j2-klm+l-1>n ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( &
                                 klu1-l+1, j1-klm+l-1 ), inca,work( mn+j1 ), work( j1 ), kb1 )
                    end do
                    if( ml>ml0 ) then
                       if( ml<=m-i+1 ) then
                          ! generate plane rotation to annihilate a(i+ml-1,i)
                          ! within the band, and apply rotation from the left
                          call stdlib${ii}$_dlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),work( mn+i+ml-1 ), &
                                    work( i+ml-1 ),ra )
                          ab( ku+ml-1, i ) = ra
                          if( i<n )call stdlib${ii}$_drot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-&
                                    1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,work( mn+i+ml-1 ), work( i+ml-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantq ) then
                       ! accumulate product of plane rotations in q
                       do j = j1, j2, kb1
                          call stdlib${ii}$_drot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,work( mn+j ), work( j &
                                    ) )
                       end do
                    end if
                    if( wantc ) then
                       ! apply plane rotations to c
                       do j = j1, j2, kb1
                          call stdlib${ii}$_drot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,work( mn+j ), &
                                    work( j ) )
                       end do
                    end if
                    if( j2+kun>n ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j-1,j+ku) above the band
                       ! and store it in work(n+1:2*n)
                       work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun )
                       ab( 1_${ik}$, j+kun ) = work( mn+j )*ab( 1_${ik}$, j+kun )
                    end do
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been generated above the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_dlargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,&
                               work( mn+j1+kun ),kb1 )
                    ! apply plane rotations from the right
                    do l = 1, kb
                       if( j2+l-1>m ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+&
                                 kun ), inca,work( mn+j1+kun ), work( j1+kun ),kb1 )
                    end do
                    if( ml==ml0 .and. mu>mu0 ) then
                       if( mu<=n-i+1 ) then
                          ! generate plane rotation to annihilate a(i,i+mu-1)
                          ! within the band, and apply rotation from the right
                          call stdlib${ii}$_dlartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),work( &
                                    mn+i+mu-1 ), work( i+mu-1 ),ra )
                          ab( ku-mu+3, i+mu-2 ) = ra
                          call stdlib${ii}$_drot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-&
                                    mu+3, i+mu-1 ), 1_${ik}$,work( mn+i+mu-1 ), work( i+mu-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantpt ) then
                       ! accumulate product of plane rotations in p**t
                       do j = j1, j2, kb1
                          call stdlib${ii}$_drot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, work( &
                                    mn+j+kun ),work( j+kun ) )
                       end do
                    end if
                    if( j2+kb>m ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j+kl+ku,j+ku-1) below the
                       ! band and store it in work(1:n)
                       work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
                       ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun )
                    end do
                    if( ml>ml0 ) then
                       ml = ml - 1_${ik}$
                    else
                       mu = mu - 1_${ik}$
                    end if
                 end do loop_80
              end do loop_90
           end if
           if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then
              ! a has been reduced to lower bidiagonal form
              ! transform lower bidiagonal form to upper bidiagonal by applying
              ! plane rotations from the left, storing diagonal elements in d
              ! and off-diagonal elements in e
              do i = 1, min( m-1, n )
                 call stdlib${ii}$_dlartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra )
                 d( i ) = ra
                 if( i<n ) then
                    e( i ) = rs*ab( 1_${ik}$, i+1 )
                    ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 )
                 end if
                 if( wantq )call stdlib${ii}$_drot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc, rs )
                 if( wantc )call stdlib${ii}$_drot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs )
                           
              end do
              if( m<=n )d( m ) = ab( 1_${ik}$, m )
           else if( ku>0_${ik}$ ) then
              ! a has been reduced to upper bidiagonal form
              if( m<n ) then
                 ! annihilate a(m,m+1) by applying plane rotations from the
                 ! right, storing diagonal elements in d and off-diagonal
                 ! elements in e
                 rb = ab( ku, m+1 )
                 do i = m, 1, -1
                    call stdlib${ii}$_dlartg( ab( ku+1, i ), rb, rc, rs, ra )
                    d( i ) = ra
                    if( i>1_${ik}$ ) then
                       rb = -rs*ab( ku, i )
                       e( i-1 ) = rc*ab( ku, i )
                    end if
                    if( wantpt )call stdlib${ii}$_drot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, rs )
                              
                 end do
              else
                 ! copy off-diagonal elements to e and diagonal elements to d
                 do i = 1, minmn - 1
                    e( i ) = ab( ku, i+1 )
                 end do
                 do i = 1, minmn
                    d( i ) = ab( ku+1, i )
                 end do
              end if
           else
              ! a is diagonal. set elements of e to zero and copy diagonal
              ! elements to d.
              do i = 1, minmn - 1
                 e( i ) = zero
              end do
              do i = 1, minmn
                 d( i ) = ab( 1_${ik}$, i )
              end do
           end if
           return
     end subroutine stdlib${ii}$_dgbbrd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$gbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, &
     !! DGBBRD: reduces a real general m-by-n band matrix A to upper
     !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
     !! The routine computes B, and optionally forms Q or P**T, or computes
     !! Q**T*C for a given matrix C.
               ldc, work, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc
           ! Array Arguments 
           real(${rk}$), intent(inout) :: ab(ldab,*), c(ldc,*)
           real(${rk}$), intent(out) :: d(*), e(*), pt(ldpt,*), q(ldq,*), work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: wantb, wantc, wantpt, wantq
           integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mn,&
                      mu, mu0, nr, nrt
           real(${rk}$) :: ra, rb, rc, rs
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           wantb = stdlib_lsame( vect, 'B' )
           wantq = stdlib_lsame( vect, 'Q' ) .or. wantb
           wantpt = stdlib_lsame( vect, 'P' ) .or. wantb
           wantc = ncc>0_${ik}$
           klu1 = kl + ku + 1_${ik}$
           info = 0_${ik}$
           if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then
              info = -1_${ik}$
           else if( m<0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kl<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ku<0_${ik}$ ) then
              info = -6_${ik}$
           else if( ldab<klu1 ) then
              info = -8_${ik}$
           else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then
              info = -12_${ik}$
           else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then
              info = -14_${ik}$
           else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGBBRD', -info )
              return
           end if
           ! initialize q and p**t to the unit matrix, if needed
           if( wantq )call stdlib${ii}$_${ri}$laset( 'FULL', m, m, zero, one, q, ldq )
           if( wantpt )call stdlib${ii}$_${ri}$laset( 'FULL', n, n, zero, one, pt, ldpt )
           ! quick return if possible.
           if( m==0 .or. n==0 )return
           minmn = min( m, n )
           if( kl+ku>1_${ik}$ ) then
              ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce
              ! first to lower bidiagonal form and then transform to upper
              ! bidiagonal
              if( ku>0_${ik}$ ) then
                 ml0 = 1_${ik}$
                 mu0 = 2_${ik}$
              else
                 ml0 = 2_${ik}$
                 mu0 = 1_${ik}$
              end if
              ! wherever possible, plane rotations are generated and applied in
              ! vector operations of length nr over the index set j1:j2:klu1.
              ! the sines of the plane rotations are stored in work(1:max(m,n))
              ! and the cosines in work(max(m,n)+1:2*max(m,n)).
              mn = max( m, n )
              klm = min( m-1, kl )
              kun = min( n-1, ku )
              kb = klm + kun
              kb1 = kb + 1_${ik}$
              inca = kb1*ldab
              nr = 0_${ik}$
              j1 = klm + 2_${ik}$
              j2 = 1_${ik}$ - kun
              loop_90: do i = 1, minmn
                 ! reduce i-th column and i-th row of matrix to bidiagonal form
                 ml = klm + 1_${ik}$
                 mu = kun + 1_${ik}$
                 loop_80: do kk = 1, kb
                    j1 = j1 + kb
                    j2 = j2 + kb
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been created below the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_${ri}$largv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, &
                              work( mn+j1 ), kb1 )
                    ! apply plane rotations from the left
                    do l = 1, kb
                       if( j2-klm+l-1>n ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_${ri}$lartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( &
                                 klu1-l+1, j1-klm+l-1 ), inca,work( mn+j1 ), work( j1 ), kb1 )
                    end do
                    if( ml>ml0 ) then
                       if( ml<=m-i+1 ) then
                          ! generate plane rotation to annihilate a(i+ml-1,i)
                          ! within the band, and apply rotation from the left
                          call stdlib${ii}$_${ri}$lartg( ab( ku+ml-1, i ), ab( ku+ml, i ),work( mn+i+ml-1 ), &
                                    work( i+ml-1 ),ra )
                          ab( ku+ml-1, i ) = ra
                          if( i<n )call stdlib${ii}$_${ri}$rot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-&
                                    1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,work( mn+i+ml-1 ), work( i+ml-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantq ) then
                       ! accumulate product of plane rotations in q
                       do j = j1, j2, kb1
                          call stdlib${ii}$_${ri}$rot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,work( mn+j ), work( j &
                                    ) )
                       end do
                    end if
                    if( wantc ) then
                       ! apply plane rotations to c
                       do j = j1, j2, kb1
                          call stdlib${ii}$_${ri}$rot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,work( mn+j ), &
                                    work( j ) )
                       end do
                    end if
                    if( j2+kun>n ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j-1,j+ku) above the band
                       ! and store it in work(n+1:2*n)
                       work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun )
                       ab( 1_${ik}$, j+kun ) = work( mn+j )*ab( 1_${ik}$, j+kun )
                    end do
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been generated above the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_${ri}$largv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,&
                               work( mn+j1+kun ),kb1 )
                    ! apply plane rotations from the right
                    do l = 1, kb
                       if( j2+l-1>m ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_${ri}$lartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+&
                                 kun ), inca,work( mn+j1+kun ), work( j1+kun ),kb1 )
                    end do
                    if( ml==ml0 .and. mu>mu0 ) then
                       if( mu<=n-i+1 ) then
                          ! generate plane rotation to annihilate a(i,i+mu-1)
                          ! within the band, and apply rotation from the right
                          call stdlib${ii}$_${ri}$lartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),work( &
                                    mn+i+mu-1 ), work( i+mu-1 ),ra )
                          ab( ku-mu+3, i+mu-2 ) = ra
                          call stdlib${ii}$_${ri}$rot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-&
                                    mu+3, i+mu-1 ), 1_${ik}$,work( mn+i+mu-1 ), work( i+mu-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantpt ) then
                       ! accumulate product of plane rotations in p**t
                       do j = j1, j2, kb1
                          call stdlib${ii}$_${ri}$rot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, work( &
                                    mn+j+kun ),work( j+kun ) )
                       end do
                    end if
                    if( j2+kb>m ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j+kl+ku,j+ku-1) below the
                       ! band and store it in work(1:n)
                       work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
                       ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun )
                    end do
                    if( ml>ml0 ) then
                       ml = ml - 1_${ik}$
                    else
                       mu = mu - 1_${ik}$
                    end if
                 end do loop_80
              end do loop_90
           end if
           if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then
              ! a has been reduced to lower bidiagonal form
              ! transform lower bidiagonal form to upper bidiagonal by applying
              ! plane rotations from the left, storing diagonal elements in d
              ! and off-diagonal elements in e
              do i = 1, min( m-1, n )
                 call stdlib${ii}$_${ri}$lartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra )
                 d( i ) = ra
                 if( i<n ) then
                    e( i ) = rs*ab( 1_${ik}$, i+1 )
                    ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 )
                 end if
                 if( wantq )call stdlib${ii}$_${ri}$rot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc, rs )
                 if( wantc )call stdlib${ii}$_${ri}$rot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs )
                           
              end do
              if( m<=n )d( m ) = ab( 1_${ik}$, m )
           else if( ku>0_${ik}$ ) then
              ! a has been reduced to upper bidiagonal form
              if( m<n ) then
                 ! annihilate a(m,m+1) by applying plane rotations from the
                 ! right, storing diagonal elements in d and off-diagonal
                 ! elements in e
                 rb = ab( ku, m+1 )
                 do i = m, 1, -1
                    call stdlib${ii}$_${ri}$lartg( ab( ku+1, i ), rb, rc, rs, ra )
                    d( i ) = ra
                    if( i>1_${ik}$ ) then
                       rb = -rs*ab( ku, i )
                       e( i-1 ) = rc*ab( ku, i )
                    end if
                    if( wantpt )call stdlib${ii}$_${ri}$rot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, rs )
                              
                 end do
              else
                 ! copy off-diagonal elements to e and diagonal elements to d
                 do i = 1, minmn - 1
                    e( i ) = ab( ku, i+1 )
                 end do
                 do i = 1, minmn
                    d( i ) = ab( ku+1, i )
                 end do
              end if
           else
              ! a is diagonal. set elements of e to zero and copy diagonal
              ! elements to d.
              do i = 1, minmn - 1
                 e( i ) = zero
              end do
              do i = 1, minmn
                 d( i ) = ab( 1_${ik}$, i )
              end do
           end if
           return
     end subroutine stdlib${ii}$_${ri}$gbbrd

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_cgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, &
     !! CGBBRD reduces a complex general m-by-n band matrix A to real upper
     !! bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! The routine computes B, and optionally forms Q or P**H, or computes
     !! Q**H*C for a given matrix C.
               ldc, work, rwork, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc
           ! Array Arguments 
           real(sp), intent(out) :: d(*), e(*), rwork(*)
           complex(sp), intent(inout) :: ab(ldab,*), c(ldc,*)
           complex(sp), intent(out) :: pt(ldpt,*), q(ldq,*), work(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           logical(lk) :: wantb, wantc, wantpt, wantq
           integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mu,&
                      mu0, nr, nrt
           real(sp) :: abst, rc
           complex(sp) :: ra, rb, rs, t
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           wantb = stdlib_lsame( vect, 'B' )
           wantq = stdlib_lsame( vect, 'Q' ) .or. wantb
           wantpt = stdlib_lsame( vect, 'P' ) .or. wantb
           wantc = ncc>0_${ik}$
           klu1 = kl + ku + 1_${ik}$
           info = 0_${ik}$
           if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then
              info = -1_${ik}$
           else if( m<0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kl<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ku<0_${ik}$ ) then
              info = -6_${ik}$
           else if( ldab<klu1 ) then
              info = -8_${ik}$
           else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then
              info = -12_${ik}$
           else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then
              info = -14_${ik}$
           else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CGBBRD', -info )
              return
           end if
           ! initialize q and p**h to the unit matrix, if needed
           if( wantq )call stdlib${ii}$_claset( 'FULL', m, m, czero, cone, q, ldq )
           if( wantpt )call stdlib${ii}$_claset( 'FULL', n, n, czero, cone, pt, ldpt )
           ! quick return if possible.
           if( m==0 .or. n==0 )return
           minmn = min( m, n )
           if( kl+ku>1_${ik}$ ) then
              ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce
              ! first to lower bidiagonal form and then transform to upper
              ! bidiagonal
              if( ku>0_${ik}$ ) then
                 ml0 = 1_${ik}$
                 mu0 = 2_${ik}$
              else
                 ml0 = 2_${ik}$
                 mu0 = 1_${ik}$
              end if
              ! wherever possible, plane rotations are generated and applied in
              ! vector operations of length nr over the index set j1:j2:klu1.
              ! the complex sines of the plane rotations are stored in work,
              ! and the real cosines in rwork.
              klm = min( m-1, kl )
              kun = min( n-1, ku )
              kb = klm + kun
              kb1 = kb + 1_${ik}$
              inca = kb1*ldab
              nr = 0_${ik}$
              j1 = klm + 2_${ik}$
              j2 = 1_${ik}$ - kun
              loop_90: do i = 1, minmn
                 ! reduce i-th column and i-th row of matrix to bidiagonal form
                 ml = klm + 1_${ik}$
                 mu = kun + 1_${ik}$
                 loop_80: do kk = 1, kb
                    j1 = j1 + kb
                    j2 = j2 + kb
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been created below the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_clargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, &
                              rwork( j1 ), kb1 )
                    ! apply plane rotations from the left
                    do l = 1, kb
                       if( j2-klm+l-1>n ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_clartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( &
                                 klu1-l+1, j1-klm+l-1 ), inca,rwork( j1 ), work( j1 ), kb1 )
                    end do
                    if( ml>ml0 ) then
                       if( ml<=m-i+1 ) then
                          ! generate plane rotation to annihilate a(i+ml-1,i)
                          ! within the band, and apply rotation from the left
                          call stdlib${ii}$_clartg( ab( ku+ml-1, i ), ab( ku+ml, i ),rwork( i+ml-1 ), &
                                    work( i+ml-1 ), ra )
                          ab( ku+ml-1, i ) = ra
                          if( i<n )call stdlib${ii}$_crot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-&
                                    1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,rwork( i+ml-1 ), work( i+ml-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantq ) then
                       ! accumulate product of plane rotations in q
                       do j = j1, j2, kb1
                          call stdlib${ii}$_crot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,rwork( j ), conjg( &
                                    work( j ) ) )
                       end do
                    end if
                    if( wantc ) then
                       ! apply plane rotations to c
                       do j = j1, j2, kb1
                          call stdlib${ii}$_crot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,rwork( j ), &
                                    work( j ) )
                       end do
                    end if
                    if( j2+kun>n ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j-1,j+ku) above the band
                       ! and store it in work(n+1:2*n)
                       work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun )
                       ab( 1_${ik}$, j+kun ) = rwork( j )*ab( 1_${ik}$, j+kun )
                    end do
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been generated above the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_clargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,&
                               rwork( j1+kun ),kb1 )
                    ! apply plane rotations from the right
                    do l = 1, kb
                       if( j2+l-1>m ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_clartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+&
                                 kun ), inca,rwork( j1+kun ), work( j1+kun ), kb1 )
                    end do
                    if( ml==ml0 .and. mu>mu0 ) then
                       if( mu<=n-i+1 ) then
                          ! generate plane rotation to annihilate a(i,i+mu-1)
                          ! within the band, and apply rotation from the right
                          call stdlib${ii}$_clartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),rwork( &
                                    i+mu-1 ), work( i+mu-1 ), ra )
                          ab( ku-mu+3, i+mu-2 ) = ra
                          call stdlib${ii}$_crot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-&
                                    mu+3, i+mu-1 ), 1_${ik}$,rwork( i+mu-1 ), work( i+mu-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantpt ) then
                       ! accumulate product of plane rotations in p**h
                       do j = j1, j2, kb1
                          call stdlib${ii}$_crot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, rwork(&
                                     j+kun ),conjg( work( j+kun ) ) )
                       end do
                    end if
                    if( j2+kb>m ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j+kl+ku,j+ku-1) below the
                       ! band and store it in work(1:n)
                       work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
                       ab( klu1, j+kun ) = rwork( j+kun )*ab( klu1, j+kun )
                    end do
                    if( ml>ml0 ) then
                       ml = ml - 1_${ik}$
                    else
                       mu = mu - 1_${ik}$
                    end if
                 end do loop_80
              end do loop_90
           end if
           if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then
              ! a has been reduced to complex lower bidiagonal form
              ! transform lower bidiagonal form to upper bidiagonal by applying
              ! plane rotations from the left, overwriting superdiagonal
              ! elements on subdiagonal elements
              do i = 1, min( m-1, n )
                 call stdlib${ii}$_clartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra )
                 ab( 1_${ik}$, i ) = ra
                 if( i<n ) then
                    ab( 2_${ik}$, i ) = rs*ab( 1_${ik}$, i+1 )
                    ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 )
                 end if
                 if( wantq )call stdlib${ii}$_crot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc,conjg( rs ) )
                           
                 if( wantc )call stdlib${ii}$_crot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs )
                           
              end do
           else
              ! a has been reduced to complex upper bidiagonal form or is
              ! diagonal
              if( ku>0_${ik}$ .and. m<n ) then
                 ! annihilate a(m,m+1) by applying plane rotations from the
                 ! right
                 rb = ab( ku, m+1 )
                 do i = m, 1, -1
                    call stdlib${ii}$_clartg( ab( ku+1, i ), rb, rc, rs, ra )
                    ab( ku+1, i ) = ra
                    if( i>1_${ik}$ ) then
                       rb = -conjg( rs )*ab( ku, i )
                       ab( ku, i ) = rc*ab( ku, i )
                    end if
                    if( wantpt )call stdlib${ii}$_crot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, &
                              conjg( rs ) )
                 end do
              end if
           end if
           ! make diagonal and superdiagonal elements real, storing them in d
           ! and e
           t = ab( ku+1, 1_${ik}$ )
           loop_120: do i = 1, minmn
              abst = abs( t )
              d( i ) = abst
              if( abst/=zero ) then
                 t = t / abst
              else
                 t = cone
              end if
              if( wantq )call stdlib${ii}$_cscal( m, t, q( 1_${ik}$, i ), 1_${ik}$ )
              if( wantc )call stdlib${ii}$_cscal( ncc, conjg( t ), c( i, 1_${ik}$ ), ldc )
              if( i<minmn ) then
                 if( ku==0_${ik}$ .and. kl==0_${ik}$ ) then
                    e( i ) = zero
                    t = ab( 1_${ik}$, i+1 )
                 else
                    if( ku==0_${ik}$ ) then
                       t = ab( 2_${ik}$, i )*conjg( t )
                    else
                       t = ab( ku, i+1 )*conjg( t )
                    end if
                    abst = abs( t )
                    e( i ) = abst
                    if( abst/=zero ) then
                       t = t / abst
                    else
                       t = cone
                    end if
                    if( wantpt )call stdlib${ii}$_cscal( n, t, pt( i+1, 1_${ik}$ ), ldpt )
                    t = ab( ku+1, i+1 )*conjg( t )
                 end if
              end if
           end do loop_120
           return
     end subroutine stdlib${ii}$_cgbbrd

     pure module subroutine stdlib${ii}$_zgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, &
     !! ZGBBRD reduces a complex general m-by-n band matrix A to real upper
     !! bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! The routine computes B, and optionally forms Q or P**H, or computes
     !! Q**H*C for a given matrix C.
               ldc, work, rwork, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc
           ! Array Arguments 
           real(dp), intent(out) :: d(*), e(*), rwork(*)
           complex(dp), intent(inout) :: ab(ldab,*), c(ldc,*)
           complex(dp), intent(out) :: pt(ldpt,*), q(ldq,*), work(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           logical(lk) :: wantb, wantc, wantpt, wantq
           integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mu,&
                      mu0, nr, nrt
           real(dp) :: abst, rc
           complex(dp) :: ra, rb, rs, t
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           wantb = stdlib_lsame( vect, 'B' )
           wantq = stdlib_lsame( vect, 'Q' ) .or. wantb
           wantpt = stdlib_lsame( vect, 'P' ) .or. wantb
           wantc = ncc>0_${ik}$
           klu1 = kl + ku + 1_${ik}$
           info = 0_${ik}$
           if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then
              info = -1_${ik}$
           else if( m<0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -4_${ik}$
           else if( kl<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ku<0_${ik}$ ) then
              info = -6_${ik}$
           else if( ldab<klu1 ) then
              info = -8_${ik}$
           else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then
              info = -12_${ik}$
           else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then
              info = -14_${ik}$
           else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGBBRD', -info )
              return
           end if
           ! initialize q and p**h to the unit matrix, if needed
           if( wantq )call stdlib${ii}$_zlaset( 'FULL', m, m, czero, cone, q, ldq )
           if( wantpt )call stdlib${ii}$_zlaset( 'FULL', n, n, czero, cone, pt, ldpt )
           ! quick return if possible.
           if( m==0 .or. n==0 )return
           minmn = min( m, n )
           if( kl+ku>1_${ik}$ ) then
              ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce
              ! first to lower bidiagonal form and then transform to upper
              ! bidiagonal
              if( ku>0_${ik}$ ) then
                 ml0 = 1_${ik}$
                 mu0 = 2_${ik}$
              else
                 ml0 = 2_${ik}$
                 mu0 = 1_${ik}$
              end if
              ! wherever possible, plane rotations are generated and applied in
              ! vector operations of length nr over the index set j1:j2:klu1.
              ! the complex sines of the plane rotations are stored in work,
              ! and the real cosines in rwork.
              klm = min( m-1, kl )
              kun = min( n-1, ku )
              kb = klm + kun
              kb1 = kb + 1_${ik}$
              inca = kb1*ldab
              nr = 0_${ik}$
              j1 = klm + 2_${ik}$
              j2 = 1_${ik}$ - kun
              loop_90: do i = 1, minmn
                 ! reduce i-th column and i-th row of matrix to bidiagonal form
                 ml = klm + 1_${ik}$
                 mu = kun + 1_${ik}$
                 loop_80: do kk = 1, kb
                    j1 = j1 + kb
                    j2 = j2 + kb
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been created below the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_zlargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, &
                              rwork( j1 ), kb1 )
                    ! apply plane rotations from the left
                    do l = 1, kb
                       if( j2-klm+l-1>n ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_zlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( &
                                 klu1-l+1, j1-klm+l-1 ), inca,rwork( j1 ), work( j1 ), kb1 )
                    end do
                    if( ml>ml0 ) then
                       if( ml<=m-i+1 ) then
                          ! generate plane rotation to annihilate a(i+ml-1,i)
                          ! within the band, and apply rotation from the left
                          call stdlib${ii}$_zlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),rwork( i+ml-1 ), &
                                    work( i+ml-1 ), ra )
                          ab( ku+ml-1, i ) = ra
                          if( i<n )call stdlib${ii}$_zrot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-&
                                    1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,rwork( i+ml-1 ), work( i+ml-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantq ) then
                       ! accumulate product of plane rotations in q
                       do j = j1, j2, kb1
                          call stdlib${ii}$_zrot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,rwork( j ), conjg( &
                                    work( j ) ) )
                       end do
                    end if
                    if( wantc ) then
                       ! apply plane rotations to c
                       do j = j1, j2, kb1
                          call stdlib${ii}$_zrot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,rwork( j ), &
                                    work( j ) )
                       end do
                    end if
                    if( j2+kun>n ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j-1,j+ku) above the band
                       ! and store it in work(n+1:2*n)
                       work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun )
                       ab( 1_${ik}$, j+kun ) = rwork( j )*ab( 1_${ik}$, j+kun )
                    end do
                    ! generate plane rotations to annihilate nonzero elements
                    ! which have been generated above the band
                    if( nr>0_${ik}$ )call stdlib${ii}$_zlargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,&
                               rwork( j1+kun ),kb1 )
                    ! apply plane rotations from the right
                    do l = 1, kb
                       if( j2+l-1>m ) then
                          nrt = nr - 1_${ik}$
                       else
                          nrt = nr
                       end if
                       if( nrt>0_${ik}$ )call stdlib${ii}$_zlartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+&
                                 kun ), inca,rwork( j1+kun ), work( j1+kun ), kb1 )
                    end do
                    if( ml==ml0 .and. mu>mu0 ) then
                       if( mu<=n-i+1 ) then
                          ! generate plane rotation to annihilate a(i,i+mu-1)
                          ! within the band, and apply rotation from the right
                          call stdlib${ii}$_zlartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),rwork( &
                                    i+mu-1 ), work( i+mu-1 ), ra )
                          ab( ku-mu+3, i+mu-2 ) = ra
                          call stdlib${ii}$_zrot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-&
                                    mu+3, i+mu-1 ), 1_${ik}$,rwork( i+mu-1 ), work( i+mu-1 ) )
                       end if
                       nr = nr + 1_${ik}$
                       j1 = j1 - kb1
                    end if
                    if( wantpt ) then
                       ! accumulate product of plane rotations in p**h
                       do j = j1, j2, kb1
                          call stdlib${ii}$_zrot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, rwork(&
                                     j+kun ),conjg( work( j+kun ) ) )
                       end do
                    end if
                    if( j2+kb>m ) then
                       ! adjust j2 to keep within the bounds of the matrix
                       nr = nr - 1_${ik}$
                       j2 = j2 - kb1
                    end if
                    do j = j1, j2, kb1
                       ! create nonzero element a(j+kl+ku,j+ku-1) below the
                       ! band and store it in work(1:n)
                       work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
                       ab( klu1, j+kun ) = rwork( j+kun )*ab( klu1, j+kun )
                    end do
                    if( ml>ml0 ) then
                       ml = ml - 1_${ik}$
                    else
                       mu = mu - 1_${ik}$
                    end if
                 end do loop_80
              end do loop_90
           end if
           if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then
              ! a has been reduced to complex lower bidiagonal form
              ! transform lower bidiagonal form to upper bidiagonal by applying
              ! plane rotations from the left, overwriting superdiagonal
              ! elements on subdiagonal elements
              do i = 1, min( m-1, n )
                 call stdlib${ii}$_zlartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra )
                 ab( 1_${ik}$, i ) = ra
                 if( i<n ) then
                    ab( 2_${ik}$, i ) = rs*ab( 1_${ik}$, i+1 )
                    ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 )
                 end if
                 if( wantq )call stdlib${ii}$_zrot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc,conjg( rs ) )
                           
                 if( wantc )call stdlib${ii}$_zrot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs )
                           
              end do
           else
              ! a has been reduced to complex upper bidiagonal form or is
              ! diagonal
              if( ku>0_${ik}$ .and. m<n ) then
                 ! annihilate a(m,m+1) by applying plane rotations from the
                 ! right
                 rb = ab( ku, m+1 )
                 do i = m, 1, -1
                    call stdlib${ii}$_zlartg( ab( ku+1, i ), rb, rc, rs, ra )
                    ab( ku+1, i ) = ra
                    if( i>1_${ik}$ ) then
                       rb = -conjg( rs )*ab( ku, i )
                       ab( ku, i ) = rc*ab( ku, i )
                    end if
                    if( wantpt )call stdlib${ii}$_zrot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, &
                              conjg( rs ) )
                 end do
              end if
           end if
           ! make diagonal and superdiagonal elements real, storing them in d
           ! and e
           t = ab( ku+1, 1_${ik}$ )
           loop_120: do i = 1, minmn
              abst = abs( t )
              d( i ) = abst
              if( abst/=zero ) then
                 t = t / abst
              else
                 t = cone
              end if
              if( wantq )call stdlib${ii}$_zscal( m, t, q( 1_${ik}$, i ), 1_${ik}$ )
              if( wantc )call stdlib${ii}$_zscal( ncc, conjg( t ), c( i, 1_${ik}$ ), ldc )
              if( i<minmn ) then
                 if( ku==0_${ik}$ .and. kl==0_${ik}$ ) then
                    e( i ) = zero
                    t = ab( 1_${ik}$, i+1 )
                 else
                    if( ku==0_${ik}$ ) then
                       t = ab( 2_${ik}$, i )*conjg( t )
                    else
                       t = ab( ku, i+1 )*conjg( t )
                    end if
                    abst = abs( t )
                    e( i ) = abst
                    if( abst/=zero ) then
                       t = t / abst
                    else
                       t = cone
                    end if
                    if( wantpt )call stdlib${ii}$_zscal( n, t, pt( i+1, 1_${ik}$ ), ldpt )
                    t = ab( ku+1, i+1 )*conjg( t )
                 end if
              end if
           end do loop_120
           return
     end subroutine stdlib${ii}$_zgbbrd

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$gbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, &
     !! ZGBBRD: reduces a complex general m-by-n band matrix A to real upper
     !! bidiagonal form B by a unitary transformation: Q**H * A * P = B.
     !! The routine computes B, and optionally forms Q or P**H, or computes
     !! Q**H*C for a given matrix C.
               ldc, work, rwork, info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: vect
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc
           ! Array Arguments 
           real(${ck}$), intent(out) :: d(*), e(*), rwork(*)
           complex(${ck}$), intent(inout) :: ab(ldab,*), c(ldc,*)
           complex(${ck}$), intent(out) :: pt(ldpt,*), q(ldq,*), work(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           logical(lk) :: wantb