Make lapy and nrm2 consistent

Avoid the introduction of NaN by divisons inf / inf
chrwarm · Aug 8, 2021 · 906adae · 906adae
1 parent 3b26987
commit 906adae
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 11 deletions.
diff --git a/SRC/dlapy2.f b/SRC/dlapy2.f
@@ -78,13 +78,16 @@ DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
       PARAMETER          ( ONE = 1.0D0 )
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION   W, XABS, YABS, Z
+      DOUBLE PRECISION   W, XABS, YABS, Z, HUGEVAL
       LOGICAL            X_IS_NAN, Y_IS_NAN
 *     ..
 *     .. External Functions ..
       LOGICAL            DISNAN
       EXTERNAL           DISNAN
 *     ..
+*     .. External Subroutines ..
+      DOUBLE PRECISION   DLAMCH
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, MIN, SQRT
 *     ..
@@ -94,13 +97,14 @@ DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
       Y_IS_NAN = DISNAN( Y )
       IF ( X_IS_NAN ) DLAPY2 = X
       IF ( Y_IS_NAN ) DLAPY2 = Y
+      HUGEVAL = DLAMCH( 'Overflow' )
 *
       IF ( .NOT.( X_IS_NAN.OR.Y_IS_NAN ) ) THEN
          XABS = ABS( X )
          YABS = ABS( Y )
          W = MAX( XABS, YABS )
          Z = MIN( XABS, YABS )
-         IF( Z.EQ.ZERO ) THEN
+         IF( Z.EQ.ZERO .OR. W.GT.HUGEVAL ) THEN
             DLAPY2 = W
          ELSE
             DLAPY2 = W*SQRT( ONE+( Z / W )**2 )

diff --git a/SRC/dlapy3.f b/SRC/dlapy3.f
@@ -81,18 +81,22 @@ DOUBLE PRECISION FUNCTION DLAPY3( X, Y, Z )
       PARAMETER          ( ZERO = 0.0D0 )
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION   W, XABS, YABS, ZABS
+      DOUBLE PRECISION   W, XABS, YABS, ZABS, HUGEVAL
+*     ..
+*     .. External Subroutines ..
+      DOUBLE PRECISION   DLAMCH
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, SQRT
 *     ..
 *     .. Executable Statements ..
 *
+      HUGEVAL = DLAMCH( 'Overflow' )
       XABS = ABS( X )
       YABS = ABS( Y )
       ZABS = ABS( Z )
       W = MAX( XABS, YABS, ZABS )
-      IF( W.EQ.ZERO ) THEN
+      IF( W.EQ.ZERO .OR. W.GT.HUGEVAL ) THEN
 *     W can be zero for max(0,nan,0)
 *     adding all three entries together will make sure
 *     NaN will not disappear.

diff --git a/SRC/slapy2.f b/SRC/slapy2.f
@@ -78,32 +78,33 @@ REAL             FUNCTION SLAPY2( X, Y )
       PARAMETER          ( ONE = 1.0E0 )
 *     ..
 *     .. Local Scalars ..
-      REAL               W, XABS, YABS, Z
+      REAL               W, XABS, YABS, Z, HUGEVAL
       LOGICAL            X_IS_NAN, Y_IS_NAN
 *     ..
 *     .. External Functions ..
       LOGICAL            SISNAN
       EXTERNAL           SISNAN
 *     ..
+*     .. External Subroutines ..
+      REAL               SLAMCH
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, MIN, SQRT
 *     ..
 *     .. Executable Statements ..
-*
-*     ..
-*     .. Executable Statements ..
 *
       X_IS_NAN = SISNAN( X )
       Y_IS_NAN = SISNAN( Y )
       IF ( X_IS_NAN ) SLAPY2 = X
       IF ( Y_IS_NAN ) SLAPY2 = Y
+      HUGEVAL = SLAMCH( 'Overflow' )
 *
       IF ( .NOT.( X_IS_NAN.OR.Y_IS_NAN ) ) THEN
          XABS = ABS( X )
          YABS = ABS( Y )
          W = MAX( XABS, YABS )
          Z = MIN( XABS, YABS )
-         IF( Z.EQ.ZERO ) THEN
+         IF( Z.EQ.ZERO .OR. W.GT.HUGEVAL ) THEN
             SLAPY2 = W
          ELSE
             SLAPY2 = W*SQRT( ONE+( Z / W )**2 )

diff --git a/SRC/slapy3.f b/SRC/slapy3.f
@@ -81,18 +81,22 @@ REAL             FUNCTION SLAPY3( X, Y, Z )
       PARAMETER          ( ZERO = 0.0E0 )
 *     ..
 *     .. Local Scalars ..
-      REAL               W, XABS, YABS, ZABS
+      REAL               W, XABS, YABS, ZABS, HUGEVAL
+*     ..
+*     .. External Subroutines ..
+      REAL               SLAMCH
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, SQRT
 *     ..
 *     .. Executable Statements ..
 *
+      HUGEVAL = SLAMCH( 'Overflow' )
       XABS = ABS( X )
       YABS = ABS( Y )
       ZABS = ABS( Z )
       W = MAX( XABS, YABS, ZABS )
-      IF( W.EQ.ZERO ) THEN
+      IF( W.EQ.ZERO .OR. W.GT.HUGEVAL ) THEN
 *     W can be zero for max(0,nan,0)
 *     adding all three entries together will make sure
 *     NaN will not disappear.