stdlib_lapack_eigv_svd_drivers3.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_sbdsqr( uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u,ldu, c, ldc, work, &
     !! SBDSQR computes the singular values and, optionally, the right and/or
     !! left singular vectors from the singular value decomposition (SVD) of
     !! a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
     !! zero-shift QR algorithm.  The SVD of B has the form
     !! B = Q * S * P**T
     !! where S is the diagonal matrix of singular values, Q is an orthogonal
     !! matrix of left singular vectors, and P is an orthogonal matrix of
     !! right singular vectors.  If left singular vectors are requested, this
     !! subroutine actually returns U*Q instead of Q, and, if right singular
     !! vectors are requested, this subroutine returns P**T*VT instead of
     !! P**T, for given real input matrices U and VT.  When U and VT are the
     !! orthogonal matrices that reduce a general matrix A to bidiagonal
     !! form:  A = U*B*VT, as computed by SGEBRD, then
     !! A = (U*Q) * S * (P**T*VT)
     !! is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
     !! for a given real input matrix C.
     !! See "Computing  Small Singular Values of Bidiagonal Matrices With
     !! Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
     !! LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
     !! no. 5, pp. 873-912, Sept 1990) and
     !! "Accurate singular values and differential qd algorithms," by
     !! B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
     !! Department, University of California at Berkeley, July 1992
     !! for a detailed description of the algorithm.
               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) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldc, ldu, ldvt, n, ncc, ncvt, nru
           ! Array Arguments 
           real(sp), intent(inout) :: c(ldc,*), d(*), e(*), u(ldu,*), vt(ldvt,*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: hndrth = 0.01_sp
           real(sp), parameter :: hndrd = 100.0_sp
           real(sp), parameter :: meigth = -0.125_sp
           integer(${ik}$), parameter :: maxitr = 6_${ik}$
           
           
           
           
           
           
           
           
           ! Local Scalars 
           logical(lk) :: lower, rotate
           integer(${ik}$) :: i, idir, isub, iter, iterdivn, j, ll, lll, m, maxitdivn, nm1, nm12, &
                     nm13, oldll, oldm
           real(sp) :: abse, abss, cosl, cosr, cs, eps, f, g, h, mu, oldcs, oldsn, r, shift, &
           sigmn, sigmx, sinl, sinr, sll, smax, smin, sminl, sminoa, sn, thresh, tol, tolmul, &
                     unfl
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           lower = stdlib_lsame( uplo, 'L' )
           if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.lower ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ncvt<0_${ik}$ ) then
              info = -3_${ik}$
           else if( nru<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ( ncvt==0_${ik}$ .and. ldvt<1_${ik}$ ) .or.( ncvt>0_${ik}$ .and. ldvt<max( 1_${ik}$, n ) ) ) then
              info = -9_${ik}$
           else if( ldu<max( 1_${ik}$, nru ) ) then
              info = -11_${ik}$
           else if( ( ncc==0_${ik}$ .and. ldc<1_${ik}$ ) .or.( ncc>0_${ik}$ .and. ldc<max( 1_${ik}$, n ) ) ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SBDSQR', -info )
              return
           end if
           if( n==0 )return
           if( n==1 )go to 160
           ! rotate is true if any singular vectors desired, false otherwise
           rotate = ( ncvt>0_${ik}$ ) .or. ( nru>0_${ik}$ ) .or. ( ncc>0_${ik}$ )
           ! if no singular vectors desired, use qd algorithm
           if( .not.rotate ) then
              call stdlib${ii}$_slasq1( n, d, e, work, info )
           ! if info equals 2, dqds didn't finish, try to finish
              if( info /= 2 ) return
              info = 0_${ik}$
           end if
           nm1 = n - 1_${ik}$
           nm12 = nm1 + nm1
           nm13 = nm12 + nm1
           idir = 0_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_slamch( 'EPSILON' )
           unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           ! if matrix lower bidiagonal, rotate to be upper bidiagonal
           ! by applying givens rotations on the left
           if( lower ) then
              do i = 1, n - 1
                 call stdlib${ii}$_slartg( d( i ), e( i ), cs, sn, r )
                 d( i ) = r
                 e( i ) = sn*d( i+1 )
                 d( i+1 ) = cs*d( i+1 )
                 work( i ) = cs
                 work( nm1+i ) = sn
              end do
              ! update singular vectors if desired
              if( nru>0_${ik}$ )call stdlib${ii}$_slasr( 'R', 'V', 'F', nru, n, work( 1_${ik}$ ), work( n ), u,ldu )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'F', n, ncc, work( 1_${ik}$ ), work( n ), c,ldc )
                        
           end if
           ! compute singular values to relative accuracy tol
           ! (by setting tol to be negative, algorithm will compute
           ! singular values to absolute accuracy abs(tol)*norm(input matrix))
           tolmul = max( ten, min( hndrd, eps**meigth ) )
           tol = tolmul*eps
           ! compute approximate maximum, minimum singular values
           smax = zero
           do i = 1, n
              smax = max( smax, abs( d( i ) ) )
           end do
           do i = 1, n - 1
              smax = max( smax, abs( e( i ) ) )
           end do
           sminl = zero
           if( tol>=zero ) then
              ! relative accuracy desired
              sminoa = abs( d( 1_${ik}$ ) )
              if( sminoa==zero )go to 50
              mu = sminoa
              do i = 2, n
                 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
                 sminoa = min( sminoa, mu )
                 if( sminoa==zero )go to 50
              end do
              50 continue
              sminoa = sminoa / sqrt( real( n,KIND=sp) )
              thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
           else
              ! absolute accuracy desired
              thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
           end if
           ! prepare for main iteration loop for the singular values
           ! (maxit is the maximum number of passes through the inner
           ! loop permitted before nonconvergence signalled.)
           maxitdivn = maxitr*n
           iterdivn = 0_${ik}$
           iter = -1_${ik}$
           oldll = -1_${ik}$
           oldm = -1_${ik}$
           ! m points to last element of unconverged part of matrix
           m = n
           ! begin main iteration loop
           60 continue
           ! check for convergence or exceeding iteration count
           if( m<=1 )go to 160
           if( iter>=n ) then
              iter = iter - n
              iterdivn = iterdivn + 1_${ik}$
              if( iterdivn>=maxitdivn )go to 200
           end if
           ! find diagonal block of matrix to work on
           if( tol<zero .and. abs( d( m ) )<=thresh )d( m ) = zero
           smax = abs( d( m ) )
           smin = smax
           do lll = 1, m - 1
              ll = m - lll
              abss = abs( d( ll ) )
              abse = abs( e( ll ) )
              if( tol<zero .and. abss<=thresh )d( ll ) = zero
              if( abse<=thresh )go to 80
              smin = min( smin, abss )
              smax = max( smax, abss, abse )
           end do
           ll = 0_${ik}$
           go to 90
           80 continue
           e( ll ) = zero
           ! matrix splits since e(ll) = 0
           if( ll==m-1 ) then
              ! convergence of bottom singular value, return to top of loop
              m = m - 1_${ik}$
              go to 60
           end if
           90 continue
           ll = ll + 1_${ik}$
           ! e(ll) through e(m-1) are nonzero, e(ll-1) is zero
           if( ll==m-1 ) then
              ! 2 by 2 block, handle separately
              call stdlib${ii}$_slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,cosr, sinl, cosl &
                        )
              d( m-1 ) = sigmx
              e( m-1 ) = zero
              d( m ) = sigmn
              ! compute singular vectors, if desired
              if( ncvt>0_${ik}$ )call stdlib${ii}$_srot( ncvt, vt( m-1, 1_${ik}$ ), ldvt, vt( m, 1_${ik}$ ), ldvt, cosr,sinr &
                        )
              if( nru>0_${ik}$ )call stdlib${ii}$_srot( nru, u( 1_${ik}$, m-1 ), 1_${ik}$, u( 1_${ik}$, m ), 1_${ik}$, cosl, sinl )
              if( ncc>0_${ik}$ )call stdlib${ii}$_srot( ncc, c( m-1, 1_${ik}$ ), ldc, c( m, 1_${ik}$ ), ldc, cosl,sinl )
                        
              m = m - 2_${ik}$
              go to 60
           end if
           ! if working on new submatrix, choose shift direction
           ! (from larger end diagonal element towards smaller)
           if( ll>oldm .or. m<oldll ) then
              if( abs( d( ll ) )>=abs( d( m ) ) ) then
                 ! chase bulge from top (big end) to bottom (small end)
                 idir = 1_${ik}$
              else
                 ! chase bulge from bottom (big end) to top (small end)
                 idir = 2_${ik}$
              end if
           end if
           ! apply convergence tests
           if( idir==1_${ik}$ ) then
              ! run convergence test in forward direction
              ! first apply standard test to bottom of matrix
              if( abs( e( m-1 ) )<=abs( tol )*abs( d( m ) ) .or.( tol<zero .and. abs( e( m-1 ) )&
                        <=thresh ) ) then
                 e( m-1 ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion forward
                 mu = abs( d( ll ) )
                 sminl = mu
                 do lll = ll, m - 1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           else
              ! run convergence test in backward direction
              ! first apply standard test to top of matrix
              if( abs( e( ll ) )<=abs( tol )*abs( d( ll ) ) .or.( tol<zero .and. abs( e( ll ) )&
                        <=thresh ) ) then
                 e( ll ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion backward
                 mu = abs( d( m ) )
                 sminl = mu
                 do lll = m - 1, ll, -1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           end if
           oldll = ll
           oldm = m
           ! compute shift.  first, test if shifting would ruin relative
           ! accuracy, and if so set the shift to zero.
           if( tol>=zero .and. n*tol*( sminl / smax )<=max( eps, hndrth*tol ) ) then
              ! use a zero shift to avoid loss of relative accuracy
              shift = zero
           else
              ! compute the shift from 2-by-2 block at end of matrix
              if( idir==1_${ik}$ ) then
                 sll = abs( d( ll ) )
                 call stdlib${ii}$_slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
              else
                 sll = abs( d( m ) )
                 call stdlib${ii}$_slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
              end if
              ! test if shift negligible, and if so set to zero
              if( sll>zero ) then
                 if( ( shift / sll )**2_${ik}$<eps )shift = zero
              end if
           end if
           ! increment iteration count
           iter = iter + m - ll
           ! if shift = 0, do simplified qr iteration
           if( shift==zero ) then
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = ll, m - 1
                    call stdlib${ii}$_slartg( d( i )*cs, e( i ), cs, sn, r )
                    if( i>ll )e( i-1 ) = oldsn*r
                    call stdlib${ii}$_slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
                    work( i-ll+1 ) = cs
                    work( i-ll+1+nm1 ) = sn
                    work( i-ll+1+nm12 ) = oldcs
                    work( i-ll+1+nm13 ) = oldsn
                 end do
                 h = d( m )*cs
                 d( m ) = h*oldcs
                 e( m-1 ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1_${ik}$ ),work( n ), &
                           vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),work( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),work( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_slartg( d( i )*cs, e( i-1 ), cs, sn, r )
                    if( i<m )e( i ) = oldsn*r
                    call stdlib${ii}$_slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
                    work( i-ll ) = cs
                    work( i-ll+nm1 ) = -sn
                    work( i-ll+nm12 ) = oldcs
                    work( i-ll+nm13 ) = -oldsn
                 end do
                 h = d( ll )*cs
                 d( ll ) = h*oldcs
                 e( ll ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),work( &
                           nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1_${ik}$ ),work( n ), u(&
                            1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1_${ik}$ ),work( n ), c(&
                            ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
              end if
           else
              ! use nonzero shift
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( ll ) )-shift )*( sign( one, d( ll ) )+shift / d( ll ) )
                 g = e( ll )
                 do i = ll, m - 1
                    call stdlib${ii}$_slartg( f, g, cosr, sinr, r )
                    if( i>ll )e( i-1 ) = r
                    f = cosr*d( i ) + sinr*e( i )
                    e( i ) = cosr*e( i ) - sinr*d( i )
                    g = sinr*d( i+1 )
                    d( i+1 ) = cosr*d( i+1 )
                    call stdlib${ii}$_slartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i ) + sinl*d( i+1 )
                    d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
                    if( i<m-1 ) then
                       g = sinl*e( i+1 )
                       e( i+1 ) = cosl*e( i+1 )
                    end if
                    work( i-ll+1 ) = cosr
                    work( i-ll+1+nm1 ) = sinr
                    work( i-ll+1+nm12 ) = cosl
                    work( i-ll+1+nm13 ) = sinl
                 end do
                 e( m-1 ) = f
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1_${ik}$ ),work( n ), &
                           vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),work( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),work( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /d( m ) )
                 g = e( m-1 )
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_slartg( f, g, cosr, sinr, r )
                    if( i<m )e( i ) = r
                    f = cosr*d( i ) + sinr*e( i-1 )
                    e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
                    g = sinr*d( i-1 )
                    d( i-1 ) = cosr*d( i-1 )
                    call stdlib${ii}$_slartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i-1 ) + sinl*d( i-1 )
                    d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
                    if( i>ll+1 ) then
                       g = sinl*e( i-2 )
                       e( i-2 ) = cosl*e( i-2 )
                    end if
                    work( i-ll ) = cosr
                    work( i-ll+nm1 ) = -sinr
                    work( i-ll+nm12 ) = cosl
                    work( i-ll+nm13 ) = -sinl
                 end do
                 e( ll ) = f
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
                 ! update singular vectors if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),work( &
                           nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1_${ik}$ ),work( n ), u(&
                            1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1_${ik}$ ),work( n ), c(&
                            ll, 1_${ik}$ ), ldc )
              end if
           end if
           ! qr iteration finished, go back and check convergence
           go to 60
           ! all singular values converged, so make them positive
           160 continue
           do i = 1, n
              if( d( i )<zero ) then
                 d( i ) = -d( i )
                 ! change sign of singular vectors, if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_sscal( ncvt, negone, vt( i, 1_${ik}$ ), ldvt )
              end if
           end do
           ! sort the singular values into decreasing order (insertion sort on
           ! singular values, but only one transposition per singular vector)
           do i = 1, n - 1
              ! scan for smallest d(i)
              isub = 1_${ik}$
              smin = d( 1_${ik}$ )
              do j = 2, n + 1 - i
                 if( d( j )<=smin ) then
                    isub = j
                    smin = d( j )
                 end if
              end do
              if( isub/=n+1-i ) then
                 ! swap singular values and vectors
                 d( isub ) = d( n+1-i )
                 d( n+1-i ) = smin
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_sswap( ncvt, vt( isub, 1_${ik}$ ), ldvt, vt( n+1-i, 1_${ik}$ ),ldvt )
                           
                 if( nru>0_${ik}$ )call stdlib${ii}$_sswap( nru, u( 1_${ik}$, isub ), 1_${ik}$, u( 1_${ik}$, n+1-i ), 1_${ik}$ )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_sswap( ncc, c( isub, 1_${ik}$ ), ldc, c( n+1-i, 1_${ik}$ ), ldc )
                           
              end if
           end do
           go to 220
           ! maximum number of iterations exceeded, failure to converge
           200 continue
           info = 0_${ik}$
           do i = 1, n - 1
              if( e( i )/=zero )info = info + 1_${ik}$
           end do
           220 continue
           return
     end subroutine stdlib${ii}$_sbdsqr

     pure module subroutine stdlib${ii}$_dbdsqr( uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u,ldu, c, ldc, work, &
     !! DBDSQR computes the singular values and, optionally, the right and/or
     !! left singular vectors from the singular value decomposition (SVD) of
     !! a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
     !! zero-shift QR algorithm.  The SVD of B has the form
     !! B = Q * S * P**T
     !! where S is the diagonal matrix of singular values, Q is an orthogonal
     !! matrix of left singular vectors, and P is an orthogonal matrix of
     !! right singular vectors.  If left singular vectors are requested, this
     !! subroutine actually returns U*Q instead of Q, and, if right singular
     !! vectors are requested, this subroutine returns P**T*VT instead of
     !! P**T, for given real input matrices U and VT.  When U and VT are the
     !! orthogonal matrices that reduce a general matrix A to bidiagonal
     !! form:  A = U*B*VT, as computed by DGEBRD, then
     !! A = (U*Q) * S * (P**T*VT)
     !! is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
     !! for a given real input matrix C.
     !! See "Computing  Small Singular Values of Bidiagonal Matrices With
     !! Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
     !! LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
     !! no. 5, pp. 873-912, Sept 1990) and
     !! "Accurate singular values and differential qd algorithms," by
     !! B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
     !! Department, University of California at Berkeley, July 1992
     !! for a detailed description of the algorithm.
               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) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldc, ldu, ldvt, n, ncc, ncvt, nru
           ! Array Arguments 
           real(dp), intent(inout) :: c(ldc,*), d(*), e(*), u(ldu,*), vt(ldvt,*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: hndrth = 0.01_dp
           real(dp), parameter :: hndrd = 100.0_dp
           real(dp), parameter :: meigth = -0.125_dp
           integer(${ik}$), parameter :: maxitr = 6_${ik}$
           
           
           
           
           
           
           
           
           ! Local Scalars 
           logical(lk) :: lower, rotate
           integer(${ik}$) :: i, idir, isub, iter, iterdivn, j, ll, lll, m, maxitdivn, nm1, nm12, &
                     nm13, oldll, oldm
           real(dp) :: abse, abss, cosl, cosr, cs, eps, f, g, h, mu, oldcs, oldsn, r, shift, &
           sigmn, sigmx, sinl, sinr, sll, smax, smin, sminl, sminoa, sn, thresh, tol, tolmul, &
                     unfl
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           lower = stdlib_lsame( uplo, 'L' )
           if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.lower ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ncvt<0_${ik}$ ) then
              info = -3_${ik}$
           else if( nru<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ( ncvt==0_${ik}$ .and. ldvt<1_${ik}$ ) .or.( ncvt>0_${ik}$ .and. ldvt<max( 1_${ik}$, n ) ) ) then
              info = -9_${ik}$
           else if( ldu<max( 1_${ik}$, nru ) ) then
              info = -11_${ik}$
           else if( ( ncc==0_${ik}$ .and. ldc<1_${ik}$ ) .or.( ncc>0_${ik}$ .and. ldc<max( 1_${ik}$, n ) ) ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DBDSQR', -info )
              return
           end if
           if( n==0 )return
           if( n==1 )go to 160
           ! rotate is true if any singular vectors desired, false otherwise
           rotate = ( ncvt>0_${ik}$ ) .or. ( nru>0_${ik}$ ) .or. ( ncc>0_${ik}$ )
           ! if no singular vectors desired, use qd algorithm
           if( .not.rotate ) then
              call stdlib${ii}$_dlasq1( n, d, e, work, info )
           ! if info equals 2, dqds didn't finish, try to finish
              if( info /= 2 ) return
              info = 0_${ik}$
           end if
           nm1 = n - 1_${ik}$
           nm12 = nm1 + nm1
           nm13 = nm12 + nm1
           idir = 0_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_dlamch( 'EPSILON' )
           unfl = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           ! if matrix lower bidiagonal, rotate to be upper bidiagonal
           ! by applying givens rotations on the left
           if( lower ) then
              do i = 1, n - 1
                 call stdlib${ii}$_dlartg( d( i ), e( i ), cs, sn, r )
                 d( i ) = r
                 e( i ) = sn*d( i+1 )
                 d( i+1 ) = cs*d( i+1 )
                 work( i ) = cs
                 work( nm1+i ) = sn
              end do
              ! update singular vectors if desired
              if( nru>0_${ik}$ )call stdlib${ii}$_dlasr( 'R', 'V', 'F', nru, n, work( 1_${ik}$ ), work( n ), u,ldu )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'F', n, ncc, work( 1_${ik}$ ), work( n ), c,ldc )
                        
           end if
           ! compute singular values to relative accuracy tol
           ! (by setting tol to be negative, algorithm will compute
           ! singular values to absolute accuracy abs(tol)*norm(input matrix))
           tolmul = max( ten, min( hndrd, eps**meigth ) )
           tol = tolmul*eps
           ! compute approximate maximum, minimum singular values
           smax = zero
           do i = 1, n
              smax = max( smax, abs( d( i ) ) )
           end do
           do i = 1, n - 1
              smax = max( smax, abs( e( i ) ) )
           end do
           sminl = zero
           if( tol>=zero ) then
              ! relative accuracy desired
              sminoa = abs( d( 1_${ik}$ ) )
              if( sminoa==zero )go to 50
              mu = sminoa
              do i = 2, n
                 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
                 sminoa = min( sminoa, mu )
                 if( sminoa==zero )go to 50
              end do
              50 continue
              sminoa = sminoa / sqrt( real( n,KIND=dp) )
              thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
           else
              ! absolute accuracy desired
              thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
           end if
           ! prepare for main iteration loop for the singular values
           ! (maxit is the maximum number of passes through the inner
           ! loop permitted before nonconvergence signalled.)
           maxitdivn = maxitr*n
           iterdivn = 0_${ik}$
           iter = -1_${ik}$
           oldll = -1_${ik}$
           oldm = -1_${ik}$
           ! m points to last element of unconverged part of matrix
           m = n
           ! begin main iteration loop
           60 continue
           ! check for convergence or exceeding iteration count
           if( m<=1 )go to 160
           if( iter>=n ) then
              iter = iter - n
              iterdivn = iterdivn + 1_${ik}$
              if( iterdivn>=maxitdivn )go to 200
           end if
           ! find diagonal block of matrix to work on
           if( tol<zero .and. abs( d( m ) )<=thresh )d( m ) = zero
           smax = abs( d( m ) )
           smin = smax
           do lll = 1, m - 1
              ll = m - lll
              abss = abs( d( ll ) )
              abse = abs( e( ll ) )
              if( tol<zero .and. abss<=thresh )d( ll ) = zero
              if( abse<=thresh )go to 80
              smin = min( smin, abss )
              smax = max( smax, abss, abse )
           end do
           ll = 0_${ik}$
           go to 90
           80 continue
           e( ll ) = zero
           ! matrix splits since e(ll) = 0
           if( ll==m-1 ) then
              ! convergence of bottom singular value, return to top of loop
              m = m - 1_${ik}$
              go to 60
           end if
           90 continue
           ll = ll + 1_${ik}$
           ! e(ll) through e(m-1) are nonzero, e(ll-1) is zero
           if( ll==m-1 ) then
              ! 2 by 2 block, handle separately
              call stdlib${ii}$_dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,cosr, sinl, cosl &
                        )
              d( m-1 ) = sigmx
              e( m-1 ) = zero
              d( m ) = sigmn
              ! compute singular vectors, if desired
              if( ncvt>0_${ik}$ )call stdlib${ii}$_drot( ncvt, vt( m-1, 1_${ik}$ ), ldvt, vt( m, 1_${ik}$ ), ldvt, cosr,sinr &
                        )
              if( nru>0_${ik}$ )call stdlib${ii}$_drot( nru, u( 1_${ik}$, m-1 ), 1_${ik}$, u( 1_${ik}$, m ), 1_${ik}$, cosl, sinl )
              if( ncc>0_${ik}$ )call stdlib${ii}$_drot( ncc, c( m-1, 1_${ik}$ ), ldc, c( m, 1_${ik}$ ), ldc, cosl,sinl )
                        
              m = m - 2_${ik}$
              go to 60
           end if
           ! if working on new submatrix, choose shift direction
           ! (from larger end diagonal element towards smaller)
           if( ll>oldm .or. m<oldll ) then
              if( abs( d( ll ) )>=abs( d( m ) ) ) then
                 ! chase bulge from top (big end) to bottom (small end)
                 idir = 1_${ik}$
              else
                 ! chase bulge from bottom (big end) to top (small end)
                 idir = 2_${ik}$
              end if
           end if
           ! apply convergence tests
           if( idir==1_${ik}$ ) then
              ! run convergence test in forward direction
              ! first apply standard test to bottom of matrix
              if( abs( e( m-1 ) )<=abs( tol )*abs( d( m ) ) .or.( tol<zero .and. abs( e( m-1 ) )&
                        <=thresh ) ) then
                 e( m-1 ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion forward
                 mu = abs( d( ll ) )
                 sminl = mu
                 do lll = ll, m - 1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           else
              ! run convergence test in backward direction
              ! first apply standard test to top of matrix
              if( abs( e( ll ) )<=abs( tol )*abs( d( ll ) ) .or.( tol<zero .and. abs( e( ll ) )&
                        <=thresh ) ) then
                 e( ll ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion backward
                 mu = abs( d( m ) )
                 sminl = mu
                 do lll = m - 1, ll, -1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           end if
           oldll = ll
           oldm = m
           ! compute shift.  first, test if shifting would ruin relative
           ! accuracy, and if so set the shift to zero.
           if( tol>=zero .and. n*tol*( sminl / smax )<=max( eps, hndrth*tol ) ) then
              ! use a zero shift to avoid loss of relative accuracy
              shift = zero
           else
              ! compute the shift from 2-by-2 block at end of matrix
              if( idir==1_${ik}$ ) then
                 sll = abs( d( ll ) )
                 call stdlib${ii}$_dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
              else
                 sll = abs( d( m ) )
                 call stdlib${ii}$_dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
              end if
              ! test if shift negligible, and if so set to zero
              if( sll>zero ) then
                 if( ( shift / sll )**2_${ik}$<eps )shift = zero
              end if
           end if
           ! increment iteration count
           iter = iter + m - ll
           ! if shift = 0, do simplified qr iteration
           if( shift==zero ) then
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = ll, m - 1
                    call stdlib${ii}$_dlartg( d( i )*cs, e( i ), cs, sn, r )
                    if( i>ll )e( i-1 ) = oldsn*r
                    call stdlib${ii}$_dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
                    work( i-ll+1 ) = cs
                    work( i-ll+1+nm1 ) = sn
                    work( i-ll+1+nm12 ) = oldcs
                    work( i-ll+1+nm13 ) = oldsn
                 end do
                 h = d( m )*cs
                 d( m ) = h*oldcs
                 e( m-1 ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1_${ik}$ ),work( n ), &
                           vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),work( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),work( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
                    if( i<m )e( i ) = oldsn*r
                    call stdlib${ii}$_dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
                    work( i-ll ) = cs
                    work( i-ll+nm1 ) = -sn
                    work( i-ll+nm12 ) = oldcs
                    work( i-ll+nm13 ) = -oldsn
                 end do
                 h = d( ll )*cs
                 d( ll ) = h*oldcs
                 e( ll ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),work( &
                           nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1_${ik}$ ),work( n ), u(&
                            1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1_${ik}$ ),work( n ), c(&
                            ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
              end if
           else
              ! use nonzero shift
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( ll ) )-shift )*( sign( one, d( ll ) )+shift / d( ll ) )
                 g = e( ll )
                 do i = ll, m - 1
                    call stdlib${ii}$_dlartg( f, g, cosr, sinr, r )
                    if( i>ll )e( i-1 ) = r
                    f = cosr*d( i ) + sinr*e( i )
                    e( i ) = cosr*e( i ) - sinr*d( i )
                    g = sinr*d( i+1 )
                    d( i+1 ) = cosr*d( i+1 )
                    call stdlib${ii}$_dlartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i ) + sinl*d( i+1 )
                    d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
                    if( i<m-1 ) then
                       g = sinl*e( i+1 )
                       e( i+1 ) = cosl*e( i+1 )
                    end if
                    work( i-ll+1 ) = cosr
                    work( i-ll+1+nm1 ) = sinr
                    work( i-ll+1+nm12 ) = cosl
                    work( i-ll+1+nm13 ) = sinl
                 end do
                 e( m-1 ) = f
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1_${ik}$ ),work( n ), &
                           vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),work( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),work( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /d( m ) )
                 g = e( m-1 )
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_dlartg( f, g, cosr, sinr, r )
                    if( i<m )e( i ) = r
                    f = cosr*d( i ) + sinr*e( i-1 )
                    e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
                    g = sinr*d( i-1 )
                    d( i-1 ) = cosr*d( i-1 )
                    call stdlib${ii}$_dlartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i-1 ) + sinl*d( i-1 )
                    d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
                    if( i>ll+1 ) then
                       g = sinl*e( i-2 )
                       e( i-2 ) = cosl*e( i-2 )
                    end if
                    work( i-ll ) = cosr
                    work( i-ll+nm1 ) = -sinr
                    work( i-ll+nm12 ) = cosl
                    work( i-ll+nm13 ) = -sinl
                 end do
                 e( ll ) = f
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
                 ! update singular vectors if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),work( &
                           nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1_${ik}$ ),work( n ), u(&
                            1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1_${ik}$ ),work( n ), c(&
                            ll, 1_${ik}$ ), ldc )
              end if
           end if
           ! qr iteration finished, go back and check convergence
           go to 60
           ! all singular values converged, so make them positive
           160 continue
           do i = 1, n
              if( d( i )<zero ) then
                 d( i ) = -d( i )
                 ! change sign of singular vectors, if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_dscal( ncvt, negone, vt( i, 1_${ik}$ ), ldvt )
              end if
           end do
           ! sort the singular values into decreasing order (insertion sort on
           ! singular values, but only one transposition per singular vector)
           do i = 1, n - 1
              ! scan for smallest d(i)
              isub = 1_${ik}$
              smin = d( 1_${ik}$ )
              do j = 2, n + 1 - i
                 if( d( j )<=smin ) then
                    isub = j
                    smin = d( j )
                 end if
              end do
              if( isub/=n+1-i ) then
                 ! swap singular values and vectors
                 d( isub ) = d( n+1-i )
                 d( n+1-i ) = smin
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_dswap( ncvt, vt( isub, 1_${ik}$ ), ldvt, vt( n+1-i, 1_${ik}$ ),ldvt )
                           
                 if( nru>0_${ik}$ )call stdlib${ii}$_dswap( nru, u( 1_${ik}$, isub ), 1_${ik}$, u( 1_${ik}$, n+1-i ), 1_${ik}$ )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_dswap( ncc, c( isub, 1_${ik}$ ), ldc, c( n+1-i, 1_${ik}$ ), ldc )
                           
              end if
           end do
           go to 220
           ! maximum number of iterations exceeded, failure to converge
           200 continue
           info = 0_${ik}$
           do i = 1, n - 1
              if( e( i )/=zero )info = info + 1_${ik}$
           end do
           220 continue
           return
     end subroutine stdlib${ii}$_dbdsqr

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$bdsqr( uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u,ldu, c, ldc, work, &
     !! DBDSQR: computes the singular values and, optionally, the right and/or
     !! left singular vectors from the singular value decomposition (SVD) of
     !! a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
     !! zero-shift QR algorithm.  The SVD of B has the form
     !! B = Q * S * P**T
     !! where S is the diagonal matrix of singular values, Q is an orthogonal
     !! matrix of left singular vectors, and P is an orthogonal matrix of
     !! right singular vectors.  If left singular vectors are requested, this
     !! subroutine actually returns U*Q instead of Q, and, if right singular
     !! vectors are requested, this subroutine returns P**T*VT instead of
     !! P**T, for given real input matrices U and VT.  When U and VT are the
     !! orthogonal matrices that reduce a general matrix A to bidiagonal
     !! form:  A = U*B*VT, as computed by DGEBRD, then
     !! A = (U*Q) * S * (P**T*VT)
     !! is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
     !! for a given real input matrix C.
     !! See "Computing  Small Singular Values of Bidiagonal Matrices With
     !! Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
     !! LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
     !! no. 5, pp. 873-912, Sept 1990) and
     !! "Accurate singular values and differential qd algorithms," by
     !! B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
     !! Department, University of California at Berkeley, July 1992
     !! for a detailed description of the algorithm.
               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) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldc, ldu, ldvt, n, ncc, ncvt, nru
           ! Array Arguments 
           real(${rk}$), intent(inout) :: c(ldc,*), d(*), e(*), u(ldu,*), vt(ldvt,*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           ! Parameters 
           real(${rk}$), parameter :: hndrth = 0.01_${rk}$
           real(${rk}$), parameter :: hndrd = 100.0_${rk}$
           real(${rk}$), parameter :: meigth = -0.125_${rk}$
           integer(${ik}$), parameter :: maxitr = 6_${ik}$
           
           
           
           
           
           
           
           
           ! Local Scalars 
           logical(lk) :: lower, rotate
           integer(${ik}$) :: i, idir, isub, iter, iterdivn, j, ll, lll, m, maxitdivn, nm1, nm12, &
                     nm13, oldll, oldm
           real(${rk}$) :: abse, abss, cosl, cosr, cs, eps, f, g, h, mu, oldcs, oldsn, r, shift, &
           sigmn, sigmx, sinl, sinr, sll, smax, smin, sminl, sminoa, sn, thresh, tol, tolmul, &
                     unfl
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           lower = stdlib_lsame( uplo, 'L' )
           if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.lower ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ncvt<0_${ik}$ ) then
              info = -3_${ik}$
           else if( nru<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ( ncvt==0_${ik}$ .and. ldvt<1_${ik}$ ) .or.( ncvt>0_${ik}$ .and. ldvt<max( 1_${ik}$, n ) ) ) then
              info = -9_${ik}$
           else if( ldu<max( 1_${ik}$, nru ) ) then
              info = -11_${ik}$
           else if( ( ncc==0_${ik}$ .and. ldc<1_${ik}$ ) .or.( ncc>0_${ik}$ .and. ldc<max( 1_${ik}$, n ) ) ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DBDSQR', -info )
              return
           end if
           if( n==0 )return
           if( n==1 )go to 160
           ! rotate is true if any singular vectors desired, false otherwise
           rotate = ( ncvt>0_${ik}$ ) .or. ( nru>0_${ik}$ ) .or. ( ncc>0_${ik}$ )
           ! if no singular vectors desired, use qd algorithm
           if( .not.rotate ) then
              call stdlib${ii}$_${ri}$lasq1( n, d, e, work, info )
           ! if info equals 2, dqds didn't finish, try to finish
              if( info /= 2 ) return
              info = 0_${ik}$
           end if
           nm1 = n - 1_${ik}$
           nm12 = nm1 + nm1
           nm13 = nm12 + nm1
           idir = 0_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_${ri}$lamch( 'EPSILON' )
           unfl = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
           ! if matrix lower bidiagonal, rotate to be upper bidiagonal
           ! by applying givens rotations on the left
           if( lower ) then
              do i = 1, n - 1
                 call stdlib${ii}$_${ri}$lartg( d( i ), e( i ), cs, sn, r )
                 d( i ) = r
                 e( i ) = sn*d( i+1 )
                 d( i+1 ) = cs*d( i+1 )
                 work( i ) = cs
                 work( nm1+i ) = sn
              end do
              ! update singular vectors if desired
              if( nru>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'R', 'V', 'F', nru, n, work( 1_${ik}$ ), work( n ), u,ldu )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'F', n, ncc, work( 1_${ik}$ ), work( n ), c,ldc )
                        
           end if
           ! compute singular values to relative accuracy tol
           ! (by setting tol to be negative, algorithm will compute
           ! singular values to absolute accuracy abs(tol)*norm(input matrix))
           tolmul = max( ten, min( hndrd, eps**meigth ) )
           tol = tolmul*eps
           ! compute approximate maximum, minimum singular values
           smax = zero
           do i = 1, n
              smax = max( smax, abs( d( i ) ) )
           end do
           do i = 1, n - 1
              smax = max( smax, abs( e( i ) ) )
           end do
           sminl = zero
           if( tol>=zero ) then
              ! relative accuracy desired
              sminoa = abs( d( 1_${ik}$ ) )
              if( sminoa==zero )go to 50
              mu = sminoa
              do i = 2, n
                 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
                 sminoa = min( sminoa, mu )
                 if( sminoa==zero )go to 50
              end do
              50 continue
              sminoa = sminoa / sqrt( real( n,KIND=${rk}$) )
              thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
           else
              ! absolute accuracy desired
              thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
           end if
           ! prepare for main iteration loop for the singular values
           ! (maxit is the maximum number of passes through the inner
           ! loop permitted before nonconvergence signalled.)
           maxitdivn = maxitr*n
           iterdivn = 0_${ik}$
           iter = -1_${ik}$
           oldll = -1_${ik}$
           oldm = -1_${ik}$
           ! m points to last element of unconverged part of matrix
           m = n
           ! begin main iteration loop
           60 continue
           ! check for convergence or exceeding iteration count
           if( m<=1 )go to 160
           if( iter>=n ) then
              iter = iter - n
              iterdivn = iterdivn + 1_${ik}$
              if( iterdivn>=maxitdivn )go to 200
           end if
           ! find diagonal block of matrix to work on
           if( tol<zero .and. abs( d( m ) )<=thresh )d( m ) = zero
           smax = abs( d( m ) )
           smin = smax
           do lll = 1, m - 1
              ll = m - lll
              abss = abs( d( ll ) )
              abse = abs( e( ll ) )
              if( tol<zero .and. abss<=thresh )d( ll ) = zero
              if( abse<=thresh )go to 80
              smin = min( smin, abss )
              smax = max( smax, abss, abse )
           end do
           ll = 0_${ik}$
           go to 90
           80 continue
           e( ll ) = zero
           ! matrix splits since e(ll) = 0
           if( ll==m-1 ) then
              ! convergence of bottom singular value, return to top of loop
              m = m - 1_${ik}$
              go to 60
           end if
           90 continue
           ll = ll + 1_${ik}$
           ! e(ll) through e(m-1) are nonzero, e(ll-1) is zero
           if( ll==m-1 ) then
              ! 2 by 2 block, handle separately
              call stdlib${ii}$_${ri}$lasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,cosr, sinl, cosl &
                        )
              d( m-1 ) = sigmx
              e( m-1 ) = zero
              d( m ) = sigmn
              ! compute singular vectors, if desired
              if( ncvt>0_${ik}$ )call stdlib${ii}$_${ri}$rot( ncvt, vt( m-1, 1_${ik}$ ), ldvt, vt( m, 1_${ik}$ ), ldvt, cosr,sinr &
                        )
              if( nru>0_${ik}$ )call stdlib${ii}$_${ri}$rot( nru, u( 1_${ik}$, m-1 ), 1_${ik}$, u( 1_${ik}$, m ), 1_${ik}$, cosl, sinl )
              if( ncc>0_${ik}$ )call stdlib${ii}$_${ri}$rot( ncc, c( m-1, 1_${ik}$ ), ldc, c( m, 1_${ik}$ ), ldc, cosl,sinl )
                        
              m = m - 2_${ik}$
              go to 60
           end if
           ! if working on new submatrix, choose shift direction
           ! (from larger end diagonal element towards smaller)
           if( ll>oldm .or. m<oldll ) then
              if( abs( d( ll ) )>=abs( d( m ) ) ) then
                 ! chase bulge from top (big end) to bottom (small end)
                 idir = 1_${ik}$
              else
                 ! chase bulge from bottom (big end) to top (small end)
                 idir = 2_${ik}$
              end if
           end if
           ! apply convergence tests
           if( idir==1_${ik}$ ) then
              ! run convergence test in forward direction
              ! first apply standard test to bottom of matrix
              if( abs( e( m-1 ) )<=abs( tol )*abs( d( m ) ) .or.( tol<zero .and. abs( e( m-1 ) )&
                        <=thresh ) ) then
                 e( m-1 ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion forward
                 mu = abs( d( ll ) )
                 sminl = mu
                 do lll = ll, m - 1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           else
              ! run convergence test in backward direction
              ! first apply standard test to top of matrix
              if( abs( e( ll ) )<=abs( tol )*abs( d( ll ) ) .or.( tol<zero .and. abs( e( ll ) )&
                        <=thresh ) ) then
                 e( ll ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion backward
                 mu = abs( d( m ) )
                 sminl = mu
                 do lll = m - 1, ll, -1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           end if
           oldll = ll
           oldm = m
           ! compute shift.  first, test if shifting would ruin relative
           ! accuracy, and if so set the shift to zero.
           if( tol>=zero .and. n*tol*( sminl / smax )<=max( eps, hndrth*tol ) ) then
              ! use a zero shift to avoid loss of relative accuracy
              shift = zero
           else
              ! compute the shift from 2-by-2 block at end of matrix
              if( idir==1_${ik}$ ) then
                 sll = abs( d( ll ) )
                 call stdlib${ii}$_${ri}$las2( d( m-1 ), e( m-1 ), d( m ), shift, r )
              else
                 sll = abs( d( m ) )
                 call stdlib${ii}$_${ri}$las2( d( ll ), e( ll ), d( ll+1 ), shift, r )
              end if
              ! test if shift negligible, and if so set to zero
              if( sll>zero ) then
                 if( ( shift / sll )**2_${ik}$<eps )shift = zero
              end if
           end if
           ! increment iteration count
           iter = iter + m - ll
           ! if shift = 0, do simplified qr iteration
           if( shift==zero ) then
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = ll, m - 1
                    call stdlib${ii}$_${ri}$lartg( d( i )*cs, e( i ), cs, sn, r )
                    if( i>ll )e( i-1 ) = oldsn*r
                    call stdlib${ii}$_${ri}$lartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
                    work( i-ll+1 ) = cs
                    work( i-ll+1+nm1 ) = sn
                    work( i-ll+1+nm12 ) = oldcs
                    work( i-ll+1+nm13 ) = oldsn
                 end do
                 h = d( m )*cs
                 d( m ) = h*oldcs
                 e( m-1 ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1_${ik}$ ),work( n ), &
                           vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),work( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),work( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_${ri}$lartg( d( i )*cs, e( i-1 ), cs, sn, r )
                    if( i<m )e( i ) = oldsn*r
                    call stdlib${ii}$_${ri}$lartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
                    work( i-ll ) = cs
                    work( i-ll+nm1 ) = -sn
                    work( i-ll+nm12 ) = oldcs
                    work( i-ll+nm13 ) = -oldsn
                 end do
                 h = d( ll )*cs
                 d( ll ) = h*oldcs
                 e( ll ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),work( &
                           nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'R', 'V', 'B', nru, m-ll+1, work( 1_${ik}$ ),work( n ), u(&
                            1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1_${ik}$ ),work( n ), c(&
                            ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
              end if
           else
              ! use nonzero shift
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( ll ) )-shift )*( sign( one, d( ll ) )+shift / d( ll ) )
                 g = e( ll )
                 do i = ll, m - 1
                    call stdlib${ii}$_${ri}$lartg( f, g, cosr, sinr, r )
                    if( i>ll )e( i-1 ) = r
                    f = cosr*d( i ) + sinr*e( i )
                    e( i ) = cosr*e( i ) - sinr*d( i )
                    g = sinr*d( i+1 )
                    d( i+1 ) = cosr*d( i+1 )
                    call stdlib${ii}$_${ri}$lartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i ) + sinl*d( i+1 )
                    d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
                    if( i<m-1 ) then
                       g = sinl*e( i+1 )
                       e( i+1 ) = cosl*e( i+1 )
                    end if
                    work( i-ll+1 ) = cosr
                    work( i-ll+1+nm1 ) = sinr
                    work( i-ll+1+nm12 ) = cosl
                    work( i-ll+1+nm13 ) = sinl
                 end do
                 e( m-1 ) = f
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1_${ik}$ ),work( n ), &
                           vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),work( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),work( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /d( m ) )
                 g = e( m-1 )
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_${ri}$lartg( f, g, cosr, sinr, r )
                    if( i<m )e( i ) = r
                    f = cosr*d( i ) + sinr*e( i-1 )
                    e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
                    g = sinr*d( i-1 )
                    d( i-1 ) = cosr*d( i-1 )
                    call stdlib${ii}$_${ri}$lartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i-1 ) + sinl*d( i-1 )
                    d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
                    if( i>ll+1 ) then
                       g = sinl*e( i-2 )
                       e( i-2 ) = cosl*e( i-2 )
                    end if
                    work( i-ll ) = cosr
                    work( i-ll+nm1 ) = -sinr
                    work( i-ll+nm12 ) = cosl
                    work( i-ll+nm13 ) = -sinl
                 end do
                 e( ll ) = f
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
                 ! update singular vectors if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),work( &
                           nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'R', 'V', 'B', nru, m-ll+1, work( 1_${ik}$ ),work( n ), u(&
                            1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_${ri}$lasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1_${ik}$ ),work( n ), c(&
                            ll, 1_${ik}$ ), ldc )
              end if
           end if
           ! qr iteration finished, go back and check convergence
           go to 60
           ! all singular values converged, so make them positive
           160 continue
           do i = 1, n
              if( d( i )<zero ) then
                 d( i ) = -d( i )
                 ! change sign of singular vectors, if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_${ri}$scal( ncvt, negone, vt( i, 1_${ik}$ ), ldvt )
              end if
           end do
           ! sort the singular values into decreasing order (insertion sort on
           ! singular values, but only one transposition per singular vector)
           do i = 1, n - 1
              ! scan for smallest d(i)
              isub = 1_${ik}$
              smin = d( 1_${ik}$ )
              do j = 2, n + 1 - i
                 if( d( j )<=smin ) then
                    isub = j
                    smin = d( j )
                 end if
              end do
              if( isub/=n+1-i ) then
                 ! swap singular values and vectors
                 d( isub ) = d( n+1-i )
                 d( n+1-i ) = smin
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_${ri}$swap( ncvt, vt( isub, 1_${ik}$ ), ldvt, vt( n+1-i, 1_${ik}$ ),ldvt )
                           
                 if( nru>0_${ik}$ )call stdlib${ii}$_${ri}$swap( nru, u( 1_${ik}$, isub ), 1_${ik}$, u( 1_${ik}$, n+1-i ), 1_${ik}$ )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_${ri}$swap( ncc, c( isub, 1_${ik}$ ), ldc, c( n+1-i, 1_${ik}$ ), ldc )
                           
              end if
           end do
           go to 220
           ! maximum number of iterations exceeded, failure to converge
           200 continue
           info = 0_${ik}$
           do i = 1, n - 1
              if( e( i )/=zero )info = info + 1_${ik}$
           end do
           220 continue
           return
     end subroutine stdlib${ii}$_${ri}$bdsqr

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_cbdsqr( uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u,ldu, c, ldc, rwork,&
     !! CBDSQR computes the singular values and, optionally, the right and/or
     !! left singular vectors from the singular value decomposition (SVD) of
     !! a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
     !! zero-shift QR algorithm.  The SVD of B has the form
     !! B = Q * S * P**H
     !! where S is the diagonal matrix of singular values, Q is an orthogonal
     !! matrix of left singular vectors, and P is an orthogonal matrix of
     !! right singular vectors.  If left singular vectors are requested, this
     !! subroutine actually returns U*Q instead of Q, and, if right singular
     !! vectors are requested, this subroutine returns P**H*VT instead of
     !! P**H, for given complex input matrices U and VT.  When U and VT are
     !! the unitary matrices that reduce a general matrix A to bidiagonal
     !! form: A = U*B*VT, as computed by CGEBRD, then
     !! A = (U*Q) * S * (P**H*VT)
     !! is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
     !! for a given complex input matrix C.
     !! See "Computing  Small Singular Values of Bidiagonal Matrices With
     !! Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
     !! LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
     !! no. 5, pp. 873-912, Sept 1990) and
     !! "Accurate singular values and differential qd algorithms," by
     !! B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
     !! Department, University of California at Berkeley, July 1992
     !! for a detailed description of the algorithm.
                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) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldc, ldu, ldvt, n, ncc, ncvt, nru
           ! Array Arguments 
           real(sp), intent(inout) :: d(*), e(*)
           real(sp), intent(out) :: rwork(*)
           complex(sp), intent(inout) :: c(ldc,*), u(ldu,*), vt(ldvt,*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: hndrth = 0.01_sp
           real(sp), parameter :: hndrd = 100.0_sp
           real(sp), parameter :: meigth = -0.125_sp
           integer(${ik}$), parameter :: maxitr = 6_${ik}$
           
           
           
           
           
           
           
           
           ! Local Scalars 
           logical(lk) :: lower, rotate
           integer(${ik}$) :: i, idir, isub, iter, j, ll, lll, m, maxit, nm1, nm12, nm13, oldll, &
                     oldm
           real(sp) :: abse, abss, cosl, cosr, cs, eps, f, g, h, mu, oldcs, oldsn, r, shift, &
           sigmn, sigmx, sinl, sinr, sll, smax, smin, sminl, sminoa, sn, thresh, tol, tolmul, &
                     unfl
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           lower = stdlib_lsame( uplo, 'L' )
           if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.lower ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ncvt<0_${ik}$ ) then
              info = -3_${ik}$
           else if( nru<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ( ncvt==0_${ik}$ .and. ldvt<1_${ik}$ ) .or.( ncvt>0_${ik}$ .and. ldvt<max( 1_${ik}$, n ) ) ) then
              info = -9_${ik}$
           else if( ldu<max( 1_${ik}$, nru ) ) then
              info = -11_${ik}$
           else if( ( ncc==0_${ik}$ .and. ldc<1_${ik}$ ) .or.( ncc>0_${ik}$ .and. ldc<max( 1_${ik}$, n ) ) ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CBDSQR', -info )
              return
           end if
           if( n==0 )return
           if( n==1 )go to 160
           ! rotate is true if any singular vectors desired, false otherwise
           rotate = ( ncvt>0_${ik}$ ) .or. ( nru>0_${ik}$ ) .or. ( ncc>0_${ik}$ )
           ! if no singular vectors desired, use qd algorithm
           if( .not.rotate ) then
              call stdlib${ii}$_slasq1( n, d, e, rwork, info )
           ! if info equals 2, dqds didn't finish, try to finish
              if( info /= 2 ) return
              info = 0_${ik}$
           end if
           nm1 = n - 1_${ik}$
           nm12 = nm1 + nm1
           nm13 = nm12 + nm1
           idir = 0_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_slamch( 'EPSILON' )
           unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
           ! if matrix lower bidiagonal, rotate to be upper bidiagonal
           ! by applying givens rotations on the left
           if( lower ) then
              do i = 1, n - 1
                 call stdlib${ii}$_slartg( d( i ), e( i ), cs, sn, r )
                 d( i ) = r
                 e( i ) = sn*d( i+1 )
                 d( i+1 ) = cs*d( i+1 )
                 rwork( i ) = cs
                 rwork( nm1+i ) = sn
              end do
              ! update singular vectors if desired
              if( nru>0_${ik}$ )call stdlib${ii}$_clasr( 'R', 'V', 'F', nru, n, rwork( 1_${ik}$ ), rwork( n ),u, ldu )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'F', n, ncc, rwork( 1_${ik}$ ), rwork( n ),c, ldc )
                        
           end if
           ! compute singular values to relative accuracy tol
           ! (by setting tol to be negative, algorithm will compute
           ! singular values to absolute accuracy abs(tol)*norm(input matrix))
           tolmul = max( ten, min( hndrd, eps**meigth ) )
           tol = tolmul*eps
           ! compute approximate maximum, minimum singular values
           smax = zero
           do i = 1, n
              smax = max( smax, abs( d( i ) ) )
           end do
           do i = 1, n - 1
              smax = max( smax, abs( e( i ) ) )
           end do
           sminl = zero
           if( tol>=zero ) then
              ! relative accuracy desired
              sminoa = abs( d( 1_${ik}$ ) )
              if( sminoa==zero )go to 50
              mu = sminoa
              do i = 2, n
                 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
                 sminoa = min( sminoa, mu )
                 if( sminoa==zero )go to 50
              end do
              50 continue
              sminoa = sminoa / sqrt( real( n,KIND=sp) )
              thresh = max( tol*sminoa, maxitr*n*n*unfl )
           else
              ! absolute accuracy desired
              thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
           end if
           ! prepare for main iteration loop for the singular values
           ! (maxit is the maximum number of passes through the inner
           ! loop permitted before nonconvergence signalled.)
           maxit = maxitr*n*n
           iter = 0_${ik}$
           oldll = -1_${ik}$
           oldm = -1_${ik}$
           ! m points to last element of unconverged part of matrix
           m = n
           ! begin main iteration loop
           60 continue
           ! check for convergence or exceeding iteration count
           if( m<=1 )go to 160
           if( iter>maxit )go to 200
           ! find diagonal block of matrix to work on
           if( tol<zero .and. abs( d( m ) )<=thresh )d( m ) = zero
           smax = abs( d( m ) )
           smin = smax
           do lll = 1, m - 1
              ll = m - lll
              abss = abs( d( ll ) )
              abse = abs( e( ll ) )
              if( tol<zero .and. abss<=thresh )d( ll ) = zero
              if( abse<=thresh )go to 80
              smin = min( smin, abss )
              smax = max( smax, abss, abse )
           end do
           ll = 0_${ik}$
           go to 90
           80 continue
           e( ll ) = zero
           ! matrix splits since e(ll) = 0
           if( ll==m-1 ) then
              ! convergence of bottom singular value, return to top of loop
              m = m - 1_${ik}$
              go to 60
           end if
           90 continue
           ll = ll + 1_${ik}$
           ! e(ll) through e(m-1) are nonzero, e(ll-1) is zero
           if( ll==m-1 ) then
              ! 2 by 2 block, handle separately
              call stdlib${ii}$_slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,cosr, sinl, cosl &
                        )
              d( m-1 ) = sigmx
              e( m-1 ) = zero
              d( m ) = sigmn
              ! compute singular vectors, if desired
              if( ncvt>0_${ik}$ )call stdlib${ii}$_csrot( ncvt, vt( m-1, 1_${ik}$ ), ldvt, vt( m, 1_${ik}$ ), ldvt,cosr, &
                        sinr )
              if( nru>0_${ik}$ )call stdlib${ii}$_csrot( nru, u( 1_${ik}$, m-1 ), 1_${ik}$, u( 1_${ik}$, m ), 1_${ik}$, cosl, sinl )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_csrot( ncc, c( m-1, 1_${ik}$ ), ldc, c( m, 1_${ik}$ ), ldc, cosl,sinl )
                        
              m = m - 2_${ik}$
              go to 60
           end if
           ! if working on new submatrix, choose shift direction
           ! (from larger end diagonal element towards smaller)
           if( ll>oldm .or. m<oldll ) then
              if( abs( d( ll ) )>=abs( d( m ) ) ) then
                 ! chase bulge from top (big end) to bottom (small end)
                 idir = 1_${ik}$
              else
                 ! chase bulge from bottom (big end) to top (small end)
                 idir = 2_${ik}$
              end if
           end if
           ! apply convergence tests
           if( idir==1_${ik}$ ) then
              ! run convergence test in forward direction
              ! first apply standard test to bottom of matrix
              if( abs( e( m-1 ) )<=abs( tol )*abs( d( m ) ) .or.( tol<zero .and. abs( e( m-1 ) )&
                        <=thresh ) ) then
                 e( m-1 ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion forward
                 mu = abs( d( ll ) )
                 sminl = mu
                 do lll = ll, m - 1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           else
              ! run convergence test in backward direction
              ! first apply standard test to top of matrix
              if( abs( e( ll ) )<=abs( tol )*abs( d( ll ) ) .or.( tol<zero .and. abs( e( ll ) )&
                        <=thresh ) ) then
                 e( ll ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion backward
                 mu = abs( d( m ) )
                 sminl = mu
                 do lll = m - 1, ll, -1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           end if
           oldll = ll
           oldm = m
           ! compute shift.  first, test if shifting would ruin relative
           ! accuracy, and if so set the shift to zero.
           if( tol>=zero .and. n*tol*( sminl / smax )<=max( eps, hndrth*tol ) ) then
              ! use a zero shift to avoid loss of relative accuracy
              shift = zero
           else
              ! compute the shift from 2-by-2 block at end of matrix
              if( idir==1_${ik}$ ) then
                 sll = abs( d( ll ) )
                 call stdlib${ii}$_slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
              else
                 sll = abs( d( m ) )
                 call stdlib${ii}$_slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
              end if
              ! test if shift negligible, and if so set to zero
              if( sll>zero ) then
                 if( ( shift / sll )**2_${ik}$<eps )shift = zero
              end if
           end if
           ! increment iteration count
           iter = iter + m - ll
           ! if shift = 0, do simplified qr iteration
           if( shift==zero ) then
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = ll, m - 1
                    call stdlib${ii}$_slartg( d( i )*cs, e( i ), cs, sn, r )
                    if( i>ll )e( i-1 ) = oldsn*r
                    call stdlib${ii}$_slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
                    rwork( i-ll+1 ) = cs
                    rwork( i-ll+1+nm1 ) = sn
                    rwork( i-ll+1+nm12 ) = oldcs
                    rwork( i-ll+1+nm13 ) = oldsn
                 end do
                 h = d( m )*cs
                 d( m ) = h*oldcs
                 e( m-1 ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1_${ik}$ ),rwork( n )&
                           , vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_clasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),rwork( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'F', m-ll+1, ncc, rwork( nm12+1 ),rwork( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_slartg( d( i )*cs, e( i-1 ), cs, sn, r )
                    if( i<m )e( i ) = oldsn*r
                    call stdlib${ii}$_slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
                    rwork( i-ll ) = cs
                    rwork( i-ll+nm1 ) = -sn
                    rwork( i-ll+nm12 ) = oldcs
                    rwork( i-ll+nm13 ) = -oldsn
                 end do
                 h = d( ll )*cs
                 d( ll ) = h*oldcs
                 e( ll ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),&
                           rwork( nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_clasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1_${ik}$ ),rwork( n ), &
                           u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1_${ik}$ ),rwork( n ), &
                           c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
              end if
           else
              ! use nonzero shift
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( ll ) )-shift )*( sign( one, d( ll ) )+shift / d( ll ) )
                 g = e( ll )
                 do i = ll, m - 1
                    call stdlib${ii}$_slartg( f, g, cosr, sinr, r )
                    if( i>ll )e( i-1 ) = r
                    f = cosr*d( i ) + sinr*e( i )
                    e( i ) = cosr*e( i ) - sinr*d( i )
                    g = sinr*d( i+1 )
                    d( i+1 ) = cosr*d( i+1 )
                    call stdlib${ii}$_slartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i ) + sinl*d( i+1 )
                    d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
                    if( i<m-1 ) then
                       g = sinl*e( i+1 )
                       e( i+1 ) = cosl*e( i+1 )
                    end if
                    rwork( i-ll+1 ) = cosr
                    rwork( i-ll+1+nm1 ) = sinr
                    rwork( i-ll+1+nm12 ) = cosl
                    rwork( i-ll+1+nm13 ) = sinl
                 end do
                 e( m-1 ) = f
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1_${ik}$ ),rwork( n )&
                           , vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_clasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),rwork( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'F', m-ll+1, ncc, rwork( nm12+1 ),rwork( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /d( m ) )
                 g = e( m-1 )
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_slartg( f, g, cosr, sinr, r )
                    if( i<m )e( i ) = r
                    f = cosr*d( i ) + sinr*e( i-1 )
                    e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
                    g = sinr*d( i-1 )
                    d( i-1 ) = cosr*d( i-1 )
                    call stdlib${ii}$_slartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i-1 ) + sinl*d( i-1 )
                    d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
                    if( i>ll+1 ) then
                       g = sinl*e( i-2 )
                       e( i-2 ) = cosl*e( i-2 )
                    end if
                    rwork( i-ll ) = cosr
                    rwork( i-ll+nm1 ) = -sinr
                    rwork( i-ll+nm12 ) = cosl
                    rwork( i-ll+nm13 ) = -sinl
                 end do
                 e( ll ) = f
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
                 ! update singular vectors if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),&
                           rwork( nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_clasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1_${ik}$ ),rwork( n ), &
                           u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_clasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1_${ik}$ ),rwork( n ), &
                           c( ll, 1_${ik}$ ), ldc )
              end if
           end if
           ! qr iteration finished, go back and check convergence
           go to 60
           ! all singular values converged, so make them positive
           160 continue
           do i = 1, n
              if( d( i )<zero ) then
                 d( i ) = -d( i )
                 ! change sign of singular vectors, if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_csscal( ncvt, negone, vt( i, 1_${ik}$ ), ldvt )
              end if
           end do
           ! sort the singular values into decreasing order (insertion sort on
           ! singular values, but only one transposition per singular vector)
           do i = 1, n - 1
              ! scan for smallest d(i)
              isub = 1_${ik}$
              smin = d( 1_${ik}$ )
              do j = 2, n + 1 - i
                 if( d( j )<=smin ) then
                    isub = j
                    smin = d( j )
                 end if
              end do
              if( isub/=n+1-i ) then
                 ! swap singular values and vectors
                 d( isub ) = d( n+1-i )
                 d( n+1-i ) = smin
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_cswap( ncvt, vt( isub, 1_${ik}$ ), ldvt, vt( n+1-i, 1_${ik}$ ),ldvt )
                           
                 if( nru>0_${ik}$ )call stdlib${ii}$_cswap( nru, u( 1_${ik}$, isub ), 1_${ik}$, u( 1_${ik}$, n+1-i ), 1_${ik}$ )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_cswap( ncc, c( isub, 1_${ik}$ ), ldc, c( n+1-i, 1_${ik}$ ), ldc )
                           
              end if
           end do
           go to 220
           ! maximum number of iterations exceeded, failure to converge
           200 continue
           info = 0_${ik}$
           do i = 1, n - 1
              if( e( i )/=zero )info = info + 1_${ik}$
           end do
           220 continue
           return
     end subroutine stdlib${ii}$_cbdsqr

     pure module subroutine stdlib${ii}$_zbdsqr( uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u,ldu, c, ldc, rwork,&
     !! ZBDSQR computes the singular values and, optionally, the right and/or
     !! left singular vectors from the singular value decomposition (SVD) of
     !! a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
     !! zero-shift QR algorithm.  The SVD of B has the form
     !! B = Q * S * P**H
     !! where S is the diagonal matrix of singular values, Q is an orthogonal
     !! matrix of left singular vectors, and P is an orthogonal matrix of
     !! right singular vectors.  If left singular vectors are requested, this
     !! subroutine actually returns U*Q instead of Q, and, if right singular
     !! vectors are requested, this subroutine returns P**H*VT instead of
     !! P**H, for given complex input matrices U and VT.  When U and VT are
     !! the unitary matrices that reduce a general matrix A to bidiagonal
     !! form: A = U*B*VT, as computed by ZGEBRD, then
     !! A = (U*Q) * S * (P**H*VT)
     !! is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
     !! for a given complex input matrix C.
     !! See "Computing  Small Singular Values of Bidiagonal Matrices With
     !! Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
     !! LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
     !! no. 5, pp. 873-912, Sept 1990) and
     !! "Accurate singular values and differential qd algorithms," by
     !! B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
     !! Department, University of California at Berkeley, July 1992
     !! for a detailed description of the algorithm.
                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) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldc, ldu, ldvt, n, ncc, ncvt, nru
           ! Array Arguments 
           real(dp), intent(inout) :: d(*), e(*)
           real(dp), intent(out) :: rwork(*)
           complex(dp), intent(inout) :: c(ldc,*), u(ldu,*), vt(ldvt,*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: hndrth = 0.01_dp
           real(dp), parameter :: hndrd = 100.0_dp
           real(dp), parameter :: meigth = -0.125_dp
           integer(${ik}$), parameter :: maxitr = 6_${ik}$
           
           
           
           
           
           
           
           
           ! Local Scalars 
           logical(lk) :: lower, rotate
           integer(${ik}$) :: i, idir, isub, iter, j, ll, lll, m, maxit, nm1, nm12, nm13, oldll, &
                     oldm
           real(dp) :: abse, abss, cosl, cosr, cs, eps, f, g, h, mu, oldcs, oldsn, r, shift, &
           sigmn, sigmx, sinl, sinr, sll, smax, smin, sminl, sminoa, sn, thresh, tol, tolmul, &
                     unfl
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           lower = stdlib_lsame( uplo, 'L' )
           if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.lower ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ncvt<0_${ik}$ ) then
              info = -3_${ik}$
           else if( nru<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ( ncvt==0_${ik}$ .and. ldvt<1_${ik}$ ) .or.( ncvt>0_${ik}$ .and. ldvt<max( 1_${ik}$, n ) ) ) then
              info = -9_${ik}$
           else if( ldu<max( 1_${ik}$, nru ) ) then
              info = -11_${ik}$
           else if( ( ncc==0_${ik}$ .and. ldc<1_${ik}$ ) .or.( ncc>0_${ik}$ .and. ldc<max( 1_${ik}$, n ) ) ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZBDSQR', -info )
              return
           end if
           if( n==0 )return
           if( n==1 )go to 160
           ! rotate is true if any singular vectors desired, false otherwise
           rotate = ( ncvt>0_${ik}$ ) .or. ( nru>0_${ik}$ ) .or. ( ncc>0_${ik}$ )
           ! if no singular vectors desired, use qd algorithm
           if( .not.rotate ) then
              call stdlib${ii}$_dlasq1( n, d, e, rwork, info )
           ! if info equals 2, dqds didn't finish, try to finish
              if( info /= 2 ) return
              info = 0_${ik}$
           end if
           nm1 = n - 1_${ik}$
           nm12 = nm1 + nm1
           nm13 = nm12 + nm1
           idir = 0_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_dlamch( 'EPSILON' )
           unfl = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
           ! if matrix lower bidiagonal, rotate to be upper bidiagonal
           ! by applying givens rotations on the left
           if( lower ) then
              do i = 1, n - 1
                 call stdlib${ii}$_dlartg( d( i ), e( i ), cs, sn, r )
                 d( i ) = r
                 e( i ) = sn*d( i+1 )
                 d( i+1 ) = cs*d( i+1 )
                 rwork( i ) = cs
                 rwork( nm1+i ) = sn
              end do
              ! update singular vectors if desired
              if( nru>0_${ik}$ )call stdlib${ii}$_zlasr( 'R', 'V', 'F', nru, n, rwork( 1_${ik}$ ), rwork( n ),u, ldu )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'F', n, ncc, rwork( 1_${ik}$ ), rwork( n ),c, ldc )
                        
           end if
           ! compute singular values to relative accuracy tol
           ! (by setting tol to be negative, algorithm will compute
           ! singular values to absolute accuracy abs(tol)*norm(input matrix))
           tolmul = max( ten, min( hndrd, eps**meigth ) )
           tol = tolmul*eps
           ! compute approximate maximum, minimum singular values
           smax = zero
           do i = 1, n
              smax = max( smax, abs( d( i ) ) )
           end do
           do i = 1, n - 1
              smax = max( smax, abs( e( i ) ) )
           end do
           sminl = zero
           if( tol>=zero ) then
              ! relative accuracy desired
              sminoa = abs( d( 1_${ik}$ ) )
              if( sminoa==zero )go to 50
              mu = sminoa
              do i = 2, n
                 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
                 sminoa = min( sminoa, mu )
                 if( sminoa==zero )go to 50
              end do
              50 continue
              sminoa = sminoa / sqrt( real( n,KIND=dp) )
              thresh = max( tol*sminoa, maxitr*n*n*unfl )
           else
              ! absolute accuracy desired
              thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
           end if
           ! prepare for main iteration loop for the singular values
           ! (maxit is the maximum number of passes through the inner
           ! loop permitted before nonconvergence signalled.)
           maxit = maxitr*n*n
           iter = 0_${ik}$
           oldll = -1_${ik}$
           oldm = -1_${ik}$
           ! m points to last element of unconverged part of matrix
           m = n
           ! begin main iteration loop
           60 continue
           ! check for convergence or exceeding iteration count
           if( m<=1 )go to 160
           if( iter>maxit )go to 200
           ! find diagonal block of matrix to work on
           if( tol<zero .and. abs( d( m ) )<=thresh )d( m ) = zero
           smax = abs( d( m ) )
           smin = smax
           do lll = 1, m - 1
              ll = m - lll
              abss = abs( d( ll ) )
              abse = abs( e( ll ) )
              if( tol<zero .and. abss<=thresh )d( ll ) = zero
              if( abse<=thresh )go to 80
              smin = min( smin, abss )
              smax = max( smax, abss, abse )
           end do
           ll = 0_${ik}$
           go to 90
           80 continue
           e( ll ) = zero
           ! matrix splits since e(ll) = 0
           if( ll==m-1 ) then
              ! convergence of bottom singular value, return to top of loop
              m = m - 1_${ik}$
              go to 60
           end if
           90 continue
           ll = ll + 1_${ik}$
           ! e(ll) through e(m-1) are nonzero, e(ll-1) is zero
           if( ll==m-1 ) then
              ! 2 by 2 block, handle separately
              call stdlib${ii}$_dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,cosr, sinl, cosl &
                        )
              d( m-1 ) = sigmx
              e( m-1 ) = zero
              d( m ) = sigmn
              ! compute singular vectors, if desired
              if( ncvt>0_${ik}$ )call stdlib${ii}$_zdrot( ncvt, vt( m-1, 1_${ik}$ ), ldvt, vt( m, 1_${ik}$ ), ldvt,cosr, &
                        sinr )
              if( nru>0_${ik}$ )call stdlib${ii}$_zdrot( nru, u( 1_${ik}$, m-1 ), 1_${ik}$, u( 1_${ik}$, m ), 1_${ik}$, cosl, sinl )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_zdrot( ncc, c( m-1, 1_${ik}$ ), ldc, c( m, 1_${ik}$ ), ldc, cosl,sinl )
                        
              m = m - 2_${ik}$
              go to 60
           end if
           ! if working on new submatrix, choose shift direction
           ! (from larger end diagonal element towards smaller)
           if( ll>oldm .or. m<oldll ) then
              if( abs( d( ll ) )>=abs( d( m ) ) ) then
                 ! chase bulge from top (big end) to bottom (small end)
                 idir = 1_${ik}$
              else
                 ! chase bulge from bottom (big end) to top (small end)
                 idir = 2_${ik}$
              end if
           end if
           ! apply convergence tests
           if( idir==1_${ik}$ ) then
              ! run convergence test in forward direction
              ! first apply standard test to bottom of matrix
              if( abs( e( m-1 ) )<=abs( tol )*abs( d( m ) ) .or.( tol<zero .and. abs( e( m-1 ) )&
                        <=thresh ) ) then
                 e( m-1 ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion forward
                 mu = abs( d( ll ) )
                 sminl = mu
                 do lll = ll, m - 1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           else
              ! run convergence test in backward direction
              ! first apply standard test to top of matrix
              if( abs( e( ll ) )<=abs( tol )*abs( d( ll ) ) .or.( tol<zero .and. abs( e( ll ) )&
                        <=thresh ) ) then
                 e( ll ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion backward
                 mu = abs( d( m ) )
                 sminl = mu
                 do lll = m - 1, ll, -1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           end if
           oldll = ll
           oldm = m
           ! compute shift.  first, test if shifting would ruin relative
           ! accuracy, and if so set the shift to zero.
           if( tol>=zero .and. n*tol*( sminl / smax )<=max( eps, hndrth*tol ) ) then
              ! use a zero shift to avoid loss of relative accuracy
              shift = zero
           else
              ! compute the shift from 2-by-2 block at end of matrix
              if( idir==1_${ik}$ ) then
                 sll = abs( d( ll ) )
                 call stdlib${ii}$_dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
              else
                 sll = abs( d( m ) )
                 call stdlib${ii}$_dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
              end if
              ! test if shift negligible, and if so set to zero
              if( sll>zero ) then
                 if( ( shift / sll )**2_${ik}$<eps )shift = zero
              end if
           end if
           ! increment iteration count
           iter = iter + m - ll
           ! if shift = 0, do simplified qr iteration
           if( shift==zero ) then
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = ll, m - 1
                    call stdlib${ii}$_dlartg( d( i )*cs, e( i ), cs, sn, r )
                    if( i>ll )e( i-1 ) = oldsn*r
                    call stdlib${ii}$_dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
                    rwork( i-ll+1 ) = cs
                    rwork( i-ll+1+nm1 ) = sn
                    rwork( i-ll+1+nm12 ) = oldcs
                    rwork( i-ll+1+nm13 ) = oldsn
                 end do
                 h = d( m )*cs
                 d( m ) = h*oldcs
                 e( m-1 ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1_${ik}$ ),rwork( n )&
                           , vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_zlasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),rwork( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'F', m-ll+1, ncc, rwork( nm12+1 ),rwork( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 cs = one
                 oldcs = one
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
                    if( i<m )e( i ) = oldsn*r
                    call stdlib${ii}$_dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
                    rwork( i-ll ) = cs
                    rwork( i-ll+nm1 ) = -sn
                    rwork( i-ll+nm12 ) = oldcs
                    rwork( i-ll+nm13 ) = -oldsn
                 end do
                 h = d( ll )*cs
                 d( ll ) = h*oldcs
                 e( ll ) = h*oldsn
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),&
                           rwork( nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1_${ik}$ ),rwork( n ), &
                           u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1_${ik}$ ),rwork( n ), &
                           c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
              end if
           else
              ! use nonzero shift
              if( idir==1_${ik}$ ) then
                 ! chase bulge from top to bottom
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( ll ) )-shift )*( sign( one, d( ll ) )+shift / d( ll ) )
                 g = e( ll )
                 do i = ll, m - 1
                    call stdlib${ii}$_dlartg( f, g, cosr, sinr, r )
                    if( i>ll )e( i-1 ) = r
                    f = cosr*d( i ) + sinr*e( i )
                    e( i ) = cosr*e( i ) - sinr*d( i )
                    g = sinr*d( i+1 )
                    d( i+1 ) = cosr*d( i+1 )
                    call stdlib${ii}$_dlartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i ) + sinl*d( i+1 )
                    d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
                    if( i<m-1 ) then
                       g = sinl*e( i+1 )
                       e( i+1 ) = cosl*e( i+1 )
                    end if
                    rwork( i-ll+1 ) = cosr
                    rwork( i-ll+1+nm1 ) = sinr
                    rwork( i-ll+1+nm12 ) = cosl
                    rwork( i-ll+1+nm13 ) = sinl
                 end do
                 e( m-1 ) = f
                 ! update singular vectors
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1_${ik}$ ),rwork( n )&
                           , vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_zlasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),rwork( &
                           nm13+1 ), u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'F', m-ll+1, ncc, rwork( nm12+1 ),rwork( &
                           nm13+1 ), c( ll, 1_${ik}$ ), ldc )
                 ! test convergence
                 if( abs( e( m-1 ) )<=thresh )e( m-1 ) = zero
              else
                 ! chase bulge from bottom to top
                 ! save cosines and sines for later singular vector updates
                 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /d( m ) )
                 g = e( m-1 )
                 do i = m, ll + 1, -1
                    call stdlib${ii}$_dlartg( f, g, cosr, sinr, r )
                    if( i<m )e( i ) = r
                    f = cosr*d( i ) + sinr*e( i-1 )
                    e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
                    g = sinr*d( i-1 )
                    d( i-1 ) = cosr*d( i-1 )
                    call stdlib${ii}$_dlartg( f, g, cosl, sinl, r )
                    d( i ) = r
                    f = cosl*e( i-1 ) + sinl*d( i-1 )
                    d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
                    if( i>ll+1 ) then
                       g = sinl*e( i-2 )
                       e( i-2 ) = cosl*e( i-2 )
                    end if
                    rwork( i-ll ) = cosr
                    rwork( i-ll+nm1 ) = -sinr
                    rwork( i-ll+nm12 ) = cosl
                    rwork( i-ll+nm13 ) = -sinl
                 end do
                 e( ll ) = f
                 ! test convergence
                 if( abs( e( ll ) )<=thresh )e( ll ) = zero
                 ! update singular vectors if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),&
                           rwork( nm13+1 ), vt( ll, 1_${ik}$ ), ldvt )
                 if( nru>0_${ik}$ )call stdlib${ii}$_zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1_${ik}$ ),rwork( n ), &
                           u( 1_${ik}$, ll ), ldu )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_zlasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1_${ik}$ ),rwork( n ), &
                           c( ll, 1_${ik}$ ), ldc )
              end if
           end if
           ! qr iteration finished, go back and check convergence
           go to 60
           ! all singular values converged, so make them positive
           160 continue
           do i = 1, n
              if( d( i )<zero ) then
                 d( i ) = -d( i )
                 ! change sign of singular vectors, if desired
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_zdscal( ncvt, negone, vt( i, 1_${ik}$ ), ldvt )
              end if
           end do
           ! sort the singular values into decreasing order (insertion sort on
           ! singular values, but only one transposition per singular vector)
           do i = 1, n - 1
              ! scan for smallest d(i)
              isub = 1_${ik}$
              smin = d( 1_${ik}$ )
              do j = 2, n + 1 - i
                 if( d( j )<=smin ) then
                    isub = j
                    smin = d( j )
                 end if
              end do
              if( isub/=n+1-i ) then
                 ! swap singular values and vectors
                 d( isub ) = d( n+1-i )
                 d( n+1-i ) = smin
                 if( ncvt>0_${ik}$ )call stdlib${ii}$_zswap( ncvt, vt( isub, 1_${ik}$ ), ldvt, vt( n+1-i, 1_${ik}$ ),ldvt )
                           
                 if( nru>0_${ik}$ )call stdlib${ii}$_zswap( nru, u( 1_${ik}$, isub ), 1_${ik}$, u( 1_${ik}$, n+1-i ), 1_${ik}$ )
                 if( ncc>0_${ik}$ )call stdlib${ii}$_zswap( ncc, c( isub, 1_${ik}$ ), ldc, c( n+1-i, 1_${ik}$ ), ldc )
                           
              end if
           end do
           go to 220
           ! maximum number of iterations exceeded, failure to converge
           200 continue
           info = 0_${ik}$
           do i = 1, n - 1
              if( e( i )/=zero )info = info + 1_${ik}$
           end do
           220 continue
           return
     end subroutine stdlib${ii}$_zbdsqr

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$bdsqr( uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u,ldu, c, ldc, rwork,&
     !! ZBDSQR: computes the singular values and, optionally, the right and/or
     !! left singular vectors from the singular value decomposition (SVD) of
     !! a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
     !! zero-shift QR algorithm.  The SVD of B has the form
     !! B = Q * S * P**H
     !! where S is the diagonal matrix of singular values, Q is an orthogonal
     !! matrix of left singular vectors, and P is an orthogonal matrix of
     !! right singular vectors.  If left singular vectors are requested, this
     !! subroutine actually returns U*Q instead of Q, and, if right singular
     !! vectors are requested, this subroutine returns P**H*VT instead of
     !! P**H, for given complex input matrices U and VT.  When U and VT are
     !! the unitary matrices that reduce a general matrix A to bidiagonal
     !! form: A = U*B*VT, as computed by ZGEBRD, then
     !! A = (U*Q) * S * (P**H*VT)
     !! is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
     !! for a given complex input matrix C.
     !! See "Computing  Small Singular Values of Bidiagonal Matrices With
     !! Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
     !! LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
     !! no. 5, pp. 873-912, Sept 1990) and
     !! "Accurate singular values and differential qd algorithms," by
     !! B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
     !! Department, University of California at Berkeley, July 1992
     !! for a detailed description of the algorithm.
                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) :: uplo
           integer(${ik}$), intent(out) :: info
           integer(${ik}$), intent(in) :: ldc, ldu, ldvt, n, ncc, ncvt, nru
           ! Array Arguments 
           real(${ck}$), intent(inout) :: d(*), e(*)
           real(${ck}$), intent(out) :: rwork(*)
           complex(${ck}$), intent(inout) :: c(ldc,*), u(ldu,*), vt(ldvt,*)
        ! =====================================================================
           ! Parameters 
           real(${ck}$), parameter :: hndrth = 0.01_${ck}$
           real(${ck}$), parameter :: hndrd = 100.0_${ck}$
           real(${ck}$), parameter :: meigth = -0.125_${ck}$
           integer(${ik}$), parameter :: maxitr = 6_${ik}$
           
           
           
           
           
           
           
           
           ! Local Scalars 
           logical(lk) :: lower, rotate
           integer(${ik}$) :: i, idir, isub, iter, j, ll, lll, m, maxit, nm1, nm12, nm13, oldll, &
                     oldm
           real(${ck}$) :: abse, abss, cosl, cosr, cs, eps, f, g, h, mu, oldcs, oldsn, r, shift, &
           sigmn, sigmx, sinl, sinr, sll, smax, smin, sminl, sminoa, sn, thresh, tol, tolmul, &
                     unfl
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           lower = stdlib_lsame( uplo, 'L' )
           if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.lower ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( ncvt<0_${ik}$ ) then
              info = -3_${ik}$
           else if( nru<0_${ik}$ ) then
              info = -4_${ik}$
           else if( ncc<0_${ik}$ ) then
              info = -5_${ik}$
           else if( ( ncvt==0_${ik}$ .and. ldvt<1_${ik}$ ) .or.( ncvt>0_${ik}$ .and. ldvt<max( 1_${ik}$, n ) ) ) then
              info = -9_${ik}$
           else if( ldu<max( 1_${ik}$, nru ) ) then
              info = -11_${ik}$
           else if( ( ncc==0_${ik}$ .and. ldc<1_${ik}$ ) .or.( ncc>0_${ik}$ .and. ldc<max( 1_${ik}$, n ) ) ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZBDSQR', -info )
              return
           end if
           if( n==0 )return
           if( n==1 )go to 160
           ! rotate is true if any singular vectors desired, false otherwise
           rotate = ( ncvt>0_${ik}$ ) .or. ( nru>0_${ik}$ ) .or. ( ncc>0_${ik}$ )
           ! if no singular vectors desired, use qd algorithm
           if( .not.rotate ) then
              call stdlib${ii}$_${c2ri(ci)}$lasq1( n, d, e, rwork, info )
           ! if info equals 2, dqds didn't finish, try to finish
              if( info /= 2 ) return
              info = 0_${ik}$
           end if
           nm1 = n - 1_${ik}$
           nm12 = nm1 + nm1
           nm13 = nm12 + nm1
           idir = 0_${ik}$
           ! get machine constants
           eps = stdlib${ii}$_${c2ri(ci)}$lamch( 'EPSILON' )
           unfl = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' )
           ! if matrix lower bidiagonal, rotate to be upper bidiagonal
           ! by applying givens rotations on the left
           if( lower ) then
              do i = 1, n - 1
                 call stdlib${ii}$_${c2ri(ci)}$lartg( d( i ), e( i ), cs, sn, r )
                 d( i ) = r
                 e( i ) = sn*d( i+1 )
                 d( i+1 ) = cs*d( i+1 )
                 rwork( i ) = cs
                 rwork( nm1+i ) = sn
              end do
              ! update singular vectors if desired
              if( nru>0_${ik}$ )call stdlib${ii}$_${ci}$lasr( 'R', 'V', 'F', nru, n, rwork( 1_${ik}$ ), rwork( n ),u, ldu )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_${ci}$lasr( 'L', 'V', 'F', n, ncc, rwork( 1_${ik}$ ), rwork( n ),c, ldc )
                        
           end if
           ! compute singular values to relative accuracy tol
           ! (by setting tol to be negative, algorithm will compute
           ! singular values to absolute accuracy abs(tol)*norm(input matrix))
           tolmul = max( ten, min( hndrd, eps**meigth ) )
           tol = tolmul*eps
           ! compute approximate maximum, minimum singular values
           smax = zero
           do i = 1, n
              smax = max( smax, abs( d( i ) ) )
           end do
           do i = 1, n - 1
              smax = max( smax, abs( e( i ) ) )
           end do
           sminl = zero
           if( tol>=zero ) then
              ! relative accuracy desired
              sminoa = abs( d( 1_${ik}$ ) )
              if( sminoa==zero )go to 50
              mu = sminoa
              do i = 2, n
                 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
                 sminoa = min( sminoa, mu )
                 if( sminoa==zero )go to 50
              end do
              50 continue
              sminoa = sminoa / sqrt( real( n,KIND=${ck}$) )
              thresh = max( tol*sminoa, maxitr*n*n*unfl )
           else
              ! absolute accuracy desired
              thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
           end if
           ! prepare for main iteration loop for the singular values
           ! (maxit is the maximum number of passes through the inner
           ! loop permitted before nonconvergence signalled.)
           maxit = maxitr*n*n
           iter = 0_${ik}$
           oldll = -1_${ik}$
           oldm = -1_${ik}$
           ! m points to last element of unconverged part of matrix
           m = n
           ! begin main iteration loop
           60 continue
           ! check for convergence or exceeding iteration count
           if( m<=1 )go to 160
           if( iter>maxit )go to 200
           ! find diagonal block of matrix to work on
           if( tol<zero .and. abs( d( m ) )<=thresh )d( m ) = zero
           smax = abs( d( m ) )
           smin = smax
           do lll = 1, m - 1
              ll = m - lll
              abss = abs( d( ll ) )
              abse = abs( e( ll ) )
              if( tol<zero .and. abss<=thresh )d( ll ) = zero
              if( abse<=thresh )go to 80
              smin = min( smin, abss )
              smax = max( smax, abss, abse )
           end do
           ll = 0_${ik}$
           go to 90
           80 continue
           e( ll ) = zero
           ! matrix splits since e(ll) = 0
           if( ll==m-1 ) then
              ! convergence of bottom singular value, return to top of loop
              m = m - 1_${ik}$
              go to 60
           end if
           90 continue
           ll = ll + 1_${ik}$
           ! e(ll) through e(m-1) are nonzero, e(ll-1) is zero
           if( ll==m-1 ) then
              ! 2 by 2 block, handle separately
              call stdlib${ii}$_${c2ri(ci)}$lasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,cosr, sinl, cosl &
                        )
              d( m-1 ) = sigmx
              e( m-1 ) = zero
              d( m ) = sigmn
              ! compute singular vectors, if desired
              if( ncvt>0_${ik}$ )call stdlib${ii}$_${ci}$drot( ncvt, vt( m-1, 1_${ik}$ ), ldvt, vt( m, 1_${ik}$ ), ldvt,cosr, &
                        sinr )
              if( nru>0_${ik}$ )call stdlib${ii}$_${ci}$drot( nru, u( 1_${ik}$, m-1 ), 1_${ik}$, u( 1_${ik}$, m ), 1_${ik}$, cosl, sinl )
                        
              if( ncc>0_${ik}$ )call stdlib${ii}$_${ci}$drot( ncc, c( m-1, 1_${ik}$ ), ldc, c( m, 1_${ik}$ ), ldc, cosl,sinl )
                        
              m = m - 2_${ik}$
              go to 60
           end if
           ! if working on new submatrix, choose shift direction
           ! (from larger end diagonal element towards smaller)
           if( ll>oldm .or. m<oldll ) then
              if( abs( d( ll ) )>=abs( d( m ) ) ) then
                 ! chase bulge from top (big end) to bottom (small end)
                 idir = 1_${ik}$
              else
                 ! chase bulge from bottom (big end) to top (small end)
                 idir = 2_${ik}$
              end if
           end if
           ! apply convergence tests
           if( idir==1_${ik}$ ) then
              ! run convergence test in forward direction
              ! first apply standard test to bottom of matrix
              if( abs( e( m-1 ) )<=abs( tol )*abs( d( m ) ) .or.( tol<zero .and. abs( e( m-1 ) )&
                        <=thresh ) ) then
                 e( m-1 ) = zero
                 go to 60
              end if
              if( tol>=zero ) then
                 ! if relative accuracy desired,
                 ! apply convergence criterion forward
                 mu = abs( d( ll ) )
                 sminl = mu
                 do lll = ll, m - 1
                    if( abs( e( lll ) )<=tol*mu ) then
                       e( lll ) = zero
                       go to 60
                    end if
                    mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
                    sminl = min( sminl, mu )
                 end do
              end if
           else
              ! run convergence test in backward direction
              ! first apply standard test to top of matrix
              if( abs( e( ll ) )<=abs( tol )*abs( d( ll ) ) .or.( tol<zero .and. abs( e( ll ) )&
                        <=thresh ) ) then