Port alpaka phi functions to HeterogeneousCore

HeterogeneousCore/AlpakaMath/BuildFile.xml

@@ -0,0 +1,4 @@
+<use name="alpaka"/>
+<export>
+  <lib name="1"/>
+</export>

HeterogeneousCore/AlpakaMath/interface/deltaPhi.h

@@ -0,0 +1,36 @@
+#ifndef HeterogeneousCore_AlpakaMath_deltaPhi_h
+#define HeterogeneousCore_AlpakaMath_deltaPhi_h
+
+#include <alpaka/alpaka.hpp>
+
+namespace cms::alpakatools {
+
+  // reduce to [-pi,pi]
+  template <typename TAcc, typename T>
+  requires std::is_floating_point_v<T>
+  ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE T reduceRange(TAcc const& acc, T x) {
+    constexpr T o2pi = 1. / (2. * M_PI);
+    if (alpaka::math::abs(acc, x) <= T(M_PI))
+      return x;
+    T n = alpaka::math::round(acc, x * o2pi);
+    return x - n * T(2. * M_PI);
+  }
+
+  template <typename TAcc, typename T>
+  ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE T phi(TAcc const& acc, T x, T y) {
+    return reduceRange(acc, M_PI + alpaka::math::atan2(acc, -y, -x));
+  }
+
+  template <typename TAcc, typename T>
+  ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE T deltaPhi(TAcc const& acc, T x1, T y1, T x2, T y2) {
+    return reduceRange(acc, alpaka::math::atan2(acc, -y2, -x2) - alpaka::math::atan2(acc, -y1, -x1));
+  }
+
+  template <typename TAcc, typename T>
+  ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE T deltaPhi(TAcc const& acc, T phi1, T phi2) {
+    return reduceRange(acc, phi1 - phi2);
+  }
+
+}  // namespace cms::alpakatools
+
+#endif

HeterogeneousCore/AlpakaMath/test/BuildFile.xml

@@ -0,0 +1,6 @@
+<bin name="alpakaTestDeltaPhi" file="alpaka/testDeltaPhi.dev.cc">
+  <use name="alpaka"/>
+  <use name="catch2"/>
+  <use name="HeterogeneousCore/AlpakaInterface"/>
+  <flags ALPAKA_BACKENDS="1"/>
+</bin>

HeterogeneousCore/AlpakaMath/test/alpaka/testDeltaPhi.dev.cc

