LAPACK: SRC/dbdsqr.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
00002      $                   LDU, C, LDC, WORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     January 2007
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          UPLO
00011       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
00015      $                   VT( LDVT, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DBDSQR computes the singular values and, optionally, the right and/or
00022 *  left singular vectors from the singular value decomposition (SVD) of
00023 *  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
00024 *  zero-shift QR algorithm.  The SVD of B has the form
00025 * 
00026 *     B = Q * S * P**T
00027 * 
00028 *  where S is the diagonal matrix of singular values, Q is an orthogonal
00029 *  matrix of left singular vectors, and P is an orthogonal matrix of
00030 *  right singular vectors.  If left singular vectors are requested, this
00031 *  subroutine actually returns U*Q instead of Q, and, if right singular
00032 *  vectors are requested, this subroutine returns P**T*VT instead of
00033 *  P**T, for given real input matrices U and VT.  When U and VT are the
00034 *  orthogonal matrices that reduce a general matrix A to bidiagonal
00035 *  form:  A = U*B*VT, as computed by DGEBRD, then
00036 *
00037 *     A = (U*Q) * S * (P**T*VT)
00038 *
00039 *  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
00040 *  for a given real input matrix C.
00041 *
00042 *  See "Computing  Small Singular Values of Bidiagonal Matrices With
00043 *  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
00044 *  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
00045 *  no. 5, pp. 873-912, Sept 1990) and
00046 *  "Accurate singular values and differential qd algorithms," by
00047 *  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
00048 *  Department, University of California at Berkeley, July 1992
00049 *  for a detailed description of the algorithm.
00050 *
00051 *  Arguments
00052 *  =========
00053 *
00054 *  UPLO    (input) CHARACTER*1
00055 *          = 'U':  B is upper bidiagonal;
00056 *          = 'L':  B is lower bidiagonal.
00057 *
00058 *  N       (input) INTEGER
00059 *          The order of the matrix B.  N >= 0.
00060 *
00061 *  NCVT    (input) INTEGER
00062 *          The number of columns of the matrix VT. NCVT >= 0.
00063 *
00064 *  NRU     (input) INTEGER
00065 *          The number of rows of the matrix U. NRU >= 0.
00066 *
00067 *  NCC     (input) INTEGER
00068 *          The number of columns of the matrix C. NCC >= 0.
00069 *
00070 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00071 *          On entry, the n diagonal elements of the bidiagonal matrix B.
00072 *          On exit, if INFO=0, the singular values of B in decreasing
00073 *          order.
00074 *
00075 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
00076 *          On entry, the N-1 offdiagonal elements of the bidiagonal
00077 *          matrix B. 
00078 *          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
00079 *          will contain the diagonal and superdiagonal elements of a
00080 *          bidiagonal matrix orthogonally equivalent to the one given
00081 *          as input.
00082 *
00083 *  VT      (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
00084 *          On entry, an N-by-NCVT matrix VT.
00085 *          On exit, VT is overwritten by P**T * VT.
00086 *          Not referenced if NCVT = 0.
00087 *
00088 *  LDVT    (input) INTEGER
00089 *          The leading dimension of the array VT.
00090 *          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
00091 *
00092 *  U       (input/output) DOUBLE PRECISION array, dimension (LDU, N)
00093 *          On entry, an NRU-by-N matrix U.
00094 *          On exit, U is overwritten by U * Q.
00095 *          Not referenced if NRU = 0.
00096 *
00097 *  LDU     (input) INTEGER
00098 *          The leading dimension of the array U.  LDU >= max(1,NRU).
00099 *
00100 *  C       (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
00101 *          On entry, an N-by-NCC matrix C.
00102 *          On exit, C is overwritten by Q**T * C.
00103 *          Not referenced if NCC = 0.
00104 *
00105 *  LDC     (input) INTEGER
00106 *          The leading dimension of the array C.
00107 *          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
00108 *
00109 *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
00110 *
00111 *  INFO    (output) INTEGER
00112 *          = 0:  successful exit
00113 *          < 0:  If INFO = -i, the i-th argument had an illegal value
00114 *          > 0:
00115 *             if NCVT = NRU = NCC = 0,
00116 *                = 1, a split was marked by a positive value in E
00117 *                = 2, current block of Z not diagonalized after 30*N
00118 *                     iterations (in inner while loop)
00119 *                = 3, termination criterion of outer while loop not met 
00120 *                     (program created more than N unreduced blocks)
00121 *             else NCVT = NRU = NCC = 0,
00122 *                   the algorithm did not converge; D and E contain the
00123 *                   elements of a bidiagonal matrix which is orthogonally
00124 *                   similar to the input matrix B;  if INFO = i, i
00125 *                   elements of E have not converged to zero.
00126 *
00127 *  Internal Parameters
00128 *  ===================
00129 *
00130 *  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
00131 *          TOLMUL controls the convergence criterion of the QR loop.
00132 *          If it is positive, TOLMUL*EPS is the desired relative
00133 *             precision in the computed singular values.
00134 *          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
00135 *             desired absolute accuracy in the computed singular
00136 *             values (corresponds to relative accuracy
00137 *             abs(TOLMUL*EPS) in the largest singular value.
00138 *          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
00139 *             between 10 (for fast convergence) and .1/EPS
00140 *             (for there to be some accuracy in the results).
00141 *          Default is to lose at either one eighth or 2 of the
00142 *             available decimal digits in each computed singular value
00143 *             (whichever is smaller).
00144 *
00145 *  MAXITR  INTEGER, default = 6
00146 *          MAXITR controls the maximum number of passes of the
00147 *          algorithm through its inner loop. The algorithms stops
00148 *          (and so fails to converge) if the number of passes
00149 *          through the inner loop exceeds MAXITR*N**2.
00150 *
00151 *  =====================================================================
00152 *
00153 *     .. Parameters ..
00154       DOUBLE PRECISION   ZERO
00155       PARAMETER          ( ZERO = 0.0D0 )
00156       DOUBLE PRECISION   ONE
00157       PARAMETER          ( ONE = 1.0D0 )
00158       DOUBLE PRECISION   NEGONE
00159       PARAMETER          ( NEGONE = -1.0D0 )
00160       DOUBLE PRECISION   HNDRTH
00161       PARAMETER          ( HNDRTH = 0.01D0 )
00162       DOUBLE PRECISION   TEN
00163       PARAMETER          ( TEN = 10.0D0 )
00164       DOUBLE PRECISION   HNDRD
00165       PARAMETER          ( HNDRD = 100.0D0 )
00166       DOUBLE PRECISION   MEIGTH
00167       PARAMETER          ( MEIGTH = -0.125D0 )
00168       INTEGER            MAXITR
00169       PARAMETER          ( MAXITR = 6 )
00170 *     ..
00171 *     .. Local Scalars ..
00172       LOGICAL            LOWER, ROTATE
00173       INTEGER            I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
00174      $                   NM12, NM13, OLDLL, OLDM
00175       DOUBLE PRECISION   ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
00176      $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
00177      $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
00178      $                   SN, THRESH, TOL, TOLMUL, UNFL
00179 *     ..
00180 *     .. External Functions ..
00181       LOGICAL            LSAME
00182       DOUBLE PRECISION   DLAMCH
00183       EXTERNAL           LSAME, DLAMCH
00184 *     ..
00185 *     .. External Subroutines ..
00186       EXTERNAL           DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT,
00187      $                   DSCAL, DSWAP, XERBLA
00188 *     ..
00189 *     .. Intrinsic Functions ..
00190       INTRINSIC          ABS, DBLE, MAX, MIN, SIGN, SQRT
00191 *     ..
00192 *     .. Executable Statements ..
00193 *
00194 *     Test the input parameters.
00195 *
00196       INFO = 0
00197       LOWER = LSAME( UPLO, 'L' )
00198       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
00199          INFO = -1
00200       ELSE IF( N.LT.0 ) THEN
00201          INFO = -2
00202       ELSE IF( NCVT.LT.0 ) THEN
00203          INFO = -3
00204       ELSE IF( NRU.LT.0 ) THEN
00205          INFO = -4
00206       ELSE IF( NCC.LT.0 ) THEN
00207          INFO = -5
00208       ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
00209      $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
00210          INFO = -9
00211       ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
00212          INFO = -11
00213       ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
00214      $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
00215          INFO = -13
00216       END IF
00217       IF( INFO.NE.0 ) THEN
00218          CALL XERBLA( 'DBDSQR', -INFO )
00219          RETURN
00220       END IF
00221       IF( N.EQ.0 )
00222      $   RETURN
00223       IF( N.EQ.1 )
00224      $   GO TO 160
00225 *
00226 *     ROTATE is true if any singular vectors desired, false otherwise
00227 *
00228       ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
00229 *
00230 *     If no singular vectors desired, use qd algorithm
00231 *
00232       IF( .NOT.ROTATE ) THEN
00233          CALL DLASQ1( N, D, E, WORK, INFO )
00234          RETURN
00235       END IF
00236 *
00237       NM1 = N - 1
00238       NM12 = NM1 + NM1
00239       NM13 = NM12 + NM1
00240       IDIR = 0
00241 *
00242 *     Get machine constants
00243 *
00244       EPS = DLAMCH( 'Epsilon' )
00245       UNFL = DLAMCH( 'Safe minimum' )
00246 *
00247 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
00248 *     by applying Givens rotations on the left
00249 *
00250       IF( LOWER ) THEN
00251          DO 10 I = 1, N - 1
00252             CALL DLARTG( D( I ), E( I ), CS, SN, R )
00253             D( I ) = R
00254             E( I ) = SN*D( I+1 )
00255             D( I+1 ) = CS*D( I+1 )
00256             WORK( I ) = CS
00257             WORK( NM1+I ) = SN
00258    10    CONTINUE
00259 *
00260 *        Update singular vectors if desired
00261 *
00262          IF( NRU.GT.0 )
00263      $      CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U,
00264      $                  LDU )
00265          IF( NCC.GT.0 )
00266      $      CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C,
00267      $                  LDC )
00268       END IF
00269 *
00270 *     Compute singular values to relative accuracy TOL
00271 *     (By setting TOL to be negative, algorithm will compute
00272 *     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
00273 *
00274       TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
00275       TOL = TOLMUL*EPS
00276 *
00277 *     Compute approximate maximum, minimum singular values
00278 *
00279       SMAX = ZERO
00280       DO 20 I = 1, N
00281          SMAX = MAX( SMAX, ABS( D( I ) ) )
00282    20 CONTINUE
00283       DO 30 I = 1, N - 1
00284          SMAX = MAX( SMAX, ABS( E( I ) ) )
00285    30 CONTINUE
00286       SMINL = ZERO
00287       IF( TOL.GE.ZERO ) THEN
00288 *
00289 *        Relative accuracy desired
00290 *
00291          SMINOA = ABS( D( 1 ) )
00292          IF( SMINOA.EQ.ZERO )
00293      $      GO TO 50
00294          MU = SMINOA
00295          DO 40 I = 2, N
00296             MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
00297             SMINOA = MIN( SMINOA, MU )
00298             IF( SMINOA.EQ.ZERO )
00299      $         GO TO 50
00300    40    CONTINUE
00301    50    CONTINUE
00302          SMINOA = SMINOA / SQRT( DBLE( N ) )
00303          THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
00304       ELSE
00305 *
00306 *        Absolute accuracy desired
00307 *
00308          THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
00309       END IF
00310 *
00311 *     Prepare for main iteration loop for the singular values
00312 *     (MAXIT is the maximum number of passes through the inner
00313 *     loop permitted before nonconvergence signalled.)
00314 *
00315       MAXIT = MAXITR*N*N
00316       ITER = 0
00317       OLDLL = -1
00318       OLDM = -1
00319 *
00320 *     M points to last element of unconverged part of matrix
00321 *
00322       M = N
00323 *
00324 *     Begin main iteration loop
00325 *
00326    60 CONTINUE
00327 *
00328 *     Check for convergence or exceeding iteration count
00329 *
00330       IF( M.LE.1 )
00331      $   GO TO 160
00332       IF( ITER.GT.MAXIT )
00333      $   GO TO 200
00334 *
00335 *     Find diagonal block of matrix to work on
00336 *
00337       IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
00338      $   D( M ) = ZERO
00339       SMAX = ABS( D( M ) )
00340       SMIN = SMAX
00341       DO 70 LLL = 1, M - 1
00342          LL = M - LLL
00343          ABSS = ABS( D( LL ) )
00344          ABSE = ABS( E( LL ) )
00345          IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
00346      $      D( LL ) = ZERO
00347          IF( ABSE.LE.THRESH )
00348      $      GO TO 80
00349          SMIN = MIN( SMIN, ABSS )
00350          SMAX = MAX( SMAX, ABSS, ABSE )
00351    70 CONTINUE
00352       LL = 0
00353       GO TO 90
00354    80 CONTINUE
00355       E( LL ) = ZERO
00356 *
00357 *     Matrix splits since E(LL) = 0
00358 *
00359       IF( LL.EQ.M-1 ) THEN
00360 *
00361 *        Convergence of bottom singular value, return to top of loop
00362 *
00363          M = M - 1
00364          GO TO 60
00365       END IF
00366    90 CONTINUE
00367       LL = LL + 1
00368 *
00369 *     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
00370 *
00371       IF( LL.EQ.M-1 ) THEN
00372 *
00373 *        2 by 2 block, handle separately
00374 *
00375          CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
00376      $                COSR, SINL, COSL )
00377          D( M-1 ) = SIGMX
00378          E( M-1 ) = ZERO
00379          D( M ) = SIGMN
00380 *
00381 *        Compute singular vectors, if desired
00382 *
00383          IF( NCVT.GT.0 )
00384      $      CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,
00385      $                 SINR )
00386          IF( NRU.GT.0 )
00387      $      CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
00388          IF( NCC.GT.0 )
00389      $      CALL DROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
00390      $                 SINL )
00391          M = M - 2
00392          GO TO 60
00393       END IF
00394 *
00395 *     If working on new submatrix, choose shift direction
00396 *     (from larger end diagonal element towards smaller)
00397 *
00398       IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
00399          IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
00400 *
00401 *           Chase bulge from top (big end) to bottom (small end)
00402 *
00403             IDIR = 1
00404          ELSE
00405 *
00406 *           Chase bulge from bottom (big end) to top (small end)
00407 *
00408             IDIR = 2
00409          END IF
00410       END IF
00411 *
00412 *     Apply convergence tests
00413 *
00414       IF( IDIR.EQ.1 ) THEN
00415 *
00416 *        Run convergence test in forward direction
00417 *        First apply standard test to bottom of matrix
00418 *
00419          IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
00420      $       ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
00421             E( M-1 ) = ZERO
00422             GO TO 60
00423          END IF
00424 *
00425          IF( TOL.GE.ZERO ) THEN
00426 *
00427 *           If relative accuracy desired,
00428 *           apply convergence criterion forward
00429 *
00430             MU = ABS( D( LL ) )
00431             SMINL = MU
00432             DO 100 LLL = LL, M - 1
00433                IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
00434                   E( LLL ) = ZERO
00435                   GO TO 60
00436                END IF
00437                MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
00438                SMINL = MIN( SMINL, MU )
00439   100       CONTINUE
00440          END IF
00441 *
00442       ELSE
00443 *
00444 *        Run convergence test in backward direction
00445 *        First apply standard test to top of matrix
00446 *
00447          IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
00448      $       ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
00449             E( LL ) = ZERO
00450             GO TO 60
00451          END IF
00452 *
00453          IF( TOL.GE.ZERO ) THEN
00454 *
00455 *           If relative accuracy desired,
00456 *           apply convergence criterion backward
00457 *
00458             MU = ABS( D( M ) )
00459             SMINL = MU
00460             DO 110 LLL = M - 1, LL, -1
00461                IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
00462                   E( LLL ) = ZERO
00463                   GO TO 60
00464                END IF
00465                MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
00466                SMINL = MIN( SMINL, MU )
00467   110       CONTINUE
00468          END IF
00469       END IF
00470       OLDLL = LL
00471       OLDM = M
00472 *
00473 *     Compute shift.  First, test if shifting would ruin relative
00474 *     accuracy, and if so set the shift to zero.
00475 *
00476       IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
00477      $    MAX( EPS, HNDRTH*TOL ) ) THEN
00478 *
00479 *        Use a zero shift to avoid loss of relative accuracy
00480 *
00481          SHIFT = ZERO
00482       ELSE
00483 *
00484 *        Compute the shift from 2-by-2 block at end of matrix
00485 *
00486          IF( IDIR.EQ.1 ) THEN
00487             SLL = ABS( D( LL ) )
00488             CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
00489          ELSE
00490             SLL = ABS( D( M ) )
00491             CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
00492          END IF
00493 *
00494 *        Test if shift negligible, and if so set to zero
00495 *
00496          IF( SLL.GT.ZERO ) THEN
00497             IF( ( SHIFT / SLL )**2.LT.EPS )
00498      $         SHIFT = ZERO
00499          END IF
00500       END IF
00501 *
00502 *     Increment iteration count
00503 *
00504       ITER = ITER + M - LL
00505 *
00506 *     If SHIFT = 0, do simplified QR iteration
00507 *
00508       IF( SHIFT.EQ.ZERO ) THEN
00509          IF( IDIR.EQ.1 ) THEN
00510 *
00511 *           Chase bulge from top to bottom
00512 *           Save cosines and sines for later singular vector updates
00513 *
00514             CS = ONE
00515             OLDCS = ONE
00516             DO 120 I = LL, M - 1
00517                CALL DLARTG( D( I )*CS, E( I ), CS, SN, R )
00518                IF( I.GT.LL )
00519      $            E( I-1 ) = OLDSN*R
00520                CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
00521                WORK( I-LL+1 ) = CS
00522                WORK( I-LL+1+NM1 ) = SN
00523                WORK( I-LL+1+NM12 ) = OLDCS
00524                WORK( I-LL+1+NM13 ) = OLDSN
00525   120       CONTINUE
00526             H = D( M )*CS
00527             D( M ) = H*OLDCS
00528             E( M-1 ) = H*OLDSN
00529 *
00530 *           Update singular vectors
00531 *
00532             IF( NCVT.GT.0 )
00533      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
00534      $                     WORK( N ), VT( LL, 1 ), LDVT )
00535             IF( NRU.GT.0 )
00536      $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
00537      $                     WORK( NM13+1 ), U( 1, LL ), LDU )
00538             IF( NCC.GT.0 )
00539      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
00540      $                     WORK( NM13+1 ), C( LL, 1 ), LDC )
00541 *
00542 *           Test convergence
00543 *
00544             IF( ABS( E( M-1 ) ).LE.THRESH )
00545      $         E( M-1 ) = ZERO
00546 *
00547          ELSE
00548 *
00549 *           Chase bulge from bottom to top
00550 *           Save cosines and sines for later singular vector updates
00551 *
00552             CS = ONE
00553             OLDCS = ONE
00554             DO 130 I = M, LL + 1, -1
00555                CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
00556                IF( I.LT.M )
00557      $            E( I ) = OLDSN*R
00558                CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
00559                WORK( I-LL ) = CS
00560                WORK( I-LL+NM1 ) = -SN
00561                WORK( I-LL+NM12 ) = OLDCS
00562                WORK( I-LL+NM13 ) = -OLDSN
00563   130       CONTINUE
00564             H = D( LL )*CS
00565             D( LL ) = H*OLDCS
00566             E( LL ) = H*OLDSN
00567 *
00568 *           Update singular vectors
00569 *
00570             IF( NCVT.GT.0 )
00571      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
00572      $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )
00573             IF( NRU.GT.0 )
00574      $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
00575      $                     WORK( N ), U( 1, LL ), LDU )
00576             IF( NCC.GT.0 )
00577      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
00578      $                     WORK( N ), C( LL, 1 ), LDC )
00579 *
00580 *           Test convergence
00581 *
00582             IF( ABS( E( LL ) ).LE.THRESH )
00583      $         E( LL ) = ZERO
00584          END IF
00585       ELSE
00586 *
00587 *        Use nonzero shift
00588 *
00589          IF( IDIR.EQ.1 ) THEN
00590 *
00591 *           Chase bulge from top to bottom
00592 *           Save cosines and sines for later singular vector updates
00593 *
00594             F = ( ABS( D( LL ) )-SHIFT )*
00595      $          ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
00596             G = E( LL )
00597             DO 140 I = LL, M - 1
00598                CALL DLARTG( F, G, COSR, SINR, R )
00599                IF( I.GT.LL )
00600      $            E( I-1 ) = R
00601                F = COSR*D( I ) + SINR*E( I )
00602                E( I ) = COSR*E( I ) - SINR*D( I )
00603                G = SINR*D( I+1 )
00604                D( I+1 ) = COSR*D( I+1 )
00605                CALL DLARTG( F, G, COSL, SINL, R )
00606                D( I ) = R
00607                F = COSL*E( I ) + SINL*D( I+1 )
00608                D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
00609                IF( I.LT.M-1 ) THEN
00610                   G = SINL*E( I+1 )
00611                   E( I+1 ) = COSL*E( I+1 )
00612                END IF
00613                WORK( I-LL+1 ) = COSR
00614                WORK( I-LL+1+NM1 ) = SINR
00615                WORK( I-LL+1+NM12 ) = COSL
00616                WORK( I-LL+1+NM13 ) = SINL
00617   140       CONTINUE
00618             E( M-1 ) = F
00619 *
00620 *           Update singular vectors
00621 *
00622             IF( NCVT.GT.0 )
00623      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
00624      $                     WORK( N ), VT( LL, 1 ), LDVT )
00625             IF( NRU.GT.0 )
00626      $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
00627      $                     WORK( NM13+1 ), U( 1, LL ), LDU )
00628             IF( NCC.GT.0 )
00629      $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
00630      $                     WORK( NM13+1 ), C( LL, 1 ), LDC )
00631 *
00632 *           Test convergence
00633 *
00634             IF( ABS( E( M-1 ) ).LE.THRESH )
00635      $         E( M-1 ) = ZERO
00636 *
00637          ELSE
00638 *
00639 *           Chase bulge from bottom to top
00640 *           Save cosines and sines for later singular vector updates
00641 *
00642             F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
00643      $          D( M ) )
00644             G = E( M-1 )
00645             DO 150 I = M, LL + 1, -1
00646                CALL DLARTG( F, G, COSR, SINR, R )
00647                IF( I.LT.M )
00648      $            E( I ) = R
00649                F = COSR*D( I ) + SINR*E( I-1 )
00650                E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
00651                G = SINR*D( I-1 )
00652                D( I-1 ) = COSR*D( I-1 )
00653                CALL DLARTG( F, G, COSL, SINL, R )
00654                D( I ) = R
00655                F = COSL*E( I-1 ) + SINL*D( I-1 )
00656                D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
00657                IF( I.GT.LL+1 ) THEN
00658                   G = SINL*E( I-2 )
00659                   E( I-2 ) = COSL*E( I-2 )
00660                END IF
00661                WORK( I-LL ) = COSR
00662                WORK( I-LL+NM1 ) = -SINR
00663                WORK( I-LL+NM12 ) = COSL
00664                WORK( I-LL+NM13 ) = -SINL
00665   150       CONTINUE
00666             E( LL ) = F
00667 *
00668 *           Test convergence
00669 *
00670             IF( ABS( E( LL ) ).LE.THRESH )
00671      $         E( LL ) = ZERO
00672 *
00673 *           Update singular vectors if desired
00674 *
00675             IF( NCVT.GT.0 )
00676      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
00677      $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )
00678             IF( NRU.GT.0 )
00679      $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
00680      $                     WORK( N ), U( 1, LL ), LDU )
00681             IF( NCC.GT.0 )
00682      $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
00683      $                     WORK( N ), C( LL, 1 ), LDC )
00684          END IF
00685       END IF
00686 *
00687 *     QR iteration finished, go back and check convergence
00688 *
00689       GO TO 60
00690 *
00691 *     All singular values converged, so make them positive
00692 *
00693   160 CONTINUE
00694       DO 170 I = 1, N
00695          IF( D( I ).LT.ZERO ) THEN
00696             D( I ) = -D( I )
00697 *
00698 *           Change sign of singular vectors, if desired
00699 *
00700             IF( NCVT.GT.0 )
00701      $         CALL DSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
00702          END IF
00703   170 CONTINUE
00704 *
00705 *     Sort the singular values into decreasing order (insertion sort on
00706 *     singular values, but only one transposition per singular vector)
00707 *
00708       DO 190 I = 1, N - 1
00709 *
00710 *        Scan for smallest D(I)
00711 *
00712          ISUB = 1
00713          SMIN = D( 1 )
00714          DO 180 J = 2, N + 1 - I
00715             IF( D( J ).LE.SMIN ) THEN
00716                ISUB = J
00717                SMIN = D( J )
00718             END IF
00719   180    CONTINUE
00720          IF( ISUB.NE.N+1-I ) THEN
00721 *
00722 *           Swap singular values and vectors
00723 *
00724             D( ISUB ) = D( N+1-I )
00725             D( N+1-I ) = SMIN
00726             IF( NCVT.GT.0 )
00727      $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
00728      $                     LDVT )
00729             IF( NRU.GT.0 )
00730      $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
00731             IF( NCC.GT.0 )
00732      $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
00733          END IF
00734   190 CONTINUE
00735       GO TO 220
00736 *
00737 *     Maximum number of iterations exceeded, failure to converge
00738 *
00739   200 CONTINUE
00740       INFO = 0
00741       DO 210 I = 1, N - 1
00742          IF( E( I ).NE.ZERO )
00743      $      INFO = INFO + 1
00744   210 CONTINUE
00745   220 CONTINUE
00746       RETURN
00747 *
00748 *     End of DBDSQR
00749 *
00750       END