Re-implement Tiles as OneToMany associator

.github/workflows/clang_format.yml

@@ -19,6 +19,5 @@ jobs:
     - name: Run clang-format style check
       uses: jidicula/[email protected]
       with:
-        clang-format-version: '17'
+        clang-format-version: '18'
         check-path: ${{ matrix.path }}
-        exclude-regex: 'CLUEstering/include/test/doctest.h'

CLUEstering/alpaka/AlpakaCore/AtomicPairCounter.h

@@ -0,0 +1,72 @@
+#ifndef AtomicPairCounter_h
+#define AtomicPairCounter_h
+
+#include <cstdint>
+
+#include <alpaka/alpaka.hpp>
+
+namespace cms::alpakatools {
+
+  class AtomicPairCounter {
+  public:
+    using DoubleWord = uint64_t;
+
+    ALPAKA_FN_HOST_ACC constexpr AtomicPairCounter() : counter_{0} {}
+    ALPAKA_FN_HOST_ACC constexpr AtomicPairCounter(uint32_t first, uint32_t second)
+        : counter_{pack(first, second)} {}
+    ALPAKA_FN_HOST_ACC constexpr AtomicPairCounter(DoubleWord values)
+        : counter_{values} {}
+
+    ALPAKA_FN_HOST_ACC constexpr AtomicPairCounter& operator=(DoubleWord values) {
+      counter_.as_doubleword = values;
+      return *this;
+    }
+
+    struct Counters {
+      uint32_t first;  // in a "One to Many" association is the number of "One"
+      uint32_t
+          second;  // in a "One to Many" association is the total number of associations
+    };
+
+    ALPAKA_FN_HOST_ACC constexpr Counters get() const { return counter_.as_counters; }
+
+    // atomically add as_counters, and return the previous value
+    template <typename TAcc>
+    ALPAKA_FN_ACC ALPAKA_FN_INLINE constexpr Counters add(const TAcc& acc, Counters c) {
+      Packer value{pack(c.first, c.second)};
+      Packer ret{0};
+      ret.as_doubleword = alpaka::atomicAdd(
+          acc, &counter_.as_doubleword, value.as_doubleword, alpaka::hierarchy::Blocks{});
+      return ret.as_counters;
+    }
+
+    // atomically increment first and add i to second, and return the previous value
+    template <typename TAcc>
+    ALPAKA_FN_ACC ALPAKA_FN_INLINE Counters constexpr inc_add(const TAcc& acc,
+                                                              uint32_t i) {
+      return add(acc, {1u, i});
+    }
+
+  private:
+    union Packer {
+      DoubleWord as_doubleword;
+      Counters as_counters;
+      constexpr Packer(DoubleWord _as_doubleword) : as_doubleword(_as_doubleword) { ; };
+      constexpr Packer(Counters _as_counters) : as_counters(_as_counters) { ; };
+    };
+
+    // pack two uint32_t values in a DoubleWord (aka uint64_t)
+    // this is needed because in c++17 a union can only be aggregate-initialised to its first type
+    // it can be probably removed with c++20, and replace with a designated initialiser
+    static constexpr DoubleWord pack(uint32_t first, uint32_t second) {
+      Packer ret{0};
+      ret.as_counters = {first, second};
+      return ret.as_doubleword;
+    }
+
+    Packer counter_;
+  };
+
+}  // namespace cms::alpakatools
+
+#endif  // HeterogeneousCore_AlpakaInterface_interface_AtomicPairCounter_h

CLUEstering/alpaka/AlpakaCore/alpakaWorkDiv.h

@@ -30,6 +30,33 @@ namespace cms::alpakatools {
     return (value + divisor - 1) / divisor;
   }
 
