diff --git a/SRC/dlaed4.f b/SRC/dlaed4.f
index 2aac7145f..4e0fe3751 100644
--- a/SRC/dlaed4.f
+++ b/SRC/dlaed4.f
@@ -168,8 +168,9 @@
       LOGICAL            ORGATI, SWTCH, SWTCH3
       INTEGER            II, IIM1, IIP1, IP1, ITER, J, NITER
       DOUBLE PRECISION   A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
-     $                   EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
-     $                   RHOINV, TAU, TEMP, TEMP1, W
+     $                   EPS, ERRETM, ETA, ETA1, ETA2, MIDPT, PHI,
+     $                   PREW, PSI, RHOINV, TAU, TEMP, TEMP1, W
+
 *     ..
 *     .. Local Arrays ..
       DOUBLE PRECISION   ZZ( 3 )
@@ -182,7 +183,7 @@
       EXTERNAL           DLAED5, DLAED6
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          ABS, MAX, MIN, SQRT
+      INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
 *     ..
 *     .. Executable Statements ..
 *
@@ -349,10 +350,23 @@
      $      ETA = -W / ( DPSI+DPHI )
          TEMP = TAU + ETA
          IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
+            ETA1 = -W / ( DPSI+DPHI )
+            TEMP = TAU + ETA1
             IF( W.LT.ZERO ) THEN
-               ETA = ( DLTUB-TAU ) / TWO
+               ETA2 = ( DLTUB-TAU ) / TWO
             ELSE
-               ETA = ( DLTLB-TAU ) / TWO
+               ETA2 = ( DLTLB-TAU ) / TWO
+            END IF
+            IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
+*
+*                 If Newton step is within boundaries,
+*                 use the geometric mean of the distance
+*                 and keep the direction (sign).
+*
+               ETA = SIGN(ONE, ETA1) *
+     $               SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
+            ELSE
+               ETA = ETA2
             END IF
          END IF
          DO 50 J = 1, N
@@ -848,10 +862,23 @@
      $         ETA = -W / DW
             TEMP = TAU + ETA
             IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
+               ETA1 = -W / DW
+               TEMP = TAU + ETA1
                IF( W.LT.ZERO ) THEN
-                  ETA = ( DLTUB-TAU ) / TWO
+                  ETA2 = ( DLTUB-TAU ) / TWO
                ELSE
-                  ETA = ( DLTLB-TAU ) / TWO
+                  ETA2 = ( DLTLB-TAU ) / TWO
+               END IF
+               IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
+*
+*                 If Newton step is within boundaries,
+*                 use the geometric mean of the distance
+*                 and keep the direction (sign).
+*
+                  ETA = SIGN( ONE, ETA1 ) *
+     $                  SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
+               ELSE
+                  ETA = ETA2
                END IF
             END IF
 *
diff --git a/SRC/slaed4.f b/SRC/slaed4.f
index 3f431d218..8fb049e59 100644
--- a/SRC/slaed4.f
+++ b/SRC/slaed4.f
@@ -168,8 +168,8 @@
       LOGICAL            ORGATI, SWTCH, SWTCH3
       INTEGER            II, IIM1, IIP1, IP1, ITER, J, NITER
       REAL               A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
-     $                   EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
-     $                   RHOINV, TAU, TEMP, TEMP1, W
+     $                   EPS, ERRETM, ETA, ETA1, ETA2, MIDPT, PHI,
+     $                   PREW, PSI, RHOINV, TAU, TEMP, TEMP1, W
 *     ..
 *     .. Local Arrays ..
       REAL               ZZ( 3 )
@@ -182,7 +182,7 @@
       EXTERNAL           SLAED5, SLAED6
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          ABS, MAX, MIN, SQRT
+      INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
 *     ..
 *     .. Executable Statements ..
 *
@@ -349,10 +349,23 @@
      $      ETA = -W / ( DPSI+DPHI )
          TEMP = TAU + ETA
          IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
+            ETA1 = -W / ( DPSI+DPHI )
+            TEMP = TAU + ETA1
             IF( W.LT.ZERO ) THEN
-               ETA = ( DLTUB-TAU ) / TWO
+               ETA2 = ( DLTUB-TAU ) / TWO
             ELSE
-               ETA = ( DLTLB-TAU ) / TWO
+               ETA2 = ( DLTLB-TAU ) / TWO
+            END IF
+            IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
+*
+*                 If Newton step is within boundaries,
+*                 use the geometric mean of the distance
+*                 and keep the direction (sign).
+*
+               ETA = SIGN(ONE, ETA1) *
+     $               SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
+            ELSE
+               ETA = ETA2
             END IF
          END IF
          DO 50 J = 1, N
@@ -848,10 +861,23 @@
      $         ETA = -W / DW
             TEMP = TAU + ETA
             IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
+               ETA1 = -W / DW
+               TEMP = TAU + ETA1
                IF( W.LT.ZERO ) THEN
-                  ETA = ( DLTUB-TAU ) / TWO
+                  ETA2 = ( DLTUB-TAU ) / TWO
                ELSE
-                  ETA = ( DLTLB-TAU ) / TWO
+                  ETA2 = ( DLTLB-TAU ) / TWO
+               END IF
+               IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
+*
+*                 If Newton step is within boundaries,
+*                 use the geometric mean of the distance
+*                 and keep the direction (sign).
+*
+                  ETA = SIGN( ONE, ETA1 ) *
+     $                  SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
+               ELSE
+                  ETA = ETA2
                END IF
             END IF
 *