"real" type of a calculation with mixed precisions

Language enhancement
1

The result of a calculation involving mixed real precisions should use the highest precision:

program precisions
  use iso_fortran env: only: sp => real32, dp => real64
  implicit none
  real(sp) :: x  = 1.23
  real(dp) :: xx = 1.23_dp
  
  print *, "single * single: x * x     = ", x*x     !  1.51289999
  print *, "double * double: xx * xx   = ", xx*xx   !  1.5128999999999999
  print *, "single * double: x * xx    = ", x*xx    !  1.5129000234603882
end program precisions

I expect the calculation using single and double precision to use the higher precision for the result.

Also, a calculation such as this:

1.23 * 1.23_dp

should automatically use the higher precision but does not.

I am porting some code from matlab that has thousands of lines and having to append the precision to each number, e.g., 1.23_dp) is a pain!

I should be able to append the precision to the first number

z = 1.23_dp * 3.56 * 7.89

and have the result be double precision without having to append _dp to each number.

2

Many compilers provide forced changes to the default floating-point precision. For example, the option -fdefault-real-8 of the GFortran compiler allows you to default to double precision for floating-point numbers instead of the single precision in the language standard.
Hopefully, this meets your needs and you no longer need to add the _dp suffix.

3

Thanks @zoziha. That should help.

Nevertheless, I think this should be the standard behaviour with mixed precision calculations without the need for compiler flags. Hence, the enhancement tag.

4

I think you could use a regex to achieve this (regex - Regular expression for floating point numbers - Stack Overflow). There may be some issues if you want to preserve the default precision in some procedures, or if you have real literals within strings.

This type of systematic change should be possible using fpt (available free of charge for academic use with open source projects): http://www.simconglobal.com/fpt_ref_systematic_changes.html. Perhaps @Jcollins can offer more advice on this.

Yup :expressionless: . There was a proposal by @FortranFan to address this: Default KINDs for constants and intrinsics · Issue #78 · j3-fortran/fortran_proposals · GitHub at the module- or procedure-scope.

This is not just a Fortran issue; particularly when porting C++ codes to GPU programming frameworks, the change from double to float literals might be desirable to improve bandwidth and the number of available arithmetic units (some hardware doesn’t even support double precision arithmetic). A couple accidental double literals and you can lose a few percent performance due to the runtime conversions.

5

These options should be used with care. The option -fdefault-real-8 will also convert double precision to 16-bytes if possible: Fortran Dialect Options (The GNU Fortran Compiler).

You could use -freal-4-real-8 to constrain the kind promotion to default reals only. Steve Kargl, one of the gfortran contributors, has warned repeatedly against the use of these options. They should be used primarily as a crutch before a full refactoring can be completed.

2 Likes
6

You code is doing precisely what it has been instructed to do.

Consider:

program test
implicit none

real(4) x
real(8) y,z

x=1.23

y=1.23d0

write(*,'(2G32.24)') x,y

write(*,'(G32.24)') x*y

z=x

write(*,'(G32.24)') z*y

z=1.23d0

write(*,'(G32.24)') z*y

end program

and its output:

   1.23000001907348632812500       1.22999999999999998223643    
   1.51290002346038821023910    
   1.51290002346038821023910    
   1.51289999999999991153743

The problem is that 1.23 is not exactly representable with floating point numbers. Thus 1.23 stored as single precision has an exact decimal expansion of 1.23000001907…

Even if I then copy x to z (which is double precision) in the code above, the multiplication z*y which is carried out as a double precision operation gives the same result as x*y. Thus we can conclude that the compiler is indeed first converting x to double precision and then performing the multiplication.

I like to use 1.23d0 and 1.23e0 in our code, which makes the intent clear to both the programmer and compiler.

2 Likes
8

Fortran converts mixed precision expressions according to well defined rules that have been in place for decades. Perhaps it would help with your understanding to write out the long form of the mixed-precision expression.

  print *, "single * single: x * x     = ", x*x           !  1.51289999
  print *, "double * double: xx * xx   = ", xx*xx         !  1.5128999999999999
  print *, "single * double: x * xx    = ", real(x,dp)*xx !  1.5129000234603882

This last expression is exactly the same as the original expression, but it uses an explicit kind conversion rather than the implicit kind conversion. It should now be clear why the last expression is not the same as either of the two previous expressions.

Perhaps you are confused why the following almost equivalent program produces this output:

program precisions
  use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64
  implicit none
  real(sp) :: x  = 1.25
  real(dp) :: xx = 1.25_dp

  print *, "single * single: x * x     = ", x*x                                                     
  print *, "double * double: xx * xx   = ", xx*xx                                                         
  print *, "single * double: x * xx    = ", x*xx                                                          
