[SYCL] Add support for SYCL math functions

added 2 commits

fbf7189b - 1 commit from branch libeigen:master
7af65e4a - Use math functions in the SYCL namespace

I understand the pragmatic reasons for this, but I am not a fan of this MR for the following reason: After this is merged, any contributor will have to reason about whether to use EIGEN_USING_STD or EIGEN_USING_MATH_STD. It is bad enough that we need EIGEN_USING_STD in the first place, because cuda/HIP does/did not implement C++11.

@cantonios @chuckyschluz what do you think?

Regardless what is decided, I think these math functions should point to numext whenever possible, and only invoke the EIGEN_USING_ macros in the numext namespace. This wouldn't be too confusing if EIGEN_USING_MATH_STD was confined to MathFunctions.h.

Good idea.

Another option is to force sycl to specialize numext::* for all their types. it might not be pretty (could be in a separate file), but at least it wouldn't pollute the rest of the code.

It is bad enough that we need EIGEN_USING_STD in the first place, because cuda/HIP does/did not implement C++11.

Could you elaborate on what you mean ?

Like any other library, CUDA is not allowed by the C++ standard to define any functions in the std namespace.

Another option is to force sycl to specialize numext::* for all their types.

SYCL does not use or define any special types, the SYCL math functions operate on float, double, etc.

I mean that most standard math functions are defined in the std namespace since C++11, see https://en.cppreference.com/w/cpp/header/cmath

OK, yes. But, as far as I know, only the C++ Standard Library can define functions in the std namespace.

So, if CUDA or SYCL need to define their own alternative implementations, they cannot define them inside std. At least, that's how I understand things...

They could certainly implement the C++ Standard Library for their device in the std namespace, in full or in parts, and I believe they have done so partially in recent years. I don't recall the reason for not doing so originally, but perhaps having to comply with the standard in terms of error handling was imposing constraints that would harm performance.

I can think of a few reasons (compatibility with the standard library used by the host compiler, licensing, avoiding ambiguity, etc.) but of course I do not know what the actual ones were.

It's also kind of beside the point... as a user of CUDA, HIP, SYCL and Eigen, and need a way for them to play well together :-)

I think these math functions should point to numext whenever possible

Is this suggestion to replace all uses (for example) of

    using std::abs;

or

    EIGEN_USING_STD(pow);

with explicit use of numext::abs and numext::pow ?

Yes, that would be the best solution, I guess, although I still regret the GPU platforms not implementing the standard library.

Yes, if we use numext::abs everywhere, we only need to use that macro once. And contributors wouldn’t have to worry about being compatible with sycl, custom scale types, and so on. That’s my hope anyway.

Looking at numext, why do you have mini and maxi instead of min and max ?

I think it's a legacy thing trying to avoid clashing with headers on older systems that define macros for "min" and "max".

I think using numext::* everywhere makes a lot of sense. SYCL could use SYCL_DEVICE_ONLY in numext:: to decide when to call mathematical functions in the sycl namespace or in std::

I tried the suggestion of replacing everything with numext:: and it worked out of the box (at least for my application) for SYCL, but broke the CUDA and HIP implementation. As soon as I have time I'll look into the error.

What's the status of this MR?

@AuroraPerego any update?

Hi, my last attempt was not successful and then I stopped working on this topic for a while. I'll look into it again and update you in the coming weeks.

Hi @cantonios, I managed to replace using std::* and EIGEN_USING_STD(*) with methods from the numext namespace. I ran into an error related to pow_impl (MathFunctions.h, line 539) because pow was called with a double for the base and an integer for the exponent (StableNorm.h, line 143). This led to an error with the result_type due to ScalarBinaryOpTraits not having a specialization for ScalarX=double and ScalarY=int. I solved it with std::common_type<ScalarX, ScalarY>::type, but let me know if you prefer another solution.

Do the other types have a specialization?

The specializations are:

ScalarBinaryOpTraits<T, T>
ScalarBinaryOpTraits<T, typename NumTraits<std::enable_if_t<NumTraits<T>::IsComplex, T>>::Real>
ScalarBinaryOpTraits<typename NumTraits<std::enable_if_t<NumTraits<T>::IsComplex, T>>::Real, T>
ScalarBinaryOpTraits<T, void>
ScalarBinaryOpTraits<void, T>
ScalarBinaryOpTraits<void, void>