zbdsqr()

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

ZBDSQR

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

Purpose:

 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.

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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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*16 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 DOUBLE PRECISION 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  DOUBLE PRECISION, 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 zbdsqr.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 ..
      DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
      COMPLEX*16         C( LDC, * ), U( LDU, * ), VT( LDVT, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d0 )
      DOUBLE PRECISION   ONE
      parameter( one = 1.0d0 )
      DOUBLE PRECISION   NEGONE
      parameter( negone = -1.0d0 )
      DOUBLE PRECISION   HNDRTH
      parameter( hndrth = 0.01d0 )
      DOUBLE PRECISION   TEN
      parameter( ten = 10.0d0 )
      DOUBLE PRECISION   HNDRD
      parameter( hndrd = 100.0d0 )
      DOUBLE PRECISION   MEIGTH
      parameter( meigth = -0.125d0 )
      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
      DOUBLE PRECISION   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
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           lsame, dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlartg, dlas2, dlasq1, dlasv2, xerbla, zdrot,
     $                   zdscal, zlasr, zswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, max, min, 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( 'ZBDSQR', -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 dlasq1( 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 = dlamch( 'Epsilon' )
      unfl = dlamch( '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 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
   10    CONTINUE
*
*        Update singular vectors if desired
*
         IF( nru.GT.0 )
     $      CALL zlasr( 'R', 'V', 'F', nru, n, rwork( 1 ), rwork( n ),
     $                  u, ldu )
         IF( ncc.GT.0 )
     $      CALL zlasr( '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( dble( 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 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.GT.0 )
     $      CALL zdrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
     $                  cosr, sinr )
         IF( nru.GT.0 )
     $      CALL zdrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
         IF( ncc.GT.0 )
     $      CALL zdrot( 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 dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
         ELSE
            sll = abs( d( m ) )
            CALL dlas2( 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 dlartg( d( i )*cs, e( i ), cs, sn, r )
               IF( i.GT.ll )
     $            e( i-1 ) = oldsn*r
               CALL 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
  120       CONTINUE
            h = d( m )*cs
            d( m ) = h*oldcs
            e( m-1 ) = h*oldsn
*
*           Update singular vectors
*
            IF( ncvt.GT.0 )
     $         CALL zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
     $                     rwork( n ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL zlasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),
     $                     rwork( nm13+1 ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL zlasr( '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 dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
               IF( i.LT.m )
     $            e( i ) = oldsn*r
               CALL 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
  130       CONTINUE
            h = d( ll )*cs
            d( ll ) = h*oldcs
            e( ll ) = h*oldsn
*
*           Update singular vectors
*
            IF( ncvt.GT.0 )
     $         CALL zlasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),
     $                     rwork( nm13+1 ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
     $                     rwork( n ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL zlasr( '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 dlartg( 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 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.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 zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
     $                     rwork( n ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL zlasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),
     $                     rwork( nm13+1 ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL zlasr( '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 dlartg( 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 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.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 zlasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),
     $                     rwork( nm13+1 ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
     $                     rwork( n ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL zlasr( '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 zdscal( 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 zswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
     $                     ldvt )
            IF( nru.GT.0 )
     $         CALL zswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
            IF( ncc.GT.0 )
     $         CALL zswap( 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 ZBDSQR
*

Here is the call graph for this function:

Here is the caller graph for this function: