From fce95e275e488d173d089ca23ab2212b0f3ef957 Mon Sep 17 00:00:00 2001
From: Angelika Schwarz <angies@cs.umu.se>
Date: Wed, 18 Sep 2024 08:11:05 +0200
Subject: [PATCH] stedc scales only if norm is tiny or huge (and not
 unconditionally)

---
 SRC/dlaed0.f |  2 +-
 SRC/dstedc.f | 44 ++++++++++++++++++++++++++++++++------------
 2 files changed, 33 insertions(+), 13 deletions(-)

diff --git a/SRC/dlaed0.f b/SRC/dlaed0.f
index 368bf3e6c4..6818584b8c 100644
--- a/SRC/dlaed0.f
+++ b/SRC/dlaed0.f
@@ -74,7 +74,7 @@
 *>         On exit, its eigenvalues.
 *> \endverbatim
 *>
-*> \param[in] E
+*> \param[in,out] E
 *> \verbatim
 *>          E is DOUBLE PRECISION array, dimension (N-1)
 *>         The off-diagonal elements of the tridiagonal matrix.
diff --git a/SRC/dstedc.f b/SRC/dstedc.f
index 1e0b44d4f5..91e3fe69cd 100644
--- a/SRC/dstedc.f
+++ b/SRC/dstedc.f
@@ -201,9 +201,10 @@ SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LQUERY
-      INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
-     $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
-      DOUBLE PRECISION   EPS, ORGNRM, P, TINY
+      INTEGER            FINISH, I, ICOMPZ, II, ISCALE, J, K, LGN,
+     $                   LIWMIN, LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
+      DOUBLE PRECISION   EPS, EPS2, ORGNRM, P, SAFMAX, SAFMIN,
+     $                   SSFMAX, SSFMIN, TINY
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -339,8 +340,15 @@ SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
          ORGNRM = DLANST( 'M', N, D, E )
          IF( ORGNRM.EQ.ZERO )
      $      GO TO 50
+*
+*        Determine the unit roundoff and over/underflow thresholds.
 *
          EPS = DLAMCH( 'Epsilon' )
+         EPS2 = EPS**2
+         SAFMIN = DLAMCH( 'S' )
+         SAFMAX = ONE / SAFMIN
+         SSFMAX = SQRT( SAFMAX ) / ( ONE + TWO )
+         SSFMIN = SQRT( SAFMIN ) / EPS2
 *
          START = 1
 *
@@ -378,12 +386,20 @@ SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
 *              Scale.
 *
                ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
-               CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ),
-     $                      M,
-     $                      INFO )
-               CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1,
-     $                      E( START ),
-     $                      M-1, INFO )
+               ISCALE = 0
+               IF( ORGNRM.GT.SSFMAX ) THEN
+                  ISCALE = 1
+                  CALL DLASCL( 'G', 0, 0, ORGNRM, SSFMAX, M, 1,
+     $                         D( START ), M, INFO )
+                  CALL DLASCL( 'G', 0, 0, ORGNRM, SSFMAX, M-1, 1,
+     $                          E( START ), M-1, INFO )
+               ELSE IF( ORGNRM.LT.SSFMIN ) THEN
+                  ISCALE = 2
+                  CALL DLASCL( 'G', 0, 0, ORGNRM, SSFMIN, M, 1,
+     $                         D( START ), M, INFO )
+                  CALL DLASCL( 'G', 0, 0, ORGNRM, SSFMIN, M-1, 1,
+     $                          E( START ), M-1, INFO )
+               END IF
 *
                IF( ICOMPZ.EQ.1 ) THEN
                   STRTRW = 1
@@ -401,9 +417,13 @@ SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
 *
 *              Scale back.
 *
-               CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ),
-     $                      M,
-     $                      INFO )
+               IF( ISCALE.EQ.1 ) THEN
+                  CALL DLASCL( 'G', 0, 0, SSFMAX, ORGNRM, M, 1,
+     $                         D( START ), M, INFO )
+               ELSE IF( ISCALE.EQ.2 ) THEN
+                  CALL DLASCL( 'G', 0, 0, SSFMIN, ORGNRM, M, 1,
+     $                         D( START ), M, INFO )
+               END IF
 *
             ELSE
                IF( ICOMPZ.EQ.1 ) THEN