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

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 (4*N)
[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, 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
      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. 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 ).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: