Incorrect division result when using double precision types on intel platform

psing51 · ‎03-02-2018

Dear All,

I am trying to run the following fortran 77 code on "Xeon E5-2690 v4 @ 2.60GHz", Suse linux 12. The result of division is incorrect when double precision dataypes are involved-

      program main
      real*8 h10
      real*8 rsdata
      real*4 div_const
      div_const=9.8
      h10=13720
      rsdata=h10/div_const
      print *,h10,9.8,rsdata
      end

user@machine:~/TESTCASE> ifort test.f
user@machine:~/TESTCASE> ./a.out
   13720.0000000000        9.800000       1399.99997275216

I tried writing fortran90 version of the same code -

program main

!for 64 bit real
!integer, parameter :: dp = selected_real_kind(15, 307)
!integer, parameter :: dp = kind(1.d0)
!real(kind=dp) :: h10
!real(kind=dp) :: rsdata


!real(kind=8) :: h10
!real(kind=8) :: rsdata


double precision :: h10
double precision :: rsdata


real(kind=4)  :: div_const
h10=13720.00
div_const=9.800
rsdata=h10/div_const
print *,h10,9.8,rsdata

end

But the result is still same

user@machine:~/TESTCASE> ifort test.f90
user@machine:~/TESTCASE> ./a.out
   13720.0000000000        9.800000       1399.99997275216

I was expecting 1400 as the output(9.8*1400=13720). I have confirmed that when i declare any of h10/rsdata as real-32 bit type, i get correct result.

Is there a way to make fortran 77 version work without involving real(32 bit) datatype?

mecej4 · ‎03-02-2018

Binary floating point arithmetic differs in many ways from decimal arithmetic. Those differences can be learned from introductory books on computer programming.

Try the following program and reason out why the answers are different.

print *,(4/3-1)*3-1
print *,(4.0/3.0-1)*3.0-1.0                                                                                                           
print *,(4d0/3d0-1d0)*3d0-1d0                                                                                                         
end

Johannes_Rieke · ‎03-02-2018

Hi,

please have a look into the following sources:

https://software.intel.com/sites/default/files/managed/de/a1/FP_Consistency_091117.pdf

https://software.intel.com/en-us/node/522979

In your code you divide a double by a single (line 7 respectively line 21) and implicit type conversion is performed by the compiler. Digits 8 to 15 of div_const are not necessarily filled with zeros. Thus, a deviation compared to exact calculus can happen. If you make div_const also double, you get the correct results. Do you need to mix single and double?

program main
  !for 64 bit real
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  real(kind=dp)      :: h10, dummy
  real(kind=dp)      :: rsdata
  real(kind=dp)      :: div_const
  h10      = 13720.00_dp
  div_const= 9.800_dp
  dummy    = 1400.0_dp * div_const
  rsdata=h10/real(div_const,dp)
  print *,h10,div_const,rsdata, dummy
end

You can override the single real definitions by compiler flag '-real_size 64' and convert literals (9.8 -> 9.8_dp) to get your code working. For this you can also use compiler flag -fpconstant.

If you stay completely in the single world (e.g. test dp = selected_real_kind(6, 30)), you get also the desired results.