Question about double precision on ARM processors

22

Following is the code run on Cortex hardware. The results are the same as for the x86 hardware. Legacy fortran code which does not incorporate the Fortran 90 features may have accuracy errors if run on Cortex ARM hardware. There is a large number of lines of legacy code in existence, eg in ACM Transactions on Mathematical Software etc.

/dbclient: Caution, skipping hostkey check for localhost

Welcome to Ubuntu in UserLAnd!
userland@localhost:~ nano testcortex.f userland@localhost:~ gfortran testcortex.f -o testcortex
userland@localhost:~ ./testcortex 8 4 0.123400003 0.12340000271797180 0.12340000000000000 userland@localhost:~ lscpu
Architecture: aarch64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
On-line CPU(s) list: 0-3
Vendor ID: ARM
Model name: Cortex-A53
Model: 4
Thread(s) per core: 1
Core(s) per cluster: 4
Socket(s): -
Cluster(s): 1
Stepping: r0p4
CPU max MHz: 2001.0000
CPU min MHz: 850.0000
BogoMIPS: 26.00
Flags: fp asimd evtstrm aes pmull sha1 sha2 crc32 cpuid
Vulnerabilities:
Itlb multihit: Not affected
L1tf: Not affected
Mds: Not affected
Meltdown: Not affected
Spec store bypass: Not affected
Spectre v1: Mitigation; __user pointer sanitization
Spectre v2: Not affected
Srbds: Not affected
Tsx async abort: Not affected
userland@localhost:~$

use, intrinsic:: iso_fortran_env, only: sp=>real32, dp=>real64
real(sp) :: s = 0.1234_sp
real(dp) :: ds = 0.1234_sp, dd = 0.1234_dp
!double precision :: dp = 0.711D25
!real :: sp = 0.711
print *,dp
print *,sp
print *,s
print *,ds
print *,dd
end program

23

Then they learn, and they stop making this error. We all did that at a moment.

1 Like
24

Sorry, but absolutely nothing that you have posted can support this assertion.

26

I accept the criticism that nothing I have posted so far can support the assertion that legacy fortran code without the Fortran 90 features may have errors in accuracy if run on Cortex hardware. Here is a code outputs the same results on Cortex ARM and Intel x86 hardware.

userland@localhost:~ nano test2.f userland@localhost:~ gfortran test2.f -o test2
userland@localhost:~ ./test2 7.1099999999999997E+024 0.711000025 userland@localhost:~ cat test2.f

! use, intrinsic:: iso_fortran_env, only: sp=>real32, dp=>real64
! real(sp) :: s = 0.1234_sp
! real(dp) :: ds = 0.1234_sp, dd = 0.1234_dp
double precision :: dp = 0.711D25
real :: sp = 0.711
print *,dp
print *,sp
! print *,s
! print *,ds
! print *,dd
end program
userland@localhost:~$ lscpu
Architecture: aarch64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
On-line CPU(s) list: 0-3
Vendor ID: ARM
Model name: Cortex-A53
Model: 4
Thread(s) per core: 1
Core(s) per cluster: 4
Socket(s): -
Cluster(s): 1
Stepping: r0p4
CPU max MHz: 2001.0000
CPU min MHz: 850.0000
BogoMIPS: 26.00

Many thanks to all the contributors to this discourse for their time and effort.

27

The above url is very interesting and discusses the -march compile and build option.

28

Google scholar on processor dependecy in Fortran Program Code
https://scholar.google.com.au/scholar?hl=en&as_sdt=0%2C5&q=processor+dependency+in+fortran&btnG=

29

Floating point arithmetic is processor dependant, that’s nothing new, and that’s not specific to Fortran.

A must read on that topic:

30

Is there any compiler where kind=8 is not double precision? One, the NAG compiler uses a different convention, but it can be made to match the rest of them using the compiler flag -kind=bytes. Given it’s a commercial compiler, most beginners probably can’t afford to use it.

Personally I prefer to just use wp=kind(1.0d0) or wp=c_double over real64, even if they are de facto all the same time on the majority of machines. Dr. Fortran (@sblionel) has explained this before:

In my view, [using the named constants REAL32, REAL64 and REAL128] is little better than the old *n extension in that it tells you that a type fits in that many bits, but nothing else about it. As an example, there’s one compiler where REAL128 is stored in 128 bits but is really the 80-bit “extended precision” real used in the old x87 floating point stack registers. If you use these you might think you’re using a portable feature, but really, you’re not and may get bitten when the kind you get doesn’t have the capabilities you need. My advice is to not use those constants unless you truly understand what they do and do not provide.

