Confusion about kinds in FORTRAN

Question

I have been in the process of writing a FORTRAN code for numerical simulations of an applied physics problem for more than two years and I've tried to follow the conventions described in Fortran Best Practices.

More specifically, I defined a parameter as

integer, parameter:: dp=kind(0.d0)

and then used it for all doubles in my code.

However, I found out (on this forum) that using KIND parameters not necessarily gives you the same precision if you compile your code using other compilers. In this question, I read that a possible solution is using the SELECTED_REAL_KIND and SELECTED_INT_KIND, which follow some convention, as far as I understand.

Later on, though, I found out about the ISO_FORTRAN_ENV module which defines the REAL32, REAL64 and REAL128 KIND parameters.

I guess that these are indeed portable and, since they belong to the FORTRAN 2008 standard (though supported by GNU), I guess that I should use these?

Therefore, I would greatly appreciate if someone with more knowledge and experience clear up the confusion.

Also, I have a follow-up question about using these KINDs, in HDF5. I was using H5T_NATIVE_DOUBLE and it was indeed working fine (as far as I know) However, in this document it is stated that this is now an obsolete feature and should not be used. In stead, they provide a function

INTEGER(HID_T) FUNCTION h5kind_to_type(kind, flag) RESULT(h5_type) . 

When I use it, and print out the exact numerical value of the HID_T integer corresponding to REAL64 gives me 50331972, whereas H5T_NATIVE_DOUBLE gives me 50331963, which is different.

If I then try to use the value calculated by H5kind_to_type, the HDF5 library runs just as fine and, using XDMF, I can plot the output in VisIt or Paraview without modifying the accompanying .xmf file.

So my second question would be (again): Is this correct usage?

Answer 1

The type double precision and the corresponding kind kind(1.d0) are perfectly well defined by the standard. But they are also not exactly fixed. Indeed there were many computers in history and use different kind of native formats for their floating point numbers and the standard must allow this!

So, a double precision is a kind of real which has a higher enough precision than the default real. The default real is also not fixed, it must correspond to what the computers can use.

Now today we have the standard for floating point numbers IEEE_754 which defines IEEE single (binary32) and IEEE double types (binary64) and some others. If the computer hardware implements this standard, as almost all computers younger than 20 years do, it is very likely that the compiler chooses these two as the real and double precision.

The Fortran 2008 standard brings the two kind constants real32 and real64 (and others). They enable you to request the real kinds which have storage size of 32 and 64 bits. It is not guaranteed it will be the IEEE types, but it is almost certain on modern computers.

To request the IEEE types (if they are available) use the intrinsic function ieee_selected_real_kind() from module ieee_arithmetic.

The IEEE types are the same an all computers (excluding endianness!), but the compiler is not required to support them, because you may have a computer which does not support these in hardware. This is only a theoretical possibility, all modern computers support them.

Now to your HDF constants, these are apparently just some indexes to some table, it does not matter if they are different or not, the important is whether they mean the same and in your case they do.

As I wrote above, it is extremely likely that on a computer which supports IEEE 754 the double precision will be identical to IEEE double. It may not be, if you use some compiler options which change this behaviour. There are compiler options which promote the default real to double and hey may also promote double precision to quad precision (128 bit) to preserve the standard semantics which requires the double precision to have more precision and storage size.

Conclusion: You can use both, or any other way to choose your kind constants (you can also use iso_c_binding's c_float and c_double), but you should be aware why are those ways different and what do they actually mean.