Move LAPACK trunk into position.

1 files changed, 742 insertions, 0 deletions

diff --git a/SRC/cbdsqr.f b/SRC/cbdsqr.f
new file mode 100644
index 00000000..cc03e132
--- /dev/null
+++ b/SRC/cbdsqr.f
@@ -0,0 +1,742 @@
+      SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
+     $                   LDU, C, LDC, RWORK, INFO )
+*
+*  -- LAPACK routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
+*     ..
+*     .. Array Arguments ..
+      REAL               D( * ), E( * ), RWORK( * )
+      COMPLEX            C( LDC, * ), U( LDU, * ), VT( LDVT, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CBDSQR computes the singular values and, optionally, the right and/or
+*  left singular vectors from the singular value decomposition (SVD) of
+*  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
+*  zero-shift QR algorithm.  The SVD of B has the form
+*  
+*     B = Q * S * P**H
+*  
+*  where S is the diagonal matrix of singular values, Q is an orthogonal
+*  matrix of left singular vectors, and P is an orthogonal matrix of
+*  right singular vectors.  If left singular vectors are requested, this
+*  subroutine actually returns U*Q instead of Q, and, if right singular
+*  vectors are requested, this subroutine returns P**H*VT instead of
+*  P**H, for given complex input matrices U and VT.  When U and VT are
+*  the unitary matrices that reduce a general matrix A to bidiagonal
+*  form: A = U*B*VT, as computed by CGEBRD, then
+*  
+*     A = (U*Q) * S * (P**H*VT)
+*  
+*  is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
+*  for a given complex input matrix C.
+*
+*  See "Computing  Small Singular Values of Bidiagonal Matrices With
+*  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
+*  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
+*  no. 5, pp. 873-912, Sept 1990) and
+*  "Accurate singular values and differential qd algorithms," by
+*  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
+*  Department, University of California at Berkeley, July 1992
+*  for a detailed description of the algorithm.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          = 'U':  B is upper bidiagonal;
+*          = 'L':  B is lower bidiagonal.
+*
+*  N       (input) INTEGER
+*          The order of the matrix B.  N >= 0.
+*
+*  NCVT    (input) INTEGER
+*          The number of columns of the matrix VT. NCVT >= 0.
+*
+*  NRU     (input) INTEGER
+*          The number of rows of the matrix U. NRU >= 0.
+*
+*  NCC     (input) INTEGER
+*          The number of columns of the matrix C. NCC >= 0.
+*
+*  D       (input/output) REAL array, dimension (N)
+*          On entry, the n diagonal elements of the bidiagonal matrix B.
+*          On exit, if INFO=0, the singular values of B in decreasing
+*          order.
+*
+*  E       (input/output) REAL array, dimension (N-1)
+*          On entry, the N-1 offdiagonal elements of the bidiagonal
+*          matrix B.
+*          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
+*          will contain the diagonal and superdiagonal elements of a
+*          bidiagonal matrix orthogonally equivalent to the one given
+*          as input.
+*
+*  VT      (input/output) COMPLEX array, dimension (LDVT, NCVT)
+*          On entry, an N-by-NCVT matrix VT.
+*          On exit, VT is overwritten by P**H * VT.
+*          Not referenced if NCVT = 0.
+*
+*  LDVT    (input) INTEGER
+*          The leading dimension of the array VT.
+*          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
+*
+*  U       (input/output) COMPLEX array, dimension (LDU, N)
+*          On entry, an NRU-by-N matrix U.
+*          On exit, U is overwritten by U * Q.
+*          Not referenced if NRU = 0.
+*
+*  LDU     (input) INTEGER
+*          The leading dimension of the array U.  LDU >= max(1,NRU).
+*
+*  C       (input/output) COMPLEX array, dimension (LDC, NCC)
+*          On entry, an N-by-NCC matrix C.
+*          On exit, C is overwritten by Q**H * C.
+*          Not referenced if NCC = 0.
+*
+*  LDC     (input) INTEGER
+*          The leading dimension of the array C.
+*          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
+*
+*  RWORK   (workspace) REAL array, dimension (2*N) 
+*          if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  If INFO = -i, the i-th argument had an illegal value
+*          > 0:  the algorithm did not converge; D and E contain the
+*                elements of a bidiagonal matrix which is orthogonally
+*                similar to the input matrix B;  if INFO = i, i
+*                elements of E have not converged to zero.
+*
+*  Internal Parameters
+*  ===================
+*
+*  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
+*          TOLMUL controls the convergence criterion of the QR loop.
+*          If it is positive, TOLMUL*EPS is the desired relative
+*             precision in the computed singular values.
+*          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
+*             desired absolute accuracy in the computed singular
+*             values (corresponds to relative accuracy
+*             abs(TOLMUL*EPS) in the largest singular value.
+*          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
+*             between 10 (for fast convergence) and .1/EPS
+*             (for there to be some accuracy in the results).
+*          Default is to lose at either one eighth or 2 of the
+*             available decimal digits in each computed singular value
+*             (whichever is smaller).
+*
+*  MAXITR  INTEGER, default = 6
+*          MAXITR controls the maximum number of passes of the
+*          algorithm through its inner loop. The algorithms stops
+*          (and so fails to converge) if the number of passes
+*          through the inner loop exceeds MAXITR*N**2.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO
+      PARAMETER          ( ZERO = 0.0E0 )
+      REAL               ONE
+      PARAMETER          ( ONE = 1.0E0 )
+      REAL               NEGONE
+      PARAMETER          ( NEGONE = -1.0E0 )
+      REAL               HNDRTH
+      PARAMETER          ( HNDRTH = 0.01E0 )
+      REAL               TEN
+      PARAMETER          ( TEN = 10.0E0 )
+      REAL               HNDRD
+      PARAMETER          ( HNDRD = 100.0E0 )
+      REAL               MEIGTH
+      PARAMETER          ( MEIGTH = -0.125E0 )
+      INTEGER            MAXITR
+      PARAMETER          ( MAXITR = 6 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            LOWER, ROTATE
+      INTEGER            I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
+     $                   NM12, NM13, OLDLL, OLDM
+      REAL               ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
+     $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
+     $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
+     $                   SN, THRESH, TOL, TOLMUL, UNFL
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      REAL               SLAMCH
+      EXTERNAL           LSAME, SLAMCH
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLASR, CSROT, CSSCAL, CSWAP, SLARTG, SLAS2,
+     $                   SLASQ1, SLASV2, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN, REAL, SIGN, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      LOWER = LSAME( UPLO, 'L' )
+      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( NCVT.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( NRU.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( NCC.LT.0 ) THEN
+         INFO = -5
+      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
+     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
+         INFO = -9
+      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
+         INFO = -11
+      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
+     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
+         INFO = -13
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CBDSQR', -INFO )
+         RETURN
+      END IF
+      IF( N.EQ.0 )
+     $   RETURN
+      IF( N.EQ.1 )
+     $   GO TO 160
+*
+*     ROTATE is true if any singular vectors desired, false otherwise
+*
+      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
+*
+*     If no singular vectors desired, use qd algorithm
+*
+      IF( .NOT.ROTATE ) THEN
+         CALL SLASQ1( N, D, E, RWORK, INFO )
+         RETURN
+      END IF
+*
+      NM1 = N - 1
+      NM12 = NM1 + NM1
+      NM13 = NM12 + NM1
+      IDIR = 0
+*
+*     Get machine constants
+*
+      EPS = SLAMCH( 'Epsilon' )
+      UNFL = SLAMCH( 'Safe minimum' )
+*
+*     If matrix lower bidiagonal, rotate to be upper bidiagonal
+*     by applying Givens rotations on the left
+*
+      IF( LOWER ) THEN
+         DO 10 I = 1, N - 1
+            CALL SLARTG( D( I ), E( I ), CS, SN, R )
+            D( I ) = R
+            E( I ) = SN*D( I+1 )
+            D( I+1 ) = CS*D( I+1 )
+            RWORK( I ) = CS
+            RWORK( NM1+I ) = SN
+   10    CONTINUE
+*
+*        Update singular vectors if desired
+*
+         IF( NRU.GT.0 )
+     $      CALL CLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ),
+     $                  U, LDU )
+         IF( NCC.GT.0 )
+     $      CALL CLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), RWORK( N ),
+     $                  C, LDC )
+      END IF
+*
+*     Compute singular values to relative accuracy TOL
+*     (By setting TOL to be negative, algorithm will compute
+*     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
+*
+      TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
+      TOL = TOLMUL*EPS
+*
+*     Compute approximate maximum, minimum singular values
+*
+      SMAX = ZERO
+      DO 20 I = 1, N
+         SMAX = MAX( SMAX, ABS( D( I ) ) )
+   20 CONTINUE
+      DO 30 I = 1, N - 1
+         SMAX = MAX( SMAX, ABS( E( I ) ) )
+   30 CONTINUE
+      SMINL = ZERO
+      IF( TOL.GE.ZERO ) THEN
+*
+*        Relative accuracy desired
+*
+         SMINOA = ABS( D( 1 ) )
+         IF( SMINOA.EQ.ZERO )
+     $      GO TO 50
+         MU = SMINOA
+         DO 40 I = 2, N
+            MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
+            SMINOA = MIN( SMINOA, MU )
+            IF( SMINOA.EQ.ZERO )
+     $         GO TO 50
+   40    CONTINUE
+   50    CONTINUE
+         SMINOA = SMINOA / SQRT( REAL( N ) )
+         THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
+      ELSE
+*
+*        Absolute accuracy desired
+*
+         THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
+      END IF
+*
+*     Prepare for main iteration loop for the singular values
+*     (MAXIT is the maximum number of passes through the inner
+*     loop permitted before nonconvergence signalled.)
+*
+      MAXIT = MAXITR*N*N
+      ITER = 0
+      OLDLL = -1
+      OLDM = -1
+*
+*     M points to last element of unconverged part of matrix
+*
+      M = N
+*
+*     Begin main iteration loop
+*
+   60 CONTINUE
+*
+*     Check for convergence or exceeding iteration count
+*
+      IF( M.LE.1 )
+     $   GO TO 160
+      IF( ITER.GT.MAXIT )
+     $   GO TO 200
+*
+*     Find diagonal block of matrix to work on
+*
+      IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
+     $   D( M ) = ZERO
+      SMAX = ABS( D( M ) )
+      SMIN = SMAX
+      DO 70 LLL = 1, M - 1
+         LL = M - LLL
+         ABSS = ABS( D( LL ) )
+         ABSE = ABS( E( LL ) )
+         IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
+     $      D( LL ) = ZERO
+         IF( ABSE.LE.THRESH )
+     $      GO TO 80
+         SMIN = MIN( SMIN, ABSS )
+         SMAX = MAX( SMAX, ABSS, ABSE )
+   70 CONTINUE
+      LL = 0
+      GO TO 90
+   80 CONTINUE
+      E( LL ) = ZERO
+*
+*     Matrix splits since E(LL) = 0
+*
+      IF( LL.EQ.M-1 ) THEN
+*
+*        Convergence of bottom singular value, return to top of loop
+*
+         M = M - 1
+         GO TO 60
+      END IF
+   90 CONTINUE
+      LL = LL + 1
+*
+*     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
+*
+      IF( LL.EQ.M-1 ) THEN
+*
+*        2 by 2 block, handle separately
+*
+         CALL SLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
+     $                COSR, SINL, COSL )
+         D( M-1 ) = SIGMX
+         E( M-1 ) = ZERO
+         D( M ) = SIGMN
+*
+*        Compute singular vectors, if desired
+*
+         IF( NCVT.GT.0 )
+     $      CALL CSROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT,
+     $                  COSR, SINR )
+         IF( NRU.GT.0 )
+     $      CALL CSROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
+         IF( NCC.GT.0 )
+     $      CALL CSROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
+     $                  SINL )
+         M = M - 2
+         GO TO 60
+      END IF
+*
+*     If working on new submatrix, choose shift direction
+*     (from larger end diagonal element towards smaller)
+*
+      IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
+         IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
+*
+*           Chase bulge from top (big end) to bottom (small end)
+*
+            IDIR = 1
+         ELSE
+*
+*           Chase bulge from bottom (big end) to top (small end)
+*
+            IDIR = 2
+         END IF
+      END IF
+*
+*     Apply convergence tests
+*
+      IF( IDIR.EQ.1 ) THEN
+*
+*        Run convergence test in forward direction
+*        First apply standard test to bottom of matrix
+*
+         IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
+     $       ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
+            E( M-1 ) = ZERO
+            GO TO 60
+         END IF
+*
+         IF( TOL.GE.ZERO ) THEN
+*
+*           If relative accuracy desired,
+*           apply convergence criterion forward
+*
+            MU = ABS( D( LL ) )
+            SMINL = MU
+            DO 100 LLL = LL, M - 1
+               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
+                  E( LLL ) = ZERO
+                  GO TO 60
+               END IF
+               MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
+               SMINL = MIN( SMINL, MU )
+  100       CONTINUE
+         END IF
+*
+      ELSE
+*
+*        Run convergence test in backward direction
+*        First apply standard test to top of matrix
+*
+         IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
+     $       ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
+            E( LL ) = ZERO
+            GO TO 60
+         END IF
+*
+         IF( TOL.GE.ZERO ) THEN
+*
+*           If relative accuracy desired,
+*           apply convergence criterion backward
+*
+            MU = ABS( D( M ) )
+            SMINL = MU
+            DO 110 LLL = M - 1, LL, -1
+               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
+                  E( LLL ) = ZERO
+                  GO TO 60
+               END IF
+               MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
+               SMINL = MIN( SMINL, MU )
+  110       CONTINUE
+         END IF
+      END IF
+      OLDLL = LL
+      OLDM = M
+*
+*     Compute shift.  First, test if shifting would ruin relative
+*     accuracy, and if so set the shift to zero.
+*
+      IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
+     $    MAX( EPS, HNDRTH*TOL ) ) THEN
+*
+*        Use a zero shift to avoid loss of relative accuracy
+*
+         SHIFT = ZERO
+      ELSE
+*
+*        Compute the shift from 2-by-2 block at end of matrix
+*
+         IF( IDIR.EQ.1 ) THEN
+            SLL = ABS( D( LL ) )
+            CALL SLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
+         ELSE
+            SLL = ABS( D( M ) )
+            CALL SLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
+         END IF
+*
+*        Test if shift negligible, and if so set to zero
+*
+         IF( SLL.GT.ZERO ) THEN
+            IF( ( SHIFT / SLL )**2.LT.EPS )
+     $         SHIFT = ZERO
+         END IF
+      END IF
+*
+*     Increment iteration count
+*
+      ITER = ITER + M - LL
+*
+*     If SHIFT = 0, do simplified QR iteration
+*
+      IF( SHIFT.EQ.ZERO ) THEN
+         IF( IDIR.EQ.1 ) THEN
+*
+*           Chase bulge from top to bottom
+*           Save cosines and sines for later singular vector updates
+*
+            CS = ONE
+            OLDCS = ONE
+            DO 120 I = LL, M - 1
+               CALL SLARTG( D( I )*CS, E( I ), CS, SN, R )
+               IF( I.GT.LL )
+     $            E( I-1 ) = OLDSN*R
+               CALL SLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
+               RWORK( I-LL+1 ) = CS
+               RWORK( I-LL+1+NM1 ) = SN
+               RWORK( I-LL+1+NM12 ) = OLDCS
+               RWORK( I-LL+1+NM13 ) = OLDSN
+  120       CONTINUE
+            H = D( M )*CS
+            D( M ) = H*OLDCS
+            E( M-1 ) = H*OLDSN
+*
+*           Update singular vectors
+*
+            IF( NCVT.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
+     $                     RWORK( N ), VT( LL, 1 ), LDVT )
+            IF( NRU.GT.0 )
+     $         CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ),
+     $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
+            IF( NCC.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ),
+     $                     RWORK( NM13+1 ), C( LL, 1 ), LDC )
+*
+*           Test convergence
+*
+            IF( ABS( E( M-1 ) ).LE.THRESH )
+     $         E( M-1 ) = ZERO
+*
+         ELSE
+*
+*           Chase bulge from bottom to top
+*           Save cosines and sines for later singular vector updates
+*
+            CS = ONE
+            OLDCS = ONE
+            DO 130 I = M, LL + 1, -1
+               CALL SLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
+               IF( I.LT.M )
+     $            E( I ) = OLDSN*R
+               CALL SLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
+               RWORK( I-LL ) = CS
+               RWORK( I-LL+NM1 ) = -SN
+               RWORK( I-LL+NM12 ) = OLDCS
+               RWORK( I-LL+NM13 ) = -OLDSN
+  130       CONTINUE
+            H = D( LL )*CS
+            D( LL ) = H*OLDCS
+            E( LL ) = H*OLDSN
+*
+*           Update singular vectors
+*
+            IF( NCVT.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ),
+     $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
+            IF( NRU.GT.0 )
+     $         CALL CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
+     $                     RWORK( N ), U( 1, LL ), LDU )
+            IF( NCC.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ),
+     $                     RWORK( N ), C( LL, 1 ), LDC )
+*
+*           Test convergence
+*
+            IF( ABS( E( LL ) ).LE.THRESH )
+     $         E( LL ) = ZERO
+         END IF
+      ELSE
+*
+*        Use nonzero shift
+*
+         IF( IDIR.EQ.1 ) THEN
+*
+*           Chase bulge from top to bottom
+*           Save cosines and sines for later singular vector updates
+*
+            F = ( ABS( D( LL ) )-SHIFT )*
+     $          ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
+            G = E( LL )
+            DO 140 I = LL, M - 1
+               CALL SLARTG( F, G, COSR, SINR, R )
+               IF( I.GT.LL )
+     $            E( I-1 ) = R
+               F = COSR*D( I ) + SINR*E( I )
+               E( I ) = COSR*E( I ) - SINR*D( I )
+               G = SINR*D( I+1 )
+               D( I+1 ) = COSR*D( I+1 )
+               CALL SLARTG( F, G, COSL, SINL, R )
+               D( I ) = R
+               F = COSL*E( I ) + SINL*D( I+1 )
+               D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
+               IF( I.LT.M-1 ) THEN
+                  G = SINL*E( I+1 )
+                  E( I+1 ) = COSL*E( I+1 )
+               END IF
+               RWORK( I-LL+1 ) = COSR
+               RWORK( I-LL+1+NM1 ) = SINR
+               RWORK( I-LL+1+NM12 ) = COSL
+               RWORK( I-LL+1+NM13 ) = SINL
+  140       CONTINUE
+            E( M-1 ) = F
+*
+*           Update singular vectors
+*
+            IF( NCVT.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
+     $                     RWORK( N ), VT( LL, 1 ), LDVT )
+            IF( NRU.GT.0 )
+     $         CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ),
+     $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
+            IF( NCC.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ),
+     $                     RWORK( NM13+1 ), C( LL, 1 ), LDC )
+*
+*           Test convergence
+*
+            IF( ABS( E( M-1 ) ).LE.THRESH )
+     $         E( M-1 ) = ZERO
+*
+         ELSE
+*
+*           Chase bulge from bottom to top
+*           Save cosines and sines for later singular vector updates
+*
+            F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
+     $          D( M ) )
+            G = E( M-1 )
+            DO 150 I = M, LL + 1, -1
+               CALL SLARTG( F, G, COSR, SINR, R )
+               IF( I.LT.M )
+     $            E( I ) = R
+               F = COSR*D( I ) + SINR*E( I-1 )
+               E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
+               G = SINR*D( I-1 )
+               D( I-1 ) = COSR*D( I-1 )
+               CALL SLARTG( F, G, COSL, SINL, R )
+               D( I ) = R
+               F = COSL*E( I-1 ) + SINL*D( I-1 )
+               D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
+               IF( I.GT.LL+1 ) THEN
+                  G = SINL*E( I-2 )
+                  E( I-2 ) = COSL*E( I-2 )
+               END IF
+               RWORK( I-LL ) = COSR
+               RWORK( I-LL+NM1 ) = -SINR
+               RWORK( I-LL+NM12 ) = COSL
+               RWORK( I-LL+NM13 ) = -SINL
+  150       CONTINUE
+            E( LL ) = F
+*
+*           Test convergence
+*
+            IF( ABS( E( LL ) ).LE.THRESH )
+     $         E( LL ) = ZERO
+*
+*           Update singular vectors if desired
+*
+            IF( NCVT.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ),
+     $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
+            IF( NRU.GT.0 )
+     $         CALL CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
+     $                     RWORK( N ), U( 1, LL ), LDU )
+            IF( NCC.GT.0 )
+     $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ),
+     $                     RWORK( N ), C( LL, 1 ), LDC )
+         END IF
+      END IF
+*
+*     QR iteration finished, go back and check convergence
+*
+      GO TO 60
+*
+*     All singular values converged, so make them positive
+*
+  160 CONTINUE
+      DO 170 I = 1, N
+         IF( D( I ).LT.ZERO ) THEN
+            D( I ) = -D( I )
+*
+*           Change sign of singular vectors, if desired
+*
+            IF( NCVT.GT.0 )
+     $         CALL CSSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
+         END IF
+  170 CONTINUE
+*
+*     Sort the singular values into decreasing order (insertion sort on
+*     singular values, but only one transposition per singular vector)
+*
+      DO 190 I = 1, N - 1
+*
+*        Scan for smallest D(I)
+*
+         ISUB = 1
+         SMIN = D( 1 )
+         DO 180 J = 2, N + 1 - I
+            IF( D( J ).LE.SMIN ) THEN
+               ISUB = J
+               SMIN = D( J )
+            END IF
+  180    CONTINUE
+         IF( ISUB.NE.N+1-I ) THEN
+*
+*           Swap singular values and vectors
+*
+            D( ISUB ) = D( N+1-I )
+            D( N+1-I ) = SMIN
+            IF( NCVT.GT.0 )
+     $         CALL CSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
+     $                     LDVT )
+            IF( NRU.GT.0 )
+     $         CALL CSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
+            IF( NCC.GT.0 )
+     $         CALL CSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
+         END IF
+  190 CONTINUE
+      GO TO 220
+*
+*     Maximum number of iterations exceeded, failure to converge
+*
+  200 CONTINUE
+      INFO = 0
+      DO 210 I = 1, N - 1
+         IF( E( I ).NE.ZERO )
+     $      INFO = INFO + 1
+  210 CONTINUE
+  220 CONTINUE
+      RETURN
+*
+*     End of CBDSQR
+*
+      END