SUBROUTINE DOQDSV(N0,A,B,WORK,INFO)

      IMPLICIT NONE

      INTEGER N0, INFO, MAXITER, IINFO
      DOUBLE PRECISION A(*), B(*), WORK(*)
      INTEGER N, M, I, J, K
      INTEGER INDRV1, INDRV2, INDRV3, INDRV4, INDRV5
      DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, TMP5
      DOUBLE PRECISION TAU, TAU1, TAU2
      DOUBLE PRECISION SIGMA, SIGMA1, SIGMA2, DESIG, T, DESIG0
      DOUBLE PRECISION C1, S1, EPS, TOL

      DOUBLE PRECISION ONE, ZERO, HALF, TWO, HUNDRD, CONST
      PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0, TWO = 2.0D0)
      PARAMETER (HUNDRD = 100.0D0, CONST=0.75D0)
*     
      EXTERNAL DLAMCH
      DOUBLE PRECISION DLAMCH
*
      IF (N0 .EQ. 1) THEN
         A(1)=ABS(A(1))
         INFO = 0
         RETURN
      ENDIF

      IF (N0 .EQ. 2) THEN
         CALL DLAS2(A(1),B(1),A(2),A(2),A(1))
         INFO = 0
         RETURN
      ENDIF
*      
      A(1:N0)=ABS(A(1:N0))
      B(1:N0-1)=ABS(B(1:N0-1))
*     
      INDRV1 = 0
      INDRV2 = INDRV1+N0
      INDRV3 = INDRV2+N0
      INDRV4 = INDRV3+N0
      INDRV5 = INDRV4+N0
*     
      EPS = DLAMCH( 'Precision' )
      TOL = HUNDRD*EPS
*     
      M = 1
      N = N0
*     
      IF (A(M) .LT. A(N)) THEN
         K=(N-M+1)/2+M-1
         DO J=M,K
            TMP1=A(J)
            A(J)=A(N+M-J)
            A(N+M-J)=TMP1
            TMP1=B(J)
            B(J)=B(N-1+M-J)
            B(N-1+M-J)=TMP1
         ENDDO
      END IF
*
      TMP1 = A(M)
      DO J = M, N-3
         IF (B(J) .LE. TOL*TMP1) THEN
            B(J) = -ZERO
            WORK(INDRV2+J) = -ZERO
            M = J+1
            TMP1 = A(J+1)
         ELSE
            TMP1 = A(J+1)*(TMP1/(TMP1+B(J)))
         ENDIF
      ENDDO
*
      IF (B(N-2) .LE. TOL*TMP1) THEN
         B(N-2) = ZERO
         TMP1 = A(N-1)
      ELSE
         TMP1 = A(N-1)*(TMP1/(TMP1+B(N-2)))
      ENDIF
*
      IF (B(N-1) .LE. TOL*TMP1) THEN
         B(N-1) = ZERO
      ENDIF
*     
      B(N) = ZERO
      WORK(INDRV2+N) = ZERO
*     
 3000 SIGMA = -B(N)
      SIGMA1 = SQRT(EPS)*SIGMA
      SIGMA2 = TOL*SIGMA
      DESIG = -WORK(INDRV2+N)

      MAXITER = 30*(N-M+1)
      DO I = 1, MAXITER
*     
 15      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
            
            CALL DLAS2(A(N-1), B(N-1), A(N), A(N), A(N-1))
            
            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 (B(N-2) .LE. SIGMA2) THEN
            
            CALL DLAS2(A(N-1), B(N-1), A(N), A(N), A(N-1))
            
            CALL DLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0)
            CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0)
            
            N = N - 2
            GO TO 15
            
         ENDIF
         
         IF (B(N-1) .LE. SIGMA2) THEN
            
            CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0)
            
            N = N - 1
            GO TO 15
            
         ENDIF
