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

      IMPLICIT NONE

      INTEGER N0, INFO, MAXITER, IINFO, METHOD
      REAL A(*), B(*), WORK(*)
      INTEGER N, M, I, J, K, OLDM, OLDN
      INTEGER INDRV1, INDRV2, INDRV3, INDRV4, INDRV5
      REAL TMP1, TMP2, TMP3, TMP4, TMP5
      REAL TAU, TAUA(4), TAUB
      REAL SIGMA, SIGMA1, SIGMA2, DESIG, T, DESIG0
      REAL C1, S1, EPS, TOL
      REAL DMIN2, DMIN1, SIGMX, SCALE, SAFMIN, RTMIN

      REAL ONE, ZERO, HALF, TWO, HUNDRD, CONST
      PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0, TWO = 2.0E0)
      PARAMETER (HUNDRD = 100.0E0, CONST=0.75E0)
*     
      INTRINSIC REAL
      EXTERNAL SLAMCH
      REAL SLAMCH
*
      IF (N0 .EQ. 1) THEN
         A(1)=ABS(A(1))
         INFO = 0
         RETURN
      ENDIF

      IF (N0 .EQ. 2) THEN
         CALL SLAS2(A(1),B(1),A(2),A(2),A(1))
         INFO = 0
         RETURN
      ENDIF
*
      TAUA(1) = ZERO
*
      A(1:N0)=ABS(A(1:N0))
      B(1:N0-1)=ABS(B(1:N0-1))
*
      EPS = SLAMCH( 'Precision' )
      TOL = HUNDRD*EPS
*     
      SAFMIN = SLAMCH( 'Safe minimum' )
      RTMIN = SQRT(SAFMIN)
      SCALE = SQRT( EPS / SAFMIN )

      SIGMX = MAX(MAXVAL(A(1:N0)),MAXVAL(B(1:N0-1)))
  
      IF (SIGMX .GT. ZERO) THEN
         CALL SLASCL( 'G', 0, 0, SIGMX, SCALE, N0, 1, 
     $        A(1), N0, IINFO )
         CALL SLASCL( 'G', 0, 0, SIGMX, SCALE, N0-1, 1, 
     $        B(1), N0-1, IINFO )
      ENDIF
*
      INDRV1 = 0
      INDRV2 = INDRV1+N0
      INDRV3 = INDRV2+N0
      INDRV4 = INDRV3+N0
      INDRV5 = INDRV4+N0
*     
      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

      OLDM = M
      OLDN = N
*     
 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 SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            
            GO TO 700
            
         ENDIF
         
         IF (N-M+1 .EQ. 2) THEN
            
            CALL SLAS2(A(N-1), B(N-1), A(N), A(N), A(N-1))
            
            CALL SLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0,RTMIN)
            CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            
            GO TO 700
         ENDIF
*
         IF (B(N-2) .LE. SIGMA2) THEN
            
            CALL SLAS2(A(N-1), B(N-1), A(N), A(N), A(N-1))
            
            CALL SLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0,RTMIN)
            CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            
            N = N - 2
            GO TO 15
            
         ENDIF
         
         IF (B(N-1) .LE. SIGMA2) THEN
            
            CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            
            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 SLAS2(A(N-1), B(N-1), A(N), TAU, TMP3)
         TAU = MIN(TAU,A(N))
         IF (TAU .LE. SIGMA1) GO TO 350

         IF (OLDM .EQ. M .AND. OLDN .EQ. N) THEN
            TAUB = DMIN1
         ELSE IF (OLDM .EQ. M .AND. OLDN-1 .EQ. N) THEN
            TAUB = DMIN2
         ELSE
            TAUB = MINVAL(A(M:N-1))
         ENDIF
         IF (TAUB .LE. SIGMA1) GO TO 350

         METHOD = 1
         IF (TAUB .LE. TAU) GO TO 140
         
         TAUB = TAU
