In dqds, set shift to zero if it is insignificant compared to the cumulative shift...

diff --git a/SRC/dlasq3.f b/SRC/dlasq3.f
index a1ce68b..343427f 100644 (file)
--- a/SRC/dlasq3.f
+++ b/SRC/dlasq3.f
@@ -331,15 +331,15 @@
 *
    70 CONTINUE
 *
-      CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
-     $             DN1, DN2, IEEE )
+      CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
+     $             DN1, DN2, IEEE, EPS )
 *
       NDIV = NDIV + ( N0-I0+2 )
       ITER = ITER + 1
 *
 *     Check status.
 *
-      IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
+      IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
 *
 *        Success.
 *

diff --git a/SRC/dlasq5.f b/SRC/dlasq5.f
index d78dd86..3b113ae 100644 (file)
--- a/SRC/dlasq5.f
+++ b/SRC/dlasq5.f
@@ -18,13 +18,13 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
-*                          DNM1, DNM2, IEEE )
+*       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
+*                          DNM1, DNM2, IEEE, EPS )
 * 
 *       .. Scalar Arguments ..
 *       LOGICAL            IEEE
 *       INTEGER            I0, N0, PP
-*       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
+*       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
 *       ..
 *       .. Array Arguments ..
 *       DOUBLE PRECISION   Z( * )
@@ -74,6 +74,12 @@
 *>        This is the shift.
 *> \endverbatim
 *>
+*> \param[in] SIGMA
+*> \verbatim
+*>          SIGMA is DOUBLE PRECISION
+*>        This is the accumulated shift up to this step.
+*> \endverbatim
+*>
 *> \param[out] DMIN
 *> \verbatim
 *>          DMIN is DOUBLE PRECISION
@@ -116,6 +122,12 @@
 *>        Flag for IEEE or non IEEE arithmetic.
 *> \endverbatim
 *
+*> \param[in] EPS
+*> \verbatim
+*>          EPS is DOUBLE PRECISION
+*>        This is the value of epsilon used.
+*> \endverbatim
+*>
 *  Authors:
 *  ========
 *
@@ -129,8 +141,8 @@
 *> \ingroup auxOTHERcomputational
 *
 *  =====================================================================
-      SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
-     $                   DNM1, DNM2, IEEE )
+      SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
+     $                   DN, DNM1, DNM2, IEEE, EPS )
 *
 *  -- LAPACK computational routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -140,7 +152,8 @@
 *     .. Scalar Arguments ..
       LOGICAL            IEEE
       INTEGER            I0, N0, PP
-      DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
+      DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
+     $                   SIGMA, EPS
 *     ..
 *     .. Array Arguments ..
       DOUBLE PRECISION   Z( * )
@@ -149,12 +162,12 @@
 *  =====================================================================
 *
 *     .. Parameter ..
-      DOUBLE PRECISION   ZERO
-      PARAMETER          ( ZERO = 0.0D0 )
+      DOUBLE PRECISION   ZERO, HALF
+      PARAMETER          ( ZERO = 0.0D0, HALF = 0.5 )
 *     ..
 *     .. Local Scalars ..
       INTEGER            J4, J4P2
-      DOUBLE PRECISION   D, EMIN, TEMP
+      DOUBLE PRECISION   D, EMIN, TEMP, DTHRESH
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MIN
@@ -164,6 +177,9 @@
       IF( ( N0-I0-1 ).LE.0 )
      $   RETURN
 *
+      DTHRESH = EPS*(SIGMA+TAU)
+      IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
+      IF( TAU.NE.ZERO ) THEN
       J4 = 4*I0 + PP - 3
       EMIN = Z( J4+4 ) 
       D = Z( J4 ) - TAU
@@ -271,7 +287,120 @@
          DMIN = MIN( DMIN, DN )
 *
       END IF