+  // Trait describing whether or not the accelerator expects the threads-per-block and elements-per-thread to be swapped
+  template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
+  struct requires_single_thread_per_block : public std::true_type {};
+
+#ifdef ALPAKA_ACC_GPU_CUDA_ENABLED
+  template <typename TDim>
+  struct requires_single_thread_per_block<alpaka::AccGpuCudaRt<TDim, Idx>>
+      : public std::false_type {};
+#endif  // ALPAKA_ACC_GPU_CUDA_ENABLED
+
+#ifdef ALPAKA_ACC_GPU_HIP_ENABLED
+  template <typename TDim>
+  struct requires_single_thread_per_block<alpaka::AccGpuHipRt<TDim, Idx>>
+      : public std::false_type {};
+#endif  // ALPAKA_ACC_GPU_HIP_ENABLED
+
+#ifdef ALPAKA_ACC_CPU_B_SEQ_T_THREADS_ENABLED
+  template <typename TDim>
+  struct requires_single_thread_per_block<alpaka::AccCpuThreads<TDim, Idx>>
+      : public std::false_type {};
+#endif  // ALPAKA_ACC_CPU_B_SEQ_T_THREADS_ENABLED
+
+  // Whether or not the accelerator expects the threads-per-block and elements-per-thread to be swapped
+  template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
+  inline constexpr bool requires_single_thread_per_block_v =
+      requires_single_thread_per_block<TAcc>::value;
+
   /*
    * Creates the accelerator-dependent workdiv for 1-dimensional operations.
    */
@@ -373,6 +400,19 @@ namespace cms::alpakatools {
            Vec<alpaka::Dim<TAcc>>::zeros();
   }
 
+  /* once_per_block
+   *
+   * `once_per_block(acc)` returns true for a single thread within the block.
+   *
+   * Usually the condition is true for thread 0, but this index should not be relied upon.
+   */
+
+  template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
+  ALPAKA_FN_ACC inline constexpr bool once_per_block(TAcc const& acc) {
+    return alpaka::getIdx<alpaka::Block, alpaka::Threads>(acc) ==
+           Vec<alpaka::Dim<TAcc>>::zeros();
+  }
+
   /*
    * Overload for elementIdxShift = 0
    */

CLUEstering/alpaka/CLUE/CLUEAlgoAlpaka.h

@@ -124,8 +124,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
                                          Queue queue_,
                                          std::size_t block_size) {
     // calculate the number of tiles and their size
-    const auto nTiles{std::ceil(h_points.n / static_cast<float>(pointsPerTile_))};
+    auto nTiles{std::ceil(h_points.n / static_cast<float>(pointsPerTile_))};
     const auto nPerDim{std::ceil(std::pow(nTiles, 1. / Ndim))};
+    nTiles = std::pow(nPerDim, Ndim);
 
     CoordinateExtremes<Ndim> min_max;
     float tile_size[Ndim];
@@ -158,12 +159,55 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         cms::alpakatools::make_host_view(h_points.m_weight.data(), h_points.n));
     alpaka::memset(queue_, *d_seeds, 0x00);
 
-    // Define the working division
+    auto tileIds = cms::alpakatools::make_device_buffer<uint32_t[]>(queue_, h_points.n);
+    // now we scan the dataset and calculate the tile of each point
+    // we can do it on the GPU
     const Idx grid_size = cms::alpakatools::divide_up_by(h_points.n, block_size);
-    const auto working_div = cms::alpakatools::make_workdiv<Acc1D>(grid_size, block_size);
+    const auto work_div = cms::alpakatools::make_workdiv<Acc1D>(grid_size, block_size);
+    alpaka::enqueue(queue_,
+                    alpaka::createTaskKernel<Acc1D>(work_div,
+                                                    KernelScanDatasetTileId{},
+                                                    m_tiles,
+                                                    d_points.view(),
+                                                    tileIds.data(),
+                                                    h_points.n));
+    alpaka::enqueue(queue_,
+                    alpaka::createTaskKernel<Acc1D>(work_div,
+                                                    KernelCalculateOffset{},
+                                                    tileIds.data(),
+                                                    m_tiles->offset(),
+                                                    h_points.n));
+    auto temp = cms::alpakatools::make_device_buffer<uint32_t[]>(queue_, h_points.n);
+    ++nTiles;
+	auto test_grid_size = cms::alpakatools::divide_up_by(nTiles, 32);
+	auto test_work_div = cms::alpakatools::make_workdiv<Acc1D>(test_grid_size, 32);
+    alpaka::enqueue(queue_,
+                    alpaka::createTaskKernel<Acc1D>(test_work_div,
+                                                    KernelOffsetAccumulate{},
+                                                    m_tiles->offset(),
+                                                    temp.data(),
+                                                    nTiles));
+    alpaka::memcpy(queue_,
+                   cms::alpakatools::make_device_view(device, m_tiles->offset(), nTiles),
+                   cms::alpakatools::make_device_view(device, temp.data(), nTiles));
+    alpaka::wait(queue_);
+    alpaka::enqueue(queue_,
+                    alpaka::createTaskKernel<Acc1D>(
+                        work_div, KernelZeroBuffer{}, temp.data(), h_points.n));
+
+    alpaka::enqueue(queue_,
+                    alpaka::createTaskKernel<Acc1D>(work_div,
+                                                    KernelFillAssociator{},
+                                                    tileIds.data(),
+                                                    m_tiles->offset(),
+                                                    temp.data(),
+                                                    m_tiles->content(),
+                                                    h_points.n));
+    alpaka::memset(queue_, (*d_seeds), 0x00);
+
     alpaka::enqueue(queue_,
                     alpaka::createTaskKernel<Acc1D>(
-                        working_div, KernelResetFollowers{}, m_followers, h_points.n));
+                        work_div, KernelResetFollowers{}, m_followers, h_points.n));
   }
 
   // Public methods
