SRC/dbdsqr.f

   1 *> \brief \b DBDSQR
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DBDSQR + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
  22 *                          LDU, C, LDC, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
  30 *      $                   VT( LDVT, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DBDSQR computes the singular values and, optionally, the right and/or
  40 *> left singular vectors from the singular value decomposition (SVD) of
  41 *> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
  42 *> zero-shift QR algorithm.  The SVD of B has the form
  43 *>
  44 *>    B = Q * S * P**T
  45 *>
  46 *> where S is the diagonal matrix of singular values, Q is an orthogonal
  47 *> matrix of left singular vectors, and P is an orthogonal matrix of
  48 *> right singular vectors.  If left singular vectors are requested, this
  49 *> subroutine actually returns U*Q instead of Q, and, if right singular
  50 *> vectors are requested, this subroutine returns P**T*VT instead of
  51 *> P**T, for given real input matrices U and VT.  When U and VT are the
  52 *> orthogonal matrices that reduce a general matrix A to bidiagonal
  53 *> form:  A = U*B*VT, as computed by DGEBRD, then
  54 *>
  55 *>    A = (U*Q) * S * (P**T*VT)
  56 *>
  57 *> is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
  58 *> for a given real input matrix C.
  59 *>
  60 *> See "Computing  Small Singular Values of Bidiagonal Matrices With
  61 *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
  62 *> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
  63 *> no. 5, pp. 873-912, Sept 1990) and
  64 *> "Accurate singular values and differential qd algorithms," by
  65 *> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
  66 *> Department, University of California at Berkeley, July 1992
  67 *> for a detailed description of the algorithm.
  68 *> \endverbatim
  69 *
  70 *  Arguments:
  71 *  ==========
  72 *
  73 *> \param[in] UPLO
  74 *> \verbatim
  75 *>          UPLO is CHARACTER*1
  76 *>          = 'U':  B is upper bidiagonal;
  77 *>          = 'L':  B is lower bidiagonal.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] N
  81 *> \verbatim
  82 *>          N is INTEGER
  83 *>          The order of the matrix B.  N >= 0.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] NCVT
  87 *> \verbatim
  88 *>          NCVT is INTEGER
  89 *>          The number of columns of the matrix VT. NCVT >= 0.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] NRU
  93 *> \verbatim
  94 *>          NRU is INTEGER
  95 *>          The number of rows of the matrix U. NRU >= 0.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] NCC
  99 *> \verbatim
 100 *>          NCC is INTEGER
 101 *>          The number of columns of the matrix C. NCC >= 0.
 102 *> \endverbatim
 103 *>
 104 *> \param[in,out] D
 105 *> \verbatim
 106 *>          D is DOUBLE PRECISION array, dimension (N)
 107 *>          On entry, the n diagonal elements of the bidiagonal matrix B.
 108 *>          On exit, if INFO=0, the singular values of B in decreasing
 109 *>          order.
 110 *> \endverbatim
 111 *>
 112 *> \param[in,out] E
 113 *> \verbatim
 114 *>          E is DOUBLE PRECISION array, dimension (N-1)
 115 *>          On entry, the N-1 offdiagonal elements of the bidiagonal
 116 *>          matrix B.
 117 *>          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
 118 *>          will contain the diagonal and superdiagonal elements of a
 119 *>          bidiagonal matrix orthogonally equivalent to the one given
 120 *>          as input.
 121 *> \endverbatim
 122 *>
 123 *> \param[in,out] VT
 124 *> \verbatim
 125 *>          VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
 126 *>          On entry, an N-by-NCVT matrix VT.
 127 *>          On exit, VT is overwritten by P**T * VT.
 128 *>          Not referenced if NCVT = 0.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] LDVT
 132 *> \verbatim
 133 *>          LDVT is INTEGER
 134 *>          The leading dimension of the array VT.
 135 *>          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
 136 *> \endverbatim
 137 *>
 138 *> \param[in,out] U
 139 *> \verbatim
 140 *>          U is DOUBLE PRECISION array, dimension (LDU, N)
 141 *>          On entry, an NRU-by-N matrix U.
 142 *>          On exit, U is overwritten by U * Q.
 143 *>          Not referenced if NRU = 0.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] LDU
 147 *> \verbatim
 148 *>          LDU is INTEGER
 149 *>          The leading dimension of the array U.  LDU >= max(1,NRU).
 150 *> \endverbatim
 151 *>
 152 *> \param[in,out] C
 153 *> \verbatim
 154 *>          C is DOUBLE PRECISION array, dimension (LDC, NCC)
 155 *>          On entry, an N-by-NCC matrix C.
 156 *>          On exit, C is overwritten by Q**T * C.
 157 *>          Not referenced if NCC = 0.
 158 *> \endverbatim
 159 *>
 160 *> \param[in] LDC
 161 *> \verbatim
 162 *>          LDC is INTEGER
 163 *>          The leading dimension of the array C.
 164 *>          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
 165 *> \endverbatim
 166 *>
 167 *> \param[out] WORK
 168 *> \verbatim
 169 *>          WORK is DOUBLE PRECISION array, dimension (4*N)
 170 *> \endverbatim
 171 *>
 172 *> \param[out] INFO
 173 *> \verbatim
 174 *>          INFO is INTEGER
 175 *>          = 0:  successful exit
 176 *>          < 0:  If INFO = -i, the i-th argument had an illegal value
 177 *>          > 0:
 178 *>             if NCVT = NRU = NCC = 0,
 179 *>                = 1, a split was marked by a positive value in E
 180 *>                = 2, current block of Z not diagonalized after 30*N
 181 *>                     iterations (in inner while loop)
 182 *>                = 3, termination criterion of outer while loop not met
 183 *>                     (program created more than N unreduced blocks)
 184 *>             else NCVT = NRU = NCC = 0,
 185 *>                   the algorithm did not converge; D and E contain the
 186 *>                   elements of a bidiagonal matrix which is orthogonally
 187 *>                   similar to the input matrix B;  if INFO = i, i
 188 *>                   elements of E have not converged to zero.
 189 *> \endverbatim
 190 *
 191 *> \par Internal Parameters:
 192 *  =========================
 193 *>
 194 *> \verbatim
 195 *>  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
 196 *>          TOLMUL controls the convergence criterion of the QR loop.
 197 *>          If it is positive, TOLMUL*EPS is the desired relative
 198 *>             precision in the computed singular values.
 199 *>          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
 200 *>             desired absolute accuracy in the computed singular
 201 *>             values (corresponds to relative accuracy
 202 *>             abs(TOLMUL*EPS) in the largest singular value.
 203 *>          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
 204 *>             between 10 (for fast convergence) and .1/EPS
 205 *>             (for there to be some accuracy in the results).
 206 *>          Default is to lose at either one eighth or 2 of the
 207 *>             available decimal digits in each computed singular value
 208 *>             (whichever is smaller).
 209 *>
 210 *>  MAXITR  INTEGER, default = 6
 211 *>          MAXITR controls the maximum number of passes of the
 212 *>          algorithm through its inner loop. The algorithms stops
 213 *>          (and so fails to converge) if the number of passes
 214 *>          through the inner loop exceeds MAXITR*N**2.
 215 *> \endverbatim
 216 *
 217 *  Authors:
 218 *  ========
 219 *
 220 *> \author Univ. of Tennessee
 221 *> \author Univ. of California Berkeley
 222 *> \author Univ. of Colorado Denver
 223 *> \author NAG Ltd.
 224 *
 225 *> \date November 2011
 226 *
 227 *> \ingroup auxOTHERcomputational
 228 *
 229 *  =====================================================================
 230       SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
 231      $                   LDU, C, LDC, WORK, INFO )
 232 *
 233 *  -- LAPACK computational routine (version 3.4.0) --
 234 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 235 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 236 *     November 2011
 237 *
 238 *     .. Scalar Arguments ..
 239       CHARACTER          UPLO
 240       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
 241 *     ..
 242 *     .. Array Arguments ..
 243       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
 244      $                   VT( LDVT, * ), WORK( * )
 245 *     ..
 246 *
 247 *  =====================================================================
 248 *
 249 *     .. Parameters ..
 250       DOUBLE PRECISION   ZERO
 251       PARAMETER          ( ZERO = 0.0D0 )
 252       DOUBLE PRECISION   ONE
 253       PARAMETER          ( ONE = 1.0D0 )
 254       DOUBLE PRECISION   NEGONE
 255       PARAMETER          ( NEGONE = -1.0D0 )
 256       DOUBLE PRECISION   HNDRTH
 257       PARAMETER          ( HNDRTH = 0.01D0 )
 258       DOUBLE PRECISION   TEN
 259       PARAMETER          ( TEN = 10.0D0 )
 260       DOUBLE PRECISION   HNDRD
 261       PARAMETER          ( HNDRD = 100.0D0 )
 262       DOUBLE PRECISION   MEIGTH
 263       PARAMETER          ( MEIGTH = -0.125D0 )
 264       INTEGER            MAXITR
 265       PARAMETER          ( MAXITR = 6 )
 266 *     ..
 267 *     .. Local Scalars ..
 268       LOGICAL            LOWER, ROTATE
 269       INTEGER            I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
 270      $                   NM12, NM13, OLDLL, OLDM
 271       DOUBLE PRECISION   ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
 272      $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
 273      $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
 274      $                   SN, THRESH, TOL, TOLMUL, UNFL
 275 *     ..
 276 *     .. External Functions ..
 277       LOGICAL            LSAME
 278       DOUBLE PRECISION   DLAMCH
 279       EXTERNAL           LSAME, DLAMCH
 280 *     ..
 281 *     .. External Subroutines ..
 282       EXTERNAL           DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT,
 283      $                   DSCAL, DSWAP, XERBLA
 284 *     ..
 285 *     .. Intrinsic Functions ..
 286       INTRINSIC          ABS, DBLE, MAX, MIN, SIGN, SQRT
 287 *     ..
 288 *     .. Executable Statements ..
 289 *
 290 *     Test the input parameters.
 291 *
 292       INFO = 0
 293       LOWER = LSAME( UPLO, 'L' )
 294       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
 295          INFO = -1
 296       ELSE IF( N.LT.0 ) THEN
 297          INFO = -2
 298       ELSE IF( NCVT.LT.0 ) THEN
 299          INFO = -3
 300       ELSE IF( NRU.LT.0 ) THEN
 301          INFO = -4
 302       ELSE IF( NCC.LT.0 ) THEN
 303          INFO = -5
 304       ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
 305      $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
 306          INFO = -9
 307       ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
 308          INFO = -11
 309       ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
 310      $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
 311          INFO = -13
 312       END IF
 313       IF( INFO.NE.0 ) THEN
 314          CALL XERBLA( 'DBDSQR', -INFO )
 315          RETURN
 316       END IF
 317       IF( N.EQ.0 )
 318      $   RETURN
 319       IF( N.EQ.1 )
 320      $   GO TO 160
 321 *
 322 *     ROTATE is true if any singular vectors desired, false otherwise
 323 *
 324       ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
 325 *
 326 *     If no singular vectors desired, use qd algorithm
 327 *
 328       IF( .NOT.ROTATE ) THEN
 329          CALL DLASQ1( N, D, E, WORK, INFO )
 330 *
 331 *     If INFO equals 2, dqds didn't finish, try to finish
 332 *
 333          IF( INFO .NE. 2 ) RETURN
 334          INFO = 0
 335       END IF
 336 *
 337       NM1 = N - 1
 338       NM12 = NM1 + NM1
 339       NM13 = NM12 + NM1
 340       IDIR = 0
 341 *
 342 *     Get machine constants
 343 *
 344       EPS = DLAMCH( 'Epsilon' )
 345       UNFL = DLAMCH( 'Safe minimum' )
 346 *
 347 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
 348 *     by applying Givens rotations on the left
 349 *
 350       IF( LOWER ) THEN
 351          DO 10 I = 1, N - 1
 352             CALL DLARTG( D( I ), E( I ), CS, SN, R )
 353             D( I ) = R
 354             E( I ) = SN*D( I+1 )
 355             D( I+1 ) = CS*D( I+1 )
 356             WORK( I ) = CS
 357             WORK( NM1+I ) = SN
 358    10    CONTINUE
 359 *
 360 *        Update singular vectors if desired
 361 *
 362          IF( NRU.GT.0 )
 363      $      CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U,
 364      $                  LDU )
 365          IF( NCC.GT.0 )
 366      $      CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C,
 367      $                  LDC )
 368       END IF
 369 *
 370 *     Compute singular values to relative accuracy TOL
 371 *     (By setting TOL to be negative, algorithm will compute
 372 *     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
 373 *
 374       TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
 375       TOL = TOLMUL*EPS
 376 *
 377 *     Compute approximate maximum, minimum singular values
 378 *
 379       SMAX = ZERO
 380       DO 20 I = 1, N
 381          SMAX = MAX( SMAX, ABS( D( I ) ) )
 382    20 CONTINUE
 383       DO 30 I = 1, N - 1
 384          SMAX = MAX( SMAX, ABS( E( I ) ) )
 385    30 CONTINUE
 386       SMINL = ZERO
 387       IF( TOL.GE.ZERO ) THEN
 388 *
 389 *        Relative accuracy desired
 390 *
 391          SMINOA = ABS( D( 1 ) )
 392          IF( SMINOA.EQ.ZERO )
 393      $      GO TO 50
 394          MU = SMINOA
 395          DO 40 I = 2, N
 396             MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
 397             SMINOA = MIN( SMINOA, MU )
 398             IF( SMINOA.EQ.ZERO )
 399      $         GO TO 50
 400    40    CONTINUE
 401    50    CONTINUE
 402          SMINOA = SMINOA / SQRT( DBLE( N ) )
 403          THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
 404       ELSE
 405 *
 406 *        Absolute accuracy desired
 407 *
 408          THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
 409       END IF
 410 *
 411 *     Prepare for main iteration loop for the singular values
 412 *     (MAXIT is the maximum number of passes through the inner
 413 *     loop permitted before nonconvergence signalled.)
 414 *
 415       MAXIT = MAXITR*N*N
 416       ITER = 0
 417       OLDLL = -1
 418       OLDM = -1
 419 *
 420 *     M points to last element of unconverged part of matrix
 421 *
 422       M = N
 423 *
 424 *     Begin main iteration loop
 425 *
 426    60 CONTINUE
 427 *
 428 *     Check for convergence or exceeding iteration count
 429 *
 430       IF( M.LE.1 )
 431      $   GO TO 160
 432       IF( ITER.GT.MAXIT )
 433      $   GO TO 200
 434 *
 435 *     Find diagonal block of matrix to work on
 436 *
 437       IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
 438      $   D( M ) = ZERO
 439       SMAX = ABS( D( M ) )
 440       SMIN = SMAX
 441       DO 70 LLL = 1, M - 1
 442          LL = M - LLL
 443          ABSS = ABS( D( LL ) )
 444          ABSE = ABS( E( LL ) )
 445          IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
 446      $      D( LL ) = ZERO
 447          IF( ABSE.LE.THRESH )
 448      $      GO TO 80
 449          SMIN = MIN( SMIN, ABSS )
 450          SMAX = MAX( SMAX, ABSS, ABSE )
 451    70 CONTINUE
 452       LL = 0
 453       GO TO 90
 454    80 CONTINUE
 455       E( LL ) = ZERO
 456 *
 457 *     Matrix splits since E(LL) = 0
 458 *
 459       IF( LL.EQ.M-1 ) THEN
 460 *
 461 *        Convergence of bottom singular value, return to top of loop
 462 *
 463          M = M - 1
 464          GO TO 60
 465       END IF
 466    90 CONTINUE
 467       LL = LL + 1
 468 *
 469 *     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
 470 *
 471       IF( LL.EQ.M-1 ) THEN
 472 *
 473 *        2 by 2 block, handle separately
 474 *
 475          CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
 476      $                COSR, SINL, COSL )
 477          D( M-1 ) = SIGMX
 478          E( M-1 ) = ZERO
 479          D( M ) = SIGMN
 480 *
 481 *        Compute singular vectors, if desired
 482 *
 483          IF( NCVT.GT.0 )
 484      $      CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,
 485      $                 SINR )
 486          IF( NRU.GT.0 )
 487      $      CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
 488          IF( NCC.GT.0 )
 489      $      CALL DROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
 490      $                 SINL )
 491          M = M - 2
 492          GO TO 60
 493       END IF
 494 *
 495 *     If working on new submatrix, choose shift direction
 496 *     (from larger end diagonal element towards smaller)
 497 *
 498       IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
 499          IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
 500 *
 501 *           Chase bulge from top (big end) to bottom (small end)
 502 *
 503             IDIR = 1
 504          ELSE
 505 *
 506 *           Chase bulge from bottom (big end) to top (small end)
 507 *
 508             IDIR = 2
 509          END IF
 510       END IF
 511 *
 512 *     Apply convergence tests
 513 *
 514       IF( IDIR.EQ.1 ) THEN
 515 *
 516 *        Run convergence test in forward direction
 517 *        First apply standard test to bottom of matrix
 518 *
 519          IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
 520      $       ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
 521             E( M-1 ) = ZERO
 522             GO TO 60
 523          END IF
 524 *
 525          IF( TOL.GE.ZERO ) THEN
 526 *
 527 *           If relative accuracy desired,
 528 *           apply convergence criterion forward
 529 *
 530             MU = ABS( D( LL ) )
 531             SMINL = MU
 532             DO 100 LLL = LL, M - 1
 533                IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
 534                   E( LLL ) = ZERO
 535                   GO TO 60
 536                END IF
 537                MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
 538                SMINL = MIN( SMINL, MU )
 539   100       CONTINUE
 540          END IF
 541 *
 542       ELSE
 543 *
 544 *        Run convergence test in backward direction
 545 *        First apply standard test to top of matrix
 546 *
 547          IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
 548      $       ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
 549             E( LL ) = ZERO
 550             GO TO 60
 551          END IF
 552 *
 553          IF( TOL.GE.ZERO ) THEN
 554 *
 555 *           If relative accuracy desired,
 556 *           apply convergence criterion backward
 557 *
 558             MU = ABS( D( M ) )
 559             SMINL = MU
 560             DO 110 LLL = M - 1, LL, -1
 561                IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
 562                   E( LLL ) = ZERO
 563                   GO TO 60
 564                END IF
 565                MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
 566                SMINL = MIN( SMINL, MU )
 567   110       CONTINUE
 568          END IF
 569       END IF
 570       OLDLL = LL
 571       OLDM = M
 572 *
 573 *     Compute shift.  First, test if shifting would ruin relative
 574 *     accuracy, and if so set the shift to zero.
 575 *
 576       IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
 577      $    MAX( EPS, HNDRTH*TOL ) ) THEN
 578 *
 579 *        Use a zero shift to avoid loss of relative accuracy
 580 *
 581          SHIFT = ZERO
 582       ELSE
 583 *
 584 *        Compute the shift from 2-by-2 block at end of matrix
 585 *
 586          IF( IDIR.EQ.1 ) THEN
 587             SLL = ABS( D( LL ) )
 588             CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
 589          ELSE
 590             SLL = ABS( D( M ) )
 591             CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
 592          END IF
 593 *
 594 *        Test if shift negligible, and if so set to zero
 595 *
 596          IF( SLL.GT.ZERO ) THEN
 597             IF( ( SHIFT / SLL )**2.LT.EPS )
 598      $         SHIFT = ZERO
 599          END IF
 600       END IF
 601 *
 602 *     Increment iteration count
 603 *
 604       ITER = ITER + M - LL
 605 *
 606 *     If SHIFT = 0, do simplified QR iteration
 607 *
 608       IF( SHIFT.EQ.ZERO ) THEN
 609          IF( IDIR.EQ.1 ) THEN
 610 *
 611 *           Chase bulge from top to bottom
 612 *           Save cosines and sines for later singular vector updates
 613 *
 614             CS = ONE
 615             OLDCS = ONE
 616             DO 120 I = LL, M - 1
 617                CALL DLARTG( D( I )*CS, E( I ), CS, SN, R )
 618                IF( I.GT.LL )
 619      $            E( I-1 ) = OLDSN*R
 620                CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
 621                WORK( I-LL+1 ) = CS
 622                WORK( I-LL+1+NM1 ) = SN
 623                WORK( I-LL+1+NM12 ) = OLDCS
 624                WORK( I-LL+1+NM13 ) = OLDSN
 625   120       CONTINUE
 626             H = D( M )*CS
 627             D( M ) = H*OLDCS
 628             E( M-1 ) = H*OLDSN
 629 *
 630 *           Update singular vectors
 631 *
 632             IF( NCVT.GT.0 )
 633      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
 634      $                     WORK( N ), VT( LL, 1 ), LDVT )
 635             IF( NRU.GT.0 )
 636      $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
 637      $                     WORK( NM13+1 ), U( 1, LL ), LDU )
 638             IF( NCC.GT.0 )
 639      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
 640      $                     WORK( NM13+1 ), C( LL, 1 ), LDC )
 641 *
 642 *           Test convergence
 643 *
 644             IF( ABS( E( M-1 ) ).LE.THRESH )
 645      $         E( M-1 ) = ZERO
 646 *
 647          ELSE
 648 *
 649 *           Chase bulge from bottom to top
 650 *           Save cosines and sines for later singular vector updates
 651 *
 652             CS = ONE
 653             OLDCS = ONE
 654             DO 130 I = M, LL + 1, -1
 655                CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
 656                IF( I.LT.M )
 657      $            E( I ) = OLDSN*R
 658                CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
 659                WORK( I-LL ) = CS
 660                WORK( I-LL+NM1 ) = -SN
 661                WORK( I-LL+NM12 ) = OLDCS
 662                WORK( I-LL+NM13 ) = -OLDSN
 663   130       CONTINUE
 664             H = D( LL )*CS
 665             D( LL ) = H*OLDCS
 666             E( LL ) = H*OLDSN
 667 *
 668 *           Update singular vectors
 669 *
 670             IF( NCVT.GT.0 )
 671      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
 672      $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )
 673             IF( NRU.GT.0 )
 674      $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
 675      $                     WORK( N ), U( 1, LL ), LDU )
 676             IF( NCC.GT.0 )
 677      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
 678      $                     WORK( N ), C( LL, 1 ), LDC )
 679 *
 680 *           Test convergence
 681 *
 682             IF( ABS( E( LL ) ).LE.THRESH )
 683      $         E( LL ) = ZERO
 684          END IF
 685       ELSE
 686 *
 687 *        Use nonzero shift
 688 *
 689          IF( IDIR.EQ.1 ) THEN
 690 *
 691 *           Chase bulge from top to bottom
 692 *           Save cosines and sines for later singular vector updates
 693 *
 694             F = ( ABS( D( LL ) )-SHIFT )*
 695      $          ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
 696             G = E( LL )
 697             DO 140 I = LL, M - 1
 698                CALL DLARTG( F, G, COSR, SINR, R )
 699                IF( I.GT.LL )
 700      $            E( I-1 ) = R
 701                F = COSR*D( I ) + SINR*E( I )
 702                E( I ) = COSR*E( I ) - SINR*D( I )
 703                G = SINR*D( I+1 )
 704                D( I+1 ) = COSR*D( I+1 )
 705                CALL DLARTG( F, G, COSL, SINL, R )
 706                D( I ) = R
 707                F = COSL*E( I ) + SINL*D( I+1 )
 708                D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
 709                IF( I.LT.M-1 ) THEN
 710                   G = SINL*E( I+1 )
 711                   E( I+1 ) = COSL*E( I+1 )
 712                END IF
 713                WORK( I-LL+1 ) = COSR
 714                WORK( I-LL+1+NM1 ) = SINR
 715                WORK( I-LL+1+NM12 ) = COSL
 716                WORK( I-LL+1+NM13 ) = SINL
 717   140       CONTINUE
 718             E( M-1 ) = F
 719 *
 720 *           Update singular vectors
 721 *
 722             IF( NCVT.GT.0 )
 723      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
 724      $                     WORK( N ), VT( LL, 1 ), LDVT )
 725             IF( NRU.GT.0 )
 726      $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
 727      $                     WORK( NM13+1 ), U( 1, LL ), LDU )
 728             IF( NCC.GT.0 )
 729      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
 730      $                     WORK( NM13+1 ), C( LL, 1 ), LDC )
 731 *
 732 *           Test convergence
 733 *
 734             IF( ABS( E( M-1 ) ).LE.THRESH )
 735      $         E( M-1 ) = ZERO
 736 *
 737          ELSE
 738 *
 739 *           Chase bulge from bottom to top
 740 *           Save cosines and sines for later singular vector updates
 741 *
 742             F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
 743      $          D( M ) )
 744             G = E( M-1 )
 745             DO 150 I = M, LL + 1, -1
 746                CALL DLARTG( F, G, COSR, SINR, R )
 747                IF( I.LT.M )
 748      $            E( I ) = R
 749                F = COSR*D( I ) + SINR*E( I-1 )
 750                E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
 751                G = SINR*D( I-1 )
 752                D( I-1 ) = COSR*D( I-1 )
 753                CALL DLARTG( F, G, COSL, SINL, R )
 754                D( I ) = R
 755                F = COSL*E( I-1 ) + SINL*D( I-1 )
 756                D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
 757                IF( I.GT.LL+1 ) THEN
 758                   G = SINL*E( I-2 )
 759                   E( I-2 ) = COSL*E( I-2 )
 760                END IF
 761                WORK( I-LL ) = COSR
 762                WORK( I-LL+NM1 ) = -SINR
 763                WORK( I-LL+NM12 ) = COSL
 764                WORK( I-LL+NM13 ) = -SINL
 765   150       CONTINUE
 766             E( LL ) = F
 767 *
 768 *           Test convergence
 769 *
 770             IF( ABS( E( LL ) ).LE.THRESH )
 771      $         E( LL ) = ZERO
 772 *
 773 *           Update singular vectors if desired
 774 *
 775             IF( NCVT.GT.0 )
 776      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
 777      $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )
 778             IF( NRU.GT.0 )
 779      $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
 780      $                     WORK( N ), U( 1, LL ), LDU )
 781             IF( NCC.GT.0 )
 782      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
 783      $                     WORK( N ), C( LL, 1 ), LDC )
 784          END IF
 785       END IF
 786 *
 787 *     QR iteration finished, go back and check convergence
 788 *
 789       GO TO 60
 790 *
 791 *     All singular values converged, so make them positive
 792 *
 793   160 CONTINUE
 794       DO 170 I = 1, N
 795          IF( D( I ).LT.ZERO ) THEN
 796             D( I ) = -D( I )
 797 *
 798 *           Change sign of singular vectors, if desired
 799 *
 800             IF( NCVT.GT.0 )
 801      $         CALL DSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
 802          END IF
 803   170 CONTINUE
 804 *
 805 *     Sort the singular values into decreasing order (insertion sort on
 806 *     singular values, but only one transposition per singular vector)
 807 *
 808       DO 190 I = 1, N - 1
 809 *
 810 *        Scan for smallest D(I)
 811 *
 812          ISUB = 1
 813          SMIN = D( 1 )
 814          DO 180 J = 2, N + 1 - I
 815             IF( D( J ).LE.SMIN ) THEN
 816                ISUB = J
 817                SMIN = D( J )
 818             END IF
 819   180    CONTINUE
 820          IF( ISUB.NE.N+1-I ) THEN
 821 *
 822 *           Swap singular values and vectors
 823 *
 824             D( ISUB ) = D( N+1-I )
 825             D( N+1-I ) = SMIN
 826             IF( NCVT.GT.0 )
 827      $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
 828      $                     LDVT )
 829             IF( NRU.GT.0 )
 830      $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
 831             IF( NCC.GT.0 )
 832      $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
 833          END IF
 834   190 CONTINUE
 835       GO TO 220
 836 *
 837 *     Maximum number of iterations exceeded, failure to converge
 838 *
 839   200 CONTINUE
 840       INFO = 0
 841       DO 210 I = 1, N - 1
 842          IF( E( I ).NE.ZERO )
 843      $      INFO = INFO + 1
 844   210 CONTINUE
 845   220 CONTINUE
 846       RETURN
 847 *
 848 *     End of DBDSQR
 849 *
 850       END