-*
+      ELSE
+*     This is the version that sets d's to zero if they are small enough
+         J4 = 4*I0 + PP - 3
+         EMIN = Z( J4+4 ) 
+         D = Z( J4 ) - TAU
+         DMIN = D
+         DMIN1 = -Z( J4 )
+         IF( IEEE ) THEN
+*     
+*     Code for IEEE arithmetic.
+*     
+            IF( PP.EQ.0 ) THEN
+               DO 50 J4 = 4*I0, 4*( N0-3 ), 4
+                  Z( J4-2 ) = D + Z( J4-1 ) 
+                  TEMP = Z( J4+1 ) / Z( J4-2 )
+                  D = D*TEMP - TAU
+                  IF( D.LT.DTHRESH ) D = ZERO
+                  DMIN = MIN( DMIN, D )
+                  Z( J4 ) = Z( J4-1 )*TEMP
+                  EMIN = MIN( Z( J4 ), EMIN )
+ 50            CONTINUE
+            ELSE
+               DO 60 J4 = 4*I0, 4*( N0-3 ), 4
+                  Z( J4-3 ) = D + Z( J4 ) 
+                  TEMP = Z( J4+2 ) / Z( J4-3 )
+                  D = D*TEMP - TAU
+                  IF( D.LT.DTHRESH ) D = ZERO
+                  DMIN = MIN( DMIN, D )
+                  Z( J4-1 ) = Z( J4 )*TEMP
+                  EMIN = MIN( Z( J4-1 ), EMIN )
+ 60            CONTINUE
+            END IF
+*     
+*     Unroll last two steps. 
+*     
+            DNM2 = D
+            DMIN2 = DMIN
+            J4 = 4*( N0-2 ) - PP
+            J4P2 = J4 + 2*PP - 1
+            Z( J4-2 ) = DNM2 + Z( J4P2 )
+            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
+            DMIN = MIN( DMIN, DNM1 )
+*     
+            DMIN1 = DMIN
+            J4 = J4 + 4
+            J4P2 = J4 + 2*PP - 1
+            Z( J4-2 ) = DNM1 + Z( J4P2 )
+            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
+            DMIN = MIN( DMIN, DN )
+*     
+         ELSE
+*     
+*     Code for non IEEE arithmetic.
+*     
+            IF( PP.EQ.0 ) THEN
+               DO 70 J4 = 4*I0, 4*( N0-3 ), 4
+                  Z( J4-2 ) = D + Z( J4-1 ) 
+                  IF( D.LT.ZERO ) THEN
+                     RETURN
+                  ELSE 
+                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
+                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
+                  END IF
+                  IF( D.LT.DTHRESH) D = ZERO
+                  DMIN = MIN( DMIN, D )
+                  EMIN = MIN( EMIN, Z( J4 ) )
+ 70            CONTINUE
+            ELSE
+               DO 80 J4 = 4*I0, 4*( N0-3 ), 4
+                  Z( J4-3 ) = D + Z( J4 ) 
+                  IF( D.LT.ZERO ) THEN
+                     RETURN
+                  ELSE 
+                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
+                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
+                  END IF
+                  IF( D.LT.DTHRESH) D = ZERO
+                  DMIN = MIN( DMIN, D )
+                  EMIN = MIN( EMIN, Z( J4-1 ) )
+ 80            CONTINUE
+            END IF
+*     
+*     Unroll last two steps. 
+*     
+            DNM2 = D
+            DMIN2 = DMIN
+            J4 = 4*( N0-2 ) - PP
+            J4P2 = J4 + 2*PP - 1
+            Z( J4-2 ) = DNM2 + Z( J4P2 )
+            IF( DNM2.LT.ZERO ) THEN
+               RETURN
+            ELSE
+               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
+            END IF
+            DMIN = MIN( DMIN, DNM1 )
+*     
+            DMIN1 = DMIN
+            J4 = J4 + 4
+            J4P2 = J4 + 2*PP - 1
+            Z( J4-2 ) = DNM1 + Z( J4P2 )
+            IF( DNM1.LT.ZERO ) THEN
+               RETURN
+            ELSE
+               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
+            END IF
+            DMIN = MIN( DMIN, DN )
+*     
+         END IF
+      END IF
+*     
       Z( J4+2 ) = DN
       Z( 4*N0-PP ) = EMIN
       RETURN