SUBROUTINE DSMALLSIGMA(FLAG,METHOD,M,N,A,B,TAUA,TAUB,TAU,
     $     WORK2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7, RTMIN, RTMAX)

      IMPLICIT NONE
      
      DOUBLE PRECISION A(*), B(*), WORK2(*)
      INTEGER N, M, J, METHOD, FLAG
      INTEGER INDRV3, INDRV4, INDRV5, INDRV6, INDRV7
      DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, TMP5
      DOUBLE PRECISION TAUA(4), TAUB, TAU
      DOUBLE PRECISION C1, S1, T, RTMIN, RTMAX
      
      DOUBLE PRECISION ONE, ZERO, HALF, TWO, CONST
      PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0, TWO = 2.0D0)
      PARAMETER (CONST = 0.75D0)
      EXTERNAL DFMA0
      DOUBLE PRECISION DFMA0

      METHOD = 1
      IF (TAUB .LE. TAU) GO TO 160
*         
      TAUB = TAU
*     
      TMP4 = A(M)-TAU
      TMP5 = A(M)+TAU
      IF (TMP4 .LE. ZERO) THEN
         GO TO 160
      ELSE
         TMP3 = SQRT(TMP4)*SQRT(TMP5)
      ENDIF
      DO J = M, N-2
         CALL DLARTG2(TMP3,B(J),C1,S1,T,RTMIN,RTMAX)
         TMP4 = DFMA0(C1,A(J+1),-TAU)
         TMP5 = DFMA0(C1,A(J+1),TAU)
         IF (TMP4 .LE. ZERO) THEN
            GO TO 160
         ELSE
            TMP3 = SQRT(TMP4)*SQRT(TMP5)
         ENDIF
      ENDDO
      CALL DLARTG2(TMP3,B(N-1),C1,S1,T,RTMIN,RTMAX)
      TMP4 = DFMA0(C1,A(N),-TAU)
      IF (TMP4 .LT. ZERO) THEN
         TMP2 = C1*A(N)
         IF (TMP2 .LT. TAU) THEN
            IF (FLAG .EQ. 0) THEN
               TAUA(METHOD+1) = TMP2
            ELSE
               TMP3 = CONST*TAU
               IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
               TAUA(METHOD+1) = MAX(TMP2,TMP3)
            ENDIF
         ELSE
            GO TO 160
         ENDIF
      ELSE
         TAUA(METHOD+1) = TAU
      ENDIF
      IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
         METHOD = METHOD+1
      ENDIF
      IF (METHOD .GE. 2) RETURN
