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, OLDM, OLDN, METHOD
      INTEGER INDRV1,INDRV2,INDRV3,INDRV4,INDRV5,INDRV6,INDRV7,INDRV8
      DOUBLE PRECISION TMP1, TMP2, TMP3, TAU, TAUA(4), TAUB, 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, DMIN1
      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

      TAUA(1) = ZERO
      
      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( EPS / 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

      OLDM = -1
      OLDN = -1

 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 .GE. TAU+SIGMA) 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(WORK(INDRV1+M:INDRV1+N-1))
         ENDIF
         IF (SIGMA .GE. TAUB+SIGMA) GO TO 350

         METHOD = 1
         
         IF (TAUB .LE. TAU) GO TO 140

         TAUB = TAU

         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) THEN
            TMP2 = TAU+DN0
            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     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
         
         TAUA(METHOD+1) = SCALE/TMP1
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF

         TAUA(METHOD+1) = SCALE/SQRT(TMP2)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF

         TMP8 = DBLE(N-M+1)
         TAUA(METHOD+1) = SCALE*(TMP8/
     $        (TMP1+SQRT(TMP8-ONE)*SQRT(TMP8*TMP2-TMP1**2)))
         IF (TAUA(METHOD+1) .EQ. TAUA(METHOD+1)) THEN
            IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
               METHOD = METHOD+1
            ENDIF
         ENDIF
         
         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
               TAUB = MIN(TMP9,SCALE/SQRT(DMIN2),TAUB)
            ELSE
               TAUB = MIN(SCALE/SQRT(DMIN2),TAUB)
            ENDIF
         ELSE
            TAUB = MIN(SCALE/SQRT(DMIN2),TAUB)
         ENDIF

         TAUA(METHOD) = MIN(TAUA(METHOD),TAUB)
         IF (METHOD .GE. 2 .AND.
     $        TAUB .LT. TWO*TAUA(METHOD)) GO TO 124
         METHOD = 1

 150     WORK(INDRV7+M:INDRV7+N) = SQRT(ONE/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
         IF (TMP1 .GT. ZERO) THEN
            TAUA(METHOD+1) = TMP3/TMP1
            IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
               METHOD = METHOD+1
            ENDIF
         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
            TAUA(METHOD+1) = 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)=SQRT(WORK(INDRV1+N))-HALF*WORK(INDRV8+N-1)
         IF (TAUA(METHOD+1) .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
            TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
         ENDDO
         TMP2=SQRT(WORK(INDRV1+M))-HALF*WORK(INDRV8+M)
         IF (TMP2 .LE. ZERO) GO TO 350
         TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)**2
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF

 124     TAU = TAUA(METHOD)
 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
            TMP2 = TAU+DMIN2
            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
         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
            TMP2 = TAU+DN1
            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
         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) THEN
            TMP2 = TAU+DN0
            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
            DN0 = ZERO
         ENDIF
         IF (DN0 .GE. WORK(INDRV1+N)) THEN
            TMP3 = DN1/WORK(INDRV3+N-1)-TAU/WORK(INDRV1+N)
            IF (TMP3 .GE. ZERO) THEN
               DN0 = MIN(WORK(INDRV1+N)*TMP3,WORK(INDRV1+N))
            ELSE
               DN0 = WORK(INDRV1+N)
            ENDIF
         ENDIF
         
         WORK(INDRV3+N)=DN0

         IF (DISNAN(DN0)) THEN
            TMP2 = CONST*TAU
            IF (TMP2 .GE. TAU) TMP2 = HALF*TAU
            TAU = TMP2
            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 .GE. 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 .GE. 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 .GE. 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 .GE. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV3+N)=DN2
         
         IF (DISNAN(DN2)) GO TO 370
         
         GO TO 400
         
 370     DN2 = WORK(INDRV1+M)
         IF (SIGMA .GE. 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 .GE. 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 .GE. 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 .GE. 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

         OLDM = M
         OLDN = N
         
         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 .GE. DN2+SIGMA) DN2 = ZERO
         DMIN1 = MIN(DMIN2,DN2)
         
         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 .GE. 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 .GE. DN2+SIGMA) DN2 = ZERO
         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
            IF (SIGMA .GE. DN2+SIGMA) DN2 = ZERO
            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 .GE. DN2+SIGMA) DN2 = ZERO
         DMIN1 = MIN(DMIN2,DN2)

         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 .GE. DN2+SIGMA) DN2 = ZERO

         WORK(INDRV1+N)=DN2

         IF (DISNAN(DN2)) GO TO 1410

         GO TO 1420

 1410    DN2 = WORK(INDRV3+M)
         IF (SIGMA .GE. DN2+SIGMA) DN2 = ZERO
         DMIN2 = DN2
         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 .GE. DN2+SIGMA) DN2 = ZERO
            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 .GE. DN2+SIGMA) DN2 = ZERO
         DMIN1 = MIN(DMIN2,DN2)

         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 .GE. 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