Type and precision in Fortran

A reccurent problem that is unfortunately ignored by Fortran developers is floating point typing mechanism in Fortran. A lazy programer, or one who relies to much on her compiler, would probably declare and define her variable this way:

real(kind=8) :: x
x = 0.1

She would probably think that everything will be alright: since x is double precision, that would be the precision she will obtain in the end. It that really true? Let's test this assertion with this simple piece of code:

real(kind=8) :: x, y
x = 0.1; y = 0.1d0

print '("x = ", F20.18, ", y = ", F20.18)', x, y

Let's compile and run it; here is what we get:

x = 0.10000000149011612, y = 0.100000000000000001

y has the precision we expect, but x is actually in simple precision! Whole programs can be corrupted by this simple error, programers thinking they are in double precision whereas they are, in effect, in simple precision.

But what just happened? Why do we get this result? The compiler just followed the Fortran norm, stipulating that assignment are done from right to left. First, the right hand side of the assignment x = 0.1 is evaluated; 0.1 is simple precision. It will then be assigned to x (x is then equal to 0.1 in simple precision), and then, x being declared as a double precision, a retyping is done, promoting it to double precision. But the promoting is not enough to get a "real" double precision, since it was not the original precision of 0.1! x is double precision, but "completed" with garbage after the first 4 bytes.

There are different ways to circumvent the issue. First of all is to use compiler options to force floating point numbers in double precision (-fdefault-real-8 with gfortran, for example). This solution can be problematic; typically, double precision numbers are then promotted in quadruple! It is also possible to force double preicision floatting points in double precision with other compiler options (as -fdefault-double-8 with gfortran), but it is not always easy to use and understand.

Another solution would be to re-read the entire code and to harmonize all the definitions and assignment... It can be awfuly fastidious and it can be easy to make other mistakes.

The last solution is to use PCF, Precision Converter for Fortran, a Python script I wrote to harmonize floating point numbers typing in Fortran sources. Feed the script with the sources to convert, and with the precision you desire (d0, q0, _dp) and everything should be fine for you then!