3

I'm using ieee_arithmetic with Fortran on a Linux machine that runs gfortran version 5.4.0.

I'm getting an error of division by zero when trying to initialize values for Inf and NaN.

There doesn't seem to be an issue with ieee_arithmetic because elsewhere in the file I can successfully call ieee_is_finite() with no issues.

I thought that ieee_arithmetic allowed division by zero to be used for these specific cases, but I must be missing something. Below is a sample of code:

module rcrlib_gnu
    use, intrinsic :: ieee_arithmetic ! requires gfortran version 5.0 or higher
    implicit none
    integer, parameter :: SP=kind(1.0), DP=selected_real_kind(9,99)
    integer, parameter :: stderr=0
    public SP, DP, is_finite, stderr, initialize

contains

subroutine initialize(infty,nan)
    real(kind=DP), intent(out) :: infty, nan
    infty = 1.0_dp/0.0_dp ! huge(1.0_dp)
    nan = 0.0_dp/0.0_dp
end subroutine initialize

elemental function is_finite(x)
    real(kind=DP), intent(in) :: x
    logical :: is_finite
    is_finite = ieee_is_finite(x) ! This call requires "ieee_arithmetic"
end function is_finite 

end module rcrlib_gnu

It seems I'm missing something basic, so I would appreciate any help.

To reproduce the error, save the above code snippet as rcrlib_gnu_example.f90 and then execute the following line: gfortran -o rcr rcrlib_gnu_example.f90

The resulting error output is

rcrlib_gnu_example.f90:12:18:

     infty = 1.0_dp/0.0_dp ! huge(1.0_dp)
                  1
Error: Division by zero at (1)
rcrlib_gnu_example.f90:13:16:

     nan = 0.0_dp/0.0_dp
                1
Error: Division by zero at (1)
Tyler R.
  • 461
  • 6
  • 15
  • Similar posts [here](https://stackoverflow.com/questions/30747019/ieee-arithmetic-intrinsic-module-in-gfortran) and [here](https://stackoverflow.com/questions/17389958/is-there-a-standard-way-to-check-for-infinite-and-nan-in-fortran-90-95) have asked about this issue, but are not a direct answer to my query. [This page](http://northstar-www.dartmouth.edu/doc/solaris-forte/manuals/fortran/prog_guide/6_floating.html) suggests that division by zero yields `+Inf` or `NaN` in `ieee_arithmetic`. – Tyler R. Aug 17 '17 at 20:23
  • 1
    It is not Fortran 90. It is Fortran 2003 or later. There was no ieee_arithmetic in the old Fortran 90. – Vladimir F Героям слава Aug 17 '17 at 20:58
  • Show us the code that causes the error or the unexpected behaviour. See [mcve]. We must be able to *reproduce* the problem. – Vladimir F Героям слава Aug 17 '17 at 20:59
  • Note that there are other ways to get `Inf` and `NaN` than division if that's what you want. Converting a decimal representation larger than the max finite float + a half ULP yields `Inf`, and once you have `Inf`, you can get `NaN` with `Inf - Inf`. (Not sure if that's useful, but I thought I should point it out in case.) – Pascal Cuoq Aug 17 '17 at 21:00
  • Loosely the IEEE module exceptions are about run-time; evaluation of a constant expression may well be compile-time. I seem to recall a question relating to this but I can't find it just at the moment. – francescalus Aug 17 '17 at 21:04
  • @francescalus You can see in my comment that `huge(1.0_dp)` is a suitable substitute for dividing by 0 to create `+Inf`. I have yet to find a workable solution for `NaN`, which is the main holdup for me. – Tyler R. Aug 17 '17 at 21:40
  • @PascalCuoq this is useful (see the comment in my code snippet), though I'm not sure if `huge(1.0_dp)-huge(1.0_dp)` will give `NaN`. Can you explain what you mean by "+ a half ULP", and how I would code that, compared with `huge(1.0_dp)`? – Tyler R. Aug 17 '17 at 21:41
  • @VladimirF I have updated the question to include information on reproducing the error, as well as the error message. Apologies for omitting this information initially. – Tyler R. Aug 17 '17 at 21:47
  • If you just want to get a NaN then you don't need to worry about doing division by zero. There are several ways, including using `ieee_value`. See also [this question](https://stackoverflow.com/q/31971836/3157076). Or are you really after an explanation about this division by zero aspect? – francescalus Aug 17 '17 at 21:49
  • @TylerR. Sorry, I am not familiar enough with FORTRAN to answer that, but I am somewhat familiar with IEEE 754. I should have said “Converting a decimal representation **significantly** larger than the max finite float” for simplicity. I just wanted to save you the surprise of using a decimal representation just above that of the max finite float. That one would be converted to the max finite float. Answering your other question, `Inf - Inf` produces `NaN` under IEEE 754 rules just like `0 / 0` does, but likely with the advantage of not being treated as a special case by your compiler. – Pascal Cuoq Aug 17 '17 at 22:35
  • Indeed, now I see the question that francescalus points to https://stackoverflow.com/q/31971836/3157076 is very closely related. Even a possible duplicate. – Vladimir F Героям слава Aug 18 '17 at 06:52
  • @VladimirF yes, thank you. I had seen that one and tried to implement it but was unsuccessful. I'll give it another try. – Tyler R. Aug 18 '17 at 21:13
  • @PascalCuoq Thank you for this; I'll give it a try. – Tyler R. Aug 18 '17 at 21:14
  • @francescalus Thank you for the pointer. I am mostly looking for a workable solution, but also curious why `0/0` doesn't work even though `ieee_arithmetic` seems to work. – Tyler R. Aug 18 '17 at 21:15
  • Because it is not alowed to delete by zero during compilation. It is all in the other question. Please read it, read also the explanation there. At least currently it is not allowed. The important point is that the operation is evaluated already during compilation. – Vladimir F Героям слава Aug 18 '17 at 21:17
  • 1
    Possible duplicate of [Having parameter (constant) variable with NaN value in FORTRAN](https://stackoverflow.com/questions/31971836/having-parameter-constant-variable-with-nan-value-in-fortran) – francescalus Aug 18 '17 at 23:20

1 Answers1

0

Thanks to Pascal Cuoq, I solved the problem.

The version of the initialize subroutine that compiles is below:

subroutine initialize(infty,nan)
    real(kind=DP), intent(out) :: infty, nan
    infty = huge(1.0_dp)+100
    nan = infty-infty
end subroutine initialize

So basically set infinity to be the largest floating point number plus 100, then set NaN to be the difference between infinity and itself.

Thanks, all, for your quick responses and patience with my lack of FORTRAN experience.

Tyler R.
  • 461
  • 6
  • 15
  • As you have access to the IEEE modules, why is this approach better than using `ieee_value()`? In particular, what happens when assigning to `infty` when the halting mode for overflow is set? – francescalus Aug 18 '17 at 23:19
  • @francescalus it may not be a better method, and I am unsure about your halting mode concern. The solution above appears to have worked (at least for compilation on my machine), but your point is well taken. – Tyler R. Aug 20 '17 at 01:33