When should I add kind suffixes like _dp to numbers in a formula?

1

Say that I have defined a central module for kind parameters as suggested in Floating Point Numbers — Fortran Programming Language . I understand that adding kind suffixes to constants is needed, such as x = 0.5_dp .But what about in parts of a formula? As a random example,

x = 0.5_dp
y = (2 * x**5) ** (1.0 / 3.0) + 3*exp(x*2)

Would the result be any different if I instead wrote

y = (2.0_dp * x**5.0_dp) ** (1.0_dp/3.0_dp) + 3.0_dp*exp(x*2.0_dp)

There’s a large difference in readability between the first and second example.

2

Long story short… always

You are allowed to don’t do it, but it is your responsibility to know that you are mixing kinds. That error might not be important at first, but it will bite you down the road, next month, next year or a few years after when you will have forgotten or it is now in someone else’s hand. Better get used to adopting good practices from the start if you can.

1 Like
3

I normally go the extra step and create parameters for commonly used constants and fractions and use them instead of literal constants. That way I can easily switch between single and double precision at compile time using some preprocessor logic. Note using parameters instead of literals is very common even in 40 year old code. The amount of extra typing (if any) is a small price to pay for greater flexibility.

Example:


module constants

use iso_fortran_env, SP=>REAL32, DP=>REAL64

#ifdef R4
   integer, parameter :: RK = SP
#else
   integer, parameter :: RK=DP
#endif
   real(RK), parameter :: one=1.0_RK
   real(RK), parameter :: two = 2.0_RK
   real(RK), parameter :: three=3.0_RK
   real(RK), parameter :: four = 4.0_RK
   real(RK), parameter :: five  = 5.0_RK
   real(RK), parameter :: half  = 0.5_RK
  real(RK), parameter  :: p5    = 0.5_RK
  real(RK), parameter :: third  = 1.0_RK/3.0_RK
  real(RK), parameter :: p33  = third

end module constants

program yourprog
   use constants

    real(rk) :: x, y

   x = half  !(or x=p5)
   y = (two*x**five)**third + three*exp(x*two)
 
end program yourprog

Now if you want to run on a GPU which only supports FP32 you just have to compile your constants.F90 file with a -DR4 compiler option. Yes a little more work now but a lot less work later

4

I suggest to run a calculation, for instance to approximate \pi

program test
use iso_fortran_env, only: dp => real64
implicit none
real :: pi_1
real(dp) :: pi_2, pi_3

pi_1 = 4 * atan(1.0)
pi_2 = 4 * atan(1.0)
pi_3 = 4 * atan(1.0_dp)

write (*, "(F13.11)") pi_1
write (*, "(F13.11)") pi_2
write (*, "(F13.11)") pi_3

end program test

Compiled with gfortran (14.2.0), this results

3.14159274101
3.14159274101
3.14159265359
1 Like
5

Add _dp to floating point literals, but leave integers alone. Your code would be

x = 0.5_dp
y = (2 * x**5) ** (1.0_dp / 3.0_dp) + 3.0_dp*exp(x*2)

The compiler may evaluate an integer power such as x**5 faster than x**5.0_dp.

2 Likes
6

Also, although the two expressions are mathematically equivalent for some values of x, they differ, for example, when x<0, and in any case, they are evaluated differently so they might not be the same computationally to the last bit. In this particular case, the x**5 expression is almost certainly evaluated with three multiplications as t*t*x where t=x*x is a compiler-generated intermediate value. The floating point exponent expression x**5.0_dp on the other hand is evaluated as something like exp(5.0_dp*log(x)), which is where the x>0 restriction arises. A compiler might recognize this identity and evaluate it in the more efficient way for an integer power when x>0, but the programmer cannot rely on that substitution.

I think the answer to the original question is really about how fortran evaluates expressions and applies its implicit type and kind conversions. Assuming that x and y are both real(dp) (and that dp is real64 and default real is real32), then the two expressions in the original post are not equivalent. The integer factors in the original expression are indeed converted to real(dp) values in that expression, but the (1.0/3.0) term is not evaluated as (1.0_dp/3.0_dp) but rather as real((1.0/3.0),dp), so the values of the two expressions differ. The initial x=0.5_dp is somewhat of a special situation because x=0.5 would be evaluated as x=real(0.5,dp), and the rhs just happens to be 0.5_dp because that value is common to all floating point kinds. If all of the implicit type and kind conversions are written explicitly, then the original expression would be evaluated as

y = ((real(2,dp) * (x**5)) ** real((1.0/3.0),dp)) + (real(3,dp) * exp(x*2))

That is not the same as the second expression in the original post. Sometimes all of the precedence rules and implicit type conversion rules are so complicated that even experienced programmers write out expressions with extra parentheses and with explicit type conversions just to make it clear what the expression means. For example, if the programmer really did want real((1.0/3.0),dp) rather than (1.0_dp/3.0_dp), then he might write it explicitly that way to make the intent clear to himself and to other programmers.

It might look like extra work is being done in this case, but in almost all situations there isn’t; the type and kind conversions are done at compile time in both cases exactly the same when possible, and both at run time in the same way otherwise.

7

I’m doing some digging and this is confusing me even further:

module kinds    
    implicit none
    
    !> Single precision real numbers, 6 digits, range 10⁻³⁷ to 10³⁷-1; 32 bits
    integer, parameter :: sp = selected_real_kind(6, 37)
    !> Double precision real numbers, 15 digits, range 10⁻³⁰⁷ to 10³⁰⁷-1; 64 bits
    integer, parameter :: dp = selected_real_kind(15, 307)
    !> Quadruple precision real numbers, 33 digits, range 10⁻⁴⁹³¹ to 10⁴⁹³¹-1; 128 bits
    integer, parameter :: qp = selected_real_kind(33, 4931)

    !> Char length for integers, range -2⁷ to 2⁷-1; 8 bits
    integer, parameter :: i1 = selected_int_kind(2)
    !> Short length for integers, range -2¹⁵ to 2¹⁵-1; 16 bits
    integer, parameter :: i2 = selected_int_kind(4)
    !> Length of default integers, range -2³¹ to 2³¹-1; 32 bits
    integer, parameter :: i4 = selected_int_kind(9)
    !> Long length for integers, range -2⁶³ to 2⁶³-1; 64 bits
    integer, parameter :: i8 = selected_int_kind(18)

    integer, parameter :: wp = dp

end module kinds



program test_kinds
    use kinds
    implicit none

    real(dp) :: x
    real(dp) :: y1, y2, y3, y4, y5, y6, y7, y8, y9, y10

    x = 0.5_dp

    ! y1 shows a loss of precision
    y1 = (2 * x**5) ** (1.0 / 3.0) + 3*exp(x*2)
    y2 = (2 * x**5) ** (1.0_dp / 3.0_dp) + 3*exp(x*2)
    y3 = (2.0_dp * x**5.0_dp) ** (1.0_dp/3.0_dp) + 3.0_dp*exp(x*2.0_dp)

    ! y4, y5 and y6 give me the same values
    y4 = 1 + y3
    y5 = 1.0 + y3
    y6 = 1.0_dp + y3

    ! y7 and y8 also give me the same values
    y7 = 0.5 + y3
    y8 = 0.5_dp + y3

    ! y9 loses some digits
    y9 = 0.325 + y3
    y10 = 0.325_dp + y3

    write(*,*) 'y1 = ', y1
    write(*,*) 'y2 = ', y2
    write(*,*) 'y3 = ', y3
    write(*,*) 'y4 = ', y4
    write(*,*) 'y5 = ', y5
    write(*,*) 'y6 = ', y6
    write(*,*) 'y7 = ', y7
    write(*,*) 'y8 = ', y8
    write(*,*) 'y9 = ', y9
    write(*,*) 'y10 = ', y10

end program

In this example, the compiler (gfortran 14.2.1) has issues with the (1.0 / 3.0) division at y1. But it seems unnecessary to specify _dp to the rest of the numbers in the formula (in summary, y2 gives me the same result as y3).

So I’ve tried some random tests. summing 1, 1.0 or 1.0_dp gives me the same result up until the 16th decimal. Summing 0.5 or 0.5_dp also gives me the same result. However, summing up a broken real such as 0.325 doesn’t give me the same result as 0.325_dp.

