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