I've been trying to find a standards-compliant way to check for Infinite and NaN values in Fortran 90/95 but it proved harder than I thought.
ieee_arithmetic
module in Fortran 2003 with the ieee_is_nan()
and ieee_is_finite()
intrinsic functions. However it's not supported by all the compilers (notably gfortran as of version 4.9).Defining infinity and NaN at the beginning like pinf = 1. / 0
and nan = 0. / 0
seems hackish to me and IMHO can raise some building problems - for example if some compilers check this in compile time one would have to provide a special flag.
Is there a way I can implement in standard Fortran 90/95?
function isinf(x) ! Returns .true. if x is infinity, .false. otherwise ... end function isinf
and isnan()
?
Not a Number (NaN) is represented by the largest value that the exponent can assume (all ones), and a nonzero fraction. Normalized REAL and DOUBLE PRECISION numbers have an implicit leading bit that provides one more bit of precision than is stored in memory.
The simple way without using the ieee_arithmatic
is to do the following.
Infinity: Define your variable infinity = HUGE(dbl_prec_var)
(or, if you have it, a quad precision variable). Then you can simply check to see if your variable is infinity by if(my_var > infinity)
.
NAN: This is even easier. By definition, NAN is not equal to anything, even itself. Simply compare the variable to itself: if(my_var /= my_var)
.
I don't have enough rep to comment so I'll "answer" regarding Rick Thompson's suggestion for testing infinity.
if (A-1 .eq. A)
This will also be true if A is a very large floating point number, and 1
is below the precision of A.
A simple test:
subroutine test_inf_1(A) real, intent(in) :: A print*, "Test (A-1 == A)" if (A-1 .eq. A) then print*, " INFINITY!!!" else print*, " NOT infinite" endif end subroutine subroutine test_inf_2(A) real, intent(in) :: A print*, "Test (A > HUGE(A))" if (A > HUGE(A)) then print*, " INFINITY!!!" else print*, " NOT infinite" endif end subroutine program test real :: A,B A=10 print*, "A = ",A call test_inf_1(A) call test_inf_2(A) print*, "" A=1e20 print*, "A = ",A call test_inf_1(A) call test_inf_2(A) print*, "" B=0.0 ! B is necessary to trick gfortran into compiling this A=1/B print*, "A = ",A call test_inf_1(A) call test_inf_2(A) print*, "" end program test
outputs:
A = 10.0000000 Test (A-1 == A) NOT infinite Test (A > HUGE(A)) NOT infinite A = 1.00000002E+20 Test (A-1 == A) INFINITY!!! Test (A > HUGE(A)) NOT infinite A = Infinity Test (A-1 == A) INFINITY!!! Test (A > HUGE(A)) INFINITY!!!
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With