@@ -179,10 +223,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
     const Idx grid_size = cms::alpakatools::divide_up_by(h_points.n, block_size);
     auto working_div = cms::alpakatools::make_workdiv<Acc1D>(grid_size, block_size);
-    alpaka::enqueue(
-        queue_,
-        alpaka::createTaskKernel<Acc1D>(
-            working_div, KernelFillTiles{}, d_points.view(), m_tiles, h_points.n));
 
     alpaka::enqueue(queue_,
                     alpaka::createTaskKernel<Acc1D>(working_div,

CLUEstering/alpaka/CLUE/CLUEAlpakaKernels.h

@@ -21,6 +21,81 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
   template <uint8_t Ndim>
   using PointsView = typename PointsAlpaka<Ndim>::PointsAlpakaView;
 
+  struct KernelScanDatasetTileId {
+    template <typename TAcc, uint8_t Ndim>
+    ALPAKA_FN_ACC void operator()(const TAcc& acc,
+                                  TilesAlpaka<Ndim>* d_tiles,
+                                  PointsView<Ndim>* d_points,
+                                  uint32_t* tileIds,
+                                  uint32_t n_points) const {
+      cms::alpakatools::for_each_element_in_grid(acc, n_points, [&](uint32_t i) -> void {
+        tileIds[i] = d_tiles->getGlobalBin(acc, d_points->coords[i]);
+      });
+    }
+  };
+
+  struct KernelCalculateOffset {
+    template <typename TAcc>
+    ALPAKA_FN_ACC void operator()(const TAcc& acc,
+                                  uint32_t* tileIds,
+                                  uint32_t* offsets,
+                                  uint32_t n_points) const {
+      if (cms::alpakatools::once_per_grid(acc)) {
+        offsets[0] = 0;
+      }
+      cms::alpakatools::for_each_element_in_grid(acc, n_points, [&](uint32_t i) -> void {
+        alpaka::atomicAdd(acc, &offsets[tileIds[i] + 1], 1u);
+      });
+    }
+  };
+
+  struct KernelZeroBuffer {
+    template <typename TAcc, typename T>
+    ALPAKA_FN_ACC void operator()(const TAcc& acc, T* buffer, uint32_t size) const {
+      cms::alpakatools::for_each_element_in_grid(
+          acc, size, [&](uint32_t i) -> void { buffer[i] = 0; });
+    }
+  };
+
+  template <typename TAcc>
+  ALPAKA_FN_ACC uint32_t accumulate(const TAcc& acc, const uint32_t* buf, uint32_t size) {
+    uint32_t sum{0};
+    for (uint32_t i{}; i <= size; ++i) {
+      sum += buf[i];
+    }
+
+    return sum;
+  };
+
+  struct KernelOffsetAccumulate {
+    template <typename TAcc>
+    ALPAKA_FN_ACC void operator()(const TAcc& acc,
+                                  const uint32_t* offset,
+                                  uint32_t* temp,
+                                  uint32_t n_tiles) const {
+      cms::alpakatools::for_each_element_in_grid(
+          acc, n_tiles, [&](uint32_t tile) -> void {
+            temp[tile] = accumulate(acc, offset, tile);
+          });
+	}
+  };
+
+  struct KernelFillAssociator {
+    template <typename TAcc>
+    ALPAKA_FN_ACC void operator()(const TAcc& acc,
+                                  const uint32_t* tileIds,
+                                  const uint32_t* offsets,
+                                  uint32_t* temp,
+                                  uint32_t* content,
+                                  uint32_t n_points) const {
+      cms::alpakatools::for_each_element_in_grid(acc, n_points, [&](uint32_t i) -> void {
+        const auto tileId{tileIds[i]};
+        const auto contentPosition{alpaka::atomicAdd(acc, &temp[tileId], 1u)};
+        content[offsets[tileId] + contentPosition] = i;
+      });
+    }
+  };
+
   struct KernelResetTiles {
     template <typename TAcc, uint8_t Ndim>
     ALPAKA_FN_ACC void operator()(TAcc const& acc,
@@ -31,7 +106,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         tiles->resizeTiles(nTiles, nPerDim);
       }
       cms::alpakatools::for_each_element_in_grid(
-          acc, nTiles, [&](uint32_t i) -> void { tiles->clear(i); });
+          acc, nTiles, [&](uint32_t i) -> void { /*tiles->clear(i);*/
+                                                 ;
+          });
     }
   };
 
