SUBROUTINE SDQDS(N0,A,B,WORK,INFO)
  
      IMPLICIT NONE
  
      REAL A(*), B(*), WORK(*)
      INTEGER N0,INFO,IINFO,MAXITER
      INTEGER N, M, I, J, K
      INTEGER INDRV1,INDRV2,INDRV3,INDRV4,INDRV5,INDRV6,INDRV7,INDRV8
      REAL TMP1, TMP2, TMP3, TAU, TAU1, TAU2, TAU3, S, T
      REAL TMP4, TMP5, TMP6, TMP7, TMP8, TMP9
      REAL EPS,SAFMIN,SCALE,SIGMX,TOL,TOL2
      REAL SIGMA,SIGMA1,DESIG,DESIG0
      REAL SCALE1
      REAL DMIN2
      REAL DN0,DN1,DN2

      REAL ONE, ZERO, HALF, TWO, HUNDRD, CONST
      PARAMETER ( ONE=1.0E0, ZERO=0.0E0, HALF=0.5E0 )
      PARAMETER ( TWO=2.0E0, HUNDRD=100.0E0, CONST=0.75E0 )

      INTRINSIC REAL
      EXTERNAL SLAMCH,SDQDS2
      REAL SLAMCH
      EXTERNAL SISNAN
      LOGICAL SISNAN

      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

      INDRV1 = 0
      INDRV2 = INDRV1+N0
      INDRV3 = INDRV2+N0
      INDRV4 = INDRV3+N0
      INDRV5 = INDRV4+N0
      INDRV6 = INDRV5+N0
      INDRV7 = INDRV6+N0
      INDRV8 = INDRV7+N0
  
      EPS = SLAMCH( 'Precision' )
      TOL = HUNDRD*EPS
      TOL2 = TOL**2
      SAFMIN = SLAMCH( 'Safe minimum' )
      SCALE = SQRT( HALF / SAFMIN )
      SCALE1 = ONE/SCALE

      WORK(INDRV1+1:INDRV1+N0)=ABS(A(1:N0))
      WORK(INDRV2+1:INDRV2+N0-1)=ABS(B(1:N0-1))

      SIGMX = MAX(MAXVAL(WORK(INDRV1+1:INDRV1+N0)),
     $     MAXVAL(WORK(INDRV2+1:INDRV2+N0-1)))

      IF (SIGMX .GT. ZERO) THEN
         CALL SLASCL( 'G', 0, 0, SIGMX, SCALE, N0, 1, 
     $        WORK(INDRV1+1), N0, IINFO )
         CALL SLASCL( 'G', 0, 0, SIGMX, SCALE, N0-1, 1, 
     $        WORK(INDRV2+1), N0-1, IINFO )
      ENDIF

      WORK(INDRV1+1:INDRV1+N0) = WORK(INDRV1+1:INDRV1+N0)**2
      WORK(INDRV2+1:INDRV2+N0-1) = WORK(INDRV2+1:INDRV2+N0-1)**2
  
      M = 1
      N = N0

      IF (WORK(INDRV1+M) .LT. WORK(INDRV1+N)) THEN
         K=(N-M+1)/2+M-1
         DO J=M,K
            TMP1=WORK(INDRV1+J)
            WORK(INDRV1+J)=WORK(INDRV1+N+M-J)
            WORK(INDRV1+N+M-J)=TMP1
            TMP1=WORK(INDRV2+J)
            WORK(INDRV2+J)=WORK(INDRV2+N-1+M-J)
            WORK(INDRV2+N-1+M-J)=TMP1
         ENDDO
      END IF

      TMP1 = WORK(INDRV1+M)
      DO J = M, N-3
         IF (WORK(INDRV2+J) .LE. TOL2*TMP1) THEN
            WORK(INDRV2+J) = -ZERO
            WORK(INDRV4+J) = -ZERO
            M = J+1
            TMP1 = WORK(INDRV1+J+1)
         ELSE
            TMP1 = WORK(INDRV1+J+1)*(TMP1/(TMP1+WORK(INDRV2+J)))
         ENDIF
      ENDDO
*     
      IF (WORK(INDRV2+N-2) .LE. TOL2*TMP1) THEN
         WORK(INDRV2+N-2) = ZERO
         TMP1 = WORK(INDRV1+N-1)
      ELSE
         TMP1 = WORK(INDRV1+N-1)*(TMP1/(TMP1+WORK(INDRV2+N-2)))
      ENDIF
*     
      IF (WORK(INDRV2+N-1) .LE. TOL2*TMP1) THEN
         WORK(INDRV2+N-1) = ZERO
      ENDIF