*
         TMP4 = A(M)-TAU
         TMP5 = A(M)+TAU
         IF (TMP4 .LE. ZERO) THEN
            GO TO 140
         ELSE
            TMP3 = SQRT(TMP4*TMP5)
         ENDIF
         DO J = M, N-2
            CALL SLARTG2(TMP3,B(J),C1,S1,TMP1,RTMIN)
            TMP4 = C1*A(J+1)-TAU
            TMP5 = C1*A(J+1)+TAU
            IF (TMP4 .LE. ZERO) THEN
               GO TO 140
            ELSE
               TMP3 = SQRT(TMP4*TMP5)
            ENDIF
         ENDDO
         CALL SLARTG2(TMP3,B(N-1),C1,S1,TMP1,RTMIN)
         TMP4 = 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
               TAUA(METHOD+1) = MAX(TMP2,TMP3)
            ELSE
               GO TO 140
            ENDIF
         ELSE
            TAUA(METHOD+1) = TAU
         ENDIF
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
         IF (METHOD .GE. 2) GO TO 124

 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
         IF (TMP1 .GT. ZERO) THEN
            TAUA(METHOD+1) = SQRT(TMP3/TMP1)
            IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
               METHOD = METHOD+1
            ENDIF
         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
            TAUA(METHOD+1) = SQRT(TMP3/TMP1)
            IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
               METHOD = METHOD+1
            ENDIF
         ENDIF

         TAUA(METHOD) = MIN(TAUA(METHOD),TAUB)
         IF (METHOD .GE. 2) GO TO 124
         
         TAUA(METHOD+1)=A(N)-HALF*B(N-1)
         IF (TAUA(METHOD+1) .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
            TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
         ENDDO
         TMP2=A(M)-HALF*B(M)
         IF (TMP2 .LE. ZERO) GO TO 350
         TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
*
 124     TAU = TAUA(METHOD)
 125     IF (TAU .LE. ZERO) GO TO 350
*     
         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*TMP5)
         ENDIF
         DO J = M, N-2
            CALL SLARTG2(TMP3,B(J),C1,S1,WORK(INDRV1+J),RTMIN)
            WORK(INDRV2+J) = S1*A(J+1)
            TMP4 = C1*A(J+1)-TAU
            TMP5 = 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*TMP5)
            ENDIF
         ENDDO
         CALL SLARTG2(TMP3,B(N-1),C1,S1,WORK(INDRV1+N-1),RTMIN)
         WORK(INDRV2+N-1) = S1*A(N)
         TMP4 = C1*A(N)-TAU
         TMP5 = 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*TMP5)
         ENDIF

         IF (TMP3 .GE. A(N)) THEN
            TMP2 = TAU/A(N)
            IF (C1 .GE. TMP2) THEN
               TMP3 = A(N)*SQRT((C1-TMP2)*(C1+TMP2))
            ENDIF
         ENDIF
         WORK(INDRV1+N) = TMP3
*     
         CALL SLARTG7(SIGMA,DESIG,TAU,T,DESIG0,RTMIN)
*
         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 SLARTG2(TMP1,B(J),C1,S1,WORK(INDRV1+J),RTMIN)
            WORK(INDRV2+J) = S1*A(J+1)
            TMP1 = C1*A(J+1)
            IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         ENDDO
         WORK(INDRV1+N) = TMP1
*
 400     OLDM = M
         OLDN = N

         TMP1 = WORK(INDRV1+M)
         DMIN2 = TMP1
         DO J = M, N-3
            CALL SLARTG2(TMP1,WORK(INDRV2+J),C1,S1,A(J),RTMIN)
            B(J) = S1*WORK(INDRV1+J+1)
            TMP1 = C1*WORK(INDRV1+J+1)
            DMIN2 = MIN(DMIN2,TMP1)
         ENDDO

         CALL SLARTG2(TMP1,WORK(INDRV2+N-2),C1,S1,A(N-2),RTMIN)
         B(N-2) = S1*WORK(INDRV1+N-1)
         TMP1 = C1*WORK(INDRV1+N-1)
         IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         DMIN1 = MIN(DMIN2,TMP1)

         CALL SLARTG2(TMP1,WORK(INDRV2+N-1),C1,S1,A(N-1),RTMIN)
         B(N-1) = S1*WORK(INDRV1+N)
         TMP1 = C1*WORK(INDRV1+N)
         IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         
         A(N) = TMP1

         IF (DMIN2 .GT. SIGMA1) GO TO 450

         TMP1 = WORK(INDRV1+M)
         IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         DMIN2 = TMP1
         DO J = M, N-3
            CALL SLARTG2(TMP1,WORK(INDRV2+J),C1,S1,A(J),RTMIN)
            B(J) = S1*WORK(INDRV1+J+1)
            TMP1 = C1*WORK(INDRV1+J+1)
            IF (TMP1 .LE. SIGMA1) TMP1=ZERO
            DMIN2 = MIN(DMIN2,TMP1)
         ENDDO

         CALL SLARTG2(TMP1,WORK(INDRV2+N-2),C1,S1,A(N-2),RTMIN)
         B(N-2) = S1*WORK(INDRV1+N-1)
         TMP1 = C1*WORK(INDRV1+N-1)
         IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         DMIN1 = MIN(DMIN2,TMP1)

         CALL SLARTG2(TMP1,WORK(INDRV2+N-1),C1,S1,A(N-1),RTMIN)
         B(N-1) = S1*WORK(INDRV1+N)
         TMP1 = C1*WORK(INDRV1+N)
         IF (TMP1 .LE. SIGMA1) TMP1=ZERO
         
         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 SLASRT( 'D', N0, A, IINFO )
*
      IF (SIGMX .GT. ZERO) THEN
         CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N0, 1, A, N0, IINFO )
      ENDIF
      INFO = 0
      RETURN
*     
      END