[CUDA] sycl::sqrt leads to IEEE754 incompatible results on NVidia cards

[anton-v-gorshkov]

Bug Description

Using of sycl::sqrt() function leads to results that are incompatible to IEEE754 on NVidia cards (sqrt.approx.f PTX instruction is generated) even for precise ffp-model. sqrtf() works fine in contrast (sqrt.rn.f PTX instruction is generated).

To Reproduce

1. Go to "sqrt" folder from archive in attachment;
2. Build test application with ./build.sh
3. Run ./sqrt
4. Expected output:

Target device: GeForce RTX 2070
Difference: 1.788139e-07

Correct output should give zero difference, can be checked with sqrtf.

Environment:

• OS: Ubuntu 18.04
• Target device and vendor: NVidia GPUs
• DPC++ version: clang version 13.0.0 (https://github.com/intel/llvm d482ac3)
• Dependencies version: NEO 21.24.20098

Additional Context

Such an issue significantly degrades user experience while CUDA to DPC++ code migration. sycl::sqrt() is a function that normally used by default to compute square roots (e.g. it's generated by DPC++ Compatibility Tool), and observed accuracy issues that this function produces may be really tricky to debug.

sqrt.zip

[zjin-lcf]

I have some questions. How do you generate PTX codes for the example ? Is -gline-tables-only required for debugging ? Thanks.

[anton-v-gorshkov]

I was using NSignt for Compute for get PTX codes. And yes, it's enough to provide "-gline-tables-only" to retrieve this info.

[zjin-lcf]

I added an example that can check the issue. I hope it will not be tricky to debug. Running the HIP, CUDA, OMP programs on AMD and Nvidia GPUs pass.
https://github.com/zjin-lcf/oneAPI-DirectProgramming/tree/master/hellinger-sycl

[npmiller]

So I've been looking into this issue, partly in #5116, but typing out the main points here:

The approx version of sqrt seems to have enough precision to match the SYCL specification requirements and is significantly faster than the regular version, therefore it makes sense to keep that as the default sqrt.

The way sqrt is implemented is by using NVidia's libdevice, which supports changing the precision of sqrt through the -prec-sqrt flag, which makes use of the __nvvm_reflect mechanism. At the moment this is not implemented for sqrt which is why it will always pick the approx variant, and it is also not exposed as an argument to clang.

I have an experimental patch to add a -Xclang -fnvvm-cuda-prec-sqrt flag to raise the precision of the sqrt function here:

• npmiller@061864a

This works as expected and with this flag both of the tests in this ticket that required higher precision pass, and the default is still the faster, lower precision version. The one problem with this approach is that this flag will also affect the native version of sqrt, so it is not possible to use both level of precision in a program.

The main remaining question is how to best expose this feature to users, as far as I understand it -ffp-model=precise is actually already the default, -ffp-model=strict would be the more "precise" one, and it seems to describe a set of more generic floating point value, I'm unsure if it makes sense to place this under that umbrella. As far as I can tell there also doesn't seem to be a SYCL flag that we could use for this, such as -cl-fp32-correctly-rounded-divide-sqrt for OpenCL. So I believe a simple target specific flag may be the best option for this.

Also this patch isn't really SYCL specific and may be better submitted directly upstream, the only SYCL specific part is making sure the nvvm-reflect metadata is also added to the module when running in SYCL mode which we should probably do regardless of the new metadata for sqrt.

[andykaylor]

The main remaining question is how to best expose this feature to users, as far as I understand it -ffp-model=precise is actually already the default, -ffp-model=strict would be the more "precise" one, and it seems to describe a set of more generic floating point value, I'm unsure if it makes sense to place this under that umbrella. As far as I can tell there also doesn't seem to be a SYCL flag that we could use for this, such as -cl-fp32-correctly-rounded-divide-sqrt for OpenCL. So I believe a simple target specific flag may be the best option for this.

We don't want -ffp-model=strict in this case. The 'strict' model is used to enable dynamic rounding mode and strict exception semantics, which aren't supported yet for SYCL offload devices. In theory, -ffp-model=precise should give IEEE-754 results, but we've kind of been cheating a bit with single precision sqrt and divide to take advantage of the SYCL specification allowance for 2.5 ulp error on these operations (SYCL 2020 section 4.17.5, referring to Open CL 1.2 section 7.4). Unfortunately, there is currently no command line option to override this behavior.

I think we need two things: (1) we need a new fp-model that will be the default for device kernels -- I've talked to several people and haven't found a good name for this, but the idea would be that it's similar to fp-model=precise but has adjustments (possibly target/language specific) like the relaxed precision for sqrt and divide, (2) we need a separate command line option to override the precision for individual operataions/functions like sqrt and divide so they can be controlled independent of the fp-model

The second item is probably the priority here so that we can get this working as soon as possible. I've seen requests to provide a mechanism to request more precision for individual math functions (i.e. not just sqrt and divide), so I'd like to have a general mechanism that can be applied to any function. I'd suggest something like this:

-ffp-max-error=[:<comma-separated function list]

So, for the problem we're talking about here, you'd use '-ffp-max-error=0.5:sqrtf'

We'll need to talk about how this can be represented in IR and SPIR-V. There is already an LLVM IR metadata value ('fpmath') that can represent this on individual instructions. I think for SPIR-V we'll need a new decorator. Since we can't assume that 0.5 ulp implementations are available for most functions, I think we'd have to say that this is "requested" accuracy and not "required" but maybe we can also have a setting to produce an error if the requested accuracy is not available.

For the new fp-model, I've mostly thought about this in terms of what we want for SYCL. OpenMP, and OpenCL. In those cases, it would mean the same as fp-model=precise except that 2.5 ulp error would be allowed for single-precision sqrt and divide plus other math library implementations would match the relative error specified by OpenCL 1.2 section 7.4. Otherwise, I know CUDA and HIP want fp-contract to default to 'fast' as opposed to 'on' (which is the default in clang trunk for fp-model=precise). I also am aware that CUDA documents different math library precision expectations than OpenCL/SYCL though I can't remember which is more precise. Host CPU math libraries aren't bound by a specification, but they tend to be accurate to within 1 ulp for most functions. Anyway, my point here is that I'd like to create just one new fp-model option that has target/language specific semantics but basically means "the expected default settings for whatever heterogeneous programming model you're using" rather than having to use "fp-model=[sycl|openmp|opencl|cuda|hip|etc.]"

How does that sound?

[zjin-lcf]

Different compiler options (-ffast-math, -ffp-model, -fp-cuda-prec-sqrt, etc) may be confusing for naive users. The general mechanism you proposed may have a few options for each math function: 'default', "x.x ulp", "lowest", and "highest".
Thanks for the proposal.

[npmiller]

That's interesting, it does make a lot of sense that the device code would need its own separate -fp-model. Although I don't think this part comes into play in this ticket, with the discrepancies you pointed out it may show up in other scenarios.

For your point (2), I've opened yesterday a PR looking into that yesterday, #5141, which essentially mimics the nvcc flag -prec-sqrt for this specific case.

And with regards to the proposal of a flag such as -ffp-max-error=0.5:sqrtf, I'm not sure if that's a great idea, because I feel like the user provided ulp will rarely be respected unless they already know which value to provide for their specific hardware. So it feels pretty confusing, with only a few values actually mattering as thresholds between the different hardware instructions.

I can only really imagine software implementations being able to make use of this, and that seems fairly restrictive for a flag like that.

[andykaylor]

@npmiller You may be right that exposing the option as if a continuous range of ulp values would be meaningful may be misleading. I suppose to some extent I let the implementation detail ("fpmath" metadata semantics) leak into the user interface. I have talked to a few expert floating point developers who suggest a control like this might be useful based on their analysis of their algorithms (without reference to the underlying hardware), but for most users a high/medium/low accuracy option would probably be easier to use.

In wanting to provide this option for the general case (any transcendental function) I'm definitely not thinking about the hardware. The backend can figure that out. In most cases, I agree that there would only be a very few cases where this would be as simple as selecting a single instruction to implement the operation. There are other cases where the compiler can combine specialized hardware instructions with a small refinement sequence to approximate a math function. And then there are cases where we may select a different library implementation.

Obviously that's getting a bit off on a tangent from the starting point here of getting correctly rounded sqrt. I just went that way because it's something that I think is needed and so I'm inclined to provide a general command line mechanism for requesting it (which would also cover sqrt).

That said, I'm not opposed to having a specialized option for requesting precise sqrt and divide. I'll comment on your PR about that.