*     
      WORK(INDRV2+N) = ZERO
      WORK(INDRV4+N) = ZERO

 3000 SIGMA = -WORK(INDRV2+N)
      SIGMA1 = TOL2*SIGMA
      DESIG = -WORK(INDRV4+N)

      MAXITER=30*(N-M+1)
      DO I=1, MAXITER, 2

 15      IF (N-M+1 .EQ. 2) THEN
            
            IF( WORK(INDRV1+N).GT.WORK(INDRV1+N-1) ) THEN
               S = WORK(INDRV1+N)
               WORK(INDRV1+N) = WORK(INDRV1+N-1)
               WORK(INDRV1+N-1) = S
            END IF
            T = HALF*( ( WORK(INDRV1+N-1)-WORK(INDRV1+N) )
     $           +WORK(INDRV2+N-1) )
            IF(WORK(INDRV2+N-1) .GT. WORK(INDRV1+N)*TOL2 
     $           .AND. T .NE. ZERO) THEN
               S = WORK(INDRV1+N)*( WORK(INDRV2+N-1) / T )
               IF( S.LE.T ) THEN
                  S = WORK(INDRV1+N)*( WORK(INDRV2+N-1) / 
     $                 ( T*( ONE+SQRT( ONE+S / T ) ) ) )
               ELSE
                  S = WORK(INDRV1+N)*( WORK(INDRV2+N-1) / 
     $                 ( T+SQRT( T )*SQRT( T+S ) ) )
               END IF
               T = WORK(INDRV1+N-1) + ( S+WORK(INDRV2+N-1) )
               WORK(INDRV1+N) = WORK(INDRV1+N)
     $              *( WORK(INDRV1+N-1) / T )
               WORK(INDRV1+N-1) = T
            END IF
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV1+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV1+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            GO TO 700
         ENDIF

         IF (N-M+1 .EQ. 1) THEN
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV1+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            GO TO 700
         ENDIF

         IF (WORK(INDRV2+N-2) .LE. SIGMA1) THEN
            
            IF( WORK(INDRV1+N).GT.WORK(INDRV1+N-1) ) THEN
               S = WORK(INDRV1+N)
               WORK(INDRV1+N) = WORK(INDRV1+N-1)
               WORK(INDRV1+N-1) = S
            END IF
            T = HALF*( ( WORK(INDRV1+N-1)-WORK(INDRV1+N) )
     $           +WORK(INDRV2+N-1) )
            IF(WORK(INDRV2+N-1) .GT. WORK(INDRV1+N)*TOL2 
     $           .AND. T .NE. ZERO) THEN
               S = WORK(INDRV1+N)*( WORK(INDRV2+N-1) / T )
               IF( S.LE.T ) THEN
                  S = WORK(INDRV1+N)*( WORK(INDRV2+N-1) / 
     $                 ( T*( ONE+SQRT( ONE+S / T ) ) ) )
               ELSE
                  S = WORK(INDRV1+N)*( WORK(INDRV2+N-1) / 
     $                 ( T+SQRT( T )*SQRT( T+S ) ) )
               END IF
               T = WORK(INDRV1+N-1) + ( S+WORK(INDRV2+N-1) )
               WORK(INDRV1+N) = WORK(INDRV1+N)
     $              *( WORK(INDRV1+N-1) / T )
               WORK(INDRV1+N-1) = T
            END IF
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV1+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV1+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            N = N - 2
            GO TO 15
         ENDIF
 
         IF (WORK(INDRV2+N-1) .LE. SIGMA1) THEN
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV1+N),T,DESIG0)
            WORK(INDRV1+N)=T

            N = N - 1
            GO TO 15
         ENDIF
     
         IF (WORK(INDRV1+M) .LT. WORK(INDRV1+N)) THEN
            K=(N-M+1)/2+M-1
            DO J=M,K
               TMP1=WORK(INDRV1+J)
               WORK(INDRV1+J)=WORK(INDRV1+N+M-J)
               WORK(INDRV1+N+M-J)=TMP1
               TMP1=WORK(INDRV2+J)
               WORK(INDRV2+J)=WORK(INDRV2+N-1+M-J)
               WORK(INDRV2+N-1+M-J)=TMP1
            ENDDO
         END IF

         TAU = WORK(INDRV1+N)
         TMP3 = WORK(INDRV1+N-1)
         IF( TAU.GT.TMP3 ) THEN
            S = TAU
            TAU = TMP3
            TMP3 = S
         END IF
         T = HALF*( ( TMP3-TAU )
     $        +WORK(INDRV2+N-1) )
         IF(WORK(INDRV2+N-1) .GT. TAU*TOL2 
     $        .AND. T .NE. ZERO) THEN
            S = TAU*( WORK(INDRV2+N-1) / T )
            IF( S.LE.T ) THEN
               S = TAU*( WORK(INDRV2+N-1) / 
     $              ( T*( ONE+SQRT( ONE+S / T ) ) ) )
            ELSE
               S = TAU*( WORK(INDRV2+N-1) / 
     $              ( T+SQRT( T )*SQRT( T+S ) ) )
            END IF
            T = TMP3 + ( S+WORK(INDRV2+N-1) )
            TAU = TAU
     $           *( TMP3 / T )
            TMP3 = T
         END IF
         TAU = MIN(TAU,WORK(INDRV1+N))
         IF (SIGMA .EQ. TAU+SIGMA) GO TO 350

         TAU2 = MINVAL(WORK(INDRV1+M:INDRV1+N-1))
         IF (SIGMA .EQ. TAU2+SIGMA) GO TO 350
         IF (TAU2 .LE. TAU) THEN
            TAU1 = TAU2
            GO TO 140
         ELSE
            TAU1 = TAU
         ENDIF

         DN2 = WORK(INDRV1+M)-TAU
         IF (DN2 .LE. ZERO) GO TO 140
         DO J=M, N-3
            DN2 = DN2*(WORK(INDRV1+J+1)/(DN2+WORK(INDRV2+J)))-TAU
            IF (DN2 .LE. ZERO) GO TO 140
         ENDDO

         TMP1 = DN2+WORK(INDRV2+N-2)
         IF (SAFMIN*TMP1 .LT. WORK(INDRV1+N-1) .AND.  
     $        SAFMIN*WORK(INDRV1+N-1) .LT. TMP1) THEN
            DN1 = DN2*(WORK(INDRV1+N-1)/TMP1)-TAU
         ELSE
            DN1 = WORK(INDRV1+N-1)*(DN2/TMP1)-TAU
         ENDIF
         IF(DN1 .LE. ZERO) GO TO 140

         TMP1 = DN1+WORK(INDRV2+N-1)
         IF (SAFMIN*TMP1 .LT. WORK(INDRV1+N) .AND.  
     $        SAFMIN*WORK(INDRV1+N) .LT. TMP1) THEN
            DN0 = DN1*(WORK(INDRV1+N)/TMP1)-TAU
         ELSE
            DN0 = WORK(INDRV1+N)*(DN1/TMP1)-TAU
         ENDIF
         IF(DN0 .LT. ZERO .AND. TAU+DN0 .EQ. TAU) DN0=ZERO
         IF(DN0 .LT. ZERO) THEN
            TAU2 = MAX(DN0+TAU,ZERO)
            IF (TAU2 .EQ. TAU) TAU2 = ZERO
            TAU3 = CONST*TAU
            IF (TAU3 .EQ. TAU) TAU3 = HALF*TAU
            TAU = MAX(TAU2,TAU3)
         ENDIF

         IF (TAU .GT. ZERO) GO TO 125

 140     TMP7=SCALE/WORK(INDRV1+M)
         TMP1=TMP7
         TMP5=TMP7**2
         TMP3=TMP5
         TMP2=TMP3
         DMIN2=TMP3
         DO J = M+1, N
            TMP4=SCALE/WORK(INDRV1+J)
            TMP8=(WORK(INDRV2+J-1)*SCALE1)*TMP4
            TMP7=TMP8*TMP7+TMP4
            TMP1=TMP1+TMP7
            TMP6=TMP7**2
            TMP3=TMP8*(TMP3+TMP5)+TMP6
            TMP2=TMP2+TMP3
            DMIN2=MAX(DMIN2,TMP3)
            TMP5=TMP6
         ENDDO

         TAU = MAX(ZERO,SCALE/TMP1,SCALE/SQRT(TMP2))
         TMP8 = REAL(N-M+1)
         TMP9 = SCALE*(TMP8/
     $        (TMP1+SQRT(TMP8-ONE)*SQRT(TMP8*TMP2-TMP1**2)))
         IF (TMP9 .EQ. TMP9) THEN
            TAU = MAX(TAU,TMP9)
         ENDIF

         TAU2 = TAU1
         TMP8 = CEILING((TMP1/TMP2)*TMP1)
         IF (TMP8 .GT. ONE) THEN
            TMP9 = SCALE*(TMP8/
     $           (TMP1+SQRT(TMP8*TMP2-TMP1**2)/SQRT(TMP8-ONE)))
            IF (TMP9 .EQ. TMP9) THEN
               TAU2 = MIN(TMP9,SCALE/SQRT(DMIN2),TAU2)
            ELSE
               TAU2 = MIN(SCALE/SQRT(DMIN2),TAU2)
            ENDIF
         ELSE
            TAU2 = MIN(SCALE/SQRT(DMIN2),TAU2)
         ENDIF

         TAU = MIN(TAU,TAU2)
         IF (TAU2 .LT. TWO*TAU) GO TO 125

 150     WORK(INDRV7+M:INDRV7+N) = ONE/SQRT(WORK(INDRV1+M:INDRV1+N))
         WORK(INDRV8+M:INDRV8+N-1) = SQRT(WORK(INDRV2+M:INDRV2+N-1))
         
         WORK(INDRV5+N) = WORK(INDRV7+N)
         TMP3 = WORK(INDRV5+N)
         DO J = N-1,M,-1
            WORK(INDRV3+J) = WORK(INDRV8+J)*WORK(INDRV7+J)
            WORK(INDRV5+J) = WORK(INDRV7+J)+
     $           WORK(INDRV3+J)*WORK(INDRV5+J+1)
            TMP3 = MAX(TMP3,WORK(INDRV5+J))
         ENDDO
         TMP3 = ONE/TMP3
         WORK(INDRV5+M) = (WORK(INDRV5+M)*TMP3)*WORK(INDRV7+M)
         TMP1 = WORK(INDRV5+M)
         DO J = M+1,N
            WORK(INDRV4+J) = WORK(INDRV8+J-1)*WORK(INDRV7+J)
            WORK(INDRV5+J) = (WORK(INDRV5+J)*TMP3)*WORK(INDRV7+J)+
     $           WORK(INDRV4+J)*WORK(INDRV5+J-1)
            TMP1 = MAX(TMP1,WORK(INDRV5+J))
         ENDDO
         TAU = ZERO
         IF (TMP1 .GT. ZERO) THEN
            TAU = MAX(TAU,TMP3/TMP1)
         ENDIF
         
         TMP1 = ONE/TMP1
         WORK(INDRV5+N) = WORK(INDRV5+N)*TMP1
         WORK(INDRV6+N) = WORK(INDRV5+N)*WORK(INDRV7+N)
         TMP3 = WORK(INDRV6+N)
         DO J = N-1,M,-1
            WORK(INDRV5+J) = WORK(INDRV5+J)*TMP1
            WORK(INDRV6+J) = WORK(INDRV5+J)*WORK(INDRV7+J)+
     $           WORK(INDRV3+J)*WORK(INDRV6+J+1)
            TMP3 = MAX(TMP3,WORK(INDRV6+J))
         ENDDO
         TMP3 = ONE/TMP3
         TMP2 = (WORK(INDRV6+M)*TMP3)*WORK(INDRV7+M)
         TMP1 = TMP2/WORK(INDRV5+M)
         DO J = M+1,N
            TMP2 = (WORK(INDRV6+J)*TMP3)*WORK(INDRV7+J)+
     $           WORK(INDRV4+J)*TMP2
            TMP1 = MAX(TMP1,TMP2/WORK(INDRV5+J))
         ENDDO
         IF (TMP1 .GT. ZERO) THEN
            TAU = MAX(TAU,TMP3/TMP1)
         ENDIF

         TAU = MIN(TAU,TAU1)
         IF (TAU .GT. ZERO) GO TO 125

         TAU=SQRT(WORK(INDRV1+N))-HALF*WORK(INDRV8+N-1)
         IF (TAU .LE. ZERO) GO TO 350
         DO J=N-1,M+1,-1
            TMP2=SQRT(WORK(INDRV1+J))
     $           -HALF*(WORK(INDRV8+J)+WORK(INDRV8+J-1))
            IF (TMP2 .LE. ZERO) GO TO 350
            TAU=MIN(TAU,TMP2)
         ENDDO
         TMP2=SQRT(WORK(INDRV1+M))-HALF*WORK(INDRV8+M)
         IF (TMP2 .LE. ZERO) GO TO 350
         TAU=MIN(TAU,TMP2)**2

 125     IF (TAU .LE. ZERO) GO TO 350

         DN2 = WORK(INDRV1+M)-TAU
         DMIN2 = DN2
         DO J=M, N-3
            WORK(INDRV3+J) = DN2+WORK(INDRV2+J)
            TMP3 = WORK(INDRV1+J+1)/WORK(INDRV3+J)
            WORK(INDRV4+J) = WORK(INDRV2+J)*TMP3
            DN2 = DN2*TMP3-TAU
            DMIN2 = MIN(DMIN2,DN2)
         ENDDO
         IF (DMIN2 .LE. ZERO) THEN 
            TAU3 = CONST*TAU
            IF (TAU3 .EQ. TAU) TAU3 = HALF*TAU
            TAU = TAU3
            GO TO 125
         ENDIF

         WORK(INDRV3+N-2) = DN2+WORK(INDRV2+N-2)
         IF (SAFMIN*WORK(INDRV3+N-2) .LT. WORK(INDRV1+N-1) .AND.  
     $        SAFMIN*WORK(INDRV1+N-1) .LT. WORK(INDRV3+N-2)) THEN
            TMP3 = WORK(INDRV1+N-1)/WORK(INDRV3+N-2)
            WORK(INDRV4+N-2) = WORK(INDRV2+N-2)*TMP3
            DN1 = DN2*TMP3-TAU
         ELSE
            WORK(INDRV4+N-2) = WORK(INDRV1+N-1)*
     $           (WORK(INDRV2+N-2)/WORK(INDRV3+N-2))
            DN1 = WORK(INDRV1+N-1)*(DN2/WORK(INDRV3+N-2))-TAU
         ENDIF
         IF(DN1 .LE. ZERO) THEN
            TAU2 = MAX(DN1+TAU,ZERO)
            IF (TAU2 .EQ. TAU) TAU2 = ZERO
            TAU3 = CONST*TAU
            IF (TAU3 .EQ. TAU) TAU3 = HALF*TAU
            TAU = MAX(TAU2,TAU3)
            GO TO 125
         ENDIF

         WORK(INDRV3+N-1) = DN1+WORK(INDRV2+N-1)
         IF (SAFMIN*WORK(INDRV3+N-1) .LT. WORK(INDRV1+N) .AND.  
     $        SAFMIN*WORK(INDRV1+N) .LT. WORK(INDRV3+N-1)) THEN
            TMP3 = WORK(INDRV1+N)/WORK(INDRV3+N-1)
            WORK(INDRV4+N-1) = WORK(INDRV2+N-1)*TMP3
            DN0 = DN1*TMP3-TAU
         ELSE
            WORK(INDRV4+N-1) = WORK(INDRV1+N)
     $           *(WORK(INDRV2+N-1)/WORK(INDRV3+N-1))
            DN0 = WORK(INDRV1+N)*(DN1/WORK(INDRV3+N-1))-TAU
         ENDIF
         IF(DN0 .LT. ZERO .AND. TAU+DN0 .EQ. TAU) DN0=ZERO
         IF(DN0 .LT. ZERO) THEN
            TAU2 = MAX(DN0+TAU,ZERO)
            IF (TAU2 .EQ. TAU) TAU2 = ZERO
            TAU3 = CONST*TAU
            IF (TAU3 .EQ. TAU) TAU3 = HALF*TAU
            TAU = MAX(TAU2,TAU3)
            GO TO 125
         ENDIF
         IF (DN0 .GE. WORK(INDRV1+N)) THEN
            TMP3 = DN1/WORK(INDRV3+N-1)-TAU/WORK(INDRV1+N)
            IF (TMP3 .GE. ZERO) THEN
               DN0 = WORK(INDRV1+N)*TMP3
            ENDIF
         ENDIF
         
         WORK(INDRV3+N)=DN0

         IF (SISNAN(DN0)) THEN
            TAU3 = CONST*TAU
            IF (TAU3 .EQ. TAU) TAU3 = HALF*TAU
            TAU = TAU3
            GO TO 125
         ENDIF

         CALL SDQDS2(SIGMA,DESIG,TAU,T,DESIG0)
         SIGMA = T
         SIGMA1 = TOL2*SIGMA
         DESIG = DESIG0

         GO TO 400
   
 350     DN2 = WORK(INDRV1+M)
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         DO J=M, N-3
            WORK(INDRV3+J) = DN2+WORK(INDRV2+J)
            TMP3 = WORK(INDRV1+J+1)/WORK(INDRV3+J)
            WORK(INDRV4+J) = WORK(INDRV2+J)*TMP3
            DN2 = DN2*TMP3
            IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         ENDDO
         
         WORK(INDRV3+N-2) = DN2+WORK(INDRV2+N-2)
         IF (SAFMIN*WORK(INDRV3+N-2) .LT. WORK(INDRV1+N-1) .AND.  
     $        SAFMIN*WORK(INDRV1+N-1) .LT. WORK(INDRV3+N-2)) THEN
            TMP2=WORK(INDRV1+N-1)/WORK(INDRV3+N-2)
            WORK(INDRV4+N-2)=WORK(INDRV2+N-2)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV4+N-2) = WORK(INDRV1+N-1)
     $           *(WORK(INDRV2+N-2)/WORK(INDRV3+N-2))
            DN2 = WORK(INDRV1+N-1)*(DN2/WORK(INDRV3+N-2))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV3+N-1) = DN2+WORK(INDRV2+N-1)
         IF (SAFMIN*WORK(INDRV3+N-1) .LT. WORK(INDRV1+N) .AND.  
     $        SAFMIN*WORK(INDRV1+N) .LT. WORK(INDRV3+N-1)) THEN
            TMP2=WORK(INDRV1+N)/WORK(INDRV3+N-1)
            WORK(INDRV4+N-1)=WORK(INDRV2+N-1)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV4+N-1) = WORK(INDRV1+N)
     $           *(WORK(INDRV2+N-1)/WORK(INDRV3+N-1))
            DN2 = WORK(INDRV1+N)*(DN2/WORK(INDRV3+N-1))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV3+N)=DN2
         
         IF (SISNAN(DN2)) GO TO 370
         
         GO TO 400
      
 370     DN2 = WORK(INDRV1+M)
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         DO J=M, N-3
            WORK(INDRV3+J) = DN2+WORK(INDRV2+J)
            IF (SAFMIN*WORK(INDRV3+J) .LT. WORK(INDRV1+J+1) .AND.  
     $           SAFMIN*WORK(INDRV1+J+1) .LT. WORK(INDRV3+J)) THEN
               TMP2=WORK(INDRV1+J+1)/WORK(INDRV3+J)
               WORK(INDRV4+J)=WORK(INDRV2+J)*TMP2
               DN2=DN2*TMP2
            ELSE
               WORK(INDRV4+J) = WORK(INDRV1+J+1)
     $              *(WORK(INDRV2+J)/WORK(INDRV3+J))
               DN2 = WORK(INDRV1+J+1)*(DN2/WORK(INDRV3+J))
            ENDIF
            IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         ENDDO

         WORK(INDRV3+N-2) = DN2+WORK(INDRV2+N-2)
         IF (SAFMIN*WORK(INDRV3+N-2) .LT. WORK(INDRV1+N-1) .AND.  
     $        SAFMIN*WORK(INDRV1+N-1) .LT. WORK(INDRV3+N-2)) THEN
            TMP2=WORK(INDRV1+N-1)/WORK(INDRV3+N-2)
            WORK(INDRV4+N-2)=WORK(INDRV2+N-2)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV4+N-2) = WORK(INDRV1+N-1)
     $           *(WORK(INDRV2+N-2)/WORK(INDRV3+N-2))
            DN2 = WORK(INDRV1+N-1)*(DN2/WORK(INDRV3+N-2))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV3+N-1) = DN2+WORK(INDRV2+N-1)
         IF (SAFMIN*WORK(INDRV3+N-1) .LT. WORK(INDRV1+N) .AND.  
     $        SAFMIN*WORK(INDRV1+N) .LT. WORK(INDRV3+N-1)) THEN
            TMP2=WORK(INDRV1+N)/WORK(INDRV3+N-1)
            WORK(INDRV4+N-1)=WORK(INDRV2+N-1)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV4+N-1) = WORK(INDRV1+N)
     $           *(WORK(INDRV2+N-1)/WORK(INDRV3+N-1))
            DN2 = WORK(INDRV1+N)*(DN2/WORK(INDRV3+N-1))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV3+N)=DN2

         GO TO 400

 390     IF (N-M+1 .EQ. 2) THEN
            
            IF( WORK(INDRV3+N).GT.WORK(INDRV3+N-1) ) THEN
               S = WORK(INDRV3+N)
               WORK(INDRV3+N) = WORK(INDRV3+N-1)
               WORK(INDRV3+N-1) = S
            END IF
            T = HALF*( ( WORK(INDRV3+N-1)-WORK(INDRV3+N) )
     $           +WORK(INDRV4+N-1) )
            IF(WORK(INDRV4+N-1) .GT. WORK(INDRV3+N)*TOL2 
     $           .AND. T .NE. ZERO) THEN
               S = WORK(INDRV3+N)*( WORK(INDRV4+N-1) / T )
               IF( S.LE.T ) THEN
                  S = WORK(INDRV3+N)*( WORK(INDRV4+N-1) / 
     $                 ( T*( ONE+SQRT( ONE+S / T ) ) ) )
               ELSE
                  S = WORK(INDRV3+N)*( WORK(INDRV4+N-1) / 
     $                 ( T+SQRT( T )*SQRT( T+S ) ) )
               END IF
               T = WORK(INDRV3+N-1) + ( S+WORK(INDRV4+N-1) )
               WORK(INDRV3+N) = WORK(INDRV3+N)
     $              *( WORK(INDRV3+N-1) / T )
               WORK(INDRV3+N-1) = T
            END IF
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV3+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV3+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            GO TO 700
         ENDIF

         IF (N-M+1 .EQ. 1) THEN
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV3+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            GO TO 700
         ENDIF

 400     IF (WORK(INDRV4+N-2) .LE. SIGMA1 .OR.
     $        WORK(INDRV2+N-2) .LE. TOL2*WORK(INDRV3+N-2)) THEN
            
            IF( WORK(INDRV3+N).GT.WORK(INDRV3+N-1) ) THEN
               S = WORK(INDRV3+N)
               WORK(INDRV3+N) = WORK(INDRV3+N-1)
               WORK(INDRV3+N-1) = S
            END IF
            T = HALF*( ( WORK(INDRV3+N-1)-WORK(INDRV3+N) )
     $           +WORK(INDRV4+N-1) )
            IF(WORK(INDRV4+N-1) .GT. WORK(INDRV3+N)*TOL2 
     $           .AND. T .NE. ZERO) THEN
               S = WORK(INDRV3+N)*( WORK(INDRV4+N-1) / T )
               IF( S.LE.T ) THEN
                  S = WORK(INDRV3+N)*( WORK(INDRV4+N-1) / 
     $                 ( T*( ONE+SQRT( ONE+S / T ) ) ) )
               ELSE
                  S = WORK(INDRV3+N)*( WORK(INDRV4+N-1) / 
     $                 ( T+SQRT( T )*SQRT( T+S ) ) )
               END IF
               T = WORK(INDRV3+N-1) + ( S+WORK(INDRV4+N-1) )
               WORK(INDRV3+N) = WORK(INDRV3+N)
     $              *( WORK(INDRV3+N-1) / T )
               WORK(INDRV3+N-1) = T
            END IF
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV3+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV3+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            N = N - 2
            GO TO 390
         ENDIF
         
         IF (WORK(INDRV4+N-1) .LE. SIGMA1 .OR.
     $        WORK(INDRV2+N-1) .LE. TOL2*WORK(INDRV3+N-1)) THEN
            
            CALL SDQDS2(SIGMA,DESIG,WORK(INDRV3+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            N = N - 1
            GO TO 390
         ENDIF
     
         DN2 = WORK(INDRV3+M)
         DMIN2 = DN2
         DO J=M, N-3
            WORK(INDRV1+J) = DN2+WORK(INDRV4+J)
            TMP2=WORK(INDRV3+J+1)/WORK(INDRV1+J)
            WORK(INDRV2+J)=WORK(INDRV4+J)*TMP2
            DN2=DN2*TMP2
            DMIN2 = MIN(DMIN2,DN2)
         ENDDO

         WORK(INDRV1+N-2) = DN2+WORK(INDRV4+N-2)
         IF (SAFMIN*WORK(INDRV1+N-2) .LT. WORK(INDRV3+N-1) .AND.  
     $           SAFMIN*WORK(INDRV3+N-1) .LT. WORK(INDRV1+N-2)) THEN
            TMP2=WORK(INDRV3+N-1)/WORK(INDRV1+N-2)
            WORK(INDRV2+N-2)=WORK(INDRV4+N-2)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV2+N-2) = WORK(INDRV3+N-1)
     $           *(WORK(INDRV4+N-2)/WORK(INDRV1+N-2))
            DN2 = WORK(INDRV3+N-1)*(DN2/WORK(INDRV1+N-2))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV1+N-1) = DN2+WORK(INDRV4+N-1)
         IF (SAFMIN*WORK(INDRV1+N-1) .LT. WORK(INDRV3+N) .AND.  
     $           SAFMIN*WORK(INDRV3+N) .LT. WORK(INDRV1+N-1)) THEN
            TMP2=WORK(INDRV3+N)/WORK(INDRV1+N-1)
            WORK(INDRV2+N-1)=WORK(INDRV4+N-1)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV2+N-1) = WORK(INDRV3+N)
     $           *(WORK(INDRV4+N-1)/WORK(INDRV1+N-1))
            DN2 = WORK(INDRV3+N)*(DN2/WORK(INDRV1+N-1))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV1+N)=DN2

         IF (SISNAN(DN2)) GO TO 1410

         IF (DMIN2+SIGMA .NE. SIGMA) GO TO 1420

         DN2 = WORK(INDRV3+M)
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         DO J=M, N-3
            WORK(INDRV1+J) = DN2+WORK(INDRV4+J)
            TMP2=WORK(INDRV3+J+1)/WORK(INDRV1+J)
            WORK(INDRV2+J)=WORK(INDRV4+J)*TMP2
            DN2=DN2*TMP2
            IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         ENDDO

         WORK(INDRV1+N-2) = DN2+WORK(INDRV4+N-2)
         IF (SAFMIN*WORK(INDRV1+N-2) .LT. WORK(INDRV3+N-1) .AND.  
     $           SAFMIN*WORK(INDRV3+N-1) .LT. WORK(INDRV1+N-2)) THEN
            TMP2=WORK(INDRV3+N-1)/WORK(INDRV1+N-2)
            WORK(INDRV2+N-2)=WORK(INDRV4+N-2)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV2+N-2) = WORK(INDRV3+N-1)
     $           *(WORK(INDRV4+N-2)/WORK(INDRV1+N-2))
            DN2 = WORK(INDRV3+N-1)*(DN2/WORK(INDRV1+N-2))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV1+N-1) = DN2+WORK(INDRV4+N-1)
         IF (SAFMIN*WORK(INDRV1+N-1) .LT. WORK(INDRV3+N) .AND.  
     $           SAFMIN*WORK(INDRV3+N) .LT. WORK(INDRV1+N-1)) THEN
            TMP2=WORK(INDRV3+N)/WORK(INDRV1+N-1)
            WORK(INDRV2+N-1)=WORK(INDRV4+N-1)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV2+N-1) = WORK(INDRV3+N)
     $           *(WORK(INDRV4+N-1)/WORK(INDRV1+N-1))
            DN2 = WORK(INDRV3+N)*(DN2/WORK(INDRV1+N-1))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV1+N)=DN2

         IF (SISNAN(DN2)) GO TO 1410

         GO TO 1420

 1410    DN2 = WORK(INDRV3+M)
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         DO J=M, N-3
            WORK(INDRV1+J) = DN2+WORK(INDRV4+J)
            IF (WORK(INDRV1+J) .EQ. ZERO) THEN
               WORK(INDRV2+J) = ZERO
               DN2 = WORK(INDRV3+J+1)
            ELSE IF (SAFMIN*WORK(INDRV1+J) .LT. WORK(INDRV3+J+1) .AND.  
     $              SAFMIN*WORK(INDRV3+J+1) .LT. WORK(INDRV1+J)) THEN
               TMP2=WORK(INDRV3+J+1)/WORK(INDRV1+J)
               WORK(INDRV2+J)=WORK(INDRV4+J)*TMP2
               DN2=DN2*TMP2
            ELSE
               WORK(INDRV2+J) = WORK(INDRV3+J+1)
     $              *(WORK(INDRV4+J)/WORK(INDRV1+J))
               DN2 = WORK(INDRV3+J+1)*(DN2/WORK(INDRV1+J))
            ENDIF
            IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO
         ENDDO

         WORK(INDRV1+N-2) = DN2+WORK(INDRV4+N-2)
         IF (SAFMIN*WORK(INDRV1+N-2) .LT. WORK(INDRV3+N-1) .AND.  
     $           SAFMIN*WORK(INDRV3+N-1) .LT. WORK(INDRV1+N-2)) THEN
            TMP2=WORK(INDRV3+N-1)/WORK(INDRV1+N-2)
            WORK(INDRV2+N-2)=WORK(INDRV4+N-2)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV2+N-2) = WORK(INDRV3+N-1)
     $           *(WORK(INDRV4+N-2)/WORK(INDRV1+N-2))
            DN2 = WORK(INDRV3+N-1)*(DN2/WORK(INDRV1+N-2))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV1+N-1) = DN2+WORK(INDRV4+N-1)
         IF (SAFMIN*WORK(INDRV1+N-1) .LT. WORK(INDRV3+N) .AND.  
     $           SAFMIN*WORK(INDRV3+N) .LT. WORK(INDRV1+N-1)) THEN
            TMP2=WORK(INDRV3+N)/WORK(INDRV1+N-1)
            WORK(INDRV2+N-1)=WORK(INDRV4+N-1)*TMP2
            DN2=DN2*TMP2
         ELSE
            WORK(INDRV2+N-1) = WORK(INDRV3+N)
     $           *(WORK(INDRV4+N-1)/WORK(INDRV1+N-1))
            DN2 = WORK(INDRV3+N)*(DN2/WORK(INDRV1+N-1))
         ENDIF
         IF (SIGMA .EQ. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV1+N)=DN2

 1420    DO J= M, N-3
            IF (WORK(INDRV2+J) .LE. SIGMA1 .OR.
     $           WORK(INDRV4+J) .LE. TOL2*WORK(INDRV1+J)) THEN
               WORK(INDRV2+J)=-SIGMA
               WORK(INDRV4+J)=-DESIG
               M=J+1
            ENDIF
         ENDDO
         IF (WORK(INDRV4+N-2) .LE. TOL2*WORK(INDRV1+N-2)) THEN
            WORK(INDRV2+N-2)=ZERO
         ENDIF
         IF (WORK(INDRV4+N-1) .LE. TOL2*WORK(INDRV1+N-1)) THEN
            WORK(INDRV2+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 (WORK(INDRV2+J) .LE. ZERO) THEN
            M=J+1
            GO TO 3000
         ENDIF
      ENDDO
      M=1

      GO TO 3000
  
 4000 A(1:N0)=SQRT(WORK(INDRV1+1:INDRV1+N0))
      CALL SLASRT( 'D', N0, A, IINFO )
      IF (SIGMX .GT. ZERO) THEN
         CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N0, 1, A, N0, IINFO )
      ENDIF
      RETURN
      END SUBROUTINE SDQDS

      SUBROUTINE SDQDS2(SIGMA,DESIG,TAU,SIGMA2,DESIG2)
      IMPLICIT NONE
      REAL SIGMA,DESIG,TAU,T,SIGMA2,DESIG2
      IF( TAU.LT.SIGMA ) THEN
         DESIG2 = DESIG + TAU
         T = SIGMA + DESIG2
         DESIG2 = DESIG2 - ( T-SIGMA )
      ELSE
         T = SIGMA + TAU
         DESIG2 = SIGMA - ( T-TAU ) + DESIG
      END IF
      SIGMA2 = T
      RETURN
      END SUBROUTINE SDQDS2