SRC/zbdsqr.f

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