I’m probably going to add _dp to anything just in case, but is there any logic to the kind conversions?

1 Like
8

1.0/3.0 and 0.325 each cannot be represented exactly in binary.

1.0/3.0 = 0.010101010101...
0.325   = 0.010100110011...

In Fortran, 1.0/3.0 and .0325 are expressions of type REAL and so stored in a 32-bit floating point number.
1.0_dp/3.0_dp and .0325_dp have type REAL(dp) and so are stored in 64 bits.

So the differences you are seeing are due to the different number of bits of precision of those expressions. You don’t have to use _dp when you know the expression can be represented exactly in the smaller type, e.g. 1.0/2.0 or 0.375. But it’s safest to always use the kind for non-default REAL types.

9

It would be good to make sure you understand how precision varies between real ( and integer! ) kinds

There are a number of different cases in the OP example.

For “(2 * x ** 5)”, using 2, 2.0 or 2.0_dp will all produce the same precision, as x is kind(dp). In this case, clarity is best and I would use “2”. For “x ** 5”, it has been preferable to use x ** 5, rather than x ** 5.0 or x ** 5.0_dp, as the compiler may recognise x ** 5 is xxxxx, while using a real constant 5.0 or 5.0_dp could invoke exp and log. Although we don’t know how smart compilers are at interpreting 5, 5.0 or 5.0_dp, which can all come with the same precision.
For “2 * x**5”, the precision/kind of x dictates the precision of the calculation.

For “3.0_dp * exp(x * 2)”, I would again go for clarity as again using 3, 3.0 or 3.0_dp will all produce the same precision, result as x is kind(dp), so I would use “3 * exp(x * 2)”

Now, “(1.0_dp / 3.0_dp)” vs “(1.0 / 3.0 )” is a clasic example of how to loose precision, as
"(1/3) = 0 ! ", “(1./3.)” and “(1.0_dp / 3.0_dp)” will all produce significantly different values when used in _dp Fortran calculations.
My preferred usage for these cases is to define a real parameter of the required precision, such as:
real(dp), parameter :: third = 1.0_dp / 3.0_dp.
This will guarantee that you have the value of what you expect 1/3 is, to 15 significant figures, while (1./3.) will have the value of 1/3 only to 6 significant figures. This is very important to recognise this, if you are aiming for double precision computation (which is typically kind(real64) or 8-byte real computation)

If your want your integer or real constants to be better than default integer or real precision, you need to be careful with annotating constants with the required kind.

This is especially important when using kind(int64) integer constants, which can often be missed when moving to the 64-bit computing environment.

1 Like
10

It is interesting to note that in all these cases, the non-annotated version below will all produce the same precision constants, such as:

real(RK), parameter :: one = 1
real(RK), parameter :: two = 2
real(RK), parameter :: three= 3
real(RK), parameter :: four = 4
real(RK), parameter :: five = 5
real(RK), parameter :: half = 0.5
real(RK), parameter :: p5 = 0.5
real(RK), parameter :: third = one / three

But significantly, the following will not
real(RK), parameter :: tenth = 0.1

11

When should I add kind suffixes like _dp to numbers in a formula?

Always for floating point (unless it is an integer).

Otherwise you need to learn a lot of details, but I prefer to focus on the actual algorithm, so I always append _dp and don’t worry about it and everything works as expected.

12

That is interesting. I’ve always used the kind annotations because I know beforehand that will give me the correct values i’m looking for. I don’t mind a little extra typing though some people appear to think that typing even a couple of extra letters is a major inconvenience. I was lucky that I took typing (on a real mechanical and electo-mechanical typewriters - yes I’m old :grin:) in high school so I was spared having to “hunt and peck”

13

No, the compiler is doing exactly what you have told it to do. This is not a compiler error. As explained in detail in my previouis post, that part of the expression is evaluated as real((1.0/3.0),dp) which does not have the same numerical value as (1.0_dp/3.0_dp). In contrast, small integer values are converted to have the same values, e.g. real(3,dp) is the same as 3.0_dp, and some other floating point values are also converted exactly, e.g. real(0.5,dp) is the same as 0.5_dp.