Surprising IFORT Complex Division Behavior

[Mentzer__Stuart]

I am getting least significant bit differences in ifort complex divisions in a couple of scenarios.

The easy one to demo (see below) gives different results depending on whether the values were set via initialization or assignment. I realize that initializer values are set at compile time but here they report as equal but yet give different division results. This anomaly is fixed with /fp:strict but is still surprising and is not seen with ifx or gfortran.

program division

! Complex division LSB ifort issue demo

! complex division gives less precise result with assigned than initialized variables
! complex variables reported as equal are apparently not equal

! Adding /fp:strict (after other /fp options) fixes this:
!  ifort /O1            division.f90 && division.exe  ->  0.28236835408287586091E-01
!  ifort /O1 /fp:strict division.f90 && division.exe  ->  0.28236835408287589561E-01

! Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.13.1 Build 20240703_000000

complex(8) :: xa, xi = ( -5.90702001956988276d0, 0.0d0 )
complex(8) :: ya, yi = ( -209.195539590677924d0, 0.0d0 )

xa = ( -5.90702001956988276d0, 0.0d0 )
ya = ( -209.195539590677924d0, 0.0d0 )

print *, xi == xa, yi == ya ! Reports that they are equal yet...
print '(G30.20)', real(xi/yi) ! 0.28236835408287589561E-01
print '(G30.20)', real(xa/ya) ! 0.28236835408287586091E-01

! Assignment from initialized variables "fixes" the division
xa = xi; ya = yi
print '(G30.20)', real(xa/ya) ! 0.28236835408287589561E-01

end program

The other scenario, which is more problematic for my actual code, is that division of identical (assigned) complex values is giving slightly different (and less precise) results with ifort than with ifx or gfortran. In this scenario /fp:strict, /fp:source, and other option changes I tried don't make a difference. Maybe there are other options to force or prevent the use of registers that could clean this up.

[JohnNichols]

The error is 2 seconds in the age of the Universe. 

[jimdempseyatthecove]

try /fp:consistent

Jim Dempsey

[Mentzer__Stuart]

For the demo code only /fp:strict eliminates the discrepancy among all the options I tried, including /fp:precise and /fp:consistent.

For the actual application code none of these or anything else I found gets ifort to give the same division result as ifx and gfortran.

And, yes @JohnNichols, this is an LSB difference. But it affects my work and it seems like at least a QoI issue for two variables to show as equal when they clearly aren't, even with options that allow register use for intermediate values.

[jimdempseyatthecove]

In addition to printing out G30.20, what does printing out the Hex values show?

IOW is the "error" related to print formatting as opposed to numerical values?

Jim Dempsey

[Mentzer__Stuart]

If formatting isn't broken then 20 digits should be sufficient to distinguish double precision floats down to the LSB.

Using TRANSFER to put the real part of the complex values into integer(8) and dumping their hex confirms that xa == xi and ya == yi, and that the result of the division differ by 1 bit.

It seems like the issue is probably related to inconsistent register use despite the various /fp options.

[jimdempseyatthecove]

I suggest you walk through the disassembly code.

BTW is your ifort build 32-bit?

ifort 32-bit may use the FPU as opposed to the AVX.

Jim Dempsey

[IanH]

(You can use Z16.16 to dump 64 bit floats.)

Curious - the "wrong" result is actually being pre-calculated by the compiler for the assigned case.

[Steve_Lionel]

That's not astonishing to me - constant folding (compile-time arithmetic) rarely honors options such as /fp and I have seen many cases over the years where constant-folded values aren't the same as those calculated in executable code.

[IanH]

I was a bit surprised that the pre-calculation wasn't done for the initialized case, given the program flow understanding required for the assigned case is really a superset of the initialised case.

[Jacob_Williams]

Gfortran also can have startling results like this. See: https://fortran-lang.discourse.group/t/odd-gfortran-result-with-atan/3590