diff --git a/HeterogeneousCore/AlpakaMath/BuildFile.xml b/HeterogeneousCore/AlpakaMath/BuildFile.xml new file mode 100644 index 0000000000000..06905828f6ea6 --- /dev/null +++ b/HeterogeneousCore/AlpakaMath/BuildFile.xml @@ -0,0 +1,4 @@ + + + + diff --git a/HeterogeneousCore/AlpakaMath/interface/deltaPhi.h b/HeterogeneousCore/AlpakaMath/interface/deltaPhi.h new file mode 100644 index 0000000000000..85c90fca79144 --- /dev/null +++ b/HeterogeneousCore/AlpakaMath/interface/deltaPhi.h @@ -0,0 +1,35 @@ +#ifndef HeterogeneousCore_AlpakaMath_deltaPhi_h +#define HeterogeneousCore_AlpakaMath_deltaPhi_h + +#include + +namespace cms::alpakatools { + + // reduce to [-pi,pi] + template + ALPAKA_FN_HOST_ACC inline T reduceRange(TAcc const& acc, T x) { + constexpr T o2pi = T{1.} / (T{2.} * std::numbers::pi_v); + if (alpaka::math::abs(acc, x) <= std::numbers::pi_v) + return x; + T n = alpaka::math::round(acc, x * o2pi); + return x - n * T{2.} * std::numbers::pi_v; + } + + template + ALPAKA_FN_HOST_ACC inline T phi(TAcc const& acc, T x, T y) { + return reduceRange(acc, std::numbers::pi_v + alpaka::math::atan2(acc, -y, -x)); + } + + template + ALPAKA_FN_HOST_ACC 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 + ALPAKA_FN_HOST_ACC inline T deltaPhi(TAcc const& acc, T phi1, T phi2) { + return reduceRange(acc, phi1 - phi2); + } + +} // namespace cms::alpakatools + +#endif diff --git a/HeterogeneousCore/AlpakaMath/test/BuildFile.xml b/HeterogeneousCore/AlpakaMath/test/BuildFile.xml new file mode 100644 index 0000000000000..2276193f8d1dd --- /dev/null +++ b/HeterogeneousCore/AlpakaMath/test/BuildFile.xml @@ -0,0 +1,6 @@ + + + + + + diff --git a/HeterogeneousCore/AlpakaMath/test/alpaka/testDeltaPhi.dev.cc b/HeterogeneousCore/AlpakaMath/test/alpaka/testDeltaPhi.dev.cc new file mode 100644 index 0000000000000..021e1d9850d7b --- /dev/null +++ b/HeterogeneousCore/AlpakaMath/test/alpaka/testDeltaPhi.dev.cc @@ -0,0 +1,95 @@ +#define CATCH_CONFIG_MAIN +#include +#include +#include + +#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; + +struct phiFuncsUnitTestsKernel { + template + ALPAKA_FN_ACC void operator()(TAcc const& acc, T* out) const { + // Unit circle typical values + out[0] = phi(acc, 1.0, 0.0); // x = 1.0, y = 0.0 => phi = 0 + out[1] = phi(acc, 0.0, 1.0); // x = 0.0, y = 1.0 => phi = pi/2 + out[2] = phi(acc, -1.0, 0.0); // x = -1.0, y = 0.0 => phi = pi + out[3] = phi(acc, 0.0, -1.0); // x = 0.0, y = -1.0 => phi = -pi/2 + out[4] = phi(acc, 0.7071, 0.7071); // x = sqrt(2)/2, y = sqrt(2)/2 => phi = pi/4 + out[5] = phi(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 + // Phi from unit circle + out[6] = deltaPhi(acc, 1.0, 0.0, 0.0, -1.0); // 3pi/2 - 0 = -pi/2 + out[7] = deltaPhi(acc, 0.0, 1.0, 0.0, -1.0); // 3pi/2 - pi/2 = pi + out[8] = deltaPhi(acc, 0.0, -1.0, 0.0, 1.0); // pi/2 - 3pi/2 = -pi + out[9] = deltaPhi(acc, 0.7071, -0.7071, -0.7071, -0.7071); // -3pi/4 - (-pi/4) = -pi/2 + + // Calculation directly from phi + out[10] = deltaPhi(acc, 3. * M_PI / 2., 0.); // 3pi/2 - 0 = -pi/2 + out[11] = deltaPhi(acc, 3. * M_PI / 2., M_PI / 2.); // 3pi/2 - pi/2 = pi + out[12] = deltaPhi(acc, M_PI / 2., 3. * M_PI / 2.); // pi/2 - 3pi/2 = -pi + out[13] = deltaPhi(acc, -3. * M_PI / 4., -M_PI / 4.); // -3pi/4 - (-pi/4) = -pi/2 + } +}; + +template +void testPhiFuncs(uint32_t size, std::vector const& res) { + // get the list of devices on the current platform + auto const& devices = cms::alpakatools::devices(); + if (devices.empty()) { + FAIL("No devices available for the " EDM_STRINGIZE(ALPAKA_ACCELERATOR_NAMESPACE) " backend, " + "the test will be skipped."); + } + + for (auto const& device : devices) { + std::cout << "...on " << alpaka::getName(device) << "\n"; + Queue queue(device); + + auto c_h = make_host_buffer(queue, size); + alpaka::memset(queue, c_h, 0.); + auto c_d = make_device_buffer(queue, size); + alpaka::memset(queue, c_d, 0.); + + alpaka::exec(queue, 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) { + CHECK_THAT(c_h.data()[i], Catch::Matchers::WithinAbs(res[i], eps)); + } + } +} + +TEST_CASE("Standard checks alpaka phi functions for the relevant data types (float and double) and for all backends") { + std::vector 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(size, res); + } + + SECTION("Tests for float data type") { + std::cout << "Testing phi functions for float data type...\n"; + testPhiFuncs(size, res); + } +}