*
 160  WORK2(INDRV7+M:INDRV7+N) = ONE/A(M:N) 

      WORK2(INDRV5+N) = WORK2(INDRV7+N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV3+J) = B(J)*WORK2(INDRV7+J)
         WORK2(INDRV5+J) = WORK2(INDRV7+J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      TMP3 = ONE/TMP3
      WORK2(INDRV5+M) = (WORK2(INDRV5+M)*TMP3)*WORK2(INDRV7+M)
      TMP1 = WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV4+J) = B(J-1)*WORK2(INDRV7+J)
         WORK2(INDRV5+J) = (WORK2(INDRV5+J)*TMP3)*WORK2(INDRV7+J)+
     $        WORK2(INDRV4+J)*WORK2(INDRV5+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV5+J))
      ENDDO
      IF (TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(SQRT(TMP3/TMP1),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      TMP1 = ONE/TMP1
      WORK2(INDRV5+N) = WORK2(INDRV5+N)*TMP1
      WORK2(INDRV6+N) = WORK2(INDRV5+N)*WORK2(INDRV7+N)
      TMP3 = WORK2(INDRV6+N)
      DO J = N-1,M,-1
         WORK2(INDRV5+J) = WORK2(INDRV5+J)*TMP1
         WORK2(INDRV6+J) = WORK2(INDRV5+J)*WORK2(INDRV7+J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV6+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV6+J))
      ENDDO
      TMP3 = ONE/TMP3
      TMP2 = (WORK2(INDRV6+M)*TMP3)*WORK2(INDRV7+M)
      TMP1 = TMP2/WORK2(INDRV5+M)
      DO J = M+1,N
         TMP2 = (WORK2(INDRV6+J)*TMP3)*WORK2(INDRV7+J)
     $        +WORK2(INDRV4+J)*TMP2
         TMP1 = MAX(TMP1,TMP2/WORK2(INDRV5+J))
      ENDDO
      IF (TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(SQRT(TMP3/TMP1),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      IF (METHOD .GE. 2) RETURN
      
 170  WORK2(INDRV5+N) = ONE/A(N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV3+J) = B(J)/A(J)
         WORK2(INDRV5+J) = ONE/A(J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      WORK2(INDRV5+M) = (WORK2(INDRV5+M)/TMP3)/A(M)
      TMP1 = WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV4+J) = B(J-1)/A(J)
         WORK2(INDRV5+J) = (WORK2(INDRV5+J)/TMP3)/A(J)+
     $        WORK2(INDRV4+J)*WORK2(INDRV5+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV5+J))
      ENDDO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(ONE/SQRT(TMP1)/SQRT(TMP3),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      WORK2(INDRV5+N) = WORK2(INDRV5+N)/TMP1
      WORK2(INDRV6+N) = WORK2(INDRV5+N)/A(N)
      TMP3 = WORK2(INDRV6+N)
      DO J = N-1,M,-1
         WORK2(INDRV5+J) = WORK2(INDRV5+J)/TMP1
         WORK2(INDRV6+J) = WORK2(INDRV5+J)/A(J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV6+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV6+J))
      ENDDO
      TMP2 = (WORK2(INDRV6+M)/TMP3)/A(M)
      TMP1 = TMP2/WORK2(INDRV5+M)
      DO J = M+1,N
         TMP2 = (WORK2(INDRV6+J)/TMP3)/A(J)+
     $        WORK2(INDRV4+J)*TMP2
         TMP1 = MAX(TMP1,TMP2/WORK2(INDRV5+J))
      ENDDO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(ONE/SQRT(TMP1)/SQRT(TMP3),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF

      IF (METHOD .GE. 2) RETURN

      TAUA(METHOD+1) = A(N)-HALF*B(N-1)
      IF (TAUA(METHOD+1) .LE. ZERO) RETURN
      DO J=N-1,M+1,-1
         TMP2 = A(J)-HALF*(B(J)+B(J-1))
         IF (TMP2 .LE. ZERO) RETURN
         TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
      ENDDO
      TMP2=A(M)-HALF*B(M)
      IF (TMP2 .LE. ZERO) RETURN
      TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
      IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
         METHOD = METHOD+1
      ENDIF
      
      RETURN
      END SUBROUTINE DSMALLSIGMA

      SUBROUTINE DOQDS(UPLO,N0,A,B,SU,LDSU,SV,LDSV,WORK,WORK2,INFO)
      
      IMPLICIT NONE
      
      CHARACTER UPLO
      INTEGER N0, INFO, LDSU, LDSV, ISUB, MAXITER, FLAG
      DOUBLE PRECISION A(*), B(*), WORK2(*)
      DOUBLE PRECISION SU(LDSU, *), SV(LDSV, *), WORK(N0,*)
      INTEGER N, M, I, J, K, OLDM, OLDN
      INTEGER INDRV1, INDRV2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7,
     $     INDRV8, METHOD
      DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, TMP5
      DOUBLE PRECISION TAU, TAUA(4), TAUB
      DOUBLE PRECISION SIGMA, SIGMA1, SIGMA2, DESIG, T, DESIG0, S
      DOUBLE PRECISION C1, S1, C2, S2, SMIN, EPS, TOL, SAFMIN
      DOUBLE PRECISION SAFMAX, RTMIN, RTMAX

      DOUBLE PRECISION ONE, ZERO, HALF, TWO, CONST
      PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0, TWO = 2.0D0)
      PARAMETER (CONST = 0.75D0)
      DOUBLE PRECISION HUNDRD 
      PARAMETER (HUNDRD = 100.0D0)
*     
      EXTERNAL DLAMCH
      DOUBLE PRECISION DLAMCH
      EXTERNAL DFMA0
      DOUBLE PRECISION DFMA0
      LOGICAL LSAME
      EXTERNAL LSAME
*
      TAUA(1) = ZERO
*     
      INDRV1 = 0
      INDRV2 = INDRV1+N0
      INDRV3 = INDRV2+N0
      INDRV4 = INDRV3+N0
      INDRV5 = INDRV4+N0
      INDRV6 = INDRV5+N0
      INDRV7 = INDRV6+N0
      INDRV8 = INDRV7+N0
*     
      EPS = DLAMCH( 'Precision' )
      TOL = HUNDRD*EPS
      SAFMIN = DLAMCH( 'S' )
      RTMIN = SQRT( SAFMIN )
      SAFMAX = ONE/SAFMIN
      RTMAX = SQRT( SAFMAX/2.0D0 )
*     
      IF ( LSAME( UPLO, 'U' ) ) THEN
         
         DO I = 1, N0-1
            IF (A(I) .LT. ZERO) THEN
               A(I) = -A(I)
               B(I) = -B(I)
               CALL DSCAL(N0,-ONE,SU(1,I),1)
            ENDIF
            IF (B(I) .LT. ZERO) THEN
               B(I) = -B(I)
               A(I+1) = -A(I+1)
               CALL DSCAL(N0,-ONE,SV(1,I+1),1)
            ENDIF
         ENDDO
         IF (A(N0) .LT. ZERO) THEN
            A(N0) = -A(N0)
            CALL DSCAL(N0,-ONE,SU(1,N0),1)
         ENDIF
         
         K=N0/2
         DO J=1,K
            TMP1=A(J)
            A(J)=A(N0+1-J)
            A(N0+1-J)=TMP1
            TMP1=B(J)
            B(J)=B(N0-J)
            B(N0-J)=TMP1
            CALL DSWAP(N0,SU(1,J),1,SU(1,N0+1-J),1)
            CALL DSWAP(N0,SV(1,J),1,SV(1,N0+1-J),1)
         ENDDO
            
      ELSE
*     
         DO I = 1, N0-1
            IF (A(I) .LT. ZERO) THEN
               A(I) = -A(I)
               B(I) = -B(I)
               CALL DSCAL(N0,-ONE,SV(1,I),1)
            ENDIF
            IF (B(I) .LT. ZERO) THEN
               B(I) = -B(I)
               A(I+1) = -A(I+1)
               CALL DSCAL(N0,-ONE,SU(1,I+1),1)
            ENDIF
         ENDDO
         IF (A(N0) .LT. ZERO) THEN
            A(N0) = -A(N0)
            CALL DSCAL(N0,-ONE,SV(1,N0),1)
         ENDIF
*
      ENDIF
*
      OLDM = -1
      OLDN = -1
*     
      DO I = 1, N0
         DO J = 1, N0
            WORK(J,I)=ZERO
         ENDDO
      ENDDO
*     
      M = 1
      N = N0
*     
      B(N) = ZERO
      WORK2(INDRV6+N) = ZERO
*     
 3000 SIGMA = -B(N)
      SIGMA1 = SQRT(EPS)*SIGMA
      SIGMA2 = TOL*SIGMA
      DESIG = -WORK2(INDRV6+N)
         
      MAXITER = 30*(N-M+1)
      DO I = 1, MAXITER
*     
 15      IF (N-M+1 .EQ. 1) THEN
            
            CALL DLARTG2(A(N),SIGMA,C1,S1,TMP1,RTMIN,RTMAX)

            CALL DROT2(N0,SU(1,N),1,WORK(1,N),1,C1,S1)
            
            CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            
            GO TO 700
            
         ENDIF

         IF (N-M+1 .EQ. 2) THEN
            CALL DLASV2(A(N-1), B(N-1), A(N), A(N), A(N-1), S1, C1, 
     $           S2, C2)
            CALL DROT2(N0,SV(1,N-1),1,SV(1,N),1,C2,S2)
            CALL DROT2(N0,SU(1,N-1),1,SU(1,N),1,C1,S1)
            CALL DROT2(N0,WORK(1,N-1),1,WORK(1,N),1,C2,S2)
            CALL DLARTG2(A(N-1),SIGMA,C1,S1,TMP1,RTMIN,RTMAX)
            CALL DROT2(N0,SU(1,N-1),1,WORK(1,N-1),1,C1,S1)
            CALL DLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0,RTMIN)
            CALL DLARTG2(A(N),SIGMA,C1,S1,TMP1,RTMIN,RTMAX)
            CALL DROT2(N0,SU(1,N),1,WORK(1,N),1,C1,S1)
            CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            GO TO 700
         ENDIF
*
         IF (M .GT. OLDN .OR. N .LT. OLDM) THEN
            IF ( A(M) .LT. A(N) ) THEN

               TMP1 = A(N)
               DO J = N-1, M, -1
                  IF (B(J) .LE. MAX(SIGMA2,EPS*TMP1)) THEN
                     B(J) = -SIGMA
                     WORK2(INDRV6+J) = -DESIG
                     M = J+1
                     GO TO 15
                  ELSE
                     CALL DLARTG2(TMP1,B(J),C1,S1,TMP2,RTMIN,RTMAX)
                     TMP1 = A(J)*C1
                  ENDIF
               ENDDO
               
               DO J=M,N
                  WORK2(INDRV5+N+M-J)=A(J)
               ENDDO
	       DO J=M,N-1
                  WORK2(INDRV6+N-1+M-J)=B(J)
               ENDDO

               K=(N-M+1)/2+M-1
               DO J=M,K
                  CALL DSWAP(N0,SU(1,J),1,SU(1,N+M-J),1)
                  CALL DSWAP(N0,WORK(1,J),1,WORK(1,N+M-J),1)
                  CALL DSWAP(N0,SV(1,J),1,SV(1,N+M-J),1)
               ENDDO

               TMP1 = WORK2(INDRV5+M)
               DO J = M, N-1
                  CALL DLARTG2(TMP1,WORK2(INDRV6+J),C1,S1,A(J),
     $                 RTMIN,RTMAX)
                  WORK2(INDRV3+J) = C1
                  WORK2(INDRV4+J) = S1
                  B(J) = S1*WORK2(INDRV5+J+1)
	          TMP1 = C1*WORK2(INDRV5+J+1)
               ENDDO
               A(N) = TMP1
*     
               DO J = M, N-1
                  CALL DROT2(N0, WORK( 1, J), 1, WORK( 1, J+1), 1,
     $                 WORK2( INDRV3+J ), WORK2( INDRV4+J ))
               ENDDO
*     
               DO J = M, N-1
                  CALL DROT2(N0, SV( 1, J), 1, SV( 1, J+1), 1,
     $                 WORK2( INDRV3+J ), WORK2( INDRV4+J ))
               ENDDO
               
               GO TO 400
            ENDIF
         ENDIF
*
         TMP1 = A(M)
         DO J = M, N-1
            IF (B(J) .LE. MAX(SIGMA2,EPS*TMP1)) THEN
               B(J) = -SIGMA
               WORK2(INDRV6+J) = -DESIG
               M = J+1
               GO TO 15
            ELSE
               CALL DLARTG2(TMP1,B(J),C1,S1,TMP2,RTMIN,RTMAX)
               TMP1 = A(J+1)*C1
            ENDIF
         ENDDO
*
         FLAG = -1
         TAU = A(M)
         IF (TAU .LE. SIGMA1) THEN
            CALL DLARTG2(SIGMA,TAU,C1,S1,T,RTMIN,RTMAX)
            FLAG = M-1
            WORK2(INDRV1+M) = C1
            WORK2(INDRV2+M) = -S1
            TAU = ZERO
         ELSE
            WORK2(INDRV1+M) = ONE
            WORK2(INDRV2+M) = ZERO
         ENDIF
         TAUB = TAU
         DO J = M, N-2
            CALL DLARTG2(TAU,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
            WORK2(INDRV7+J) = C1
            WORK2(INDRV8+J) = S1
            WORK2(INDRV6+J) = S1*A(J+1)
            TAU = C1*A(J+1)
            IF (TAU .LE. SIGMA1) THEN
               CALL DLARTG2(SIGMA,TAU,C1,S1,T,RTMIN,RTMAX)
               IF (FLAG .EQ. -1) FLAG = J
               WORK2(INDRV1+J+1) = C1
               WORK2(INDRV2+J+1) = -S1
               TAU = ZERO
            ELSE
               WORK2(INDRV1+J+1) = ONE
               WORK2(INDRV2+J+1) = ZERO
            ENDIF
            TAUB = MIN(TAUB,TAU)
         ENDDO
         CALL DLARTG2(TAU,B(N-1),C1,S1,WORK2(INDRV5+N-1),RTMIN,RTMAX)
         WORK2(INDRV7+N-1) = C1
         WORK2(INDRV8+N-1) = S1
         WORK2(INDRV6+N-1) = S1*A(N)
         TAU = C1*A(N)
         IF (TAU .LE. SIGMA1) THEN
            CALL DLARTG2(SIGMA,TAU,C1,S1,T,RTMIN,RTMAX)
            IF (FLAG .EQ. -1) FLAG = N-1
            WORK2(INDRV1+N) = C1
            WORK2(INDRV2+N) = -S1
            TAU = ZERO
         ELSE
            WORK2(INDRV1+N) = ONE
            WORK2(INDRV2+N) = ZERO
         ENDIF
         WORK2(INDRV5+N) = TAU
         IF (FLAG .GE. 0) GO TO 360

         CALL DLAS2(A(N-1), B(N-1), A(N), TMP2, TMP3)
         TAU = MIN(TMP2,TAU)
         CALL DSMALLSIGMA(1,METHOD,M,N,A,B,TAUA,TAUB,TAU,
     $        WORK2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7, 
     $        RTMIN, RTMAX)
         TAU = TAUA(METHOD)
*     
 125     IF (TAU .EQ. ZERO) GO TO 350
*     
         IF (SIGMA .EQ. ZERO) THEN
*     
            T = TAU
            DESIG0 = ZERO
*     
            TMP4 = A(M)-TAU
            TMP5 = A(M)+TAU
            IF (TMP4 .LE. ZERO) THEN
               TMP2 = A(M)
               IF (TMP2 .GE. TAU) TMP2 = ZERO
               TMP3 = CONST*TAU
               IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
               TMP2 = MAX(TMP2,TMP3)
               IF (METHOD .GE. 2) THEN
                  IF (TAUA(METHOD-1) .GE. TMP2) THEN
                     METHOD = METHOD-1
                     TAU = TAUA(METHOD)
                  ELSE
                     TAU = TMP2
                  ENDIF
               ELSE
                  TAU = TMP2
               ENDIF
               GO TO 125
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            CALL DLARTG2(TMP3,TAU,C1,S1,S,RTMIN,RTMAX)
            
            WORK2(INDRV1+M) = C1
            WORK2(INDRV2+M) = -S1
            DO J = M, N-2
               CALL DLARTG2(TMP3,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
               WORK2(INDRV7+J) = C1
               WORK2(INDRV8+J) = S1
               WORK2(INDRV6+J) = S1*A(J+1)
               TMP4 = DFMA0(C1,A(J+1),-TAU)
               TMP5 = DFMA0(C1,A(J+1),TAU)
               IF (TMP4 .LE. ZERO) THEN
                  TMP2 = C1*A(J+1)
                  IF (TMP2 .GE. TAU) TMP2 = ZERO
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ELSE
                  TMP3 = SQRT(TMP4)*SQRT(TMP5)
               ENDIF
               
               CALL DLARTG2(TMP3,TAU,C1,S1,S,RTMIN,RTMAX)
               WORK2(INDRV1+J+1) = C1
               WORK2(INDRV2+J+1) = -S1
            ENDDO
            CALL DLARTG2(TMP3,B(N-1),C1,S1,WORK2(INDRV5+N-1),
     $           RTMIN,RTMAX)
            WORK2(INDRV7+N-1) = C1
            WORK2(INDRV8+N-1) = S1
            WORK2(INDRV6+N-1) = S1*A(N)
            TMP4 = DFMA0(C1,A(N),-TAU)
            TMP5 = DFMA0(C1,A(N),TAU)
            IF (TMP4 .LT. ZERO) THEN
               TMP2 = C1*A(N)
               IF (TMP2 .LT. TAU) THEN
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ENDIF
               TMP3 = ZERO
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            IF (TMP3 .GE. A(N)) THEN
               TMP2 = TAU/A(N)
               IF (C1 .GE. TMP2) THEN
                  TMP3 = MIN(A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2),A(N))
               ELSE
                  TMP3 = A(N)
               ENDIF
            ENDIF
            
            CALL DLARTG2(TMP3,TAU,C1,S1,S,RTMIN,RTMAX)
            WORK2(INDRV1+N) = C1
            WORK2(INDRV2+N) = -S1
            WORK2(INDRV5+N) = TMP3
*     
         ELSE
*     
            CALL DLARTG7(SIGMA,DESIG,TAU,T,DESIG0,RTMIN)
*     
            TMP4 = A(M)-TAU
            TMP5 = A(M)+TAU
            IF (TMP4 .LE. ZERO) THEN
               TMP2 = A(M)
               IF (TMP2 .GE. TAU) TMP2 = ZERO
               TMP3 = CONST*TAU
               IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
               TMP2 = MAX(TMP2,TMP3)
               IF (METHOD .GE. 2) THEN
                  IF (TAUA(METHOD-1) .GE. TMP2) THEN
                     METHOD = METHOD-1
                     TAU = TAUA(METHOD)
                  ELSE
                     TAU = TMP2
                  ENDIF
               ELSE
                  TAU = TMP2
               ENDIF
               GO TO 125
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            CALL DLARTG6(A(M),TMP3,SIGMA,T,C1,S1,rtmin,rtmax)
            WORK2(INDRV1+M) = C1
            WORK2(INDRV2+M) = S1
            DO J = M, N-2
               CALL DLARTG2(TMP3,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
               WORK2(INDRV7+J) = C1
               WORK2(INDRV8+J) = S1
               WORK2(INDRV6+J) = S1*A(J+1)
               
               TMP4 = DFMA0(C1,A(J+1),-TAU)
               TMP5 = DFMA0(C1,A(J+1),TAU)
               IF (TMP4 .LE. ZERO) THEN
                  TMP2 = C1*A(J+1)
                  IF (TMP2 .GE. TAU) TMP2 = ZERO
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ELSE
                  TMP3 = SQRT(TMP4)*SQRT(TMP5)
               ENDIF
               
               CALL DLARTG6(C1*A(J+1),TMP3,SIGMA,T,C1,S1,rtmin,rtmax)
               WORK2(INDRV1+J+1) = C1
               WORK2(INDRV2+J+1) = S1
            ENDDO
            CALL DLARTG2(TMP3,B(N-1),C1,S1,WORK2(INDRV5+N-1),
     $           RTMIN,RTMAX)
            WORK2(INDRV7+N-1) = C1
            WORK2(INDRV8+N-1) = S1
            WORK2(INDRV6+N-1) = S1*A(N)
            
            TMP4 = DFMA0(C1,A(N),-TAU)
            TMP5 = DFMA0(C1,A(N),TAU)
            IF (TMP4 .LT. ZERO) THEN
               TMP2 = C1*A(N)
               IF (TMP2 .LT. TAU) THEN
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ENDIF
               TMP3 = ZERO
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            IF (TMP3 .GE. A(N)) THEN
               TMP2 = TAU/A(N)
               IF (C1 .GE. TMP2) THEN
                  TMP3 = MIN(A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2),A(N))
               ELSE
                  TMP3 = A(N)
               ENDIF
            ENDIF
            
            IF (TMP3 .EQ. ZERO) THEN
               CALL DLARTG2(SIGMA,C1*A(N),C1,S1,S,RTMIN,RTMAX)
               S1 = -S1
            ELSE
               CALL DLARTG6(C1*A(N),TMP3,SIGMA,T,C1,S1,rtmin,rtmax)
            ENDIF
            WORK2(INDRV1+N) = C1
            WORK2(INDRV2+N) = S1
            WORK2(INDRV5+N) = TMP3

         ENDIF
*     
         SIGMA = T
         SIGMA1 = SQRT(EPS)*SIGMA
         SIGMA2 = TOL*SIGMA
         DESIG = DESIG0
*     
         TMP1 = WORK2(INDRV5+M)
         DO J = M, N-1
            CALL DLARTG2(TMP1,WORK2(INDRV6+J),C1,S1,A(J),RTMIN,RTMAX)
            WORK2(INDRV3+J) = C1
            WORK2(INDRV4+J) = S1
            B(J) = S1*WORK2(INDRV5+J+1)
            TMP1 = C1*WORK2(INDRV5+J+1)
         ENDDO
         A(N) = TMP1
*     
         CALL DROT2(N0,SU(1,M),1,WORK(1,M),1,
     $        WORK2(INDRV1+M),WORK2(INDRV2+M))
         DO J = M, N-1
            CALL DROT2(N0,SU(1,J),1,SU(1,J+1),1,WORK2(INDRV7+J),
     $           WORK2(INDRV8+J))
            CALL DROT2(N0,SU(1,J+1),1,WORK(1,J+1),1,
     $           WORK2(INDRV1+J+1),WORK2(INDRV2+J+1))
            CALL DROT2(N0,WORK(1,J),1,WORK(1,J+1),1,
     $           WORK2(INDRV3+J),WORK2(INDRV4+J))
         ENDDO
*     
         DO J = M, N-1
            CALL DROT2(N0, SV( 1, J), 1, SV( 1, J+1), 1,
     $           WORK2( INDRV3+J ), WORK2( INDRV4+J ))
         ENDDO
*
         GO TO 400
*     
 350     FLAG = -1
         TMP1 = A(M)
         IF (TMP1 .LE. SIGMA1) THEN
            CALL DLARTG2(SIGMA,TMP1,C1,S1,T,RTMIN,RTMAX)
            FLAG = M-1
            WORK2(INDRV1+M) = C1
            WORK2(INDRV2+M) = -S1
            TMP1 = ZERO
         ELSE
            WORK2(INDRV1+M) = ONE
            WORK2(INDRV2+M) = ZERO
         ENDIF
         DO J = M, N-1
            CALL DLARTG2(TMP1,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
            WORK2(INDRV7+J) = C1
            WORK2(INDRV8+J) = S1
            WORK2(INDRV6+J) = S1*A(J+1)
            TMP1 = C1*A(J+1)
            IF (TMP1 .LE. SIGMA1) THEN
               CALL DLARTG2(SIGMA,TMP1,C1,S1,T,RTMIN,RTMAX)
               IF (FLAG .EQ. -1) FLAG = J
               WORK2(INDRV1+J+1) = C1
               WORK2(INDRV2+J+1) = -S1
               TMP1 = ZERO
            ELSE
               WORK2(INDRV1+J+1) = ONE
               WORK2(INDRV2+J+1) = ZERO
            ENDIF
         ENDDO
         WORK2(INDRV5+N) = TMP1
*
 360     TMP1 = WORK2(INDRV5+M)
         DO J = M, N-1
            CALL DLARTG2(TMP1,WORK2(INDRV6+J),C1,S1,A(J),RTMIN,RTMAX)
            WORK2(INDRV3+J) = C1
            WORK2(INDRV4+J) = S1
            B(J) = S1*WORK2(INDRV5+J+1)
            TMP1 = C1*WORK2(INDRV5+J+1)
         ENDDO
         A(N) = TMP1
*     
         IF (FLAG .EQ. -1) THEN

            DO J = M, N-1
               CALL DROT2(N0, SU( 1, J), 1, SU( 1, J+1), 1,
     $              WORK2( INDRV7+J ), WORK2( INDRV8+J ))
            ENDDO

         ELSE

            DO J = M, FLAG
               CALL DROT2(N0, SU( 1, J), 1, SU( 1, J+1), 1,
     $              WORK2( INDRV7+J ), WORK2( INDRV8+J ))
            ENDDO

            CALL DROT2(N0,SU(1,FLAG+1),1,WORK(1,FLAG+1),1,
     $           WORK2(INDRV1+FLAG+1),WORK2(INDRV2+FLAG+1))

            DO J = FLAG+1, N-1
               CALL DROT2(N0, SU( 1, J), 1, SU( 1, J+1), 1,
     $              WORK2( INDRV7+J ), WORK2( INDRV8+J ))
            ENDDO

         ENDIF
*
         DO J = M, N-1
            CALL DROT2(N0, WORK( 1, J), 1, WORK( 1, J+1), 1,
     $           WORK2( INDRV3+J ), WORK2( INDRV4+J ))
         ENDDO
*     
         DO J = M, N-1
            CALL DROT2(N0, SV( 1, J), 1, SV( 1, J+1), 1,
     $           WORK2( INDRV3+J ), WORK2( INDRV4+J ))
         ENDDO
*     
 400     OLDM = M
         OLDN = N
*     
      ENDDO
      INFO = 2
      RETURN
      
 700  N = M-1
      IF (N .LE. 0) THEN
         GO TO 4000
      ENDIF
*     
      DO J = M-2,1,-1
         IF (B(J) .LE. ZERO) THEN
            M = J+1
            GO TO 3000
         ENDIF
      ENDDO
      M = 1
      GO TO 3000
*     
 4000 DO I = 1, N0 - 1
*     
*     Scan for smallest D(I)
*     
         ISUB = 1
         SMIN = A( 1 )
         DO J = 2, N0 + 1 - I
            IF( A( J ).LE.SMIN ) THEN
               ISUB = J
               SMIN = A( J )
            END IF
         ENDDO
         IF( ISUB.NE.N0+1-I ) THEN
*     
*     Swap singular values and vectors
*     
            A( ISUB ) = A( N0+1-I )
            A( N0+1-I ) = SMIN
            CALL DSWAP( N0, SV( 1, ISUB ), 1, SV( 1, N0+1-I ), 1 ) 
            CALL DSWAP( N0, SU( 1, ISUB ), 1, SU( 1, N0+1-I ), 1 )
         END IF
      ENDDO
*
 30   INFO = 0
      RETURN
*     
      END