dqds, equivalent changes to single precision versions

diff --git a/SRC/slasq3.f b/SRC/slasq3.f
index 458d798..4b7d825 100644 (file)
--- a/SRC/slasq3.f
+++ b/SRC/slasq3.f
@@ -331,15 +331,15 @@
 *
    70 CONTINUE
 *
-      CALL SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
-     $             DN1, DN2, IEEE )
+      CALL SLASQ5( 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/slasq5.f b/SRC/slasq5.f
index 85dab34..b7db2d9 100644 (file)
--- a/SRC/slasq5.f
+++ b/SRC/slasq5.f
@@ -129,8 +129,8 @@
 *> \ingroup auxOTHERcomputational
 *
 *  =====================================================================
-      SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
-     $                   DNM1, DNM2, IEEE )
+      SUBROUTINE SLASQ5( 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 +140,8 @@
 *     .. Scalar Arguments ..
       LOGICAL            IEEE
       INTEGER            I0, N0, PP
-      REAL               DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
+      REAL               DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
+     $                   SIGMA, EPS
 *     ..
 *     .. Array Arguments ..
       REAL               Z( * )
@@ -149,12 +150,12 @@
 *  =====================================================================
 *
 *     .. Parameter ..
-      REAL               ZERO
-      PARAMETER          ( ZERO = 0.0E0 )
+      REAL               ZERO, HALF
+      PARAMETER          ( ZERO = 0.0E0, HALF = 0.5 )
 *     ..
 *     .. Local Scalars ..
       INTEGER            J4, J4P2
-      REAL               D, EMIN, TEMP
+      REAL               D, EMIN, TEMP, DTHRESH
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MIN
@@ -164,114 +165,231 @@
       IF( ( N0-I0-1 ).LE.0 )
      $   RETURN
 *
-      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 10 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
-               DMIN = MIN( DMIN, D )
-               Z( J4 ) = Z( J4-1 )*TEMP
-               EMIN = MIN( Z( J4 ), EMIN )
-   10       CONTINUE
+      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
+         DMIN = D
+         DMIN1 = -Z( J4 )
+*     
+         IF( IEEE ) THEN
+*     
+*     Code for IEEE arithmetic.
+*     
+            IF( PP.EQ.0 ) THEN
+               DO 10 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
+                  DMIN = MIN( DMIN, D )
+                  Z( J4 ) = Z( J4-1 )*TEMP
+                  EMIN = MIN( Z( J4 ), EMIN )
+ 10            CONTINUE
+            ELSE
+               DO 20 J4 = 4*I0, 4*( N0-3 ), 4
+                  Z( J4-3 ) = D + Z( J4 )
+                  TEMP = Z( J4+2 ) / Z( J4-3 )
+                  D = D*TEMP - TAU
+                  DMIN = MIN( DMIN, D )
+                  Z( J4-1 ) = Z( J4 )*TEMP
+                  EMIN = MIN( Z( J4-1 ), EMIN )
+ 20            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
-            DO 20 J4 = 4*I0, 4*( N0-3 ), 4
-               Z( J4-3 ) = D + Z( J4 )
-               TEMP = Z( J4+2 ) / Z( J4-3 )
-               D = D*TEMP - TAU
-               DMIN = MIN( DMIN, D )
-               Z( J4-1 ) = Z( J4 )*TEMP
-               EMIN = MIN( Z( J4-1 ), EMIN )
-   20       CONTINUE
+*     
+*     Code for non IEEE arithmetic.
+*     
+            IF( PP.EQ.0 ) THEN
+               DO 30 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
+                  DMIN = MIN( DMIN, D )
+                  EMIN = MIN( EMIN, Z( J4 ) )
+ 30            CONTINUE
+            ELSE
+               DO 40 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
+                  DMIN = MIN( DMIN, D )
+                  EMIN = MIN( EMIN, Z( J4-1 ) )
+ 40            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
 *
-*        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 30 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
-               DMIN = MIN( DMIN, D )
-               EMIN = MIN( EMIN, Z( J4 ) )
-   30       CONTINUE
-         ELSE
-            DO 40 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
-               DMIN = MIN( DMIN, D )
-               EMIN = MIN( EMIN, Z( J4-1 ) )
-   40       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
+*     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
-         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
+            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
-         DMIN = MIN( DMIN, DN )
-*
+*     
       END IF
-*
       Z( J4+2 ) = DN
       Z( 4*N0-PP ) = EMIN
       RETURN