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
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?
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
please have a look into the following sources:
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.