Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
compute_histogram.hpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
10#pragma once
11
18
19#include "shambase/exception.hpp"
20#include "shambase/overloaded.hpp"
21#include "shambase/string.hpp"
22#include "shamalgs/ImplVariant.hpp"
23#include "shambackends/Device.hpp"
24#include "shambackends/DeviceBuffer.hpp"
25#include "shambackends/DeviceScheduler.hpp"
26#include "shambackends/kernel_call.hpp"
27#include "shamcomm/logs.hpp"
28#include <shambackends/sycl.hpp>
29#include <optional>
30#include <stdexcept>
31#include <tuple>
32#include <utility>
33#include <vector>
34
35namespace shamalgs::primitives {
36
37 namespace impl {
38
40 struct Reference {
41 static constexpr std::string_view variant_type_name = "reference";
42 };
43
44 struct NaiveGpu {
45 static constexpr std::string_view variant_type_name = "naive_gpu";
46 };
47
48 struct GpuTeamFetching {
49 static constexpr std::string_view variant_type_name = "gpu_team_fetching";
50 };
51
52 struct GpuOversubscribe {
53 static constexpr std::string_view variant_type_name = "gpu_oversubscribe";
54 };
55
56 inline shamalgs::ImplVariantGlobal<Reference, NaiveGpu, GpuTeamFetching, GpuOversubscribe>
57 compute_histogram_impl{[](const sham::DeviceScheduler_ptr &dev_sched, auto &self) {
58 if (dev_sched->ctx->device->prop.type == sham::DeviceType::GPU) {
59 self.set(GpuOversubscribe{});
60 } else {
61 self.set(NaiveGpu{}); // it is portable and fast everywhere
62 }
63 }};
64
66 inline std::vector<std::string> get_default_impl_list_compute_histogram() {
67 return compute_histogram_impl.get_default_config_list();
68 }
69
71 inline std::string get_current_impl_compute_histogram() {
72 return compute_histogram_impl.get_current_config();
73 }
74
76 inline bool is_impl_set_compute_histogram() { return compute_histogram_impl.is_set(); }
77
79 inline void set_impl_compute_histogram(const std::string &impl) {
80 shamlog_info_ln("algs", "setting compute_histogram implementation to impl :", impl);
81 compute_histogram_impl.set(impl);
82 }
83
85 inline void autoselect_impl_compute_histogram(const sham::DeviceScheduler_ptr &dev_sched) {
86 compute_histogram_impl.autoselect(dev_sched);
87 shamlog_info_ln(
88 "algs",
89 "defaulting compute_histogram implementation to impl :",
90 get_current_impl_compute_histogram());
91 }
92
93 template<class T, class Tbins, class... Targs, class Tfunctor>
94 inline void compute_histogram_reference(
95 const sham::DeviceBuffer<Tbins> &bin_edge_inf,
96 const sham::DeviceBuffer<Tbins> &bin_edge_sup,
97 size_t nbins,
98 size_t element_count,
99 Tfunctor &&functor,
100 sham::DeviceBuffer<T> &result,
101 const sham::DeviceBuffer<Targs> &...input_data) {
102
103 auto result_vec = result.copy_to_stdvec();
104
105 auto cpu_basic_impl = [&](const std::vector<Tbins> &bin_edge_inf,
106 const std::vector<Tbins> &bin_edge_sup,
107 const std::vector<Targs> &...in_data,
108 std::vector<T> &result) {
109 for (size_t ibin = 0; ibin < nbins; ibin++) {
110 Tbins edge_inf = bin_edge_inf[ibin];
111 Tbins edge_sup = bin_edge_sup[ibin];
112
113 T accumulator = 0;
114
115 for (size_t i = 0; i < element_count; i++) {
116 bool has_value = false;
117 auto tmp = functor(edge_inf, edge_sup, in_data[i]..., has_value);
118 if (has_value) {
119 accumulator += tmp;
120 }
121 }
122
123 result[ibin] = accumulator;
124 }
125 };
126
127 cpu_basic_impl(
128 bin_edge_inf.copy_to_stdvec(),
129 bin_edge_sup.copy_to_stdvec(),
130 input_data.copy_to_stdvec()...,
131 result_vec);
132
133 result.copy_from_stdvec(result_vec);
134 }
135
136 template<class T, class Tbins, class... Targs, class Tfunctor>
137 inline void compute_histogram_naive_gpu(
138 const sham::DeviceScheduler_ptr &dev_sched,
139 const sham::DeviceBuffer<Tbins> &bin_edge_inf,
140 const sham::DeviceBuffer<Tbins> &bin_edge_sup,
141 size_t nbins,
142 size_t element_count,
143 Tfunctor &&functor,
144 sham::DeviceBuffer<T> &result,
145 const sham::DeviceBuffer<Targs> &...input_data) {
146
147 sham::kernel_call(
148 dev_sched->get_queue(),
149 sham::MultiRef{bin_edge_inf, bin_edge_sup, input_data...},
150 sham::MultiRef{result},
151 nbins,
152 [element_count, functor](
153 u32 ibin,
154 const Tbins *__restrict bin_edge_inf,
155 const Tbins *__restrict bin_edge_sup,
156 const Targs *__restrict... in_data,
157 T *__restrict result) {
158 Tbins edge_inf = bin_edge_inf[ibin];
159 Tbins edge_sup = bin_edge_sup[ibin];
160
161 T accumulator = 0;
162
163 for (size_t i = 0; i < element_count; i++) {
164 bool has_value = false;
165 T tmp = functor(edge_inf, edge_sup, in_data[i]..., has_value);
166 if (has_value) {
167 accumulator += tmp;
168 }
169 }
170
171 result[ibin] = accumulator;
172 });
173 }
174
175 template<class T, class Tbins, class... Targs, class Tfunctor>
176 inline void compute_histogram_gpu_team_fetching(
177 const sham::DeviceScheduler_ptr &dev_sched,
178 const sham::DeviceBuffer<Tbins> &bin_edge_inf,
179 const sham::DeviceBuffer<Tbins> &bin_edge_sup,
180 size_t nbins,
181 size_t element_count,
182 Tfunctor &&functor,
183 sham::DeviceBuffer<T> &result,
184 const sham::DeviceBuffer<Targs> &...input_data) {
185
186 sham::kernel_call_hndl(
187 dev_sched->get_queue(),
188 sham::MultiRef{bin_edge_inf, bin_edge_sup, input_data...},
189 sham::MultiRef{result},
190 nbins,
191 [element_count, functor](
192 u32 nbins,
193 const Tbins *__restrict bin_edge_inf,
194 const Tbins *__restrict bin_edge_sup,
195 const Targs *__restrict... in_data,
196 T *__restrict result) {
197 return [=, in_data = std::tuple{in_data...}](sycl::handler &cgh) {
198 u32 group_size = 128;
199 u32 group_cnt = shambase::group_count(nbins, group_size);
200
201 // roundup to next multiple of 4
202 group_cnt = (group_cnt + 3) / 4 * 4;
203 u32 corrected_len = group_cnt * group_size;
204
205 auto locals
206 = sycl::local_accessor<std::tuple<Targs...>, 1>(group_size, cgh);
207
208 cgh.parallel_for(
209 sycl::nd_range<1>{corrected_len, group_size},
210 [=](sycl::nd_item<1> item) {
211 u32 local_id = item.get_local_id(0);
212 u32 group_tile_id = item.get_group_linear_id();
213 u32 ibin = group_tile_id * group_size + local_id;
214
215 bool is_valid_point = (ibin < nbins);
216 Tbins edge_inf = is_valid_point ? bin_edge_inf[ibin] : Tbins{};
217 Tbins edge_sup = is_valid_point ? bin_edge_sup[ibin] : Tbins{};
218
219 T local_sum = 0;
220
221 for (size_t i = 0; i < element_count; i += group_size) {
222
223 item.barrier(sycl::access::fence_space::local_space);
224
225 if (i + local_id < element_count) {
226 std::apply(
227 [&](auto &...in_data) {
228 locals[local_id]
229 = std::tuple{in_data[i + local_id]...};
230 },
231 in_data);
232 }
233
234 item.barrier(sycl::access::fence_space::local_space);
235
236 if (is_valid_point) {
237 for (size_t lane = 0; lane < group_size; lane++) {
238 if (i + lane >= element_count) {
239 continue;
240 }
241 bool has_value = false;
242 T tmp = std::apply(
243 [&](auto &...local_accs) {
244 return functor(
245 edge_inf,
246 edge_sup,
247 local_accs...,
248 has_value);
249 },
250 locals[lane]);
251 if (has_value) {
252 local_sum += tmp;
253 }
254 }
255 }
256
257 item.barrier(sycl::access::fence_space::local_space);
258 }
259
260 if (is_valid_point) {
261 result[ibin] = local_sum;
262 }
263 });
264 };
265 });
266 }
267
268 template<class T, class Tbins, class... Targs, class Tfunctor>
269 inline void compute_histogram_gpu_oversubscribe(
270 const sham::DeviceScheduler_ptr &dev_sched,
271 u32 group_size,
272 const sham::DeviceBuffer<Tbins> &bin_edge_inf,
273 const sham::DeviceBuffer<Tbins> &bin_edge_sup,
274 size_t nbins,
275 size_t element_count,
276 Tfunctor &&functor,
277 sham::DeviceBuffer<T> &result,
278 const sham::DeviceBuffer<Targs> &...input_data) {
279
280 sham::kernel_call_hndl(
281 dev_sched->get_queue(),
282 sham::MultiRef{bin_edge_inf, bin_edge_sup, input_data...},
283 sham::MultiRef{result},
284 nbins * group_size,
285 [element_count, functor, group_size, nbins](
286 u32 nbins_oversubscribed,
287 const Tbins *__restrict bin_edge_inf,
288 const Tbins *__restrict bin_edge_sup,
289 const Targs *__restrict... in_data,
290 T *__restrict result) {
291 return [=, in_data = std::tuple{in_data...}](sycl::handler &cgh) {
292 u32 group_cnt = shambase::group_count(nbins_oversubscribed, group_size);
293
294 // roundup to next multiple of 4
295 group_cnt = (group_cnt + 3) / 4 * 4;
296
297 u32 corrected_len = group_cnt * group_size;
298
299 cgh.parallel_for(
300 sycl::nd_range<1>{corrected_len, group_size},
301 [=](sycl::nd_item<1> item) {
302 u32 local_id = item.get_local_id(0);
303 u32 ibin = item.get_group_linear_id();
304
305 bool is_valid_point = (ibin < nbins);
306 Tbins edge_inf = is_valid_point ? bin_edge_inf[ibin] : Tbins{};
307 Tbins edge_sup = is_valid_point ? bin_edge_sup[ibin] : Tbins{};
308
309 // for each thread this will the sum of all the
310 // "func(in_data[group_size*i + local_id]) for all i"
311 T local_sum = 0;
312
313 for (size_t i = 0; i < element_count; i += group_size) {
314
315 if (i + local_id < element_count) {
316
317 bool has_value = false;
318
319 // coalesced read of the data and then
320 // compute the value to accumulate
321 T tmp = std::apply(
322 [&](auto &...in_data) {
323 return functor(
324 edge_inf,
325 edge_sup,
326 in_data[i + local_id]...,
327 has_value);
328 },
329 in_data);
330
331 if (has_value) {
332 // add it to the local sum of this thread
333 local_sum += tmp;
334 }
335 }
336 }
337
338 // we have all the terms scattered across the threads of the group,
339 // we can just accumulate the result
340 auto group_sum = sycl::reduce_over_group(
341 item.get_group(), local_sum, sycl::plus<T>{});
342
343 if (is_valid_point && local_id == 0) {
344 result[ibin] = group_sum;
345 }
346 });
347 };
348 });
349 }
350
351 } // namespace impl
352
353 template<class T, class Tbins, class... Targs, class Tfunctor>
354 inline sham::DeviceBuffer<T> compute_histogram(
355 const sham::DeviceScheduler_ptr &dev_sched,
356 const sham::DeviceBuffer<Tbins> &bin_edge_inf,
357 const sham::DeviceBuffer<Tbins> &bin_edge_sup,
358 size_t element_count,
359 Tfunctor &&functor,
360 const sham::DeviceBuffer<Targs> &...input_data) {
361
362 using namespace impl;
363
364 size_t nbins = bin_edge_inf.get_size();
365
366 if (nbins != bin_edge_sup.get_size()) {
367 shambase::make_except_with_loc<std::invalid_argument>(
368 "bin_edge_inf and bin_edge_sup must have the same size");
369 }
370
371 sham::DeviceBuffer<T> result(nbins, dev_sched);
372
373 if (!impl::compute_histogram_impl.is_set()) {
374 impl::autoselect_impl_compute_histogram(dev_sched);
375 }
376
377 std::visit(
378 shambase::overloaded{
379 [&](impl::Reference) {
380 compute_histogram_reference(
381 bin_edge_inf,
382 bin_edge_sup,
383 nbins,
384 element_count,
385 std::forward<Tfunctor>(functor),
386 result,
387 input_data...);
388 },
389 [&](impl::NaiveGpu) {
390 compute_histogram_naive_gpu(
391 dev_sched,
392 bin_edge_inf,
393 bin_edge_sup,
394 nbins,
395 element_count,
396 std::forward<Tfunctor>(functor),
397 result,
398 input_data...);
399 },
400 [&](impl::GpuTeamFetching) {
401 compute_histogram_gpu_team_fetching(
402 dev_sched,
403 bin_edge_inf,
404 bin_edge_sup,
405 nbins,
406 element_count,
407 std::forward<Tfunctor>(functor),
408 result,
409 input_data...);
410 },
411 [&](impl::GpuOversubscribe) {
412 compute_histogram_gpu_oversubscribe(
413 dev_sched,
414 256,
415 bin_edge_inf,
416 bin_edge_sup,
417 nbins,
418 element_count,
419 std::forward<Tfunctor>(functor),
420 result,
421 input_data...);
422 },
423 },
424 impl::compute_histogram_impl.get());
425
426 return result;
427 }
428
429 template<class T>
430 inline sham::DeviceBuffer<T> compute_histogram_basic(
431 const sham::DeviceScheduler_ptr &dev_sched,
432 const sham::DeviceBuffer<T> &bin_edge_inf,
433 const sham::DeviceBuffer<T> &bin_edge_sup,
434 const sham::DeviceBuffer<T> &positions) {
435
436 size_t element_count = positions.get_size();
437
438 return compute_histogram<T>(
439 dev_sched,
440 bin_edge_inf,
441 bin_edge_sup,
442 element_count,
443 [](const T &bin_edge_inf, const T &bin_edge_sup, const T &position, bool &has_value) {
444 has_value = position >= bin_edge_inf && position < bin_edge_sup;
445 return has_value ? 1 : 0;
446 },
447 positions);
448 }
449
450} // namespace shamalgs::primitives
ImplVariant.hpp
Generic std::variant-based implementation selector.
u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27
sham::DeviceBuffer
A buffer allocated in USM (Unified Shared Memory).
Definition DeviceBuffer.hpp:106
sham::DeviceBuffer::copy_from_stdvec
void copy_from_stdvec(const std::vector< T > &vec)
Copy the content of a std::vector into the buffer.
Definition DeviceBuffer.hpp:666
sham::DeviceBuffer::copy_to_stdvec
std::vector< T > copy_to_stdvec() const
Copy the content of the buffer to a std::vector.
Definition DeviceBuffer.hpp:510
sham::DeviceBuffer::get_size
size_t get_size() const
Gets the number of elements in the buffer.
Definition DeviceBuffer.hpp:476
shamalgs::ImplVariantGlobal
Drop-in replacement for the hand-rolled "global variable + enum + name mapping.
Definition ImplVariant.hpp:265
exception.hpp
This header file contains utility functions related to exception handling in the code.
kernel_call.hpp
sham::kernel_call
void kernel_call(sham::DeviceQueue &q, RefIn in, RefOut in_out, u32 n, Functor &&func, SourceLocation &&callsite=SourceLocation{})
Submit a kernel to a SYCL queue.
Definition kernel_call.hpp:344
shamalgs::primitives::impl
namespace to control implementation behavior
Definition compute_histogram.hpp:37
shamalgs::primitives::impl::autoselect_impl_compute_histogram
void autoselect_impl_compute_histogram(const sham::DeviceScheduler_ptr &dev_sched)
Select the default implementation for compute_histogram.
Definition compute_histogram.hpp:85
shamalgs::primitives::impl::set_impl_compute_histogram
void set_impl_compute_histogram(const std::string &impl)
Set the implementation for compute_histogram.
Definition compute_histogram.hpp:79
shamalgs::primitives::impl::get_current_impl_compute_histogram
std::string get_current_impl_compute_histogram()
Get the current implementation for compute_histogram.
Definition compute_histogram.hpp:71
shamalgs::primitives::impl::get_default_impl_list_compute_histogram
std::vector< std::string > get_default_impl_list_compute_histogram()
Get list of available compute_histogram implementations.
Definition compute_histogram.hpp:66
shamalgs::primitives::impl::is_impl_set_compute_histogram
bool is_impl_set_compute_histogram()
Check if an implementation has been selected for compute_histogram.
Definition compute_histogram.hpp:76
shamalgs::primitives
namespace for primitive algorithm (e.g. sort, scan, reductions, ...)
Definition append_subset_to.hpp:22
shambase::make_except_with_loc
ExcptTypes make_except_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Create an exception with a message and a location.
Definition exception.hpp:105
sham::MultiRef
A class that references multiple buffers or similar objects.
Definition MultiRef.hpp:33
shamalgs::primitives::impl::GpuOversubscribe
GPU kernel oversubscribing a work-group per bin, reducing locally.
Definition compute_histogram.hpp:52
shamalgs::primitives::impl::GpuTeamFetching
GPU kernel using a team-local cache to coalesce input reads.
Definition compute_histogram.hpp:48
shamalgs::primitives::impl::NaiveGpu
Portable GPU kernel, one work-item per bin.
Definition compute_histogram.hpp:44
shamalgs::primitives::impl::Reference
CPU reference implementation, computed on a host copy of the buffers.
Definition compute_histogram.hpp:40
shambase::overloaded
Build an overload set out of several callables, for use with std::visit.
Definition overloaded.hpp:35
Generated on for Shamrock by doxygen 1.15.0