Hi,
there is several ways to declare a double precision real in Fortran:
real*8 :: w
real(8) :: x
real(kind(0d0)) :: y
double precision :: z
My habits are to use double precision
. Is it a good practice ?
Hi,
there is several ways to declare a double precision real in Fortran:
real*8 :: w
real(8) :: x
real(kind(0d0)) :: y
double precision :: z
My habits are to use double precision
. Is it a good practice ?
The following is more or less, my go-to module that I use it everywhere:
module mod_types
use iso_fortran_env, only: int8, int16, int32, int64, real32, real64, &
input_unit, output_unit, error_unit
implicit none
integer, parameter :: sp = real32
integer, parameter :: dp = real64
integer, parameter :: stdin = input_unit
integer, parameter :: stdout = output_unit
integer, parameter :: stderr = error_unit
end module
An example use is:
module foo
use mod_types, only: wp=>dp
implicit none
real(wp)::bar
end module
wp
stands for working precision and if for some reason I want to change precision, I just change this wp=>sp
so I do not need to change manually all real
variables. I believe something in that line is considered the â€śstandard wayâ€ť now.
I also use the constants defined in iso_fortran_env
, though I have recently learned (from this discussion in stdlib
) that there are still concerns over the portability of these constants since they only guarantee a byte width but nothing about precision.
In his post â€śIt takes all kindsâ€ť, Steve Lionel recommends use of the intrinsic function selected_real_kind
above all other methods.
Itâ€™s also worth mentioning the constants provided by iso_c_binding
which are required for c interoperability:
use iso_c_binding, only: sp=>c_float, dp=>c_double
Note that the real*8
syntax is an extension and not standard. Also with real(8)
you are assuming that â€śreal kind 8â€ť is an 8-byte real (which conceivably isnâ€™t even double precision), but thereâ€™s nothing in the standard that says kind values correspond to the byte size. In fact with the NAG compiler this isnâ€™t the case. They number their kinds consecutively from 1 (though there is a compiler flag you can use to change that).
If one truly knows what decimal precision and exponent range they need for their problem, they can use selected_real_kind
and let the system decide the kind (and hope it has one). But if youâ€™re like me, all I really know is that the hardware supports 32-bit and 64-bit floating point and that I want to use 64-bit. So in the past Iâ€™d figure out a precision and range that would ensure that selected_real_kind
gave me the 64-bit kind. Thankfully that silliness isnâ€™t needed anymore with the constants from iso_fortran_env
. I hate to disagree with Steve, but if you know you want a specific fp size, just use the constant and donâ€™t play games with selected_real_kind
.
You can disagree with me, thatâ€™s fine. All I wanted to say in that post was that use of real32, real64 etc. was little better than real(4) and real(8). If youâ€™re aware that all youâ€™re doing is asking for a physical storage size (and no other properties), and document that, then OK. But I would be happier to see a module that defines constants using selected_real_kind()
- you need do this only once.
Of course, if youâ€™re interfacing with C, then use the kinds from ISO_C_BINDING
.
P.S. Thank you for reading my blog posts!
Thank you all for your valuable answers (and â€śIt Takes All KINDsâ€ť was so interesting and pleasant to read). This is a very interesting discussion, I am learning or re-learning many things.
Perhaps using iso_fortran_env
or selected_real_kind
depends on the kind of computer scientist you are:
iso_fortran_env
can probably do the job. That kind of machine is relatively stable in time (concerning the reals in the CPU). Although, if GPU acceleration is democratized things can become different. And in 20 years, I donâ€™t know what kind of CPU we will use.selected_real_kind
can surely be very useful. Moreover, you probably work in a bigger team and your code will perhaps still used by other people in 20 years, so keeping the memory of the needed precision will be helpful.Now, what I am interested in is getting the best compromise between precision and speed. On my PC, I can use 32 bits real, 64 bits, 80 bits, 128 bits:
My question would now be: is it possible with Fortran to know the biggest â€śhardwareâ€ť reals offered by the CPU ?
Another question: I thought 80 bits were hard wired, no ? But in that case why the computing time is 50 % higher than with 64 bits reals ?
80-bit floating is only possible with the old x87 instructions. You canâ€™t take advantage of any of the vector SSE or AVX instructions that are generally much faster.
Since you mention speed, it should be pointed out that there are performance advantages to using 32bit, even though the fundamental operations (+ / *) are computed in the same time as 64bit:
Functions implemented with polynomial series will be quicker since fewer terms are needed for the reduced precision;
If your program can be vectorized, then throughput will be higher with 32bit than 64bit since more values can be placed in the vector registers;
Many intensive programs are memory-bound these days, so if your program operates on a large-enough dataset (i.e. one that doesnâ€™t fit in the cache), and it involves a lot of memory accesses per flop, then you will benefit from ~2x speedup when using 32bit simply by being able to load more values into the processor.
Obviously, these arenâ€™t options if you need double precision (though there are so-called mixed-precision solvers that attempt to exploit this).
@sblionel
Thank you for you explanation. I suspected it was that, but I was surprised that those 80 bits reals works also with gfortran -mfpmath=sse
. I guess that even with this flag the compiler choose to use x87 when needed.
@lkedward
Thank you very much for you answer, you are right, I was convinced there was no difference, based on arithmetical measurements.
When computing 2e9 cos(), the 32 bits version is 2.66x times faster than the 64 bits.
Finally, itâ€™s another good reason to adopt a practice like writing real(wp)::bar
. By just changing the wp
value I could study the difference in computing time.
Hi everyone,
About the nag compiler, youâ€™re right, kind=8 is not working with the default option. However, when the â€ś-kind=byteâ€ť option is used, old codes (and new ones) can be compile with real*8 or real (kind=8).
Why does this have to be so complicated? (itâ€™s a rhetorical questionâ€¦I know how we got here).
But the fact is, this is all very confusing for beginners and people coming from other languages. This is why every day on stack overflow we see questions about why isnâ€™t this working right:
real*8 pi = 3.14159265358979323
Making people have to use iso_fortran_env
or selected_real_kind
just to declare a real number is just user-hostile. Nobody wants to type all that (hence real*8
).
Is there nothing that can be done about it?
In most other languages that people might come to Fortran from, double precision is the default. Can we just make double precision the default in Fortran? I feel like that would fix 80% of the problems people have.
It was said many times. If one hates iso-s and kinds, one can still write
DOUBLE PRECISION x=3.14159265358979323d0
On comp.lang.fortran I asked if we can Standardize real*8 as real(kind=real64)? There was not much support for doing so. I donâ€™t think making double precision the default is feasible, because it would break the codes calling procedures with single precision dummy arguments where the actual argument is a literal such as 3.14159265358979323.
Note, also, that the C world is not that-better. Indeed, the literal FP constants in C default to double
but it is only because the pre-ANSI C did all FP arithmetic in double precision. But then if you want to define a long double
value, you have to add L
suffix, otherwise it will get truncated to double
precision before being assigned to long double
variable.
PS. for some reason, this is not the case with integer C literal constants, for which the compiler assigns the first (â€śshortestâ€ť, starting from int
) type which contains the value.
I am pretty sure that because there are multiple ways of defining types in Fortran and people are unsure as to which on is â€śthe correct oneâ€ť, my research group eventually gave up and are using the way more satanic flags -fdefault-real-8
and -fdefault-double-8
That is great comment
Intel Fortran can add flags to define the size of default integer, real, double precision,
gfortran can easily do such things too, as @gnikit pointed out.
I like this style.
real(8) :: pi = 3.14159265358979323d0
I do the following,
integer, parameter :: r8=selected_real_kind(15,9)
real(kind=r8), parameter :: pi=4.0_r8*atan(1.0_r8)