Enable fast math library in CUDA compilation

One of the biggest performance bottlenecks of our code right now are the
math functions we use from the standard library. To illustrate the
problem, these functions have the problem that they contain a _lot_ of
edge case handling code for NaN values, infinities, etc. which are not
necessary for our use case. This causes an incredible amount of code to
be added and causes a huge amount of instruction cache misses when
though CUDA contains native intrinsics of high accurary which model the
same behaviour. This PR serves to illustrate the massive performance
gain that we can achieve by using intrinsic math rather than the math
provided by the C++ standard library.

Author: stephenswat

cmake/traccc-compiler-options-cuda.cmake

@@ -25,7 +25,7 @@ endif()
 
 # Allow to use functions in device code that are constexpr, even if they are
 # not marked with __device__.
-traccc_add_flag( CMAKE_CUDA_FLAGS "--expt-relaxed-constexpr" )
+traccc_add_flag( CMAKE_CUDA_FLAGS "--expt-relaxed-constexpr --use_fast_math" )
 
 # Make CUDA generate debug symbols for the device code as well in a debug
 # build.

core/include/traccc/definitions/math.hpp

@@ -7,6 +7,8 @@
 
 #pragma once
 
+#include "traccc/definitions/qualifiers.hpp"
+
 // SYCL include(s).
 #if defined(CL_SYCL_LANGUAGE_VERSION) || defined(SYCL_LANGUAGE_VERSION)
 #include <sycl/sycl.hpp>
@@ -15,13 +17,63 @@
 // System include(s).
 #include <cmath>
 
-namespace traccc {
+namespace traccc::math {
 
-/// Namespace to pick up math functions from
 #if defined(CL_SYCL_LANGUAGE_VERSION) || defined(SYCL_LANGUAGE_VERSION)
-namespace math = ::sycl;
+using ::sycl::abs;
+using ::sycl::acos;
+using ::sycl::asin;
+using ::sycl::atan;
+using ::sycl::atan2;
+using ::sycl::cos;
+using ::sycl::exp;
+using ::sycl::fabs;
+using ::sycl::floor;
+using ::sycl::fmod;
+using ::sycl::log;
+using ::sycl::max;
+using ::sycl::min;
+using ::sycl::pow;
+using ::sycl::sin;
+using ::sycl::sqrt;
+using ::sycl::tan;
 #else
-namespace math = std;
-#endif  // SYCL
+using std::abs;
+using std::acos;
+using std::asin;
+using std::atan;
+using std::atan2;
+using std::cos;
+using std::exp;
+using std::fabs;
+using std::floor;
+using std::fmod;
+using std::log;
+using std::max;
+using std::min;
+using std::pow;
+using std::sin;
+using std::sqrt;
+using std::tan;
+#endif
 
-}  // namespace traccc
+/**
+ * @brief Perform IEEE-754 division, even if fast math is enabled.
+ *
+ * @returns x divided by y
+ */
+template <typename T>
+TRACCC_HOST_DEVICE inline __attribute__((always_inline)) T div_ieee754(T x,
+                                                                       T y) {
+    static_assert(std::is_same_v<T, double> || std::is_same_v<T, float>);
+#ifdef __CUDA_ARCH__
+    if constexpr (std::is_same_v<T, double>) {
+        return __ddiv_rn(x, y);
+    } else {
+        return __fdiv_rn(x, y);
+    }
+#else
+    return x / y;
+#endif
+}
+}  // namespace traccc::math

device/cuda/src/ambiguity_resolution/greedy_ambiguity_resolution_algorithm.cu

@@ -25,6 +25,7 @@
 #include "./kernels/sort_tracks_per_measurement.cuh"
 #include "./kernels/sort_updated_tracks.cuh"
 #include "traccc/cuda/ambiguity_resolution/greedy_ambiguity_resolution_algorithm.hpp"
+#include "traccc/definitions/math.hpp"
 
 // Thrust include(s).
 #include <thrust/execution_policy.h>
@@ -42,7 +43,8 @@ namespace traccc::cuda {
 struct devide_op {
     TRACCC_HOST_DEVICE
     traccc::scalar operator()(unsigned int a, unsigned int b) const {
-        return static_cast<traccc::scalar>(a) / static_cast<traccc::scalar>(b);
+        return math::div_ieee754(static_cast<traccc::scalar>(a),
+                                 static_cast<traccc::scalar>(b));
     }
 };
 

device/cuda/src/ambiguity_resolution/kernels/remove_tracks.cu

@@ -6,6 +6,7 @@
  */
 
 // Project include(s).
+#include "traccc/definitions/math.hpp"
 #include "traccc/utils/pair.hpp"
 
 // Local include(s).
@@ -433,8 +434,9 @@ __launch_bounds__(512) __global__
             updated_tracks[pos2] = tid;
             is_updated[tid] = 1;
 
-            rel_shared.at(tid) = static_cast<traccc::scalar>(n_shared.at(tid)) /
-                                 static_cast<traccc::scalar>(n_meas.at(tid));
+            rel_shared.at(tid) =
+                math::div_ieee754(static_cast<traccc::scalar>(n_shared.at(tid)),
+                                  static_cast<traccc::scalar>(n_meas.at(tid)));
         }
     }
 }