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

      DOUBLE PRECISION ONE, ZERO, HALF, TWO, HUNDRD, CONST
      PARAMETER ( ONE=1.0D0, ZERO=0.0D0, HALF=0.5D0 )
      PARAMETER ( TWO=2.0D0, HUNDRD=100.0D0, CONST=0.75D0 )

      EXTERNAL DLAMCH,DDQDS2
      DOUBLE PRECISION DLAMCH
      EXTERNAL DISNAN
      LOGICAL DISNAN

      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

      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
      TOL2 = TOL**2
      SAFMIN = DLAMCH( '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 DLASCL( 'G', 0, 0, SIGMX, SCALE, N0, 1, 
     $        WORK(INDRV1+1), N0, IINFO )
         CALL DLASCL( '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 DDQDS2(SIGMA,DESIG,WORK(INDRV1+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL DDQDS2(SIGMA,DESIG,WORK(INDRV1+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            GO TO 700
         ENDIF

         IF (N-M+1 .EQ. 1) THEN
            
            CALL DDQDS2(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 DDQDS2(SIGMA,DESIG,WORK(INDRV1+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL DDQDS2(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 DDQDS2(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 = DBLE(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 (DISNAN(DN0)) THEN
            TAU3 = CONST*TAU
            IF (TAU3 .EQ. TAU) TAU3 = HALF*TAU
            TAU = TAU3
            GO TO 125
         ENDIF

         CALL DDQDS2(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 (DISNAN(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 DDQDS2(SIGMA,DESIG,WORK(INDRV3+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL DDQDS2(SIGMA,DESIG,WORK(INDRV3+N),T,DESIG0)
            WORK(INDRV1+N)=T
            
            GO TO 700
         ENDIF

         IF (N-M+1 .EQ. 1) THEN
            
            CALL DDQDS2(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 DDQDS2(SIGMA,DESIG,WORK(INDRV3+N-1),T,DESIG0)
            WORK(INDRV1+N-1)=T
            CALL DDQDS2(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 DDQDS2(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 (DISNAN(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 (DISNAN(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 DLASRT( 'D', N0, A, IINFO )
      IF (SIGMX .GT. ZERO) THEN
         CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N0, 1, A, N0, IINFO )
      ENDIF
      RETURN
      END SUBROUTINE DDQDS

      SUBROUTINE DDQDS2(SIGMA,DESIG,TAU,SIGMA2,DESIG2)
      IMPLICIT NONE
      DOUBLE PRECISION 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 DDQDS2