*     
         IF (A(M) .LT. A(N)) THEN
            K=(N-M+1)/2+M-1
            DO J=M,K
               TMP1=A(J)
               A(J)=A(N+M-J)
               A(N+M-J)=TMP1
               TMP1=B(J)
               B(J)=B(N-1+M-J)
               B(N-1+M-J)=TMP1
            ENDDO
         END IF
*
         CALL DLAS2(A(N-1), B(N-1), A(N), TAU, TMP3)
         TAU = MIN(TAU,A(N))
         IF (TAU .LE. SIGMA1) GO TO 350

         TAU2 = MINVAL(A(M:N-1))
         IF (TAU2 .LE. SIGMA1) GO TO 350
         
         IF (TAU2 .LE. TAU) THEN
            TAU1 = TAU2
            GO TO 140
         ELSE
            TAU1 = TAU
         ENDIF
*
         TMP4 = A(M)-TAU
         TMP5 = A(M)+TAU
         IF (TMP4 .LE. ZERO) THEN
            GO TO 140
         ELSE
            TMP3 = SQRT(TMP4)*SQRT(TMP5)
         ENDIF
         DO J = M, N-2
            C1 = ONE/SQRT(ONE+(B(J)/TMP3)**2)
            TMP4 = C1*A(J+1)-TAU
            TMP5 = C1*A(J+1)+TAU
            IF (TMP4 .LE. ZERO) THEN
               GO TO 140
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
         ENDDO
         TMP4 = A(N)/SQRT(ONE+(B(N-1)/TMP3)**2)-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) GO TO 125
         
 140     WORK(INDRV5+M:INDRV5+N)=ONE/A(M:N)
         WORK(INDRV1+N) = WORK(INDRV5+N)
         TMP3 = WORK(INDRV1+N)
         DO J = N-1,M,-1
            WORK(INDRV3+J) = B(J)*WORK(INDRV5+J)
            WORK(INDRV1+J) = WORK(INDRV5+J)+
     $           WORK(INDRV3+J)*WORK(INDRV1+J+1)
            TMP3 = MAX(TMP3,WORK(INDRV1+J))
         ENDDO
         TMP3 = ONE/TMP3
         WORK(INDRV1+M) = (WORK(INDRV1+M)*TMP3)*WORK(INDRV5+M)
         TMP1 = WORK(INDRV1+M)
         DO J = M+1,N
            WORK(INDRV4+J) = B(J-1)*WORK(INDRV5+J)
            WORK(INDRV1+J) = (WORK(INDRV1+J)*TMP3)*WORK(INDRV5+J)+
     $           WORK(INDRV4+J)*WORK(INDRV1+J-1)
            TMP1 = MAX(TMP1,WORK(INDRV1+J))
         ENDDO
         TAU = ZERO
         IF (TMP1 .GT. ZERO) THEN
            TAU = MAX(TAU,SQRT(TMP3)/SQRT(TMP1))
         ENDIF
         
         TMP1 = ONE/TMP1
         WORK(INDRV1+N) = WORK(INDRV1+N)*TMP1
         WORK(INDRV2+N) = WORK(INDRV1+N)*WORK(INDRV5+N)
         TMP3 = WORK(INDRV2+N)
         DO J = N-1,M,-1
            WORK(INDRV1+J) = WORK(INDRV1+J)*TMP1
            WORK(INDRV2+J) = WORK(INDRV1+J)*WORK(INDRV5+J)+
     $           WORK(INDRV3+J)*WORK(INDRV2+J+1)
            TMP3 = MAX(TMP3,WORK(INDRV2+J))
         ENDDO
         TMP3 = ONE/TMP3
         TMP2 = (WORK(INDRV2+M)*TMP3)*WORK(INDRV5+M)
         TMP1 = TMP2/WORK(INDRV1+M)
         DO J = M+1,N
            TMP2 = (WORK(INDRV2+J)*TMP3)*WORK(INDRV5+J)
     $           +WORK(INDRV4+J)*TMP2
            TMP1 = MAX(TMP1,TMP2/WORK(INDRV1+J))
         ENDDO
         IF (TMP1 .GT. ZERO) THEN
            TAU = MAX(TAU,SQRT(TMP3)/SQRT(TMP1))
         ENDIF

         TAU = MIN(TAU,TAU1)
         IF (TAU .GT. ZERO) GO TO 125
         
         TAU=A(N)-HALF*B(N-1)
         IF (TAU .LE. ZERO) GO TO 350
         DO J=N-1,M+1,-1
            TMP2=A(J)-HALF*(B(J)+B(J-1))
            IF (TMP2 .LE. ZERO) GO TO 350
            TAU=MIN(TAU,TMP2)
         ENDDO
         TMP2=A(M)-HALF*B(M)
         IF (TMP2 .LE. ZERO) GO TO 350
         TAU=MIN(TAU,TMP2)
