cbdsqr()

subroutine cbdsqr	(	character	uplo,
		integer	n,
		integer	ncvt,
		integer	nru,
		integer	ncc,
		real, dimension( * )	d,
		real, dimension( * )	e,
		complex, dimension( ldvt, * )	vt,
		integer	ldvt,
		complex, dimension( ldu, * )	u,
		integer	ldu,
		complex, dimension( ldc, * )	c,
		integer	ldc,
		real, dimension( * )	rwork,
		integer	info )

CBDSQR

Download CBDSQR + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> 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  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
!>  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.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> = 'U': B is upper bidiagonal; !> = 'L': B is lower bidiagonal. !>
[in]	N	!> N is INTEGER !> The order of the matrix B. N >= 0. !>
[in]	NCVT	!> NCVT is INTEGER !> The number of columns of the matrix VT. NCVT >= 0. !>
[in]	NRU	!> NRU is INTEGER !> The number of rows of the matrix U. NRU >= 0. !>
[in]	NCC	!> NCC is INTEGER !> The number of columns of the matrix C. NCC >= 0. !>
[in,out]	D	!> D is REAL array, dimension (N) !> On entry, the n diagonal elements of the bidiagonal matrix B. !> On exit, if INFO=0, the singular values of B in decreasing !> order. !>
[in,out]	E	!> E is REAL array, dimension (N-1) !> On entry, the N-1 offdiagonal elements of the bidiagonal !> matrix B. !> On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E !> will contain the diagonal and superdiagonal elements of a !> bidiagonal matrix orthogonally equivalent to the one given !> as input. !>
[in,out]	VT	!> VT is COMPLEX array, dimension (LDVT, NCVT) !> On entry, an N-by-NCVT matrix VT. !> On exit, VT is overwritten by P**H * VT. !> Not referenced if NCVT = 0. !>
[in]	LDVT	!> LDVT is INTEGER !> The leading dimension of the array VT. !> LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. !>
[in,out]	U	!> U is COMPLEX array, dimension (LDU, N) !> On entry, an NRU-by-N matrix U. !> On exit, U is overwritten by U * Q. !> Not referenced if NRU = 0. !>
[in]	LDU	!> LDU is INTEGER !> The leading dimension of the array U. LDU >= max(1,NRU). !>
[in,out]	C	!> C is COMPLEX array, dimension (LDC, NCC) !> On entry, an N-by-NCC matrix C. !> On exit, C is overwritten by Q**H * C. !> Not referenced if NCC = 0. !>
[in]	LDC	!> LDC is INTEGER !> The leading dimension of the array C. !> LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. !>
[out]	RWORK	!> RWORK is REAL array, dimension (LRWORK) !> LRWORK = 4*N, if NCVT = NRU = NCC = 0, and !> LRWORK = 4*(N-1), otherwise !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: If INFO = -i, the i-th argument had an illegal value !> > 0: the algorithm did not converge; D and E contain the !> elements of a bidiagonal matrix which is orthogonally !> similar to the input matrix B; if INFO = i, i !> elements of E have not converged to zero. !>

Internal Parameters:

