 |
ORB5
4.00
|
|
ORB5
4.00
|
Guidelines for setting the floating point number precision.
A computer stores a number in a finite number of bits, so it can only represent finitely many values. Writing a program therefore means choosing, for every variable and every constant, how many bits to spend on it. That choice is a trade-off: more bits give more accurate results, but cost more memory and more memory traffic.
ORB5 used to make this choice in the build system, by asking the compiler to widen every real to 64 bits. It no longer does. The precision is now stated in the source code, using Fortran kind parameters. The responsibility has therefore moved to you: if you do not say what you want, the compiler gives you the default, and the default is single precision.
The short version:
The rest of this page explains what those rules buy you, and why they are worth the noise.
The two floating point widths in use are 32 bits (single precision) and 64 bits (double precision). What the extra bits buy is best stated in decimal digits, which is what Fortran's own precision and range intrinsics report:
| Width | Significant decimal digits | Decimal exponent range | Largest integer |
|---|---|---|---|
| 32 bits | 6 | \(10^{\pm 37}\) | 2'147'483'647 |
| 64 bits | 15 | \(10^{\pm 307}\) | 9'223'372'036'854'775'807 |
Six significant digits is the number to remember, and it is not a lot. Worse, it is a starting point rather than a guarantee: rounding errors accumulate over the millions of operations a time step performs, so a quantity that begins with six good digits can end with noticeably fewer. This is why ORB5 works in double precision throughout.
The integer column matters too. The default integer stops at \(-2'147'483'648 \leq x \leq 2'147'483'647\), that is about 2.1 billion. This is usually plenty, but not always: ORB5 can run with more than 2 billion Larmor points, and a Larmor point count or a Larmor point identifier will overflow a 32-bit integer silently.
Here is the whole problem in one program. All five assignments write to the same double precision variable, printed with the same f19.17 format:
Read this carefully, because it is not intuitive.
Declaring x as r8 did not protect it. In x = 1.0 / 3.0, the literals 1.0 and 3.0 are default reals, so the division is carried out in single precision, and only then widened to fill x. The extra digits were lost before the assignment ever happened, and widening afterwards cannot invent them back: the result is a 64-bit variable holding a value that is only correct to six digits. The same happens to 0.1, which was never 0.1 to begin with.
This is exactly the failure mode the old promotion flags used to hide, and it is why moving away from them meant rewriting the literals throughout the code, not merely the declarations.
The unifying rule: Fortran decides the precision of an expression from its operands, never from what you assign the result to. The assignment happens last, when it is already too late.
Since Fortran 90, every numeric type carries a kind parameter that selects the representation the compiler should use. Kind parameters are positive integers, but the values themselves carry no meaning: a kind is an opaque key identifying a representation, chosen by the compiler vendor.
That has a practical consequence, because there are three ways to spell "give me a 64-bit real" and only one of them is right:
| Spelling | Verdict |
|---|---|
| real*8 :: x | Non-standard. A vendor extension from the 70s that has never appeared in any Fortran standard. It was used because the FORTRAN language at the time had no mechanism for specifying the size of a variable. Compilers accept it out of habit, but gfortran -std=f2018 rejects it outright. Do not use it. |
| real(8) :: x | Standard, but not portable. For the compilers we use, 8 means 64 bits but this is guaranteed by nothing: kind values are keys, not byte counts. Do not use it. |
| real(kind=r8) :: x | Correct. It names the precision we want, and lets a single module decide how to obtain it. |
Nothing in the standard says that kind 4 is 32 bits, or that the default real is 32 bits, or that the default integer is 32 bits. What the standard actually requires is much weaker: that double precision be more precise than real, and that the compiler document the kinds it supports. Every compiler ORB5 is built with happens to agree on the usual widths, but that agreement is a convention rather than a guarantee, and it is invisible to someone reading the code. Naming the kinds explicitly fixes both problems at once.
You will meet the double precision keyword in older code. It is a declaration form rather than a type of its own: it declares a real, of the kind kind(0.0d0), and the standard requires only that this kind be more precise than the default real, not by how much. On our platforms it is exactly real(kind=r8), and the two are so interchangeable that a compiler rejects a generic interface overloaded on both as ambiguous. In other words, prefer using real(kind=r8) over double precision.
The related literal form is 1.0d0, meaning "the number 1, in double precision". It is correct on our platforms and it is not a precision bug, but prefer 1.0_r8, so that literals and declarations are tied to the same parameter and would follow it if it ever changed.
The kind parameters live in src/precision_kinds.F90, which defines four of them as thin aliases over the standard iso_fortran_env module. That module is the portable answer to the problem above: it names kinds by the property we care about and lets the compiler supply the value.
| Name | Value | Meaning |
|---|---|---|
| r4 | real32 | 32 bits float |
| r8 | real64 | 64 bits float |
| i4 | int32 | 32 bits integer |
| i8 | int64 | 64 bits integer |
The kind parameter must appear everywhere the default is not what you want: in declarations, in literals, and in conversion functions.
In Fortran the real intrinsic can be used for converting any number to a real, or to take the real part of a complex number. In the latter case, the resulting real part is of the same kind as the complex number. For this reason, we choose not to repeat the kind in this case.
A last pitfall concerning conversions is the cmplx intrinsic, used to create a complex number. If we wish to create a complex number from two literals, we must proceed as follows:
Note the three kind mentions. The first two target the literals, while the last specifies the kind of the resulting complex number. If it is omitted, a single precision complex number is generated by default.
By default, we work in double precision, that is r8, for reals and complex numbers, and we use the default integer for integers. Those two rules cover almost all the code. There are deliberate exceptions to each.
Single precision, r4, is used where halving the memory footprint of a large array matters more than the extra digits: the diagnostic marker attributes, which are only ever written out for diagnostics, the deposition weights (W_TYPE in src/deposition.inc), and the sorting work arrays in src/sorting.F90. Such values are computed in r8 and demoted explicitly when stored.
Eight-byte integers, i8, are used for the marker counts and marker identifiers, which can exceed the 2.1 billion ceiling of a 32-bit integer. The same reasoning already appears in Inputs for the input_integer8 input type.
The i4 kind is not used in declarations: the default integer is already 32 bits on every platform we target, so there is nothing to state, and writing integer(kind=i4) everywhere would be noise. It appears instead in narrowing conversions, where an expression computed in i8 has to come back to a default integer, as in src/sorting.F90 and src/parmove.F90:
The intermediate product would overflow a 32-bit integer, so it is computed in i8. The result is a loop bound and is known to fit, so it is converted back explicitly.
A forgotten _r8 is silent. It is not a compiler error, and it will not raise a warning either: evaluating in single precision and widening on assignment is perfectly legal Fortran, and the compiler cannot know it is not what you meant. Nothing in the build will tell you.
What catches these in practice is code review, and the regression tests, which compare against reference outputs and drift once a hot expression quietly loses nine digits. That feedback loop is slow, so the cheapest place to catch a missing suffix is while writing it.
We could also test by turning the old promotion flags into a diagnostic. Rebuild with -fdefault-real-8 -fdefault-double-8 (or the vendor equivalent listed in Background: why we do not use compiler flags anymore) and compare the results against a normal build. Those flags promote default reals and leave explicitly kinded declarations alone, so in code that names every kind they have nothing left to act on and the two builds should agree. Any difference means a default real was doing arithmetic somewhere, and the routines whose output moved are where to look.
Used this way, the flags find unstated kinds only. They say nothing about integer division, and they cannot see a kind that is explicit but wrong: an r4 where you meant r8 involves no default real at all, so both builds agree and the bug survives.
In summary:
ORB5 used to obtain double precision a different way: the build promoted every default real to 64 bits with compiler flags. Depending on the vendor, that meant -r8 (Intel, NVHPC), -fdefault-real-8 -fdefault-double-8 (GNU), or -s real64 (Cray). The source code was then written with bare real declarations and bare literals, and the flag silently made them double precision.
This is convenient, and it is dangerous:
For these reasons, the promotion flags were removed and the precision is now stated in the source, where it belongs.
1.15.0 on