LAPACK: SRC/cbdsqr.f Source File

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