*     
 125     IF (TAU .LE. ZERO) GO TO 350
*     
         TMP4 = A(M)-TAU
         TMP5 = A(M)+TAU
         IF (TMP4 .LE. 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
         DO J = M, N-2
            CALL DLARTG(TMP3,B(J),C1,S1,WORK(INDRV1+J))
            WORK(INDRV2+J) = S1*A(J+1)
            TMP4 = C1*A(J+1)-TAU
            TMP5 = C1*A(J+1)+TAU
            IF (TMP4 .LE. 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
         ENDDO
         CALL DLARTG(TMP3,B(N-1),C1,S1,WORK(INDRV1+N-1))
         WORK(INDRV2+N-1) = S1*A(N)
         TMP4 = C1*A(N)-TAU
         TMP5 = 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
         WORK(INDRV1+N) = TMP3
*     
         CALL DLARTG7(SIGMA,DESIG,TAU,T,DESIG0)
*
         SIGMA = T
         SIGMA1 = SQRT(EPS)*SIGMA
         SIGMA2 = TOL*SIGMA
         DESIG = DESIG0
*     
         GO TO 400
*     
 350     TMP1 = A(M)
         IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         DO J = M, N-1
            CALL DLARTG(TMP1,B(J),C1,S1,WORK(INDRV1+J))
            WORK(INDRV2+J) = S1*A(J+1)
            TMP1 = C1*A(J+1)
            IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         ENDDO
         WORK(INDRV1+N) = TMP1
*     
 400     TMP1 = WORK(INDRV1+M)
         TMP2 = TMP1
         DO J = M, N-1
            CALL DLARTG(TMP1,WORK(INDRV2+J),C1,S1,A(J))
            B(J) = S1*WORK(INDRV1+J+1)
            TMP1 = C1*WORK(INDRV1+J+1)
            TMP2 = MIN(TMP2,TMP1)
         ENDDO
         A(N) = TMP1

         IF (TMP2 .GT. SIGMA1) GO TO 450

         TMP1 = WORK(INDRV1+M)
         IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         DO J = M, N-1
            CALL DLARTG(TMP1,WORK(INDRV2+J),C1,S1,A(J))
            B(J) = S1*WORK(INDRV1+J+1)
            TMP1 = C1*WORK(INDRV1+J+1)
            IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         ENDDO
         A(N) = TMP1
*  
 450     DO J = M, N-3
            IF (B(J) .LE. SIGMA2 .OR.
     $           WORK(INDRV2+J) .LE. TOL*A(J)) THEN
               B(J) = -SIGMA
               WORK(INDRV2+J) = -DESIG
               M = J+1
            ENDIF
         ENDDO
         IF (WORK(INDRV2+N-2) .LE. TOL*A(N-2)) THEN
            B(N-2) = ZERO
         ENDIF
         IF (WORK(INDRV2+N-1) .LE. TOL*A(N-1)) THEN
            B(N-1) = ZERO
         ENDIF
*     
      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 CALL DLASRT( 'D', N0, A, IINFO )
*
      INFO = 0
      RETURN
*     
      END