SUBROUTINE DSMALLSIGMA(M,N,A,B,TAU,WORK2,
     $     INDRV3, INDRV4, INDRV5, INDRV6)

      IMPLICIT NONE
      
      DOUBLE PRECISION A(*), B(*), WORK2(*)
      INTEGER N, M, J
      INTEGER INDRV3, INDRV4, INDRV5, INDRV6
      DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, TMP5
      DOUBLE PRECISION TAU, TAU1, TAU2
      DOUBLE PRECISION C1, S1, T
      
      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

      CALL DLAS2U(A(N-1), B(N-1), A(N), TAU, TMP3)
      TAU = MIN(TAU,A(N))
         
      TAU2 = MINVAL(A(M:N-1))
      IF (TAU2 .LE. TAU) THEN
         TAU1 = TAU2
         GO TO 160
      ELSE
         TAU1 = TAU
      ENDIF
*     
      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 DLARTG(TMP3,B(J),C1,S1,T)
         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 DLARTG(TMP3,B(N-1),C1,S1,T)
      TMP4 = DFMA0(C1,A(N),-TAU)
      IF (TMP4 .LT. ZERO .AND. TAU+TMP4 .EQ. TAU) TMP4 = ZERO
      IF (TMP4 .LT. ZERO) THEN
         TAU2 = MAX(TMP4+TAU,ZERO)
         IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU
         IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU
         TAU = TAU2
      ENDIF
      IF (TAU .GT. ZERO) RETURN
*     
 160  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
      TAU = ZERO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAU = MAX(TAU,ONE/SQRT(TMP1)/SQRT(TMP3))
      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
         TAU = MAX(TAU,ONE/SQRT(TMP1)/SQRT(TMP3))
      ENDIF
      
      TAU = MIN(TAU,TAU1)
      IF (TAU .GT. ZERO) RETURN
      
      TAU=A(N)-HALF*B(N-1)
      IF (TAU .LE. ZERO) THEN
         TAU = ZERO
         RETURN
      ENDIF
      DO J=N-1,M+1,-1
         TMP2=A(J)-HALF*(B(J)+B(J-1))
         IF (TMP2 .LE. ZERO) THEN
            TAU = ZERO
            RETURN
         ENDIF
         TAU=MIN(TAU,TMP2)
      ENDDO
      TMP2=A(M)-HALF*B(M)
      IF (TMP2 .LE. ZERO) THEN
         TAU = ZERO
         RETURN
      ENDIF
      TAU=MIN(TAU,TMP2)

      RETURN
      END SUBROUTINE DSMALLSIGMA

      SUBROUTINE DOQDSQ(L,N0,A,B,SU,LDSU,WORK2,INFO)

      IMPLICIT NONE
      
      INTEGER N0, INFO, LDSU, MAXITER, FLAG
      DOUBLE PRECISION A(*), B(*), WORK2(*)
      DOUBLE PRECISION SU(LDSU, *)
      INTEGER N, M, I, J, K, L, OLDM, OLDN
      INTEGER INDRV3, INDRV4, INDRV5, INDRV6
      DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, TMP5
      DOUBLE PRECISION TAU, TAU2
      DOUBLE PRECISION SIGMA, SIGMA1, DESIG, T, DESIG0
      DOUBLE PRECISION C1, S1, C2, S2, EPS, TOL
      
      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
*     
      INDRV3 = 0
      INDRV4 = INDRV3+N0
      INDRV5 = INDRV4+N0
      INDRV6 = INDRV5+N0
*     
      EPS = DLAMCH( 'Precision' )
      TOL = HUNDRD*EPS
*     
      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)
         ENDIF
      ENDDO
      IF (A(N0) .LT. ZERO) THEN
         A(N0) = -A(N0)
         CALL DSCAL(N0,-ONE,SU(1,N0),1)
      ENDIF
*
      OLDM = -1
      OLDN = -1
*     
      M = 1
      N = N0
*
      B(N) = ZERO
      WORK2(INDRV6+N) = ZERO
*     
 3000 SIGMA = -B(N)
      SIGMA1 = SQRT(EPS)*SIGMA
      DESIG = -WORK2(INDRV6+N)
         
      MAXITER = 30*(N-M+1)
      DO I = 1, MAXITER
*
 15      IF (N .LE. N0-L) GO TO 4000
