1

I am playing around with selected_real_kind() against real(kind=8) to check how much extra processor time I require to improve numerical precision. The set of parameters I use are:

module param
implicit none
!extended double precision
INTEGER, PARAMETER :: xtd = SELECTED_REAL_KIND(16, 40)
save

        !List of Input paramters and constants
real(xtd), parameter :: PI =4.0_xtd*atan(1.0_xtd)
integer, parameter :: L = 512
integer, parameter :: steps = 160000
real(xtd), parameter :: dt = 0.001_xtd

real(xtd), parameter :: lambda = 1.0_xtd
real(xtd), parameter :: mu = 0.0_xtd
real(xtd), parameter :: a = dt*(lambda+mu), b = dt*(lambda-mu)

end module param

The output when I print these parameters is:

Pi =    3.14159265358979323851      
 L =          512
 steps =       160000
 dt =    9.99999999999999999958E-0004
 lambda =    1.00000000000000000000      
 mu =    0.00000000000000000000      
 a = dt*(Lambda+Mu) =    9.99999999999999999958E-0004
 b = dt*(Lambda-Mu) =    9.99999999999999999958E-0004

The code for double precision real(kind=8),

module param
implicit none

save

        !List of Input paramters and constants
real(kind=8), parameter :: PI =4.d0*atan(1.d0)
integer, parameter :: L = 512
integer, parameter :: steps = 160000
real(kind=8), parameter :: dt = 0.001d0

real(kind=8), parameter :: lambda = 1.0d0
real(kind=8), parameter :: mu = 0.0d0
real(kind=8), parameter :: a = dt*(lambda+mu), b = dt*(lambda-mu)

end module param

the output for which is

Pi =    3.1415926535897931     
 L =         512
 steps =        160000
 dt =    1.0000000000000000E-003
 lambda =    1.0000000000000000     
 mu =    0.0000000000000000     
 a = dt*(Lambda+Mu) =    1.0000000000000000E-003
 b = dt*(Lambda-Mu) =    1.0000000000000000E-003

I fear this difference in numerical constants may be contributing to the slower run-time for a marginally precise constant. Is there a better way to write these parameters using selected_real_kind()?

ferro11001
  • 13
  • 4
  • 1
    Where exactly do you use `kind=8`? Normally, you would do `INTEGER, PARAMETER :: xtd = 8` instead. It is much better. Please show the **full** code so we see how exactly you declare your numeric literals and constants. – Vladimir F Героям слава Jun 02 '22 at 16:33
  • Namely, if you do `real(kind=8), parameter :: dt = 0.001`, you are doing it wrong. It should be `real(8), parameter :: dt = 0.001_8`, but a named constant `xtd` is, of course, better. Regarding your fears, you would have to explain them more thoroughly, but they are likely unsubstantiated. – Vladimir F Героям слава Jun 02 '22 at 16:34
  • @VladimirFГероямслава my concern is that the numerical constant dt when written with the extended precision is 9.99999999999999999958E-0004 instead of 1.0000000000000000E-003. Could that be contributing to a slower run time for the same parameter values? – ferro11001 Jun 02 '22 at 17:03
  • 1
    I'll humbly suggest that you need to read Goldberg's paper (perhaps, several times). https://www.itu.dk/~sestoft/bachelor/IEEE754_article.pdf – steve Jun 02 '22 at 17:10
  • BTW `1d0` is not equivalent to `kind=8`. 1d0 is for double precision, but kind=8 might be something else. The correct literal for kind 8 is `1.0_8`. See https://stackoverflow.com/questions/838310/fortran-90-kind-parameter – Vladimir F Героям слава Jun 02 '22 at 18:59

1 Answers1

1

Many options for real/integer datatypes have been stacking over the decades in the Fortran standard and selected_real_kind is definitely not the best one.

First you should print its output to se what's being used:

print *, 'my compiler chose real kind = ',xtd

I suspect it's going to be 16 and not 8 (typically chosen with selected_real_kind(15,307)). In other words, with selected_real_kind you usually cherry-pick the digit numbers to be sure your compiler will pick the right datatype.

This is ugly, bad, unsafe. It may work in scattered code, but will almost certainly cause headaches in anything that's not elementary. My suggestion is to never use it: use the iso_fortran_env provided datatypes instead:

use iso_fortran_env, only: real32,real64,real128
integer, parameter :: x64 = real64 
integer, parameter :: x32 = real32

real(x64) :: pi64 = acos(-1.0_x64) ! this is a double pi
real(x32) :: pi32 = acos(-1.0_x32) ! this is a float  pi

print *, '32-bit pi = ',pi32
print *, '64-bit pi = ',pi64

In your case, you're using quad-precision (128bit floats), which is far more computationally intensive than double precision, just to get one more digit.

Federico Perini
  • 1,414
  • 8
  • 13
  • 3
    I think the answer could be improved by a short discussing about using kinds for constants - the op is using the very Fortran77 `d0` which doesn't mix well with proper use of kinds. I also disagree with the vehemence you show against `selected_real_kind`, if used correctly it can be exactly thee required tool. But I'm not going to argue with it, in most cases `iso_fortran_env` is simpler and just what is required. – Ian Bush Jun 03 '22 at 07:14
  • The is fortran env constants speak about storage size, not about the numerical precision or a C equivalent (float/double). – Vladimir F Героям слава Jun 03 '22 at 08:20
  • But foremost, in my reading, the question is not at all about the source of the kind number, but about actually using it and the differences in numerical results it causes. For the compiler in question, there should not be any output difference if the value 8 is entered directly or using `real64` or as a result of an intrisic function. I can see even why someone would set the kind number to 8 directly, e.g. in some compiler specific tests. What I hate is directly putting the magic number all over the code instead of one named constant. If the programmer knows why they chose value 8, so be it. – Vladimir F Героям слава Jun 03 '22 at 08:23
  • ok, I'll drop using d0 for the constants, and use selected_real_kind with the appropriate integer parameter for double precision or its modification thereof. As iso_fortran_env is compatible with fortran-2003 or later, while I am working with fortran-90, I'll remember it for future use. – ferro11001 Jun 03 '22 at 10:41