--- rpl/lapack/lapack/dlapy3.f 2018/05/29 07:17:58 1.17 +++ rpl/lapack/lapack/dlapy3.f 2023/08/07 08:38:55 1.18 @@ -31,7 +31,7 @@ *> \verbatim *> *> DLAPY3 returns sqrt(x**2+y**2+z**2), taking care not to cause -*> unnecessary overflow. +*> unnecessary overflow and unnecessary underflow. *> \endverbatim * * Arguments: @@ -61,17 +61,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 -* *> \ingroup OTHERauxiliary * * ===================================================================== DOUBLE PRECISION FUNCTION DLAPY3( X, Y, Z ) * -* -- LAPACK auxiliary routine (version 3.7.0) -- +* -- LAPACK auxiliary routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 * * .. Scalar Arguments .. DOUBLE PRECISION X, Y, Z @@ -84,18 +81,22 @@ 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.