Correct bug sent by NAG people

Got a new bug for you! Mick Pont found this problem in DLAED6. The code in question is the 40 loop - when DSCALE(I)=TAU you get a divide by zero (rare in practice). This can cause some compilers to immediately stop, e.g. the Sun compiler.

Mick proposed solution is below:

        DO 40 I = 1, 3
           IF (DSCALE( I ).NE.TAU) THEN
              TEMP = ONE / ( DSCALE( I )-TAU )
              TEMP1 = ZSCALE( I )*TEMP
              TEMP2 = TEMP1*TEMP
              TEMP3 = TEMP2*TEMP
              TEMP4 = TEMP1 / DSCALE( I )
              FC = FC + TEMP4
              ERRETM = ERRETM + ABS( TEMP4 )
              DF = DF + TEMP2
              DDF = DDF + TEMP3
           ELSE
*              On rare occasions dscale(i) can be exactly equal to
*              tau, leading to division by zero. If no trap occurs,
*              there is no problem; the quantities above all overflow
*              and the test on abs(f) below sends you to the end
*              with good results. If a trap occurs, though, the
*              user program will stop. Avoid that happening by
*              jumping directly out.
              GO TO 60
           END IF
  40    CONTINUE

This seems to work OK in our testing.

diff --git a/SRC/dlaed6.f b/SRC/dlaed6.f
index 1b55c1a..e8e9709 100644 (file)
--- a/SRC/dlaed6.f
+++ b/SRC/dlaed6.f
@@ -115,7 +115,7 @@
 *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd. 
 *
-*> \date November 2011
+*> \date April 2012
 *
 *> \ingroup auxOTHERcomputational
 *
@@ -140,10 +140,10 @@
 *  =====================================================================
       SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
 *
-*  -- LAPACK computational routine (version 3.4.0) --
+*  -- LAPACK computational routine (version 3.4.1) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2011
+*     April 2012
 *
 *     .. Scalar Arguments ..
       LOGICAL            ORGATI
@@ -370,15 +370,19 @@
          DF = ZERO
          DDF = ZERO
          DO 40 I = 1, 3
-            TEMP = ONE / ( DSCALE( I )-TAU )
-            TEMP1 = ZSCALE( I )*TEMP
-            TEMP2 = TEMP1*TEMP
-            TEMP3 = TEMP2*TEMP
-            TEMP4 = TEMP1 / DSCALE( I )
-            FC = FC + TEMP4
-            ERRETM = ERRETM + ABS( TEMP4 )
-            DF = DF + TEMP2
-            DDF = DDF + TEMP3
+            IF (DSCALE( I ).NE.TAU) THEN
+               TEMP = ONE / ( DSCALE( I )-TAU )
+               TEMP1 = ZSCALE( I )*TEMP
+               TEMP2 = TEMP1*TEMP
+               TEMP3 = TEMP2*TEMP
+               TEMP4 = TEMP1 / DSCALE( I )
+               FC = FC + TEMP4
+               ERRETM = ERRETM + ABS( TEMP4 )
+               DF = DF + TEMP2
+               DDF = DDF + TEMP3
+            ELSE
+              GO TO 60
+            END IF
    40    CONTINUE
          F = FINIT + TAU*FC
          ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +

diff --git a/SRC/slaed6.f b/SRC/slaed6.f
index 515d5f0..9419419 100644 (file)
--- a/SRC/slaed6.f
+++ b/SRC/slaed6.f
@@ -115,7 +115,7 @@
 *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd. 
 *
-*> \date November 2011
+*> \date April 2012
 *
 *> \ingroup auxOTHERcomputational
 *
@@ -140,10 +140,10 @@
 *  =====================================================================
       SUBROUTINE SLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
 *
-*  -- LAPACK computational routine (version 3.4.0) --
+*  -- LAPACK computational routine (version 3.4.1) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2011
+*     April 2012
 *
 *     .. Scalar Arguments ..
       LOGICAL            ORGATI
@@ -370,15 +370,19 @@
          DF = ZERO
          DDF = ZERO
          DO 40 I = 1, 3
-            TEMP = ONE / ( DSCALE( I )-TAU )
-            TEMP1 = ZSCALE( I )*TEMP
-            TEMP2 = TEMP1*TEMP
-            TEMP3 = TEMP2*TEMP
-            TEMP4 = TEMP1 / DSCALE( I )
-            FC = FC + TEMP4
-            ERRETM = ERRETM + ABS( TEMP4 )
-            DF = DF + TEMP2
-            DDF = DDF + TEMP3
+            IF (DSCALE( I ).NE.TAU) THEN
+               TEMP = ONE / ( DSCALE( I )-TAU )
+               TEMP1 = ZSCALE( I )*TEMP
+               TEMP2 = TEMP1*TEMP
+               TEMP3 = TEMP2*TEMP
+               TEMP4 = TEMP1 / DSCALE( I )
+               FC = FC + TEMP4
+               ERRETM = ERRETM + ABS( TEMP4 )
+               DF = DF + TEMP2
+               DDF = DDF + TEMP3
+            ELSE
+               GO TO 60
+            END IF
    40    CONTINUE
          F = FINIT + TAU*FC
          ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +