introduce inf_values to fix the failure with gfortran-13 -Ofast #23

inf.f90

@@ -1,12 +1,22 @@
 module inf_mod
-! infnan_mod together with inf_mod provides functions for checking Inf/NaN. They aim to work even when 
+! infnan_mod together with inf_mod provides functions for checking Inf/NaN. They aim to work even when
 ! compilers are invoked with aggressive optimization flags, e.g., `gfortran -Ofast`.
-! See infnan.f90 for more comments. 
-
+! See infnan.f90 for more comments.
+use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_negative_inf, ieee_positive_inf
+use, non_intrinsic :: consts_mod, only : SP, DP
+use, non_intrinsic :: inf_values_mod, only : posinf, neginf
 implicit none
 private
 public :: is_finite, is_inf, is_posinf, is_neginf
 
+!interface posinf
+!    module procedure posinf_sp, posinf_dp
+!end interface posinf
+
+!interface neginf
+!    module procedure neginf_sp, neginf_dp
+!end interface neginf
+
 interface is_finite
     module procedure is_finite_sp, is_finite_dp
 end interface is_finite
@@ -27,20 +37,58 @@ module inf_mod
 contains
 
 
+!pure function posinf_sp(x) result(y)
+!use, non_intrinsic :: consts_mod, only : SP
+!implicit none
+!real(SP), intent(in) :: x
+!real(QP) :: y
+!y = ieee_value(x, ieee_positive_inf)
+!end function posinf_sp
+
+!pure function posinf_dp(x) result(y)
+!use, non_intrinsic :: consts_mod, only : DP
+!implicit none
+!real(DP), intent(in) :: x
+!real(QP) :: y
+!y = ieee_value(x, ieee_positive_inf)
+!end function posinf_dp
+
+
+!pure function neginf_sp(x) result(y)
+!use, non_intrinsic :: consts_mod, only : SP
+!implicit none
+!real(SP), intent(in) :: x
+!real(QP) :: y
+!y = ieee_value(x, ieee_negative_inf)
+!end function neginf_sp
+
+!pure function neginf_dp(x) result(y)
+!use, non_intrinsic :: consts_mod, only : DP
+!implicit none
+!real(DP), intent(in) :: x
+!real(QP) :: y
+!y = ieee_value(x, ieee_negative_inf)
+!end function neginf_dp
+
+
 elemental pure function is_finite_sp(x) result(y)
 use consts_mod, only : SP
 implicit none
 real(SP), intent(in) :: x
 logical :: y
-y = (x <= huge(x) .and. x >= -huge(x))
+!y = (x <= huge(x) .and. x >= -huge(x))
+y = (x < posinf(x) .and. x > neginf(x))
+!y = abs(x) < posinf(x)
 end function is_finite_sp
 
 elemental pure function is_finite_dp(x) result(y)
 use consts_mod, only : DP
 implicit none
 real(DP), intent(in) :: x
 logical :: y
-y = (x <= huge(x) .and. x >= -huge(x))
+!y = (x <= huge(x) .and. x >= -huge(x))
+y = (x < posinf(x) .and. x > neginf(x))
+!y = abs(x) < posinf(x)
 end function is_finite_dp
 
 
@@ -49,15 +97,17 @@ elemental pure function is_inf_sp(x) result(y)
 implicit none
 real(SP), intent(in) :: x
 logical :: y
-y = (abs(x) > huge(x))
+!y = (abs(x) > huge(x))
+y = (abs(x) >= posinf(x))
 end function is_inf_sp
 
 elemental pure function is_inf_dp(x) result(y)
 use consts_mod, only : DP
 implicit none
 real(DP), intent(in) :: x
 logical :: y
-y = (abs(x) > huge(x))
+!y = (abs(x) > huge(x))
+y = (abs(x) >= posinf(x))
 end function is_inf_dp
 
 
@@ -66,15 +116,17 @@ elemental pure function is_posinf_sp(x) result(y)
 implicit none
 real(SP), intent(in) :: x
 logical :: y
-y = (abs(x) > huge(x)) .and. (x > 0)
+!y = (abs(x) > huge(x)) .and. (x > 0)
+y = ((abs(x) >= posinf(x)) .and. (x > 0))
 end function is_posinf_sp
 
 elemental pure function is_posinf_dp(x) result(y)
 use consts_mod, only : DP
 implicit none
 real(DP), intent(in) :: x
 logical :: y
-y = (abs(x) > huge(x)) .and. (x > 0)
+!y = (abs(x) > huge(x)) .and. (x > 0)
+y = ((abs(x) >= posinf(x)) .and. (x > 0))
 end function is_posinf_dp
 
 
