Chimie Théorique » OpenPath

      SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,

     $                   LDU, C, LDC, WORK, INFO )

*

*  -- LAPACK routine (version 3.2) --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*     January 2007

*

*     .. Scalar Arguments ..

      CHARACTER          UPLO

      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU

*     ..

*     .. Array Arguments ..

      DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),

     $                   VT( LDVT, * ), WORK( * )

*     ..

*

*  Purpose

*  =======

*

*  DBDSQR 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**T

*

*  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**T*VT instead of

*  P**T, for given real input matrices U and VT.  When U and VT are the

*  orthogonal matrices that reduce a general matrix A to bidiagonal

*  form:  A = U*B*VT, as computed by DGEBRD, then

*

*     A = (U*Q) * S * (P**T*VT)

*

*  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C

*  for a given real 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) DOUBLE PRECISION array, dimension (LDVT, NCVT)

*          On entry, an N-by-NCVT matrix VT.

*          On exit, VT is overwritten by P**T * 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDC, NCC)

*          On entry, an N-by-NCC matrix C.

*          On exit, C is overwritten by Q**T * 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.

*

*  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)

*

*  INFO    (output) INTEGER

*          = 0:  successful exit

*          < 0:  If INFO = -i, the i-th argument had an illegal value

*          > 0:

*             if NCVT = NRU = NCC = 0,

*                = 1, a split was marked by a positive value in E

*                = 2, current block of Z not diagonalized after 30*N

*                     iterations (in inner while loop)

*                = 3, termination criterion of outer while loop not met

*                     (program created more than N unreduced blocks)

*             else NCVT = NRU = NCC = 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, DLASR, DLASV2, DROT,

     $                   DSCAL, DSWAP, XERBLA

*     ..

*     .. 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( 'DBDSQR', -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, WORK, 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 )

            WORK( I ) = CS

            WORK( NM1+I ) = SN

   10    CONTINUE

*

*        Update singular vectors if desired

*

         IF( NRU.GT.0 )

     $      CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U,

     $                  LDU )

         IF( NCC.GT.0 )

     $      CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( 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 DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,

     $                 SINR )

         IF( NRU.GT.0 )

     $      CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )

         IF( NCC.GT.0 )

     $      CALL DROT( 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 ) )

               WORK( I-LL+1 ) = CS

               WORK( I-LL+1+NM1 ) = SN

               WORK( I-LL+1+NM12 ) = OLDCS

               WORK( 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 DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),

     $                     WORK( N ), VT( LL, 1 ), LDVT )

            IF( NRU.GT.0 )

     $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),

     $                     WORK( NM13+1 ), U( 1, LL ), LDU )

            IF( NCC.GT.0 )

     $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),

     $                     WORK( 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 ) )

               WORK( I-LL ) = CS

               WORK( I-LL+NM1 ) = -SN

               WORK( I-LL+NM12 ) = OLDCS

               WORK( 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 DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),

     $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )

            IF( NRU.GT.0 )

     $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),

     $                     WORK( N ), U( 1, LL ), LDU )

            IF( NCC.GT.0 )

     $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),

     $                     WORK( 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

               WORK( I-LL+1 ) = COSR

               WORK( I-LL+1+NM1 ) = SINR

               WORK( I-LL+1+NM12 ) = COSL

               WORK( I-LL+1+NM13 ) = SINL

  140       CONTINUE

            E( M-1 ) = F

*

*           Update singular vectors

*

            IF( NCVT.GT.0 )

     $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),

     $                     WORK( N ), VT( LL, 1 ), LDVT )

            IF( NRU.GT.0 )

     $         CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),

     $                     WORK( NM13+1 ), U( 1, LL ), LDU )

            IF( NCC.GT.0 )

     $         CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),

     $                     WORK( 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

               WORK( I-LL ) = COSR

               WORK( I-LL+NM1 ) = -SINR

               WORK( I-LL+NM12 ) = COSL

               WORK( 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 DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),

     $                     WORK( NM13+1 ), VT( LL, 1 ), LDVT )

            IF( NRU.GT.0 )

     $         CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),

     $                     WORK( N ), U( 1, LL ), LDU )

            IF( NCC.GT.0 )

     $         CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),

     $                     WORK( 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 DSCAL( 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 DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),

     $                     LDVT )

            IF( NRU.GT.0 )

     $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )

            IF( NCC.GT.0 )

     $         CALL DSWAP( 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 DBDSQR

*

      END