SUBROUTINE DSMALLSIGMA(FLAG,METHOD,M,N,A,B,TAUA,TAUB,TAU,
     $     WORK2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7, RTMIN, RTMAX)

      IMPLICIT NONE
      
      DOUBLE PRECISION A(*), B(*), WORK2(*)
      INTEGER N, M, J, METHOD, FLAG
      INTEGER INDRV3, INDRV4, INDRV5, INDRV6, INDRV7
      DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, TMP5
      DOUBLE PRECISION TAUA(4), TAUB, TAU
      DOUBLE PRECISION C1, S1, T, RTMIN, RTMAX
      
      DOUBLE PRECISION ONE, ZERO, HALF, TWO, CONST
      PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0, TWO = 2.0D0)
      PARAMETER (CONST = 0.75D0)
      EXTERNAL DFMA0
      DOUBLE PRECISION DFMA0

      METHOD = 1
      IF (TAUB .LE. TAU) GO TO 160
*         
      TAUB = TAU
*     
      TMP4 = A(M)-TAU
      TMP5 = A(M)+TAU
      IF (TMP4 .LE. ZERO) THEN
         GO TO 160
      ELSE
         TMP3 = SQRT(TMP4)*SQRT(TMP5)
      ENDIF
      DO J = M, N-2
         CALL DLARTG2(TMP3,B(J),C1,S1,T,RTMIN,RTMAX)
         TMP4 = DFMA0(C1,A(J+1),-TAU)
         TMP5 = DFMA0(C1,A(J+1),TAU)
         IF (TMP4 .LE. ZERO) THEN
            GO TO 160
         ELSE
            TMP3 = SQRT(TMP4)*SQRT(TMP5)
         ENDIF
      ENDDO
      CALL DLARTG2(TMP3,B(N-1),C1,S1,T,RTMIN,RTMAX)
      TMP4 = DFMA0(C1,A(N),-TAU)
      IF (TMP4 .LT. ZERO) THEN
         TMP2 = C1*A(N)
         IF (TMP2 .LT. TAU) THEN
            IF (FLAG .EQ. 0) THEN
               TAUA(METHOD+1) = TMP2
            ELSE
               TMP3 = CONST*TAU
               IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
               TAUA(METHOD+1) = MAX(TMP2,TMP3)
            ENDIF
         ELSE
            GO TO 160
         ENDIF
      ELSE
         TAUA(METHOD+1) = TAU
      ENDIF
      IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
         METHOD = METHOD+1
      ENDIF
      IF (METHOD .GE. 2) RETURN