@@ -83,15 +135,19 @@ elemental pure function is_neginf_sp(x) result(y)
 implicit none
 real(SP), intent(in) :: x
 logical :: y
-y = (abs(x) > huge(x)) .and. (x < 0)
+!y = (abs(x) > huge(x)) .and. (x < 0)
+y = ((abs(x) >= posinf(x)) .and. (x < 0))
+!y = (x <= neginf(x))
 end function is_neginf_sp
 
 elemental pure function is_neginf_dp(x) result(y)
 use consts_mod, only : DP
 implicit none
 real(DP), intent(in) :: x
 logical :: y
-y = (abs(x) > huge(x)) .and. (x < 0)
+!y = (abs(x) > huge(x)) .and. (x < 0)
+y = ((abs(x) >= posinf(x)) .and. (x < 0))
+!y = (x <= neginf(x))
 end function is_neginf_dp
 
 end module inf_mod

inf_values.f90

@@ -0,0 +1,52 @@
+module inf_values_mod
+
+use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_negative_inf, ieee_positive_inf
+use, non_intrinsic :: consts_mod, only : SP, DP
+implicit none
+private
+public :: posinf, neginf
+interface posinf
+    module procedure posinf_sp, posinf_dp
+end interface posinf
+
+interface neginf
+    module procedure neginf_sp, neginf_dp
+end interface neginf
+
+contains
+
+pure function posinf_sp(x) result(y)
+use, non_intrinsic :: consts_mod, only : SP
+implicit none
+real(SP), intent(in) :: x
+real(SP) :: y
+y = ieee_value(x, ieee_positive_inf)
+end function posinf_sp
+
+pure function posinf_dp(x) result(y)
+use, non_intrinsic :: consts_mod, only : DP
+implicit none
+real(DP), intent(in) :: x
+real(DP) :: y
+y = ieee_value(x, ieee_positive_inf)
+end function posinf_dp
+
+
+pure function neginf_sp(x) result(y)
+use, non_intrinsic :: consts_mod, only : SP
+implicit none
+real(SP), intent(in) :: x
+real(SP) :: y
+y = ieee_value(x, ieee_negative_inf)
+end function neginf_sp
+
+pure function neginf_dp(x) result(y)
+use, non_intrinsic :: consts_mod, only : DP
+implicit none
+real(DP), intent(in) :: x
+real(DP) :: y
+y = ieee_value(x, ieee_negative_inf)
+end function neginf_dp
+
+
+end module inf_values_mod

infnan.f90

@@ -25,7 +25,8 @@ module infnan_mod
 ! 6. Even though the functions involve invocation of ABS and HUGE, their performance (in terms of
 ! CPU time) turns out comparable to or even better than the functions in `ieee_arithmetic`.
 
-use inf_mod, only : is_finite, is_inf, is_posinf, is_neginf
+use, non_intrinsic :: inf_values_mod, only : posinf, neginf
+use, non_intrinsic :: inf_mod, only : is_finite, is_inf, is_posinf, is_neginf
 implicit none
 private
 public :: is_finite, is_inf, is_posinf, is_neginf, is_nan
@@ -45,7 +46,8 @@ elemental pure function is_nan_sp(x) result(y)
 real(SP), intent(in) :: x
 logical :: y
 !y = (.not. (x <= huge(x) .and. x >= -huge(x))) .and. (.not. abs(x) > huge(x))  ! Does not always work
-y = (.not. is_finite(x)) .and. (.not. is_inf(x))
+y = (.not. (x < posinf(x) .and. x >= neginf(x))) .and. (.not. abs(x) >= posinf(x))  ! Does not always work
+y = y .or. ((.not. is_finite(x)) .and. (.not. is_inf(x)))
 end function is_nan_sp
 
 elemental pure function is_nan_dp(x) result(y)
@@ -55,7 +57,8 @@ elemental pure function is_nan_dp(x) result(y)
 real(DP), intent(in) :: x
 logical :: y
 !y = (.not. (x <= huge(x) .and. x >= -huge(x))) .and. (.not. abs(x) > huge(x))  ! Does not always work
-y = (.not. is_finite(x)) .and. (.not. is_inf(x))
+y = (.not. (x < posinf(x) .and. x >= neginf(x))) .and. (.not. abs(x) >= posinf(x))  ! Does not always work
+y = y .or. ((.not. is_finite(x)) .and. (.not. is_inf(x)))
 end function is_nan_dp
 
 end module infnan_mod