One of our code bases still uses double precision to define most floating point variables. We are about to change that to real(dp) but as everyone we’re ending up in a discussion on selecting the right definition for the real kind dp. Based on various posts and a discussion by Dr Fortran we have arrived at a shortlist of two/three options:
use stdlib, only: dp … focusing on interoperability with stdlib which internally uses selected_real_kind(15) to define dp
use iso_c_binding, only: dp => c_double … focusing on interoperability with C
use ieee_arithmetic, only: ieee_selected_real_kind and dp = ieee_selected_real_kind(15)
Since most hardware seems to be using IEEE arithmetic by default, the benefit of the third option seems to be limited. So, let’s focus on the first two.
For pure Fortran code, following the first option seems most logical: the precision is well defined (well at least concerning the precision, not yet concerning exponent or radix) and the code would be interoperable with all functionality of stdlib.
However, we interoperate at a number of places C code, e.g. for writing netCDF files, solving equations using PetSc, and exchanging many large arrays via API with codes in C and Python. Therefore, I’m thinking that second option would be the safer bet. However, that would add a series of conversions when calling stdlib and that doesn’t feel right.
So, how likely is it that we go wrong by just assuming that the stdlib dp kind and the iso_c_binding c_double kind are in fact equal? Stated differently: how likely is it that a stdlib real(dp) array is interpreted incorrectly when being accessed from C code as double? Surely there will be platforms where this assumption doesn’t hold, but what are those hardware platforms and how common are they?
I’m not sure why you think there would be conversions. (If anything, you should get compile-time errors if there are explicit interfaces for the library routines and the kind doesn’t match.)
Just going on what you wrote here, I’d go for option 2.
Yes, unfortunately there is no clear advice yet on this. I once (15 years ago) asked this very question at a Fortran mailinglist, and a long debate ensured where there was no consensus.
Have you and your colleagues found any such platform? Or, to think differently, do you know of any platform in use now where the 3 options you list in your original post are not the same?
I didn’t mean to suggest that there would be automatic conversion. Obviously, you’re right to state that compilers should notice the real kind mismatch and report an error. So, to avoid that we would have to write some wrapper routines for all stdlib routines that we use to explicitly convert all real precision values/arrays passed across the interface from c_double to dp and back … even though the kinds may in most cases be effectively the same. … The alternative would be to redefine the stdlib dp precision in our local clone of the source code to c_double but that feels like a very slippery path since we would deviate from the standard code base and may use the code in ways not intended or tested.
@FortranFan,
No we don’t have examples yet. My “surely” was more intended to reflect that anything that can go wrong, will go wrong … so, even though if I don’t know of the system, the wider community likely does. Once we know on what kind of platform this assumption does not hold, we can decide how relevant those are for us and either address the issue or ignore it with a statement that we explicitly don’t support those platforms.
Thank you for your reply. Yes, it’s my first question on the forum. Usually, I follow the advise of @Arjen, but this time he pointed to the forum.
I agree that if the real kinds are defined in one place, it’s easy to adjust and one of our codes has followed the Fortran Best Practices — Fortran90 1.0 documentation for a similar period of 15 years. However, with other code catching up and new colleagues joining the team, it’s time to reflect whether that approach is still the best. Indeed the selected_real_kind is now preferred, and following stdlib is a pragmatic choice in that context. However, all of that discussion doesn’t address the potential issues with language interoperability which I consider an critical part of our code.
According at the change log of GCC, REAL16 is now provided in ISO_FORTRAN_ENV.
(I do not know, to what extent there is support of actual half-precision types).
Has that name REAL16 been adopted by the standard committee, or are the GCC people anticipating future support in the stanard? And are there corresponding entries in the GCC version of iso_c_binding?
16.10.2.27 REAL16, REAL32, REAL64, and REAL128
The values of these default integer scalar named constants shall be those of the kind type parameters that specify a REAL type whose storage size expressed in bits is 16, 32, 64, and 128 respectively. If, for any of these constants, the processor supports more than one kind of that size, it is processor dependent which kind value is provided. If the processor supports no kind of a particular size, that constant shall be equal to −2 if the processor supports kinds of a larger size and −1 otherwise.
That’s the first time I see a max operator in a real kind definition. Neat approach to implement a fall back if real16 is not implemented … but what if real16 and real32 are both implemented. Which precision would x then be? It seems easy enough to write a program that tells which precision was selected, but does the standard state that real16 should be less or more than real32? If not, it’s compiler specific at best … although common sense bias might suggest that real32 will typically be larger than real16 and hence x would always be defined as real32 whether real16 is defined or not (assuming that the more common type real32 is defined), which would make the max construct utterly useless.
I would say that a library writer’s task is both simple and complicated. It is simple because “all” that is necessary is to write versions of the code that work with every element of the real_kinds(:) array for a given compiler. (I’m ignoring the complementary integer and logical kinds problems, and I’m assuming the complex kinds code will follow the real kinds codes.) Sometimes that will involve only declaring the right kind values in each specific library routine, while other times it might involve more significant changes to the underlying algorithms. However, the complicated part is that the fortran language provides no support at all for a programmer to do that. The real_kinds(:) array was introduced in f2003, so it has been some 20 years, and still this is a difficult task. In my opinion, this is something that should have been provided at the time the real_kinds(:) array was introduced in the language. It isn’t like this was some unanticipated problem that came out of nowhere, it is really THE central problem with writing library codes in a portable way that can be compiled by every compiler.