Consistent floating point arithmetic in fortran with different compilers on different computers?

floating-point, math

Solution

We'll turn to floating-point arithmetic below, but let me deal first with a misconception you have. Your use of list-directed formatting in the write statement means that the format chosen by the compiler to write out a variable is, well, the compiler's choice; it's not mandated by the language standard. Your use of the second `*` in `write(*,*)` tells the compiler to write out the value of a variable as it wishes. So what you have is not evidence that there is something wrong with the arithmetic, but that there are differences between `gfortran` and `ifort`. If I modify your write statement to

write(*,'(f21.18)') a

then my Intel Fortran program writes

0.900000000000000022

to the console.

SO is littered with questions arising from unfamiliarity with the details of floating-point arithmetic so I'm not going to write a treatise, just a few observations relevant to your question.

IEEE-754 64-bit floating-point numbers (which is probably what you get with the declaration `real(kind=8)`) only provide approximately 16 decimal digits of useful information. Actually, since they're binary and there is no easy correspondence between the number of digits in one base and in another, it's actually 15.95 decimal digits and many users round this down and never look at anything after the 15th significant figure in a floating-point number's decimal representation. So both `ifort` and `gfortran` are misleading you with their trailing `2`s.

IEEE-754 defines not only the formats for f-p numbers, but also some rules for rounding and for arithmetic operations. A carefully-written program which uses only those arithmetic operations (I think square root is also specified) and which takes care about rounding modes and rounding operations should produce the same results on two different processors. Of course, not many useful numerical programs confine themselves to the basic arithmetic operations.

Since the 2003 standard Fortran has included an intrinsic module called `ieee_arithmetic` which gives the programmer direct access to the underlying IEEE-754 capabilities of the hardware on which it is running -- but note that it does not require that the hardware have any such capabilities. Using `ieee_arithmetic` and another intrinsic module called `ieee_exceptions` you should, if your hardware provides the necessary support, be able to write programs which compile under both `gfortran` and `ifort` and which produce, when executed, the same results to the last bit of every f-p number.

You'll also need to familiarise yourself with the optimisation options and numerical arithmetic options of the compiler(s) you're using. Most compilers have an option whose meaning is `sacrifice IEEE compliance for speed` (for `ifort` I think it's `fp-model`). In general adherence to IEEE-754 for operations will slow your programs.

Problem

Original source