Moving lapack to right place

1 files changed, 0 insertions, 742 deletions

diff --git a/src/lib/lapack/zbdsqr.f b/src/lib/lapack/zbdsqr.f
deleted file mode 100644
index f9086be5..00000000
--- a/src/lib/lapack/zbdsqr.f
+++ /dev/null
@@ -1,742 +0,0 @@
-      SUBROUTINE ZBDSQR( 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 ..
-      DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
-      COMPLEX*16         C( LDC, * ), U( LDU, * ), VT( LDVT, * )
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  ZBDSQR 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 ZGEBRD, 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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*16 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*16 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*16 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) DOUBLE PRECISION 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  DOUBLE PRECISION, 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 ..
-      DOUBLE PRECISION   ZERO
-      PARAMETER          ( ZERO = 0.0D0 )
-      DOUBLE PRECISION   ONE
-      PARAMETER          ( ONE = 1.0D0 )
-      DOUBLE PRECISION   NEGONE
-      PARAMETER          ( NEGONE = -1.0D0 )
-      DOUBLE PRECISION   HNDRTH
-      PARAMETER          ( HNDRTH = 0.01D0 )
-      DOUBLE PRECISION   TEN
-      PARAMETER          ( TEN = 10.0D0 )
-      DOUBLE PRECISION   HNDRD
-      PARAMETER          ( HNDRD = 100.0D0 )
-      DOUBLE PRECISION   MEIGTH
-      PARAMETER          ( MEIGTH = -0.125D0 )
-      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
-      DOUBLE PRECISION   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
-      DOUBLE PRECISION   DLAMCH
-      EXTERNAL           LSAME, DLAMCH
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL           DLARTG, DLAS2, DLASQ1, DLASV2, XERBLA, ZDROT,
-     $                   ZDSCAL, ZLASR, ZSWAP
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          ABS, DBLE, MAX, MIN, 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( 'ZBDSQR', -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 DLASQ1( N, D, E, RWORK, INFO )
-         RETURN
-      END IF
-*
-      NM1 = N - 1
-      NM12 = NM1 + NM1
-      NM13 = NM12 + NM1
-      IDIR = 0
-*
-*     Get machine constants
-*
-      EPS = DLAMCH( 'Epsilon' )
-      UNFL = DLAMCH( '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 DLARTG( 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 ZLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ),
-     $                  U, LDU )
-         IF( NCC.GT.0 )
-     $      CALL ZLASR( '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( DBLE( 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 DLASV2( 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 ZDROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT,
-     $                  COSR, SINR )
-         IF( NRU.GT.0 )
-     $      CALL ZDROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
-         IF( NCC.GT.0 )
-     $      CALL ZDROT( 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 DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
-         ELSE
-            SLL = ABS( D( M ) )
-            CALL DLAS2( 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 DLARTG( D( I )*CS, E( I ), CS, SN, R )
-               IF( I.GT.LL )
-     $            E( I-1 ) = OLDSN*R
-               CALL DLARTG( 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 ZLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
-     $                     RWORK( N ), VT( LL, 1 ), LDVT )
-            IF( NRU.GT.0 )
-     $         CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ),
-     $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
-            IF( NCC.GT.0 )
-     $         CALL ZLASR( '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 DLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
-               IF( I.LT.M )
-     $            E( I ) = OLDSN*R
-               CALL DLARTG( 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 ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ),
-     $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
-            IF( NRU.GT.0 )
-     $         CALL ZLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
-     $                     RWORK( N ), U( 1, LL ), LDU )
-            IF( NCC.GT.0 )
-     $         CALL ZLASR( '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 DLARTG( 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 DLARTG( 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 ZLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
-     $                     RWORK( N ), VT( LL, 1 ), LDVT )
-            IF( NRU.GT.0 )
-     $         CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ),
-     $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
-            IF( NCC.GT.0 )
-     $         CALL ZLASR( '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 DLARTG( 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 DLARTG( 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 ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ),
-     $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
-            IF( NRU.GT.0 )
-     $         CALL ZLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
-     $                     RWORK( N ), U( 1, LL ), LDU )
-            IF( NCC.GT.0 )
-     $         CALL ZLASR( '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 ZDSCAL( 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 ZSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
-     $                     LDVT )
-            IF( NRU.GT.0 )
-     $         CALL ZSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
-            IF( NCC.GT.0 )
-     $         CALL ZSWAP( 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 ZBDSQR
-*
-      END