In an huge numerical simulation code, written by someone else, in many code lines 'implicit type conversion' between integer, real32 and real64 is used. During some tests I encountered a precision loss in intermediate data (a rigid body rotation of an input mesh) and I traced it down to a possible loss due to implicit type conversion.
To have a closer look on what happens I built a simple generic example, which uses some calculations similar to the simulation code:
program Implicit_Type_Conversion use, intrinsic :: ISO_FORTRAN_ENV, only : i8 => int64, r4 => real32, r8 => real64, r16 => real128, output_unit implicit none ! Variables integer(i8) :: i, j integer(i8), parameter :: leng = 10000 real(r4 ), dimension(leng) :: vec_r4 real(r8 ), dimension(leng) :: vec_r8 real(r16), dimension(leng) :: vec_r16 real(r16) :: sum_of_values_i8 real(r16) :: sum_of_values_r4 real(r16) :: sum_of_values_r8 real(r16) :: sum_of_values_r16 ! Body of Implicit_Type_Conversion write(output_unit,*) 'Implicit Type Conversion -> Impact on Precision' ! i8 -> r16 sum_of_values_i8 = 0.0_r16 do i = 1, leng ! implicit type conversion i8 to r16 sum_of_values_i8 = sum_of_values_i8 + sind((i*i)/2.0_r16) end do ! r4 -> r16 sum_of_values_r4 = 0.0_r16 do i = 1, leng ! implicit type conversion r4 to r16 sum_of_values_r4 = sum_of_values_r4 + sind((real(i,r4)*real(i,r4))/2.0_r16) end do ! r8 -> r16 sum_of_values_r8 = 0.0_r16 do i = 1, leng ! implicit type conversion r8 to r16 sum_of_values_r8 = sum_of_values_r8 + sind((real(i,r8)*real(i,r8))/2.0_r16) end do ! r16 = r16 sum_of_values_r16 = 0.0_r16 do i = 1, leng ! no implicit type conversion sum_of_values_r16 = sum_of_values_r16 + sind((real(i,r16)*real(i,r16))/2.0_r16) end do write(output_unit,'(2ES30.20)') sum_of_values_i8 , sum_of_values_r16-sum_of_values_i8 write(output_unit,'(2ES30.20)') sum_of_values_r4 , sum_of_values_r16-sum_of_values_r4 write(output_unit,'(2ES30.20)') sum_of_values_r8 , sum_of_values_r16-sum_of_values_r8 write(output_unit,'(2ES30.20)') sum_of_values_r16, sum_of_values_r16-sum_of_values_r16 end program Implicit_Type_Conversion
The results are:
Implicit Type Conversion -> Impact on Precision 3.69082119314211831229E+02 0.00000000000000000000E+00 3.69113651696722043610E+02 -3.15323825102123810941E-02 3.69082119314211831229E+02 0.00000000000000000000E+00 3.69082119314211831229E+02 0.00000000000000000000E+00
I used PSXE 16.0.2, options: /nologo /debug:full /Od /warn:interfaces /Qinit:snan /fp:source /module:"x64\Debug\\" /object:"x64\Debug\\" /Fd"x64\Debug\vc140.pdb" /traceback /check:bounds /check:stack /libs:dll /threads /dbglibs /c
I would personally not use this mixture of different types, but it would be very interesting to know, why especially the second case (r4 -> r16) differs from the other that much.
gfortran supports the '-Wconversion'/'-Wconversion-extra', which help to find out, where a type conversion is happening. Is it possible to get a 'precision loss' warning? A warning where you might loose precision would be more useful IMHO. Is there some flag for the Intel compiler? Is it possible to do this at all?
Best regards, Johannes
Here is a partial explanation for the second result: 32-bit reals have a 23+1 bit mantissa. When i = 10000, the expression i*i has the value 100000000. However, with only 24 bits available, this number as well as all the integers from 99999996 to 100000004 have the same internal representation as 32-bit reals.
Furthermore, the non-standard function SIND takes a 32-bit real as argument and returns a 32-bit real. In effect, then, you are evaluating the sine function with an argument that may be in error by + or - 2 degrees, and I hope that you will agree that 2 degrees out of 90 is plenty of error.
Depending on what you wish to expose with your test program, you may need to use the functions DSIND and QSIND for the other cases.
you're absolutely right. The test case wasn't adequate to test the precision impact of implicit type conversion, due to reaching the precision limit of real32. The difference comes from the sensitivity of sin (or the above used non-standard sind) to the last digits. As you wrote in real32 99999996 to 100000004 cannot be distinguished...
Reducing leng to 1000 all calculations result in the same result. --> implicit type conversion does not seem itself to have an impact on the precision. Other not shown cases could be covered by /fpconstant or by adjusting the code... hopefully.
Nevertheless, I intentionally used the sin (respectively sind) and did not use sin/dsin/qsin. Since I divide through a quadruple precision two firstly (2.0_r16, e.g. line 23) qsin should be used in every case.
Many thanks again.
@Steve and the Intel team: If there is a feature request list for something similar to gfortran's -Wconversion/ Wconversion-extra, please put my name on that list.
Best regards, Johannes