@@ -0,0 +1,86 @@
+#define CATCH_CONFIG_MAIN
+#include <catch.hpp>
+#include <vector>
+
+#include "HeterogeneousCore/AlpakaInterface/interface/config.h"
+#include "HeterogeneousCore/AlpakaInterface/interface/memory.h"
+#include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
+#include "HeterogeneousCore/AlpakaMath/interface/deltaPhi.h"
+
+using namespace cms::alpakatools;
+using namespace ALPAKA_ACCELERATOR_NAMESPACE;
+
+template <typename T>
+void ok(T a, T b, T eps) { assert(std::abs(a - b) < eps); }
+
+struct phiFuncsUnitTestsKernel {
+  template <typename TAcc, typename T>
+  ALPAKA_FN_ACC void operator()(TAcc const& acc, T* out) const {
+
+    // Unit circle typical values
+    out[0] = phi<TAcc, T>(acc, 1.0, 0.0); // x = 1.0, y = 0.0 => phi = 0
+    out[1] = phi<TAcc, T>(acc, 0.0, 1.0); // x = 0.0, y = 1.0 => phi = pi/2
+    out[2] = phi<TAcc, T>(acc, -1.0, 0.0); // x = -1.0, y = 0.0 => phi = pi
+    out[3] = phi<TAcc, T>(acc, 0.0, -1.0); // x = 0.0, y = -1.0 => phi = -pi/2
+    out[4] = phi<TAcc, T>(acc, 0.7071, 0.7071); // x = sqrt(2)/2, y = sqrt(2)/2 => phi = pi/4
+    out[5] = phi<TAcc, T>(acc, -0.7071, -0.7071); // x = sqrt(2)/2, y = sqrt(2)/2 => phi = -3pi/4
+
+    // Making sure that delta phi is within [-pi, pi] range
+    out[6] = deltaPhi<TAcc, T>(acc, 1.0, 0.0, 0.0, -1.0); // phi from unit circle, 3pi/2 - 0 = -pi/2
+    out[7] = deltaPhi<TAcc, T>(acc, 0.0, 1.0, 0.0, -1.0); // phi from unit circle, 3pi/2 - pi/2 = pi
+    out[8] = deltaPhi<TAcc, T>(acc, 0.0, -1.0, 0.0, 1.0); // phi from unit circle, pi/2 - 3pi/2 = -pi
+    out[9] = deltaPhi<TAcc, T>(acc, 0.7071, -0.7071, -0.7071, -0.7071); // phi from unit circle, -3pi/4 - (-pi/4) = -pi/2
+
+    out[10] = deltaPhi<TAcc, T>(acc, 3.*M_PI/2., 0.); // calculation directly from phi, 3pi/2 - 0 = -pi/2
+    out[11] = deltaPhi<TAcc, T>(acc, 3.*M_PI/2., M_PI/2.); // calculation directly from phi, 3pi/2 - pi/2 = pi
+    out[12] = deltaPhi<TAcc, T>(acc, M_PI/2., 3.*M_PI/2.); // calculation directly from phi, pi/2 - 3pi/2 = -pi
+    out[13] = deltaPhi<TAcc, T>(acc, -3.*M_PI/4., -M_PI/4.); // calculation directly from phi, -3pi/4 - (-pi/4) = -pi/2
+  }
+};
+
+template <typename T>
+void testPhiFuncs(uint32_t size, std::vector<double> const& res) {
+  // get the list of devices on the current platform
+  auto const& devices = cms::alpakatools::devices<Platform>();
+  if (devices.empty()) {
+    std::cerr << "No devices available for the " EDM_STRINGIZE(ALPAKA_ACCELERATOR_NAMESPACE) " backend, "
+      "the test will be skipped.\n";
+    exit(EXIT_FAILURE);
+  }
+
+  for (auto const& device : devices) {
+    std::cout << "...on " << alpaka::getName(device) << "\n";
+    Queue queue(device);
+
+    auto c_h = make_host_buffer<T[]>(queue, size);
+    alpaka::memset(queue, c_h, 0.);
+    auto c_d = make_device_buffer<T[]>(queue, size);
+    alpaka::memset(queue, c_d, 0.);
+
+    alpaka::enqueue(queue, alpaka::createTaskKernel<Acc1D>(WorkDiv1D{1u, 1u, 1u}, phiFuncsUnitTestsKernel(), c_d.data()));
+    alpaka::memcpy(queue, c_h, c_d);
+    alpaka::wait(queue);
+
+    constexpr T eps = 1.e-5;
+    for (size_t i = 0; i < size; ++i) {
+      ok(c_h.data()[i], static_cast<T>(res[i]), eps);
+    }
+  }
+}
+
+TEST_CASE("Standard checks alpaka phi functions for the relevant data types (float and double) and for all backends") {
+  std::vector<double> res = {0.0, M_PI/2., M_PI, -M_PI/2., M_PI/4., -3.*M_PI/4.,
+                            -M_PI/2., M_PI, -M_PI, -M_PI/2.,
+                            -M_PI/2., M_PI, -M_PI, -M_PI/2.}; // Expected results
+  uint32_t size = res.size(); // Number of tests
+
+  SECTION("Tests for double data type") {
+    std::cout << "Testing phi functions for double data type...\n";
+    testPhiFuncs<double>(size, res);
+  }
+
+  SECTION("Tests for float data type") {
+    std::cout << "Testing phi functions for float data type...\n";
+    testPhiFuncs<float>(size, res);
+  }
+}