I agree it’s complicated for beginners but I don’t see it as much easier in other languages. When using NumPy it can also be tricky to tell if an array is np.float32, np.float64, or np.float128 (not quad!) and miss the fact conversions are happening. To make things more confusing Pythons default float type is actually C’s double under the hood.

In Julia, the numeric type system is completely general, and you can use even special types like MPFR, GMP, unums, etc. The same genericity can leads to a problem known as type instability, where the AoT compiler can’t predict the type in advance, and is forced to move back to interpreted execution, destroying performance.

In C++, literals such as 3.0 are double by default, which leads to problems when people want to use float’s and forget to write them as 3.0f. In mixed-precision expressions, the compiler will cast the variables to the longer type. On a consumer GPU, where double precision units are scarce, this can kill the performance.

At the end of the day, on an x86-64 machine, your numbers are packed into the same registers, and just sign-extended with zeros… Don’t get me wrong, the differences are huge when it comes to precision and bandwidth. I just want to highlight that the types and kind parameters are already an abstraction in the first place.

smult_:
        vmovss  xmm0, DWORD PTR [rdi]
        vmulss  xmm0, xmm0, DWORD PTR [rsi]
        ret
dmult_:
        vmovsd  xmm0, QWORD PTR [rdi]
        vmulsd  xmm0, xmm0, QWORD PTR [rsi]
        ret

function smult(a,b)
real(kind(1.0)) :: a, b, smult
smult = a*b
end function

function dmult(a,b)
real(kind(1.0d0)) :: a, b, dmult
dmult = a*b
end function
31

Anyways, I agree with your meme, I’d only change the second-to-last bullet to:

  • Define wp as wp = kind(1.0d0) to get double precision. Put it in a module if you want to use it in multiple files.
32

I think a compiler should warn about this by default: LFortran mode to warn when there is a conversion loss (double to single precision) ¡ Issue #316 ¡ lfortran/lfortran ¡ GitHub (it can be turned off if you know what you are doing). It will catch many bugs.

3 Likes
33

If the f77 code was written in an ambiguous way, then yes there can be accuracy errors and portability issues. But it was always possible to write the code in an unambiguous way, in which case there would be no accuracy issues and no surprising behavior. Portability was always an issue with pre f90 fortran because the precision of real and double precision was not specified by the standard, it was purposefully left as compiler (processor) dependent. For example, REAL on an IBM mainframe was 32-bit floating point with radix 16, while REAL on a Cray was 64-bit floating point with radix 2, and there were many may other combinations of radix and word size supported by fortran. The fortran KIND system introduced in f90 is/was a major improvement to this prior situation.

34

The following was run on Android 13 on an ARM processor on the Userland app from the Google play store. Note the difference in output values between the print and formatted write statements

7.1099999999999997E+024
0.71100000000000D+25
0.711000025

     double precision dp /.711D25/
      real   sp  /.711/
       print *,dp
       write(6,9000) dp

9000 format(" ",D20.14)
print *,sp
end program

35

Do you think this is an error? One shows two more decimal digits than the other, so the second one is rounded. Also, you keep initializing the variables with nonstandard syntax. Why not use the right syntax?

2 Likes
36

I honestly don’t get what you are trying to show or prove with the multiple tests you have posted… Keep in mind that:

  • floating point numbers do not represent the mathematical real numbers
  • floating point arithmetic is an approximation of mathematical real number arithmetic
  • this approximation is processor dependent
  • don’t confuse what is printed/written with the actual internal representation of a floating point number
  • it’s up to the developer (you) to check the precision of the types they are using, and to select the right type one according to the task
  • this precision can be processor dependent
  • literal floating point constants such as 0.711 or 0.711e25 have the default REAL type, with the associated precision (formally 6 digits on most machines).
  • literal floating point constants such as 0.711d25 have the DOUBLE PRECISION type. The standard requires this type to have a precision of at least 10 digits, and it is rather 15 on most machines.
2 Likes
37

If compiler options are allowed, I think the correct answer is yes.

program kinds
   use, intrinsic :: iso_fortran_env, only: real32, real64, real128
   write(*,*) real32, real64, real128
   write(*,*) 'default real=', kind(1.0), ' double precision=', kind(1.D0)
end program kinds
$ gfortran kinds.f90 && a.out
           4           8          16
 default real=           4  double precision=           8
