Bug fix from Osni Marques, Beresford Parlett and Jim Demmel.

This fixes bugs 032 and 056.
Edits in cstein.f, dstein.f, sstein.f and zstein.f

From Osni on Monday, May 4:

As we discussed in our last conference call, I am attaching a new version of
_STEIN, which fixes bugs 032 and 056.  (The bugs were reported for DSTEIN, but
I have propagated the fix to the other versions.) Also, in the process of
testing the fix for those two bugs with more difficult cases, we stumbled upon
a matrix that led SSTEIN to return NaNs in the eigenvectors ... We fixed that
too and, in summary, these are the changes in _STEIN:

1) The assignment after the GO TO 60 needs to be GPIND = J1 (instead of GPIND =
B1).

2) In 'Normalize and scale the righthand side vector Pb' (after GO TO 100 or GO
TO 120) we have replaced _ASUM with I_AMAX.

diff --git a/SRC/cstein.f b/SRC/cstein.f
index 2f2ae26..ea934ef 100644 (file)
--- a/SRC/cstein.f
+++ b/SRC/cstein.f
@@ -308,7 +308,7 @@
          BLKSIZ = BN - B1 + 1
          IF( BLKSIZ.EQ.1 )
      $      GO TO 60
-         GPIND = B1
+         GPIND = J1
 *
 *        Compute reorthogonalization criterion and stopping criterion.
 *
@@ -381,9 +381,10 @@
 *
 *           Normalize and scale the righthand side vector Pb.
 *
+            JMAX = ISAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
             SCL = BLKSIZ*ONENRM*MAX( EPS,
      $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
-     $            SASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
+     $            ABS( WORK( INDRV1+JMAX ) )
             CALL SSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
 *
 *           Solve the system LU = Pb.

diff --git a/SRC/dstein.f b/SRC/dstein.f
index 7cc372b..ddc84fe 100644 (file)
--- a/SRC/dstein.f
+++ b/SRC/dstein.f
@@ -297,7 +297,7 @@
          BLKSIZ = BN - B1 + 1
          IF( BLKSIZ.EQ.1 )
      $      GO TO 60
-         GPIND = B1
+         GPIND = J1
 *
 *        Compute reorthogonalization criterion and stopping criterion.
 *
@@ -370,9 +370,10 @@
 *
 *           Normalize and scale the righthand side vector Pb.
 *
+            JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
             SCL = BLKSIZ*ONENRM*MAX( EPS,
      $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
-     $            DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
+     $            ABS( WORK( INDRV1+JMAX ) )
             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
 *
 *           Solve the system LU = Pb.

diff --git a/SRC/sstein.f b/SRC/sstein.f
index 0c2ab02..0e3cd24 100644 (file)
--- a/SRC/sstein.f
+++ b/SRC/sstein.f
@@ -297,7 +297,7 @@
          BLKSIZ = BN - B1 + 1
          IF( BLKSIZ.EQ.1 )
      $      GO TO 60
-         GPIND = B1
+         GPIND = J1
 *
 *        Compute reorthogonalization criterion and stopping criterion.
 *
@@ -370,9 +370,10 @@
 *
 *           Normalize and scale the righthand side vector Pb.
 *
+            JMAX = ISAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
             SCL = BLKSIZ*ONENRM*MAX( EPS,
      $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
-     $            SASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
+     $            ABS( WORK( INDRV1+JMAX ) )
             CALL SSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
 *
 *           Solve the system LU = Pb.

diff --git a/SRC/zstein.f b/SRC/zstein.f
index 1f6a5fd..d134884 100644 (file)
--- a/SRC/zstein.f
+++ b/SRC/zstein.f
@@ -308,7 +308,7 @@
          BLKSIZ = BN - B1 + 1
          IF( BLKSIZ.EQ.1 )
      $      GO TO 60
-         GPIND = B1
+         GPIND = J1
 *
 *        Compute reorthogonalization criterion and stopping criterion.
 *
@@ -381,9 +381,10 @@
 *
 *           Normalize and scale the righthand side vector Pb.
 *
+            JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
             SCL = BLKSIZ*ONENRM*MAX( EPS,
      $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
-     $            DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
+     $            ABS( WORK( INDRV1+JMAX ) )
             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
 *
 *           Solve the system LU = Pb.