*         
         IF (N-M+1 .EQ. 1) THEN
            
            CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0)
            
            GO TO 700
            
         ENDIF

         IF (N-M+1 .EQ. 2) THEN

            IF (A(M) .GE. A(N)) THEN

               IF (A(M) .EQ. B(N-1)+A(M)) THEN
                  B(N-1) = -SIGMA
                  WORK2(INDRV6+N-1) = -DESIG
                  M = N
                  CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0)
                  GO TO 700
               ENDIF

            ELSE

               IF (A(N) .EQ. B(N-1)+A(N)) THEN
                  B(N-1) = -SIGMA
                  WORK2(INDRV6+N-1) = -DESIG
                  M = N
                  CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0)
                  GO TO 700
               ENDIF

            ENDIF

         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,SU(1,N-1),1,SU(1,N),1,C2,S2)
            
            CALL DLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0)
            
            CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0)
            
            GO TO 700
         ENDIF
*
         IF (M .GT. OLDN .OR. N .LT. OLDM) THEN
            IF ( A(M) .LT. A(N) ) THEN

               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)
               ENDDO

               TMP1 = WORK2(INDRV5+M)
               DO J = M, N-1
                  CALL DLARTG(TMP1,WORK2(INDRV6+J),C1,S1,A(J))
                  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, SU( 1, J), 1, SU( 1, J+1), 1,
     $                 WORK2( INDRV3+J ), WORK2( INDRV4+J ))
               ENDDO
               
               GO TO 400
            ENDIF
         ENDIF
*
         TMP1 = A(M)
         DO J = M, N-2
            TMP1 = A(J+1)*(TMP1/(TMP1+B(J)))
         ENDDO

         IF (B(N-1) .LE. TOL*SIGMA .OR.
     $        TMP1 .EQ. TMP1+B(N-1)) THEN

            IF(A(N) .EQ. ZERO) THEN
               B(N-1) = -SIGMA
               WORK2(INDRV6+N-1) = -DESIG
               M = N
               A(N) = SIGMA
               GO TO 700
            ENDIF

            CALL DLARTG7(SIGMA,DESIG,A(N),T,DESIG0)
            IF(T .LE. SIGMA) THEN
               B(N-1) = -SIGMA
               WORK2(INDRV6+N-1) = -DESIG
               M = N
               A(N) = T
               GO TO 700
            ENDIF
            
            CALL DSMALLSIGMA(M,N-1,A,B,TMP2,WORK2,
     $           INDRV3, INDRV4, INDRV5, INDRV6)

            IF(A(N) .LE. TMP2) THEN
               B(N-1) = -SIGMA
               WORK2(INDRV6+N-1) = -DESIG
               M = N
               A(N) = T
               GO TO 700
            ENDIF
         ENDIF
*
         FLAG = 0
         TMP1 = A(M)
         IF (TMP1 .LE. SIGMA1) THEN
            TMP1 = ZERO
            FLAG = 1
         ENDIF
         DO J = M, N-1
            IF (TMP1 .EQ. ZERO .AND. B(J) .EQ. ZERO) THEN
               C1 = ZERO
               S1 = ONE
               WORK2(INDRV5+J) = ZERO
            ELSE
               CALL DLARTG(TMP1,B(J),C1,S1,WORK2(INDRV5+J))
            ENDIF
            WORK2(INDRV6+J) = S1*A(J+1)
            TMP1 = C1*A(J+1)
            IF (TMP1 .LE. SIGMA1) THEN
               TMP1 = ZERO
               FLAG = 1
            ENDIF
         ENDDO
         WORK2(INDRV5+N) = TMP1
         IF (FLAG .EQ. 1) GO TO 360
*         
         CALL DSMALLSIGMA(M,N,A,B,TAU,WORK2,
     $        INDRV3, INDRV4, INDRV5, INDRV6)
*     
 125     IF (TAU .LE. 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
               IF (TMP4 .EQ. ZERO .AND. B(M) .EQ. ZERO) THEN
                  TMP3 = ZERO
               ELSE
                  TAU2 = MAX(TMP4+TAU,ZERO)
                  IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU
                  IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU
                  TAU = TAU2
                  GO TO 125
               ENDIF
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            DO J = M, N-2
               CALL DLARTG(TMP3,B(J),C1,S1,WORK2(INDRV5+J))
               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
                  IF (TMP4 .EQ. ZERO .AND. B(J+1) .EQ. ZERO) THEN
                     TMP3 = ZERO
                  ELSE
                     TAU2 = MAX(TMP4+TAU,ZERO)
                     IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU
                     IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU
                     TAU = TAU2
                     GO TO 125
                  ENDIF
               ELSE
                  TMP3 = SQRT(TMP4)*SQRT(TMP5)
               ENDIF
               
            ENDDO
            CALL DLARTG(TMP3,B(N-1),C1,S1,WORK2(INDRV5+N-1))
            WORK2(INDRV6+N-1) = S1*A(N)
            
            TMP4 = DFMA0(C1,A(N),-TAU)
            TMP5 = DFMA0(C1,A(N),TAU)
            IF (TMP4 .LT. ZERO .AND. TAU+TMP4 .EQ. TAU) TMP4 = ZERO
            IF (TMP4 .LT. ZERO) THEN
               TAU2 = MAX(TMP4+TAU,ZERO)
               IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU
               IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU
               TAU = TAU2
               GO TO 125
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            IF (TMP3 .GE. A(N)) THEN
               TMP2 = TAU/A(N)
               IF (C1 .GE. TMP2) THEN
                  TMP3 = A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2)
               ENDIF
            ENDIF
            WORK2(INDRV5+N) = TMP3
*     
         ELSE
*     
            CALL DLARTG7(SIGMA,DESIG,TAU,T,DESIG0)
*     
            TMP4 = A(M)-TAU
            TMP5 = A(M)+TAU
            IF (TMP4 .LE. ZERO) THEN
               IF (TMP4 .EQ. ZERO .AND. B(M) .EQ. ZERO) THEN
                  TMP3 = ZERO
               ELSE
                  TAU2 = MAX(TMP4+TAU,ZERO)
                  IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU
                  IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU
                  TAU = TAU2
                  GO TO 125
               ENDIF
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            DO J = M, N-2
               CALL DLARTG(TMP3,B(J),C1,S1,WORK2(INDRV5+J))
               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
                  IF (TMP4 .EQ. ZERO .AND. B(J+1) .EQ. ZERO) THEN
                     TMP3 = ZERO
                  ELSE
                     TAU2 = MAX(TMP4+TAU,ZERO)
                     IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU
                     IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU
                     TAU = TAU2
                     GO TO 125
                  ENDIF
               ELSE
                  TMP3 = SQRT(TMP4)*SQRT(TMP5)
               ENDIF
               
            ENDDO
            CALL DLARTG(TMP3,B(N-1),C1,S1,WORK2(INDRV5+N-1))
            WORK2(INDRV6+N-1) = S1*A(N)
            
            TMP4 = DFMA0(C1,A(N),-TAU)
            TMP5 = DFMA0(C1,A(N),TAU)
            IF (TMP4 .LT. ZERO .AND. TAU+TMP4 .EQ. TAU) TMP4 = ZERO
            IF (TMP4 .LT. ZERO) THEN
               TAU2 = MAX(TMP4+TAU,ZERO)
               IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU
               IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU
               TAU = TAU2
               GO TO 125
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            IF (TMP3 .GE. A(N)) THEN
               TMP2 = TAU/A(N)
               IF (C1 .GE. TMP2) THEN
                  TMP3 = A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2)
               ENDIF
            ENDIF
            WORK2(INDRV5+N) = TMP3
*     
         ENDIF
*     
         SIGMA = T
         SIGMA1 = SQRT(EPS)*SIGMA
         DESIG = DESIG0
*     
         TMP1 = WORK2(INDRV5+M)
         DO J = M, N-1
            CALL DLARTG(TMP1,WORK2(INDRV6+J),C1,S1,A(J))
            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, SU( 1, J), 1, SU( 1, J+1), 1,
     $           WORK2( INDRV3+J ), WORK2( INDRV4+J ))
         ENDDO
*
         GO TO 400
*     
 350     TMP1 = A(M)
         IF (TMP1 .LE. SIGMA1) THEN
            TMP1 = ZERO
         ENDIF
         DO J = M, N-1
            IF (TMP1 .EQ. ZERO .AND. B(J) .EQ. ZERO) THEN
               C1 = ZERO
               S1 = ONE
               WORK2(INDRV5+J) = ZERO
            ELSE
               CALL DLARTG(TMP1,B(J),C1,S1,WORK2(INDRV5+J))
            ENDIF
            WORK2(INDRV6+J) = S1*A(J+1)
            TMP1 = C1*A(J+1)
            IF (TMP1 .LE. SIGMA1) THEN
               TMP1 = ZERO
            ENDIF
         ENDDO
         WORK2(INDRV5+N) = TMP1

 360     TMP1 = WORK2(INDRV5+M)
         DO J = M, N-1
            CALL DLARTG(TMP1,WORK2(INDRV6+J),C1,S1,A(J))
            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, SU( 1, J), 1, SU( 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
      M = 1
      GO TO 3000
*     
 4000 INFO = 0
      RETURN
*     
      END SUBROUTINE DOQDSQ