end program precisions
$ gfortran xxx.f90 && a.out
 single * single: x * x     =    1.56250000    
 double * double: xx * xx   =    1.5625000000000000     
 single * double: x * xx    =    1.5625000000000000

There is an important but perhaps subtle, difference in the two programs. Do you understand this output?

I don’t know how Matlab does mixed precision arithmetic, but if it does something different from fortran, I would think it would be difficult for programmers to use. As far as appending the kind values for constants of nondefault kind, there are tools that will do that automatically, but in this case you would need to be careful to avoid converting the values that are intended to be default kind. That is, the converted expression x=1.23_dp might not be what is intended, and it would certainly be confusing for a human to look at that expression and to determine what was the original intent. Using compiler options to automatically promote kind values of variables and/or constants has similar issues, they can obscure or even conflict with the original programmer’s intent.

9

Maybe I’m missing something, but this doesn’t appear to be related to mixed precision arithmetic at all. This is simply a result of Fortran’s default real kind being “single precision”.

When you provide a hard coded constant of 1.23 it is a “single precision” real number, which is rounded before it’s stored because it’s not exactly representable in binary. If you then upconvert it to a “double precision” real kind (whether implicitly or explicitly) you don’t somehow gain back the rounding error. Whatever you do with the number afterward has reduced precision.

As far as I know, Matlab and C treat hard coded constants a double precision by default. That’s the entire difference you’re seeing.

10

I’m replying to my own previous post because there is one other interesting feature demonstrated with these values. Here is the code with an extra print statement:

program precisions
  use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64
  implicit none
  real(sp) :: x  = 1.23
  real(dp) :: xx = 1.23_dp
  
  print *, "single * single: x * x     = ", x*x     !  1.51289999
  print *, "double * double: xx * xx   = ", xx*xx   !  1.5128999999999999
  print *, "single * double: x * xx    = ", x*xx    !  1.5129000234603882
  print '(b36.32/b65.64)', x, xx
end program precisions
$ gfortran xxx.f90 && a.out
 single * single: x * x     =    1.51289999    
 double * double: xx * xx   =    1.5128999999999999     
 single * double: x * xx    =    1.5129000234603882     
    00111111100111010111000010100100
 0011111111110011101011100001010001111010111000010100011110101110

I shifted the binary representation of the two floating point values to line up the high-order bits of the mantissas. From this, you can see that the sp value is the correct rounded value from the dp value and that the sp value is larger than the dp value. However, when printed with list-directed i/o, it appears as if the sp value is smaller than the dp value. When the products are printed, it then appears as if the small*big product is larger than the big*big product, an apparent paradox. It is only when the binary representations, rather than the rounded decimal values, are printed out that this paradox is understood. It is just an artifact of the rounded decimal values.

11

I find that using integer values in REAL variables can make the situation clearer

program precisions
  use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64
  real(sp) :: x  = 2**24+1._sp ! RHS sp, can't represent that odd number
  real(dp) :: xx = 2**24+1._sp ! RHS sp, still can't
  real(dp) :: x2 = 2**24+1._dp ! RHS dp, finally, an odd number

  print 999, "x     = ", x,   MOD(INT(x),2)
  print 999, "xx    = ", xx , MOD(INT(xx),2)
  print 999, "x2    = ", x2,  MOD(INT(x2),2)
  print '(b36.32/b65.64/b65.64)', x, xx,x2
  999 Format (A30,F17.0,5X,I1)
end

gives output

                      x     =         16777216.     0
                      xx    =         16777216.     0
                      x2    =         16777217.     1
    01001011100000000000000000000000
 0100000101110000000000000000000000000000000000000000000000000000
 0100000101110000000000000000000000010000000000000000000000000000
1 Like
12

Here is another followup to my previous post. This one shows the full set of numerical comparisons.

program precisions
  use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64
  implicit none
  real(sp) :: big   = 1.23_sp
  real(dp) :: small = 1.23_dp

  print '(a,b36.32/a,b65.64)', '   big=', big, ' small=', small
  
  print *, "big * big     =", big * big
  print *, "small * small =", small * small
  print *, "big * small   =", big * small
  print *, 'big > small             =', big > small
  print *, 'big*small > small*small =', big*small > small*small
  print *, 'big*big > small*small   =', big*big > small*small
  print *, 'big*big > big*small     =', big*big > big*small
end program precisions
$ gfortran xxx.f90 && a.out
   big=    00111111100111010111000010100100
 small= 0011111111110011101011100001010001111010111000010100011110101110
 big * big     =   1.51289999    
 small * small =   1.5128999999999999     
 big * small   =   1.5129000234603882     
 big > small             = T
 big*small > small*small = T
 big*big > small*small   = F
 big*big > big*small     = F

In exact arithmetic, all four of those last expressions would result in T. It is roundoff errors and decimal conversions that result in the apparent paradox.

1 Like