EBGeometry
Compact, header-only C++ library for fast evaluation of signed distance functions
Loading...
Searching...
No Matches
EBGeometry_SIMD.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2026 Robert Marskar <robert.marskar@sintef.no>
2//
3// SPDX-License-Identifier: GPL-3.0-or-later
4
16#ifndef EBGEOMETRY_SIMD_HPP
17#define EBGEOMETRY_SIMD_HPP
18
19// Std includes
20#include <cmath>
21#include <cstddef>
22#include <type_traits>
23
24#if (defined(__SSE4_1__) || defined(__AVX__) || defined(__AVX512F__)) && !defined(EBGEOMETRY_DEVICE_COMPILE)
25#include <immintrin.h>
26#endif
27
28// Our includes
29#include "EBGeometry_GPU.hpp"
30
31namespace EBGeometry {
32
36namespace SIMD {
37
46template <class T>
47[[nodiscard]] EBGEOMETRY_HOST_DEVICE
48inline T
49select(const bool a_mask, const T a_ifTrue, const T a_ifFalse) noexcept
50{
51 return a_mask ? a_ifTrue : a_ifFalse;
52}
53
60[[nodiscard]] EBGEOMETRY_HOST_DEVICE
61inline bool
62both(const bool a_lhs, const bool a_rhs) noexcept
63{
64 return a_lhs && a_rhs;
65}
66
72[[nodiscard]] EBGEOMETRY_HOST_DEVICE
73inline float
74sqrt(const float a_x) noexcept
75{
76 return std::sqrt(a_x);
77}
78
84[[nodiscard]] EBGEOMETRY_HOST_DEVICE
85inline double
86sqrt(const double a_x) noexcept
87{
88 return std::sqrt(a_x);
89}
90
96[[nodiscard]] EBGEOMETRY_HOST_DEVICE
97inline float
98abs(const float a_x) noexcept
99{
100 return std::abs(a_x);
101}
102
108[[nodiscard]] EBGEOMETRY_HOST_DEVICE
109inline double
110abs(const double a_x) noexcept
111{
112 return std::abs(a_x);
113}
114
122template <class T>
123[[nodiscard]] EBGEOMETRY_HOST_DEVICE
124constexpr size_t
125WidestPack() noexcept
126{
127 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SIMD: T must be float or double");
128
129#if defined(EBGEOMETRY_DEVICE_COMPILE)
130 return 1;
131#elif defined(__AVX512F__)
132 return 64 / sizeof(T);
133#elif defined(__AVX__)
134 return 32 / sizeof(T);
135#elif defined(__SSE4_1__)
136 return 16 / sizeof(T);
137#else
138 return 1;
139#endif
140}
141
149template <class T, size_t W>
150[[nodiscard]] EBGEOMETRY_HOST_DEVICE
151constexpr size_t
152NativeWidth() noexcept
153{
154 // The narrowest register is 16 bytes (SSE), so no pack is narrower than that.
155 for (size_t n = WidestPack<T>(); n * sizeof(T) >= 16; n /= 2) {
156 if (W % n == 0) {
157 return n;
158 }
159 }
160
161 return 1;
162}
163
164template <class T, size_t N>
165struct Pack;
166
171template <class V>
172struct Scalar
173{
177 using type = V;
178};
179
185template <class T, size_t N>
186struct Scalar<Pack<T, N>>
187{
191 using type = T;
192};
193
205template <class T, size_t N>
206struct Pack;
207
208#if !defined(EBGEOMETRY_DEVICE_COMPILE)
209
210#if defined(__SSE4_1__)
211
215struct MaskF4
216{
220 __m128 m;
221
223 friend inline MaskF4
224 both(const MaskF4 a_lhs, const MaskF4 a_rhs) noexcept
225 {
226 return MaskF4{_mm_and_ps(a_lhs.m, a_rhs.m)};
227 }
229};
230
234struct MaskD2
235{
239 __m128d m;
240
242 friend inline MaskD2
243 both(const MaskD2 a_lhs, const MaskD2 a_rhs) noexcept
244 {
245 return MaskD2{_mm_and_pd(a_lhs.m, a_rhs.m)};
246 }
248};
249
253template <>
254struct Pack<float, 4>
255{
259 __m128 v;
260
264 Pack() noexcept = default;
265
270 explicit Pack(const __m128 a_v) noexcept : v(a_v)
271 {}
272
277 explicit Pack(const float a_x) noexcept : v(_mm_set1_ps(a_x))
278 {}
279
285 static Pack
286 load(const float* a_p) noexcept
287 {
288 return Pack(_mm_load_ps(a_p));
289 }
290
295 void
296 store(float* a_p) const noexcept
297 {
298 _mm_store_ps(a_p, v);
299 }
300
301 // clang-format off
303 friend Pack operator+(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_add_ps(a_lhs.v, a_rhs.v)); }
304 friend Pack operator-(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_sub_ps(a_lhs.v, a_rhs.v)); }
305 friend Pack operator*(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_mul_ps(a_lhs.v, a_rhs.v)); }
306 friend Pack operator/(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_div_ps(a_lhs.v, a_rhs.v)); }
307 friend MaskF4 operator<(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF4{_mm_cmplt_ps(a_lhs.v, a_rhs.v)}; }
308 friend MaskF4 operator>(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF4{_mm_cmpgt_ps(a_lhs.v, a_rhs.v)}; }
309 friend MaskF4 operator>=(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF4{_mm_cmpge_ps(a_lhs.v, a_rhs.v)}; }
310 friend Pack select(const MaskF4 a_mask, const Pack a_ifTrue, const Pack a_ifFalse) noexcept { return Pack(_mm_blendv_ps(a_ifFalse.v, a_ifTrue.v, a_mask.m)); }
311 friend Pack sqrt(const Pack a_x) noexcept { return Pack(_mm_sqrt_ps(a_x.v)); }
312 friend Pack abs(const Pack a_x) noexcept { return Pack(_mm_andnot_ps(_mm_set1_ps(-0.0F), a_x.v)); }
314 // clang-format on
315};
316
320template <>
321struct Pack<double, 2>
322{
326 __m128d v;
327
331 Pack() noexcept = default;
332
337 explicit Pack(const __m128d a_v) noexcept : v(a_v)
338 {}
339
344 explicit Pack(const double a_x) noexcept : v(_mm_set1_pd(a_x))
345 {}
346
352 static Pack
353 load(const double* a_p) noexcept
354 {
355 return Pack(_mm_load_pd(a_p));
356 }
357
362 void
363 store(double* a_p) const noexcept
364 {
365 _mm_store_pd(a_p, v);
366 }
367
368 // clang-format off
370 friend Pack operator+(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_add_pd(a_lhs.v, a_rhs.v)); }
371 friend Pack operator-(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_sub_pd(a_lhs.v, a_rhs.v)); }
372 friend Pack operator*(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_mul_pd(a_lhs.v, a_rhs.v)); }
373 friend Pack operator/(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm_div_pd(a_lhs.v, a_rhs.v)); }
374 friend MaskD2 operator<(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD2{_mm_cmplt_pd(a_lhs.v, a_rhs.v)}; }
375 friend MaskD2 operator>(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD2{_mm_cmpgt_pd(a_lhs.v, a_rhs.v)}; }
376 friend MaskD2 operator>=(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD2{_mm_cmpge_pd(a_lhs.v, a_rhs.v)}; }
377 friend Pack select(const MaskD2 a_mask, const Pack a_ifTrue, const Pack a_ifFalse) noexcept { return Pack(_mm_blendv_pd(a_ifFalse.v, a_ifTrue.v, a_mask.m)); }
378 friend Pack sqrt(const Pack a_x) noexcept { return Pack(_mm_sqrt_pd(a_x.v)); }
379 friend Pack abs(const Pack a_x) noexcept { return Pack(_mm_andnot_pd(_mm_set1_pd(-0.0), a_x.v)); }
381 // clang-format on
382};
383
384#endif
385
386#if defined(__AVX__)
387
391struct MaskF8
392{
396 __m256 m;
397
399 friend inline MaskF8
400 both(const MaskF8 a_lhs, const MaskF8 a_rhs) noexcept
401 {
402 return MaskF8{_mm256_and_ps(a_lhs.m, a_rhs.m)};
403 }
405};
406
410struct MaskD4
411{
415 __m256d m;
416
418 friend inline MaskD4
419 both(const MaskD4 a_lhs, const MaskD4 a_rhs) noexcept
420 {
421 return MaskD4{_mm256_and_pd(a_lhs.m, a_rhs.m)};
422 }
424};
425
429template <>
430struct Pack<float, 8>
431{
435 __m256 v;
436
440 Pack() noexcept = default;
441
446 explicit Pack(const __m256 a_v) noexcept : v(a_v)
447 {}
448
453 explicit Pack(const float a_x) noexcept : v(_mm256_set1_ps(a_x))
454 {}
455
461 static Pack
462 load(const float* a_p) noexcept
463 {
464 return Pack(_mm256_load_ps(a_p));
465 }
466
471 void
472 store(float* a_p) const noexcept
473 {
474 _mm256_store_ps(a_p, v);
475 }
476
477 // clang-format off
479 friend Pack operator+(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_add_ps(a_lhs.v, a_rhs.v)); }
480 friend Pack operator-(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_sub_ps(a_lhs.v, a_rhs.v)); }
481 friend Pack operator*(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_mul_ps(a_lhs.v, a_rhs.v)); }
482 friend Pack operator/(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_div_ps(a_lhs.v, a_rhs.v)); }
483 friend MaskF8 operator<(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF8{_mm256_cmp_ps(a_lhs.v, a_rhs.v, _CMP_LT_OQ)}; }
484 friend MaskF8 operator>(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF8{_mm256_cmp_ps(a_lhs.v, a_rhs.v, _CMP_GT_OQ)}; }
485 friend MaskF8 operator>=(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF8{_mm256_cmp_ps(a_lhs.v, a_rhs.v, _CMP_GE_OQ)}; }
486 friend Pack select(const MaskF8 a_mask, const Pack a_ifTrue, const Pack a_ifFalse) noexcept { return Pack(_mm256_blendv_ps(a_ifFalse.v, a_ifTrue.v, a_mask.m)); }
487 friend Pack sqrt(const Pack a_x) noexcept { return Pack(_mm256_sqrt_ps(a_x.v)); }
488 friend Pack abs(const Pack a_x) noexcept { return Pack(_mm256_andnot_ps(_mm256_set1_ps(-0.0F), a_x.v)); }
490 // clang-format on
491};
492
496template <>
497struct Pack<double, 4>
498{
502 __m256d v;
503
507 Pack() noexcept = default;
508
513 explicit Pack(const __m256d a_v) noexcept : v(a_v)
514 {}
515
520 explicit Pack(const double a_x) noexcept : v(_mm256_set1_pd(a_x))
521 {}
522
528 static Pack
529 load(const double* a_p) noexcept
530 {
531 return Pack(_mm256_load_pd(a_p));
532 }
533
538 void
539 store(double* a_p) const noexcept
540 {
541 _mm256_store_pd(a_p, v);
542 }
543
544 // clang-format off
546 friend Pack operator+(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_add_pd(a_lhs.v, a_rhs.v)); }
547 friend Pack operator-(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_sub_pd(a_lhs.v, a_rhs.v)); }
548 friend Pack operator*(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_mul_pd(a_lhs.v, a_rhs.v)); }
549 friend Pack operator/(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm256_div_pd(a_lhs.v, a_rhs.v)); }
550 friend MaskD4 operator<(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD4{_mm256_cmp_pd(a_lhs.v, a_rhs.v, _CMP_LT_OQ)}; }
551 friend MaskD4 operator>(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD4{_mm256_cmp_pd(a_lhs.v, a_rhs.v, _CMP_GT_OQ)}; }
552 friend MaskD4 operator>=(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD4{_mm256_cmp_pd(a_lhs.v, a_rhs.v, _CMP_GE_OQ)}; }
553 friend Pack select(const MaskD4 a_mask, const Pack a_ifTrue, const Pack a_ifFalse) noexcept { return Pack(_mm256_blendv_pd(a_ifFalse.v, a_ifTrue.v, a_mask.m)); }
554 friend Pack sqrt(const Pack a_x) noexcept { return Pack(_mm256_sqrt_pd(a_x.v)); }
555 friend Pack abs(const Pack a_x) noexcept { return Pack(_mm256_andnot_pd(_mm256_set1_pd(-0.0), a_x.v)); }
557 // clang-format on
558};
559
560#endif
561
562#if defined(__AVX512F__)
563
567struct MaskF16
568{
572 __mmask16 m;
573
575 friend inline MaskF16
576 both(const MaskF16 a_lhs, const MaskF16 a_rhs) noexcept
577 {
578 return MaskF16{static_cast<__mmask16>(a_lhs.m & a_rhs.m)};
579 }
581};
582
586struct MaskD8
587{
591 __mmask8 m;
592
594 friend inline MaskD8
595 both(const MaskD8 a_lhs, const MaskD8 a_rhs) noexcept
596 {
597 return MaskD8{static_cast<__mmask8>(a_lhs.m & a_rhs.m)};
598 }
600};
601
607template <>
608struct Pack<float, 16>
609{
613 __m512 v;
614
618 Pack() noexcept = default;
619
624 explicit Pack(const __m512 a_v) noexcept : v(a_v)
625 {}
626
631 explicit Pack(const float a_x) noexcept : v(_mm512_set1_ps(a_x))
632 {}
633
639 static Pack
640 load(const float* a_p) noexcept
641 {
642 return Pack(_mm512_load_ps(a_p));
643 }
644
649 void
650 store(float* a_p) const noexcept
651 {
652 _mm512_store_ps(a_p, v);
653 }
654
655 // clang-format off
657 friend Pack operator+(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_add_ps(a_lhs.v, a_rhs.v)); }
658 friend Pack operator-(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_sub_ps(a_lhs.v, a_rhs.v)); }
659 friend Pack operator*(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_mul_ps(a_lhs.v, a_rhs.v)); }
660 friend Pack operator/(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_div_ps(a_lhs.v, a_rhs.v)); }
661 friend MaskF16 operator<(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF16{_mm512_cmp_ps_mask(a_lhs.v, a_rhs.v, _CMP_LT_OQ)}; }
662 friend MaskF16 operator>(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF16{_mm512_cmp_ps_mask(a_lhs.v, a_rhs.v, _CMP_GT_OQ)}; }
663 friend MaskF16 operator>=(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskF16{_mm512_cmp_ps_mask(a_lhs.v, a_rhs.v, _CMP_GE_OQ)}; }
664 friend Pack select(const MaskF16 a_mask, const Pack a_ifTrue, const Pack a_ifFalse) noexcept { return Pack(_mm512_mask_blend_ps(a_mask.m, a_ifFalse.v, a_ifTrue.v)); }
665 friend Pack sqrt(const Pack a_x) noexcept { return Pack(_mm512_sqrt_ps(a_x.v)); }
666 friend Pack abs(const Pack a_x) noexcept
667 {
668 return Pack(_mm512_castsi512_ps(_mm512_andnot_si512(_mm512_castps_si512(_mm512_set1_ps(-0.0F)), _mm512_castps_si512(a_x.v))));
669 }
671 // clang-format on
672};
673
678template <>
679struct Pack<double, 8>
680{
684 __m512d v;
685
689 Pack() noexcept = default;
690
695 explicit Pack(const __m512d a_v) noexcept : v(a_v)
696 {}
697
702 explicit Pack(const double a_x) noexcept : v(_mm512_set1_pd(a_x))
703 {}
704
710 static Pack
711 load(const double* a_p) noexcept
712 {
713 return Pack(_mm512_load_pd(a_p));
714 }
715
720 void
721 store(double* a_p) const noexcept
722 {
723 _mm512_store_pd(a_p, v);
724 }
725
726 // clang-format off
728 friend Pack operator+(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_add_pd(a_lhs.v, a_rhs.v)); }
729 friend Pack operator-(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_sub_pd(a_lhs.v, a_rhs.v)); }
730 friend Pack operator*(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_mul_pd(a_lhs.v, a_rhs.v)); }
731 friend Pack operator/(const Pack a_lhs, const Pack a_rhs) noexcept { return Pack(_mm512_div_pd(a_lhs.v, a_rhs.v)); }
732 friend MaskD8 operator<(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD8{_mm512_cmp_pd_mask(a_lhs.v, a_rhs.v, _CMP_LT_OQ)}; }
733 friend MaskD8 operator>(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD8{_mm512_cmp_pd_mask(a_lhs.v, a_rhs.v, _CMP_GT_OQ)}; }
734 friend MaskD8 operator>=(const Pack a_lhs, const Pack a_rhs) noexcept { return MaskD8{_mm512_cmp_pd_mask(a_lhs.v, a_rhs.v, _CMP_GE_OQ)}; }
735 friend Pack select(const MaskD8 a_mask, const Pack a_ifTrue, const Pack a_ifFalse) noexcept { return Pack(_mm512_mask_blend_pd(a_mask.m, a_ifFalse.v, a_ifTrue.v)); }
736 friend Pack sqrt(const Pack a_x) noexcept { return Pack(_mm512_sqrt_pd(a_x.v)); }
737 friend Pack abs(const Pack a_x) noexcept
738 {
739 return Pack(_mm512_castsi512_pd(_mm512_andnot_si512(_mm512_castpd_si512(_mm512_set1_pd(-0.0)), _mm512_castpd_si512(a_x.v))));
740 }
742 // clang-format on
743};
744
745#endif
746
747#endif
748
749} // namespace SIMD
750
751} // namespace EBGeometry
752
753#endif
EBGeometry_GPU.hpp
GPU portability layer: host/device decoration macros and backend detection.
EBGEOMETRY_HOST_DEVICE
#define EBGEOMETRY_HOST_DEVICE
Mark a function as callable from both host and device code (no-op on host compilers).
Definition EBGeometry_GPU.hpp:80
EBGeometry::SIMD::WidestPack
constexpr size_t WidestPack() noexcept
The widest SIMD register, in lanes of T, that this compilation pass has: 16 floats or 8 doubles under...
Definition EBGeometry_SIMD.hpp:125
EBGeometry::SIMD::NativeWidth
constexpr size_t NativeWidth() noexcept
The lane count to evaluate a group of W values of T with: the widest register this compilation pass h...
Definition EBGeometry_SIMD.hpp:152
EBGeometry::SIMD::both
bool both(const bool a_lhs, const bool a_rhs) noexcept
Whether both conditions hold (scalar form).
Definition EBGeometry_SIMD.hpp:62
EBGeometry::SIMD::select
T select(const bool a_mask, const T a_ifTrue, const T a_ifFalse) noexcept
Pick a_ifTrue where a_mask holds, else a_ifFalse (scalar form).
Definition EBGeometry_SIMD.hpp:49
EBGeometry::SIMD::sqrt
float sqrt(const float a_x) noexcept
Square root (scalar form).
Definition EBGeometry_SIMD.hpp:74
EBGeometry::SIMD::abs
float abs(const float a_x) noexcept
Absolute value (scalar form).
Definition EBGeometry_SIMD.hpp:98
EBGeometry
Namespace containing all of EBGeometry's functionality.
Definition EBGeometry_Affine.hpp:26
EBGeometry::operator*
constexpr Vec2T< T > operator*(const R &s, const Vec2T< T > &a_other) noexcept
Scalar-left multiplication: s * vec.
EBGeometry::operator/
constexpr Vec2T< T > operator/(const R &s, const Vec2T< T > &a_other) noexcept
Component-wise scalar-left division: s / vec.
EBGeometry::SIMD::Pack
N lanes of T in one SIMD register. Defined only for the (T, N) pairs the compiler's flags provide (se...
Definition EBGeometry_SIMD.hpp:165
EBGeometry::SIMD::Scalar< Pack< T, N > >::type
T type
The scalar type.
Definition EBGeometry_SIMD.hpp:191
EBGeometry::SIMD::Scalar
The scalar type of a scalar or pack type: T for T, and T for Pack<T, N>.
Definition EBGeometry_SIMD.hpp:173
EBGeometry::SIMD::Scalar::type
V type
The scalar type.
Definition EBGeometry_SIMD.hpp:177
Generated by doxygen 1.9.8