dbdsqr.f

Go to the documentation of this file.

*> \brief \b DBDSQR
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download DBDSQR + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*       SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
*                          LDU, C, LDC, WORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
*      $                   VT( LDVT, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  B is upper bidiagonal;
*>          = 'L':  B is lower bidiagonal.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix B.  N >= 0.
*> \endverbatim
*>
*> \param[in] NCVT
*> \verbatim
*>          NCVT is INTEGER
*>          The number of columns of the matrix VT. NCVT >= 0.
*> \endverbatim
*>
*> \param[in] NRU
*> \verbatim
*>          NRU is INTEGER
*>          The number of rows of the matrix U. NRU >= 0.
*> \endverbatim
*>
*> \param[in] NCC
*> \verbatim
*>          NCC is INTEGER
*>          The number of columns of the matrix C. NCC >= 0.
*> \endverbatim
*>
*> \param[in,out] D
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in,out] E
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in,out] VT
*> \verbatim
*>          VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
*>          On entry, an N-by-NCVT matrix VT.
*>          On exit, VT is overwritten by P**T * VT.
*>          Not referenced if NCVT = 0.
*> \endverbatim
*>
*> \param[in] LDVT
*> \verbatim
*>          LDVT is INTEGER
*>          The leading dimension of the array VT.
*>          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
*> \endverbatim
*>
*> \param[in,out] U
*> \verbatim
*>          U is DOUBLE PRECISION 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.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of the array U.  LDU >= max(1,NRU).
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*>          C is DOUBLE PRECISION array, dimension (LDC, NCC)
*>          On entry, an N-by-NCC matrix C.
*>          On exit, C is overwritten by Q**T * C.
*>          Not referenced if NCC = 0.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>          The leading dimension of the array C.
*>          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (LWORK)
*>          LWORK = 4*N, if NCVT = NRU = NCC = 0, and
*>          LWORK = 4*(N-1), otherwise
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  If INFO = -i, the i-th argument had an illegal value
*>          > 0:
*>             if NCVT = NRU = NCC = 0,
*>                = 1, a split was marked by a positive value in E
*>                = 2, current block of Z not diagonalized after 30*N
*>                     iterations (in inner while loop)
*>                = 3, termination criterion of outer while loop not met
*>                     (program created more than N unreduced blocks)
*>             else NCVT = NRU = NCC = 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.
*> \endverbatim
*
*> \par Internal Parameters:
*  =========================
*>
*> \verbatim
*>  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.
*>
*> \endverbatim
*
*> \par Note:
*  ===========
*>
*> \verbatim
*>  Bug report from Cezary Dendek.
*>  On March 23rd 2017, 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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup bdsqr
*
*  =====================================================================
      SUBROUTINE dbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
     $                   LDU, C, LDC, WORK, 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..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
     $                   vt( ldvt, * ), work( * )
*     ..
*
*  =====================================================================
*
*     .. 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, dlasr, dlasv2,
     $                   drot,
     $                   dscal, dswap, xerbla
*     ..
*     .. 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( 'DBDSQR', -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, work, 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 )
            work( i ) = cs
            work( nm1+i ) = sn
   10    CONTINUE
*
*        Update singular vectors if desired
*
         IF( nru.GT.0 )
     $      CALL dlasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ),
     $                  u,
     $                  ldu )
         IF( ncc.GT.0 )
     $      CALL dlasr( 'L', 'V', 'F', n, ncc, work( 1 ), 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 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 drot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
     $                 cosr,
     $                 sinr )
         IF( nru.GT.0 )
     $      CALL drot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl,
     $                 sinl )
         IF( ncc.GT.0 )
     $      CALL drot( 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 ) )
               work( i-ll+1 ) = cs
               work( i-ll+1+nm1 ) = sn
               work( i-ll+1+nm12 ) = oldcs
               work( 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 dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
     $                     work( n ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL dlasr( 'R', 'V', 'F', nru, m-ll+1,
     $                     work( nm12+1 ),
     $                     work( nm13+1 ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc,
     $                     work( nm12+1 ),
     $                     work( 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 ) )
               work( i-ll ) = cs
               work( i-ll+nm1 ) = -sn
               work( i-ll+nm12 ) = oldcs
               work( 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 dlasr( 'L', 'V', 'B', m-ll+1, ncvt,
     $                     work( nm12+1 ),
     $                     work( nm13+1 ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
     $                     work( n ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
     $                     work( 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
               work( i-ll+1 ) = cosr
               work( i-ll+1+nm1 ) = sinr
               work( i-ll+1+nm12 ) = cosl
               work( i-ll+1+nm13 ) = sinl
  140       CONTINUE
            e( m-1 ) = f
*
*           Update singular vectors
*
            IF( ncvt.GT.0 )
     $         CALL dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
     $                     work( n ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL dlasr( 'R', 'V', 'F', nru, m-ll+1,
     $                     work( nm12+1 ),
     $                     work( nm13+1 ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc,
     $                     work( nm12+1 ),
     $                     work( 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
               work( i-ll ) = cosr
               work( i-ll+nm1 ) = -sinr
               work( i-ll+nm12 ) = cosl
               work( 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 dlasr( 'L', 'V', 'B', m-ll+1, ncvt,
     $                     work( nm12+1 ),
     $                     work( nm13+1 ), vt( ll, 1 ), ldvt )
            IF( nru.GT.0 )
     $         CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
     $                     work( n ), u( 1, ll ), ldu )
            IF( ncc.GT.0 )
     $         CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
     $                     work( 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 dscal( 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 dswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
     $                     ldvt )
            IF( nru.GT.0 )
     $         CALL dswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
            IF( ncc.GT.0 )
     $         CALL dswap( 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 DBDSQR
*
      END