$ gfortran -fdefault-real-8 kinds.f90 && a.out
           4           8          16
 default real=           8  double precision=          16

$ nagfor kinds.f90 && a.out
NAG Fortran Compiler Release 7.1(Hanzomon) Build 7114
[NAG Fortran Compiler normal termination]
 1 2 3
 default real= 1  double precision= 2
$ nagfor -r8 kinds.f90 && a.out
NAG Fortran Compiler Release 7.1(Hanzomon) Build 7114
[NAG Fortran Compiler normal termination]
 1 2 3
 default real= 2  double precision= 3

I read an article a while back that said some modern IBM hardware supports nine different floating point formats. These included the old IBM 360 style radix 16 formats for backwards compatibility, along with more modern binary and decimal formats. When I read that article, I thought how well fortran is situated to handle that hardware. All that would be required is to assign a different KIND value to all of those different floating point formats. Of course, that is a nontrivial programming task for the compiler writers because they also need to support all possible interconversions and mixed-kind expression evaluations, but at least the language itself is already up to speed to that task. And there is also the new 16-bit floating point kinds that are becoming more popular, due to gpu hardware I think; the language is ready and waiting for those.

I also thought back to the 1980s about the VAX hardware. I’m pretty sure that among many fields of science and engineering, the DEC VAX was the most popular computer architecture. It supported one real type, two double precision types, and one quad precision type, but with their f77 compiler, only one of those double precision types could be supported within a compilation step. When I was reading about the upcoming fortran 8x language standard, I thought the DEC VAX would immediately have a really nice feature of allowing the programmer to select, variable by variable, which of those double precision types to use (one had more precision but a smaller exponent range than the other). But my optimism was never realized. DEC as a company aggressively opposed the adoption of fortran 8x, and the VAX hardware would never support a full-fledged f90 compiler before the hardware was obsolete and the company evaporated.

So this raises the question, why aren’t there any fortran compilers, right now, that do this type of spectacular KIND gymnastics? As far as I know, there is no IBM fortran compiler that fully supports all of those floating point formats. The NAG compiler is the only one I know of that supports the REAL16 KIND (it does it by converting back and forth with REAL32). The KIND convention was part of fortran 8x in the mid 1980s, and it has been part of the standard since f90, over three decades ago. No other language is positioned so well to take advantage of all those different formats. It seems like it would be a big win for everyone, hardware vendors, compiler writers, programmers, and end users. What is missing, what is it that has been, and is, preventing this from happening?

1 Like
38

You mentioned historical hardware that was not supported, but such hardware is not in wide use anymore. For modern hardware, besides the f16 and bf16, which other floating point format would you like to see supported besides the usual f32, f64 and f128?

39

well Nvidia has tf32 which is in between float16 and float32. also ieee currently has a committee talking about float8 formats.

40

The real kind I would like in compilers that don’t offer it is gfprtran’s precision=18
one. Last time I looked my kinds.f90 program gave


 Properties of real kinds:

                        single     double   extended       quad precision
          kind . . . . . .   4          8         10         16
          minexponent     -125      -1021     -16381     -16381
          maxexponent  . . 128       1024      16384      16384
          range             37        307       4931       4931
          digits . . . . .  24         53         64        113
          precision          6         15         18         33
          radix  . . . . .   2          2          2          2
          bits (iolength)   32         64        128        128
          bits (IEEE 754)   32         64         80        128
          epsilon   1.19E-0007 2.22E-0016 1.08E-0019 1.93E-0034
                  = 2.0**  -23        -52        -63       -112
          huge      3.40E+0038 1.80E+0308 1.19E+4932 1.19E+4932
                  ~ 2.0**  128       1024      16384      16384
          tiny      1.18E-0038 2.23E-0308 3.36E-4932 3.36E-4932
                  = 2.0** -126      -1022     -16382     -16382
41

It is not a high priority, but when I heard about the IEEE decimal formats, I thought I would like to experiment with those. There is some hardware that supports those formats, but I don’t know of any fortran compilers anywhere that supports them. One might think that fortran would be the first language that anyone used to program with these data types, but for some unknown reason, that is not the case.

42

It would certainly be an interesting development if f128 was hardware supported, including AVX512.
Given the huge changes in memory capability, running large models with reduced round-off errors should be available one day ?
I had experimented with f80 for my structural finite element problems, but it appeared that modelling approaches were more significant than round-off errors when looking at accuracy being achieved. Floating point precision does not account for all the errors.
Present f128 performance is not a viable option for my problems.

1 Like