Implicit Type Conversion -> Precision

[Johannes_Rieke]

Hi all,

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

[Steven_L_Intel1]

Intel Fortran doesn't currently have such a warning.

[mecej4]

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.

[Johannes_Rieke]

Hi mecej4,

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