!>  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
!>          TOLMUL controls the convergence criterion of the QR loop.
!>          If it is positive, TOLMUL*EPS is the desired relative
!>             precision in the computed singular values.
!>          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
!>             desired absolute accuracy in the computed singular
!>             values (corresponds to relative accuracy
!>             abs(TOLMUL*EPS) in the largest singular value.
!>          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
!>             between 10 (for fast convergence) and .1/EPS
!>             (for there to be some accuracy in the results).
!>          Default is to lose at either one eighth or 2 of the
!>             available decimal digits in each computed singular value
!>             (whichever is smaller).
!>
!>  MAXITR  INTEGER, default = 6
!>          MAXITR controls the maximum number of passes of the
!>          algorithm through its inner loop. The algorithms stops
!>          (and so fails to converge) if the number of passes
!>          through the inner loop exceeds MAXITR*N**2.
!>
!> 

Note:

!>  Bug report from Cezary Dendek.
!>  On November 3rd 2023, the INTEGER variable MAXIT = MAXITR*N**2 is
!>  removed since it can overflow pretty easily (for N larger or equal
!>  than 18,919). We instead use MAXITDIVN = MAXITR*N.
!> 

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 231 of file cbdsqr.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
*     ..
*     .. Array Arguments ..
      REAL               D( * ), E( * ), RWORK( * )
      COMPLEX            C( LDC, * ), U( LDU, * ), VT( LDVT, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      parameter( zero = 0.0e0 )
      REAL               ONE
      parameter( one = 1.0e0 )
      REAL               NEGONE
      parameter( negone = -1.0e0 )
      REAL               HNDRTH
      parameter( hndrth = 0.01e0 )
      REAL               TEN
      parameter( ten = 10.0e0 )
      REAL               HNDRD
      parameter( hndrd = 100.0e0 )
      REAL               MEIGTH
      parameter( meigth = -0.125e0 )
      INTEGER            MAXITR
      parameter( maxitr = 6 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LOWER, ROTATE
      INTEGER            I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
     $                   MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM
      REAL               ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
     $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
     $                   SINR, SLL, SMAX, SMIN, SMINOA,
     $                   SN, THRESH, TOL, TOLMUL, UNFL
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH
      EXTERNAL           lsame, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           clasr, csrot, csscal, cswap, slartg,
     $                   slas2,
     $                   slasq1, slasv2, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, real, sign, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      lower = lsame( uplo, 'L' )
      IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( ncvt.LT.0 ) THEN
         info = -3
      ELSE IF( nru.LT.0 ) THEN
         info = -4
      ELSE IF( ncc.LT.0 ) THEN
         info = -5
      ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
     $         ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
         info = -9
      ELSE IF( ldu.LT.max( 1, nru ) ) THEN
         info = -11
      ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
     $         ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
         info = -13
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CBDSQR', -info )
         RETURN
      END IF
      IF( n.EQ.0 )
     $   RETURN
      IF( n.EQ.1 )
     $   GO TO 160
*
*     ROTATE is true if any singular vectors desired, false otherwise
*
      rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
*
*     If no singular vectors desired, use qd algorithm
*
      IF( .NOT.rotate ) THEN
         CALL slasq1( n, d, e, rwork, info )
*
*     If INFO equals 2, dqds didn't finish, try to finish
*
         IF( info .NE. 2 ) RETURN
         info = 0
      END IF
*
      nm1 = n - 1
      nm12 = nm1 + nm1
      nm13 = nm12 + nm1
      idir = 0
*
*     Get machine constants
*
      eps = slamch( 'Epsilon' )
      unfl = slamch( 'Safe minimum' )
*
*     If matrix lower bidiagonal, rotate to be upper bidiagonal
*     by applying Givens rotations on the left
*
      IF( lower ) THEN
         DO 10 i = 1, n - 1
            CALL 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
   10    CONTINUE
*
*        Update singular vectors if desired
*
         IF( nru.GT.0 )
     $      CALL clasr( 'R', 'V', 'F', nru, n, rwork( 1 ),
     $                  rwork( n ),
     $                  u, ldu )
         IF( ncc.GT.0 )
     $      CALL clasr( 'L', 'V', 'F', n, ncc, rwork( 1 ),
     $                  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 20 i = 1, n
         smax = max( smax, abs( d( i ) ) )
   20 CONTINUE
      DO 30 i = 1, n - 1
         smax = max( smax, abs( e( i ) ) )
   30 CONTINUE
      smin = zero
      IF( tol.GE.zero ) THEN
*
*        Relative accuracy desired
*
         sminoa = abs( d( 1 ) )
         IF( sminoa.EQ.zero )
     $      GO TO 50
         mu = sminoa
         DO 40 i = 2, n
            mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
            sminoa = min( sminoa, mu )
            IF( sminoa.EQ.zero )
     $         GO TO 50
   40    CONTINUE
   50    CONTINUE
         sminoa = sminoa / sqrt( real( n ) )
         thresh = max( tol*sminoa,
     $                 real(maxitr)*(real(n)*(real(n)*unfl)) )
      ELSE
*
*        Absolute accuracy desired
*
         thresh = max( abs( tol )*smax,
     $                 real(maxitr)*(real(n)*(real(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
      iter = -1
      oldll = -1
      oldm = -1
*
*     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.LE.1 )
     $   GO TO 160
      IF( iter.GE.n ) THEN
         iter = iter - n
         iterdivn = iterdivn + 1
         IF( iterdivn.GE.maxitdivn )
     $      GO TO 200
      END IF
*
*     Find diagonal block of matrix to work on
*
      IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
     $   d( m ) = zero
      smax = abs( d( m ) )
      DO 70 lll = 1, m - 1
         ll = m - lll
         abss = abs( d( ll ) )
         abse = abs( e( ll ) )
         IF( tol.LT.zero .AND. abss.LE.thresh )
     $      d( ll ) = zero
         IF( abse.LE.thresh )
     $      GO TO 80
         smax = max( smax, abss, abse )
   70 CONTINUE
      ll = 0
      GO TO 90
   80 CONTINUE
      e( ll ) = zero
*
*     Matrix splits since E(LL) = 0
*
      IF( ll.EQ.m-1 ) THEN
*
*        Convergence of bottom singular value, return to top of loop
*
         m = m - 1
         GO TO 60
      END IF
   90 CONTINUE
      ll = ll + 1
*
*     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
*
      IF( ll.EQ.m-1 ) THEN
*
*        2 by 2 block, handle separately
*
         CALL 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.GT.0 )
     $      CALL csrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
     $                  cosr, sinr )
         IF( nru.GT.0 )
     $      CALL csrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl,
     $                  sinl )
         IF( ncc.GT.0 )
     $      CALL csrot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
     $                  sinl )
         m = m - 2
         GO TO 60
      END IF
*
*     If working on new submatrix, choose shift direction
*     (from larger end diagonal element towards smaller)
*
      IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
         IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
*
*           Chase bulge from top (big end) to bottom (small end)
*
            idir = 1
         ELSE
*
*           Chase bulge from bottom (big end) to top (small end)
*
            idir = 2
         END IF
      END IF
*
*     Apply convergence tests
*
      IF( idir.EQ.1 ) THEN
*
*        Run convergence test in forward direction
*        First apply standard test to bottom of matrix
*
         IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
     $       ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
            e( m-1 ) = zero
            GO TO 60
         END IF
*
         IF( tol.GE.zero ) THEN
*
*           If relative accuracy desired,
*           apply convergence criterion forward
*
            mu = abs( d( ll ) )
            smin = mu
            DO 100 lll = ll, m - 1
               IF( abs( e( lll ) ).LE.tol*mu ) THEN
                  e( lll ) = zero
                  GO TO 60
               END IF
               mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
               smin = min( smin, mu )
  100       CONTINUE
         END IF
*
      ELSE
*
*        Run convergence test in backward direction
*        First apply standard test to top of matrix
*
         IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
     $       ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
            e( ll ) = zero
            GO TO 60
         END IF
*
         IF( tol.GE.zero ) THEN
*
*           If relative accuracy desired,
*           apply convergence criterion backward
*
            mu = abs( d( m ) )
            smin = mu
            DO 110 lll = m - 1, ll, -1
               IF( abs( e( lll ) ).LE.tol*mu ) THEN
                  e( lll ) = zero
                  GO TO 60
               END IF
               mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
               smin = min( smin, mu )
  110       CONTINUE
         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.GE.zero .AND. real( n )*tol*( smin / smax ).LE.
     $    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.EQ.1 ) THEN
            sll = abs( d( ll ) )
            CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
         ELSE
            sll = abs( d( m ) )
            CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
         END IF
*
*        Test if shift negligible, and if so set to zero
*
         IF( sll.GT.zero ) THEN
            IF( ( shift / sll )**2.LT.eps )
     $         shift = zero
         END IF
      END IF
*
*     Increment iteration count
*
      iter = iter + m - ll
*
*     If SHIFT = 0, do simplified QR iteration
*
      IF( shift.EQ.zero ) THEN
         IF( idir.EQ.1 ) THEN
*
*           Chase bulge from top to bottom
*           Save cosines and sines for later singular vector updates
*
            cs = one
            oldcs = one
            DO 120 i = ll, m - 1
               CALL slartg( d( i )*cs, e( i ), cs, sn, r )
               IF( i.GT.ll )
     $            e( i-1 ) = oldsn*r
               CALL 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
  120       CONTINUE
            h = d( m )*cs
            d( m ) = h*oldcs
            e( m-1 ) = h*oldsn
*
*           Update singular vectors
*
            IF( ncvt.GT.0 )
     $         CALL clasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
     $                     rwork( n ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL clasr( 'R', 'V', 'F', nru, m-ll+1,
     $                     rwork( nm12+1 ),
     $                     rwork( nm13+1 ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL clasr( 'L', 'V', 'F', m-ll+1, ncc,
     $                     rwork( nm12+1 ),
     $                     rwork( nm13+1 ), c( ll, 1 ), ldc )
*
*           Test convergence
*
            IF( abs( e( m-1 ) ).LE.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 130 i = m, ll + 1, -1
               CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
               IF( i.LT.m )
     $            e( i ) = oldsn*r
               CALL 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
  130       CONTINUE
            h = d( ll )*cs
            d( ll ) = h*oldcs
            e( ll ) = h*oldsn
*
*           Update singular vectors
*
            IF( ncvt.GT.0 )
     $         CALL clasr( 'L', 'V', 'B', m-ll+1, ncvt,
     $                     rwork( nm12+1 ),
     $                     rwork( nm13+1 ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL clasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
     $                     rwork( n ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL clasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1 ),
     $                     rwork( n ), c( ll, 1 ), ldc )
*
*           Test convergence
*
            IF( abs( e( ll ) ).LE.thresh )
     $         e( ll ) = zero
         END IF
      ELSE
*
*        Use nonzero shift
*
         IF( idir.EQ.1 ) 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 140 i = ll, m - 1
               CALL slartg( f, g, cosr, sinr, r )
               IF( i.GT.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 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.LT.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
  140       CONTINUE
            e( m-1 ) = f
*
*           Update singular vectors
*
            IF( ncvt.GT.0 )
     $         CALL clasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
     $                     rwork( n ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL clasr( 'R', 'V', 'F', nru, m-ll+1,
     $                     rwork( nm12+1 ),
     $                     rwork( nm13+1 ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL clasr( 'L', 'V', 'F', m-ll+1, ncc,
     $                     rwork( nm12+1 ),
     $                     rwork( nm13+1 ), c( ll, 1 ), ldc )
*
*           Test convergence
*
            IF( abs( e( m-1 ) ).LE.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 150 i = m, ll + 1, -1
               CALL slartg( f, g, cosr, sinr, r )
               IF( i.LT.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 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.GT.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
  150       CONTINUE
            e( ll ) = f
*
*           Test convergence
*
            IF( abs( e( ll ) ).LE.thresh )
     $         e( ll ) = zero
*
*           Update singular vectors if desired
*
            IF( ncvt.GT.0 )
     $         CALL clasr( 'L', 'V', 'B', m-ll+1, ncvt,
     $                     rwork( nm12+1 ),
     $                     rwork( nm13+1 ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL clasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
     $                     rwork( n ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL clasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1 ),
     $                     rwork( n ), c( ll, 1 ), 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 170 i = 1, n
         IF( d( i ).EQ.zero ) THEN
*
*           Avoid -ZERO
*
            d( i ) = zero
         END IF
         IF( d( i ).LT.zero ) THEN
            d( i ) = -d( i )
*
*           Change sign of singular vectors, if desired
*
            IF( ncvt.GT.0 )
     $         CALL csscal( ncvt, negone, vt( i, 1 ), ldvt )
         END IF
  170 CONTINUE
*
*     Sort the singular values into decreasing order (insertion sort on
*     singular values, but only one transposition per singular vector)
*
      DO 190 i = 1, n - 1
*
*        Scan for smallest D(I)
*
         isub = 1
         smin = d( 1 )
         DO 180 j = 2, n + 1 - i
            IF( d( j ).LE.smin ) THEN
               isub = j
               smin = d( j )
            END IF
  180    CONTINUE
         IF( isub.NE.n+1-i ) THEN
*
*           Swap singular values and vectors
*
            d( isub ) = d( n+1-i )
            d( n+1-i ) = smin
            IF( ncvt.GT.0 )
     $         CALL cswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
     $                     ldvt )
            IF( nru.GT.0 )
     $         CALL cswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
            IF( ncc.GT.0 )
     $         CALL cswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ),
     $                     ldc )
         END IF
  190 CONTINUE
      GO TO 220
*
*     Maximum number of iterations exceeded, failure to converge
*
  200 CONTINUE
      info = 0
      DO 210 i = 1, n - 1
         IF( e( i ).NE.zero )
     $      info = info + 1
  210 CONTINUE
  220 CONTINUE
      RETURN
*
*     End of CBDSQR
*

Here is the call graph for this function:

Here is the caller graph for this function: