#: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