@@ -45,17 +122,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     }
   };
 
-  struct KernelFillTiles {
-    template <typename TAcc, uint8_t Ndim>
-    ALPAKA_FN_ACC void operator()(const TAcc& acc,
-                                  PointsView<Ndim>* points,
-                                  TilesAlpaka<Ndim>* tiles,
-                                  uint32_t n_points) const {
-      cms::alpakatools::for_each_element_in_grid(
-          acc, n_points, [&](uint32_t i) { tiles->fill(acc, points->coords[i], i); });
-    }
-  };
-
   template <typename TAcc, uint8_t Ndim, uint8_t N_, typename KernelType>
   ALPAKA_FN_HOST_ACC void for_recursion(
       const TAcc& acc,
@@ -70,9 +136,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
       float dc,
       uint32_t point_id) {
     if constexpr (N_ == 0) {
-      int binId{tiles->getGlobalBinByBin(acc, base_vec)};
+      auto binId{tiles->getGlobalBinByBin(acc, base_vec)};
       // get the size of this bin
-      int binSize{static_cast<int>((*tiles)[binId].size())};
+      auto binSize{static_cast<int>((*tiles)[binId].size())};
 
       // iterate inside this bin
       for (int binIter{}; binIter < binSize; ++binIter) {
@@ -172,9 +238,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
       float dm_sq,
       uint32_t point_id) {
     if constexpr (N_ == 0) {
-      int binId{tiles->getGlobalBinByBin(acc, base_vec)};
+      auto binId{tiles->getGlobalBinByBin(acc, base_vec)};
       // get the size of this bin
-      int binSize{(*tiles)[binId].size()};
+      auto binSize{(*tiles)[binId].size()};
 
       // iterate inside this bin
       for (int binIter{}; binIter < binSize; ++binIter) {

CLUEstering/alpaka/DataFormats/Points.h

@@ -11,15 +11,15 @@ struct Points {
   Points() = default;
   Points(const std::vector<VecArray<float, Ndim>>& coords,
          const std::vector<float>& weight)
-      : m_coords{coords}, m_weight{weight}, n{weight.size()} {
+      : m_coords{coords}, m_weight{weight}, n{static_cast<uint32_t>(weight.size())} {
     m_rho.resize(n);
     m_delta.resize(n);
     m_nearestHigher.resize(n);
     m_clusterIndex.resize(n);
     m_isSeed.resize(n);
   }
   Points(const std::vector<std::vector<float>>& coords, const std::vector<float>& weight)
-      : m_weight{weight}, n{weight.size()} {
+      : m_weight{weight}, n{static_cast<uint32_t>(weight.size())} {
     for (const auto& x : coords) {
       VecArray<float, Ndim> temp_vecarray;
       for (auto value : x) {
@@ -43,7 +43,7 @@ struct Points {
   std::vector<int> m_clusterIndex;
   std::vector<int> m_isSeed;
 
-  size_t n;
+  uint32_t n;
 };
 
 #endif