*
 160  WORK2(INDRV7+M:INDRV7+N) = ONE/A(M:N) 

      WORK2(INDRV5+N) = WORK2(INDRV7+N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV3+J) = B(J)*WORK2(INDRV7+J)
         WORK2(INDRV5+J) = WORK2(INDRV7+J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      TMP3 = ONE/TMP3
      WORK2(INDRV5+M) = (WORK2(INDRV5+M)*TMP3)*WORK2(INDRV7+M)
      TMP1 = WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV4+J) = B(J-1)*WORK2(INDRV7+J)
         WORK2(INDRV5+J) = (WORK2(INDRV5+J)*TMP3)*WORK2(INDRV7+J)+
     $        WORK2(INDRV4+J)*WORK2(INDRV5+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV5+J))
      ENDDO
      IF (TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(SQRT(TMP3/TMP1),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      TMP1 = ONE/TMP1
      WORK2(INDRV5+N) = WORK2(INDRV5+N)*TMP1
      WORK2(INDRV6+N) = WORK2(INDRV5+N)*WORK2(INDRV7+N)
      TMP3 = WORK2(INDRV6+N)
      DO J = N-1,M,-1
         WORK2(INDRV5+J) = WORK2(INDRV5+J)*TMP1
         WORK2(INDRV6+J) = WORK2(INDRV5+J)*WORK2(INDRV7+J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV6+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV6+J))
      ENDDO
      TMP3 = ONE/TMP3
      WORK2(INDRV6+M) = (WORK2(INDRV6+M)*TMP3)*WORK2(INDRV7+M)
      TMP1 = WORK2(INDRV6+M)/WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV6+J) = (WORK2(INDRV6+J)*TMP3)*WORK2(INDRV7+J)
     $        +WORK2(INDRV4+J)*WORK2(INDRV6+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV6+J)/WORK2(INDRV5+J))
      ENDDO
      IF (TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(SQRT(TMP3/TMP1),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      TMP1 = ONE/TMP1
      WORK2(INDRV6+N) = WORK2(INDRV6+N)*TMP1
      WORK2(INDRV5+N) = WORK2(INDRV6+N)*WORK2(INDRV7+N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV6+J) = WORK2(INDRV6+J)*TMP1
         WORK2(INDRV5+J) = WORK2(INDRV6+J)*WORK2(INDRV7+J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      TMP3 = ONE/TMP3
      TMP2 = (WORK2(INDRV5+M)*TMP3)*WORK2(INDRV7+M)
      TMP1 = TMP2/WORK2(INDRV6+M)
      DO J = M+1,N
         TMP2 = (WORK2(INDRV5+J)*TMP3)*WORK2(INDRV7+J)
     $        +WORK2(INDRV4+J)*TMP2
         TMP1 = MAX(TMP1,TMP2/WORK2(INDRV6+J))
      ENDDO
      IF (TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(SQRT(TMP3/TMP1),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF

      IF (METHOD .GE. 2) RETURN
      
 170  WORK2(INDRV5+N) = ONE/A(N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV3+J) = B(J)/A(J)
         WORK2(INDRV5+J) = ONE/A(J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      WORK2(INDRV5+M) = (WORK2(INDRV5+M)/TMP3)/A(M)
      TMP1 = WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV4+J) = B(J-1)/A(J)
         WORK2(INDRV5+J) = (WORK2(INDRV5+J)/TMP3)/A(J)+
     $        WORK2(INDRV4+J)*WORK2(INDRV5+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV5+J))
      ENDDO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(ONE/SQRT(TMP1)/SQRT(TMP3),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      WORK2(INDRV5+N) = WORK2(INDRV5+N)/TMP1
      WORK2(INDRV6+N) = WORK2(INDRV5+N)/A(N)
      TMP3 = WORK2(INDRV6+N)
      DO J = N-1,M,-1
         WORK2(INDRV5+J) = WORK2(INDRV5+J)/TMP1
         WORK2(INDRV6+J) = WORK2(INDRV5+J)/A(J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV6+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV6+J))
      ENDDO
      WORK2(INDRV6+M) = (WORK2(INDRV6+M)/TMP3)/A(M)
      TMP1 = WORK2(INDRV6+M)/WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV6+J) = (WORK2(INDRV6+J)/TMP3)/A(J)+
     $        WORK2(INDRV4+J)*WORK2(INDRV6+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV6+J)/WORK2(INDRV5+J))
      ENDDO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(ONE/SQRT(TMP1)/SQRT(TMP3),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF

      WORK2(INDRV6+N) = WORK2(INDRV6+N)/TMP1
      WORK2(INDRV5+N) = WORK2(INDRV6+N)/A(N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV6+J) = WORK2(INDRV6+J)/TMP1
         WORK2(INDRV5+J) = WORK2(INDRV6+J)/A(J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      TMP2 = (WORK2(INDRV5+M)/TMP3)/A(M)
      TMP1 = TMP2/WORK2(INDRV6+M)
      DO J = M+1,N
         TMP2 = (WORK2(INDRV5+J)/TMP3)/A(J)+
     $        WORK2(INDRV4+J)*TMP2
         TMP1 = MAX(TMP1,TMP2/WORK2(INDRV6+J))
      ENDDO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(ONE/SQRT(TMP1)/SQRT(TMP3),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF

      IF (METHOD .GE. 2) RETURN

      TAUA(METHOD+1) = A(N)-HALF*B(N-1)
      IF (TAUA(METHOD+1) .LE. ZERO) RETURN
      DO J=N-1,M+1,-1
         TMP2 = A(J)-HALF*(B(J)+B(J-1))
         IF (TMP2 .LE. ZERO) RETURN
         TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
      ENDDO
      TMP2=A(M)-HALF*B(M)
      IF (TMP2 .LE. ZERO) RETURN
      TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
      IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
         METHOD = METHOD+1
      ENDIF
      
      RETURN
      END SUBROUTINE DSMALLSIGMA

      SUBROUTINE DOQDQ(N, A, B, WORK2, INDRV5, INDRV6, SIGMA1,
     $     EPS, RTMIN, RTMAX, FLAG)

      IMPLICIT NONE
      INTEGER N, INDRV5, INDRV6, J, FLAG
      DOUBLE PRECISION A(*), B(*), WORK2(*)
      DOUBLE PRECISION TMP1, C1, S1
      DOUBLE PRECISION EPS, RTMIN, RTMAX, SIGMA1
      DOUBLE PRECISION ONE, ZERO, HALF, TWO
      PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0, TWO = 2.0D0)

      FLAG = 0
      
      TMP1 = A(1)
      IF (TMP1 .LE. SIGMA1) THEN
         TMP1 = ZERO
         FLAG = 1
      ENDIF
      DO J = 1, N-1
         IF (TMP1 .EQ. ZERO .AND. B(J) .EQ. ZERO) THEN
            C1 = ZERO
            S1 = ONE
            WORK2(INDRV5+J) = ZERO
         ELSE
            CALL DLARTG2(TMP1,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
         ENDIF
         WORK2(INDRV6+J) = S1*A(J+1)
         TMP1 = C1*A(J+1)
         IF (TMP1 .LE. SIGMA1) THEN
            TMP1 = ZERO
            FLAG = 1
         ENDIF
      ENDDO
      WORK2(INDRV5+N) = TMP1

      END SUBROUTINE DOQDQ
      
      SUBROUTINE DOQDSQ(L,N0,M0,A,B,SU,LDSU,WORK2,WORK3,
     $     IWORK,IWORK2,INFO)

      IMPLICIT NONE
      
      INTEGER N0, N1, M0, INFO, LDSU, FLAG, OLDN, IWORK(N0),
     $     IWORK2(M0)
      DOUBLE PRECISION A(*), B(*), WORK2(*), WORK3(2*N0,M0)
      DOUBLE PRECISION SU(LDSU, *)
      INTEGER N, I, I0, J1, J, L
      INTEGER INDRV3, INDRV4, INDRV5, INDRV6, INDRV7, METHOD
      DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, TMP5
      DOUBLE PRECISION TAU, TAUA(4), TAUB
      DOUBLE PRECISION SIGMA, SIGMA1, DESIG, SIGMA2, T, DESIG0
      DOUBLE PRECISION C1, S1, DMIN1
      DOUBLE PRECISION EPS, TOL, SAFMIN, SAFMAX, RTMIN, RTMAX
      
      DOUBLE PRECISION ONE, ZERO, HALF, TWO, CONST
      PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0, TWO = 2.0D0)
      PARAMETER (CONST = 0.75D0)
      DOUBLE PRECISION HUNDRD 
      PARAMETER (HUNDRD = 100.0D0)
*     
      EXTERNAL DLAMCH
      DOUBLE PRECISION DLAMCH
      EXTERNAL DFMA0
      DOUBLE PRECISION DFMA0
      LOGICAL LSAME
      EXTERNAL LSAME
*
!$omp parallel
!$omp single
      TAUA(1) = ZERO
      
      INDRV3 = 0
      INDRV4 = INDRV3+N0
      INDRV5 = INDRV4+N0
      INDRV6 = INDRV5+N0
      INDRV7 = INDRV6+N0
*     
      EPS = DLAMCH( 'Precision' )
      SAFMIN = DLAMCH( 'Safe minimum' )
      RTMIN = SQRT( SAFMIN )
      SAFMAX = ONE/SAFMIN
      RTMAX = SQRT( SAFMAX/2.0D0 )
      TOL = HUNDRD*EPS
*     
      DO I = 1, N0-1
         IF (A(I) .LT. ZERO) THEN
            A(I) = -A(I)
            B(I) = -B(I)
!            CALL DSCAL(N0,-ONE,SU(1,I),1)
            IWORK(I) = -1
         ELSE
            IWORK(I) = 1
         ENDIF
         IF (B(I) .LT. ZERO) THEN
            B(I) = -B(I)
            A(I+1) = -A(I+1)
         ENDIF
      ENDDO
      IF (A(N0) .LT. ZERO) THEN
         A(N0) = -A(N0)
!         CALL DSCAL(N0,-ONE,SU(1,N0),1)
         IWORK(N0) = -1
      ELSE
         IWORK(N0) = 1
      ENDIF
*
      OLDN = -1
      N = N0
*
 3000 SIGMA = ZERO
      SIGMA1 = SAFMIN
      DESIG = ZERO
      SIGMA2 = TOL*SIGMA
      
      DO I = 1, M0+1
*
 15      IF (N .LE. N0-L) GO TO 4000
*         
         IF (N .EQ. 1) THEN
            
            CALL DLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)

!$omp task firstprivate(N) private(N1,J1,J,I0,TMP1,TMP2,TMP3,TMP4,C1,S1)
            N1 = N0+1-N
            DO J = 1,N0
               SU(J,N1)=ZERO
            ENDDO
            SU(N,N1)=ONE
            DO I0=I-1,1,-1
               J1 = IWORK2(I0)
               TMP2 = SU( J1, N1)
               DO J = J1-1,1,-1
                  TMP4 = WORK3(2*J-1,I0)
                  TMP3 = WORK3(2*J,I0)
                  TMP1 = SU( J, N1)
                  IF (SIGN(ONE,TMP4) .GT. ZERO) THEN
                     SU( J+1, N1) = TMP4*(TMP3*TMP2+TMP1)+TMP2
                     TMP2 = TMP4*(TMP3*TMP1-TMP2)+TMP1
                  ELSE
                     TMP2 = -TMP2
                     SU( J+1, N1) = TMP4*(TMP3*TMP1+TMP2)+TMP1
                     TMP2 = TMP4*(TMP3*TMP2-TMP1)+TMP2
                  ENDIF
               ENDDO
               SU( 1, N1) = TMP2
            ENDDO
            DO J = N0,1,-1
               IF (IWORK(J) .LT. 0) THEN
                  SU(J,N1) = - SU(J,N1)
               ENDIF
            ENDDO
!$omp end task
            
            GO TO 4000
            
         ENDIF
*
         IF (B(N-1) .LE. SIGMA2) THEN

            IF(A(N) .EQ. ZERO) THEN
               A(N) = SIGMA

!$omp task firstprivate(N) private(N1,J1,J,I0,TMP1,TMP2,TMP3,TMP4,C1,S1)
               N1 = N0+1-N
               DO J = 1,N0
                  SU(J,N1)=ZERO
               ENDDO
               SU(N,N1)=ONE
               DO I0=I-1,1,-1
                  J1 = IWORK2(I0)
                  TMP2 = SU( J1, N1)
                  DO J = J1-1,1,-1
                     TMP4 = WORK3(2*J-1,I0)
                     TMP3 = WORK3(2*J,I0)
                     TMP1 = SU( J, N1)
                     IF (SIGN(ONE,TMP4) .GT. ZERO) THEN
                        SU( J+1, N1) = TMP4*(TMP3*TMP2+TMP1)+TMP2
                        TMP2 = TMP4*(TMP3*TMP1-TMP2)+TMP1
                     ELSE
                        TMP2 = -TMP2
                        SU( J+1, N1) = TMP4*(TMP3*TMP1+TMP2)+TMP1
                        TMP2 = TMP4*(TMP3*TMP2-TMP1)+TMP2
                     ENDIF
                  ENDDO
                  SU( 1, N1) = TMP2
               ENDDO
               DO J = N0,1,-1
                  IF (IWORK(J) .LT. 0) THEN
                     SU(J,N1) = - SU(J,N1)
                  ENDIF
               ENDDO
!$omp end task
               
               N = N-1
               IF (N .LE. 0) THEN
                  GO TO 4000
               ENDIF

               GO TO 15
               
            ENDIF

            CALL DLARTG7(SIGMA,DESIG,A(N),T,DESIG0,RTMIN)
            IF(T .LE. SIGMA) THEN
               A(N) = T

!$omp task firstprivate(N) private(N1,J1,J,I0,TMP1,TMP2,TMP3,TMP4,C1,S1)
               N1 = N0+1-N
               DO J = 1,N0
                  SU(J,N1)=ZERO
               ENDDO
               SU(N,N1)=ONE
               DO I0=I-1,1,-1
                  J1 = IWORK2(I0)
                  TMP2 = SU( J1, N1)
                  DO J = J1-1,1,-1
                     TMP4 = WORK3(2*J-1,I0)
                     TMP3 = WORK3(2*J,I0)
                     TMP1 = SU( J, N1)
                     IF (SIGN(ONE,TMP4) .GT. ZERO) THEN
                        SU( J+1, N1) = TMP4*(TMP3*TMP2+TMP1)+TMP2
                        TMP2 = TMP4*(TMP3*TMP1-TMP2)+TMP1
                     ELSE
                        TMP2 = -TMP2
                        SU( J+1, N1) = TMP4*(TMP3*TMP1+TMP2)+TMP1
                        TMP2 = TMP4*(TMP3*TMP2-TMP1)+TMP2
                     ENDIF
                  ENDDO
                  SU( 1, N1) = TMP2
               ENDDO
               DO J = N0,1,-1
                  IF (IWORK(J) .LT. 0) THEN
                     SU(J,N1) = - SU(J,N1)
                  ENDIF
               ENDDO
!$omp end task
               
               N = N-1
               IF (N .LE. 0) THEN
                  GO TO 4000
               ENDIF

               GO TO 15
               
            ENDIF

            IF (N .EQ. 2) THEN
               TAU = A(N-1)
            ELSE
               TAUB = MINVAL(A(1:N-2))
               CALL DLAS2(A(N-2), B(N-2), A(N-1), TMP2, TMP3)
               TAU = MIN(TMP2,A(N-1))
               CALL DSMALLSIGMA(0,METHOD,1,N-1,A,B,TAUA,TAUB,TAU,
     $              WORK2,INDRV3, INDRV4, INDRV5, INDRV6, INDRV7,
     $              RTMIN, RTMAX)
               TAU = TAUA(METHOD)
            ENDIF
               
            IF(A(N) .LE. TAU) THEN
               A(N) = T

!$omp task firstprivate(N) private(N1,J1,J,I0,TMP1,TMP2,TMP3,TMP4,C1,S1)
               N1 = N0+1-N
               DO J = 1,N0
                  SU(J,N1)=ZERO
               ENDDO
               SU(N,N1)=ONE
               DO I0=I-1,1,-1
                  J1 = IWORK2(I0)
                  TMP2 = SU( J1, N1)
                  DO J = J1-1,1,-1
                     TMP4 = WORK3(2*J-1,I0)
                     TMP3 = WORK3(2*J,I0)
                     TMP1 = SU( J, N1)
                     IF (SIGN(ONE,TMP4) .GT. ZERO) THEN
                        SU( J+1, N1) = TMP4*(TMP3*TMP2+TMP1)+TMP2
                        TMP2 = TMP4*(TMP3*TMP1-TMP2)+TMP1
                     ELSE
                        TMP2 = -TMP2
                        SU( J+1, N1) = TMP4*(TMP3*TMP1+TMP2)+TMP1
                        TMP2 = TMP4*(TMP3*TMP2-TMP1)+TMP2
                     ENDIF
                  ENDDO
                  SU( 1, N1) = TMP2
               ENDDO
               DO J = N0,1,-1
                  IF (IWORK(J) .LT. 0) THEN
                     SU(J,N1) = - SU(J,N1)
                  ENDIF
               ENDDO
!$omp end task
                
               N = N-1
               IF (N .LE. 0) THEN
                  GO TO 4000
               ENDIF

               GO TO 15
               
            ENDIF
         ENDIF
*
         IF (I .EQ. M0+1) EXIT
*
         CALL DLAS2(A(N-1), B(N-1), A(N), TAU, TMP3)
         TAU = MIN(TAU,A(N))
         IF (TAU .LE. SIGMA1) THEN
            CALL DOQDQ(N, A, B, WORK2, INDRV5, INDRV6, SIGMA1,
     $           EPS, RTMIN, RTMAX, FLAG)
            IF (FLAG .EQ. 1) GO TO 360
         ENDIF

         IF (OLDN .EQ. N) THEN
            TAUB = DMIN1
         ELSE
            TAUB = MINVAL(A(1:N-1))
         ENDIF
         IF (TAUB .LE. SIGMA1) THEN
            CALL DOQDQ(N, A, B, WORK2, INDRV5, INDRV6, SIGMA1,
     $           EPS, RTMIN, RTMAX, FLAG)
            IF (FLAG .EQ. 1) GO TO 360
         ENDIF
            
         CALL DSMALLSIGMA(1,METHOD,1,N,A,B,TAUA,TAUB,TAU,
     $        WORK2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7,
     $        RTMIN, RTMAX)
         TAU = TAUA(METHOD)
*     
 125     IF (TAU .LE. ZERO) THEN
            CALL DOQDQ(N, A, B, WORK2, INDRV5, INDRV6, SIGMA1,
     $           EPS, RTMIN, RTMAX, FLAG)
            GO TO 360
         ENDIF
*     
         TMP4 = A(1)-TAU
         TMP5 = A(1)+TAU
         IF (TMP4 .LE. ZERO) THEN
            IF (TMP4 .EQ. ZERO .AND. B(1) .EQ. ZERO) THEN
               TMP3 = ZERO
            ELSE
               TMP2 = A(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
            ENDIF
         ELSE
            TMP3 = SQRT(TMP4)*SQRT(TMP5)
         ENDIF
         
         DO J = 1, N-2
            CALL DLARTG2(TMP3,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
            WORK2(INDRV6+J) = S1*A(J+1)
            
            TMP4 = DFMA0(C1,A(J+1),-TAU)
            TMP5 = DFMA0(C1,A(J+1),TAU)
            IF (TMP4 .LE. ZERO) THEN
               IF (TMP4 .EQ. ZERO .AND. B(J+1) .EQ. ZERO) THEN
                  TMP3 = ZERO
               ELSE
                  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
               ENDIF
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
         ENDDO
         CALL DLARTG2(TMP3,B(N-1),C1,S1,WORK2(INDRV5+N-1),RTMIN,RTMAX)
         WORK2(INDRV6+N-1) = S1*A(N)
         
         TMP4 = DFMA0(C1,A(N),-TAU)
         TMP5 = DFMA0(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)*SQRT(TMP5)
         ENDIF
         
         IF (TMP3 .GE. A(N)) THEN
            TMP2 = TAU/A(N)
            IF (C1 .GE. TMP2) THEN
               TMP3 = MIN(A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2),A(N))
            ELSE
               TMP3 = A(N)
            ENDIF
         ENDIF
         WORK2(INDRV5+N) = TMP3
*
         CALL DLARTG7(SIGMA,DESIG,TAU,T,DESIG0,RTMIN)
*     
         SIGMA = T
         SIGMA1 = MAX(SQRT(EPS)*SIGMA,SAFMIN)
         DESIG = DESIG0
         SIGMA2 = TOL*SIGMA
*     
 360     OLDN = N

         TMP1 = WORK2(INDRV5+1)
         IF (TMP1 .LE. SAFMIN) TMP1 = ZERO
         DMIN1 = TMP1
         DO J = 1, N-2
            IF (WORK2(INDRV6+J) .LE. EPS*TMP1) WORK2(INDRV6+J) = ZERO
            CALL DLARTG2(TMP1,WORK2(INDRV6+J),C1,S1,A(J),RTMIN,RTMAX)
            IF (C1 .GE. S1) THEN
               WORK3(2*J-1,I)=SIGN(S1,ONE)
               WORK3(2*J,I)=-S1/(ONE+C1)
            ELSE
               WORK3(2*J-1,I)=SIGN(C1,-ONE)
               WORK3(2*J,I)=C1/(ONE+S1)
            ENDIF
            B(J) = S1*WORK2(INDRV5+J+1)
            TMP1 = C1*WORK2(INDRV5+J+1)
            IF (TMP1 .LE. SAFMIN) TMP1 = ZERO
            DMIN1 = MIN(DMIN1,TMP1)
         ENDDO
         IF (WORK2(INDRV6+N-1) .LE. EPS*TMP1) WORK2(INDRV6+N-1) = ZERO
         CALL DLARTG2(TMP1,WORK2(INDRV6+N-1),C1,S1,A(N-1),RTMIN,RTMAX)
         IF (C1 .GE. S1) THEN
            WORK3(2*N-3,I)=SIGN(S1,ONE)
            WORK3(2*N-2,I)=-S1/(ONE+C1)
         ELSE
            WORK3(2*N-3,I)=SIGN(C1,-ONE)
            WORK3(2*N-2,I)=C1/(ONE+S1)
         ENDIF
         B(N-1) = S1*WORK2(INDRV5+N)
         TMP1 = C1*WORK2(INDRV5+N)
         IF (TMP1 .LE. SAFMIN) TMP1 = ZERO
         A(N) = TMP1
*
         DO J = 1, N-1
            IF (B(J) .LE. SIGMA2) THEN
               B(J) = ZERO
            ENDIF
         ENDDO
*         
 400     IWORK2(I)=N
      ENDDO
      INFO = 2
      GO TO 4100
*     
 4000 INFO = 0
 4100 CONTINUE
!$omp taskwait
!$omp end single
!$omp end parallel
      RETURN
*     
      END SUBROUTINE DOQDSQ