--- rpl/lapack/lapack/dlapy2.f	2018/05/29 07:17:58	1.18
+++ rpl/lapack/lapack/dlapy2.f	2023/08/07 08:38:55	1.19
@@ -31,7 +31,7 @@
 *> \verbatim
 *>
 *> DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary
-*> overflow.
+*> overflow and unnecessary underflow.
 *> \endverbatim
 *
 *  Arguments:
@@ -56,17 +56,14 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \date June 2017
-*
 *> \ingroup OTHERauxiliary
 *
 *  =====================================================================
       DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
 *
-*  -- LAPACK auxiliary routine (version 3.7.1) --
+*  -- LAPACK auxiliary routine --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     June 2017
 *
 *     .. Scalar Arguments ..
       DOUBLE PRECISION   X, Y
@@ -81,13 +78,16 @@
       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
 *     ..
@@ -97,13 +97,14 @@
       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 )