I tried another approach where I enabled fma for the host code as well and specified the same fma() instructions in the host code as well.
I still notice that upto 700 iterations the residue value matches upto 14 decimals and instantly after it starts matching only to 9 decimals.
Disabling fma in both cases helps with the “output matching” (although it is a wrong practice like you specified), I just wanted to put forth my observations.
I also made a change to the epsi calculation using the herbie tool
Keep in kind that numerical discrepancies can creep in computation other than the num, den, temp computations that I looked at. When examining numerical differences in the final results, you will have to do an end-to-end analysis. The most likely sources of larger numerical differences are cases of subtractive cancellation.
With rare exceptions involving corner cases, two platforms adhering to the IEEE-754 standard will produce bit-identical results for the same sequence of floating-point operations (and from personal experience, this applies when comparing modern GPUs with SSE/AVX), however most programming languages and their tool chains do not offer sufficiently tight controls to ensure such an identical sequence of operations for non-trivial computation.
While the CUDA compiler is quite conservative in the handling of floating-point arithmetic with the exception of a default setting of -fmad=true, host compilers typically require a “big hammer” flag like -fp-model=strict or /fp:strict to achieve close compliance with IEEE-754, and some may require a whole set of flags to turn off all the optimizations that can interfere with that (gcc comes to mind).
The problem of computing on two different platforms is often described using this old adage: A person with one clock always knows what time it is. A person with two clocks can never be sure.
A practical solution in either case is to use a higher-precision reference (NIST’s atomic clock, in the case of two clocks) to assess the accuracy of any particular device or computation. I highly recommend that approach.
As for use of the Herbie tool, you will need to carefully double check any solution it proposes. It is based on sampling the numerical space using the constraints on variables supplied by the user, and this (frequently, in my experience) leads to sub-optimal and even impractical or non-sensical proposals. While recent versions are aware of FMA, whether Herbie uses it or not seems to be a toss-up, so sometimes it does not use it where it seems obvious that it should. This may be due to phase-ordering issues.
You can call fma functions directly, they are all defined in the cudadevice module.
So with:
use cudadevice
you will have all the fma variants at your disposal (__dfma_r[n|z|u|d], __fma_r[n|z|u|d] ).
Full list is in this pdf (https://www.pgroup.com/lit/literature/pgi-cuf-qrg-2019.pdf)