Added nbody benchmark: behaves as expected

Author: leobeltra

benchmark.h

@@ -380,4 +380,70 @@ void BM_CPUHardRW(benchmark::State &state, T t)
     state.counters["n_elem"] = n;
 }
 
+inline float rand_float() {
+    return static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
+}
+
+template <typename T>
+void BM_nbody(benchmark::State &state, T t)
+{
+    auto n = state.range(0);
+    float dt = 0.01f;
+    const float softening = 1e-9f;
+
+    std::vector<float> Fx(n, 0.0f);
+    std::vector<float> Fy(n, 0.0f);
+    std::vector<float> Fz(n, 0.0f);
+
+    // Inizializza le posizioni e velocità
+    for (int i = 0; i < n; ++i) {
+        MEMBER_ACCESS(t, x, i) = rand_float();
+        MEMBER_ACCESS(t, y, i) = rand_float();
+        MEMBER_ACCESS(t, z, i) = rand_float();
+        MEMBER_ACCESS(t, vx, i) = rand_float();
+        MEMBER_ACCESS(t, vy, i) = rand_float();
+        MEMBER_ACCESS(t, vz, i) = rand_float();
+    }
+
+    for (auto _ : state) {
+        // Calcolo delle forze
+        for (int i = 0; i < n; ++i) {
+            Fx[i] = 0.0f;
+            Fy[i] = 0.0f;
+            Fz[i] = 0.0f;
+
+            for (int j = 0; j < n; ++j) {
+                if (i != j) {
+                    float dx = MEMBER_ACCESS(t, x, j) - MEMBER_ACCESS(t, x, i);
+                    float dy = MEMBER_ACCESS(t, y, j) - MEMBER_ACCESS(t, y, i);
+                    float dz = MEMBER_ACCESS(t, z, j) - MEMBER_ACCESS(t, z, i);
+                    float distSqr = dx * dx + dy * dy + dz * dz + softening;
+                    float invDist = 1.0f / std::sqrt(distSqr);
+                    float invDist3 = invDist * invDist * invDist;
+
+                    Fx[i] += dx * invDist3;
+                    Fy[i] += dy * invDist3;
+                    Fz[i] += dz * invDist3;
+                }
+            }
+
+            MEMBER_ACCESS(t, vx, i) += dt * Fx[i];
+            MEMBER_ACCESS(t, vy, i) += dt * Fy[i];
+            MEMBER_ACCESS(t, vz, i) += dt * Fz[i];
+        }
+
+        // Integrazione posizioni
+        for (int i = 0; i < n; ++i) {
+            MEMBER_ACCESS(t, x, i) += MEMBER_ACCESS(t, vx, i) * dt;
+            MEMBER_ACCESS(t, y, i) += MEMBER_ACCESS(t, vy, i) * dt;
+            MEMBER_ACCESS(t, z, i) += MEMBER_ACCESS(t, vz, i) * dt;
+        }
+    }
+
+    state.counters["n_elem"] = n;
+    state.counters["N^2_interactions"] = benchmark::Counter(
+        static_cast<double>(n) * static_cast<double>(n),
+        benchmark::Counter::kIsRate);
+}
+
 #endif // BENCHMARK_H

soa_boost.cpp

@@ -97,6 +97,17 @@ GENERATE_SOA_LAYOUT(BigSoALayout,
 using BigSoA = BigSoALayout<>;
 using BigSoAView = BigSoA::View;
 
+GENERATE_SOA_LAYOUT(SoANbodyLayout,
+    SOA_COLUMN(float, x),
+    SOA_COLUMN(float, y),
+    SOA_COLUMN(float, z),
+    SOA_COLUMN(float, vx),
+    SOA_COLUMN(float, vy),
+    SOA_COLUMN(float, vz))
+
+using SoANbody = SoANbodyLayout<>;
+using SoANbodyView = SoANbody::View;
+
 int main(int argc, char** argv) {
     std::vector<void *> free_list;
 
@@ -133,6 +144,14 @@ int main(int argc, char** argv) {
         free_list.push_back(buffer);
     }
 
+    for (auto n : N) {
+        auto buffer = reinterpret_cast<std::byte *>(aligned_alloc(SoANbody::alignment, SoANbody::computeDataSize(n)));
+        SoANbody nbodySoA(buffer, n);
+        SoANbodyView nbodySoAView{nbodySoA};
+        benchmark::RegisterBenchmark("BM_nbody", BM_nbody<SoANbodyView>, nbodySoAView)->Arg(n)->Unit(benchmark::kMillisecond);
+        free_list.push_back(buffer);
+    }
+
     benchmark::Initialize(&argc, argv);
     benchmark::RunSpecifiedBenchmarks();
     benchmark::Shutdown();

soa_manual.cpp

@@ -226,6 +226,31 @@ struct S64 {
     }
 };
 
+struct Snbody {
+    float* __restrict__ x, *__restrict__ y, *__restrict__ z;
+    float* __restrict__ vx, *__restrict__ vy, *__restrict__ vz;
+
+    Snbody(std::byte* buf, size_t n) {
+        size_t offset = 0;
+
+        x = reinterpret_cast<float* __restrict__>(buf);
+        offset += align_size(n * sizeof(float));
+        y = reinterpret_cast<float* __restrict__>(buf + offset);
+        offset += align_size(n * sizeof(float));
+        z = reinterpret_cast<float* __restrict__>(buf + offset);
+        offset += align_size(n * sizeof(float));
+        vx = reinterpret_cast<float* __restrict__>(buf + offset);
+        offset += align_size(n * sizeof(float));
+        vy = reinterpret_cast<float* __restrict__>(buf + offset);
+        offset += align_size(n * sizeof(float));
+        vz = reinterpret_cast<float* __restrict__>(buf + offset);
+    }
+
+    static size_t size_bytes(size_t n) {
+        return align_size(sizeof(float[n])) * 6;
+    }
+};
+
 int main(int argc, char** argv) {
     benchmark::Initialize(&argc, argv);
     std::vector<void *> free_list;
@@ -258,6 +283,13 @@ int main(int argc, char** argv) {
         free_list.push_back(buffer);
     }
 
+    for (auto n : N) {
+        auto buffer = reinterpret_cast<std::byte * __restrict__>(std::aligned_alloc(Alignment, Snbody::size_bytes(n)));
+        Snbody t64(buffer, n);
+        benchmark::RegisterBenchmark("BM_nbody", BM_nbody<Snbody>, t64)->Arg(n)->Unit(benchmark::kMillisecond);
+        free_list.push_back(buffer);
+    }
+
     benchmark::RunSpecifiedBenchmarks();
     benchmark::Shutdown();
 

soa_wrapper.cpp

@@ -41,6 +41,13 @@ struct S64 {
     F<Eigen::Matrix3d> x51, x52, x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63;
 };
 
+template <template <class> class F>
+struct Snbody {
+    template<template <class> class F_new>
+    operator Snbody<F_new>() { return {x, y, z, vx, vy, vz}; }
+    F<float> x, y, z, vx, vy, vz;
+};
+
 int main(int argc, char** argv) {
     constexpr wrapper::layout L = wrapper::layout::soa;
 
@@ -86,6 +93,16 @@ int main(int argc, char** argv) {
         benchmark::RegisterBenchmark("BM_CPUHardRW", BM_CPUHardRW<wrapper_type>, t_span)->Arg(n)->Unit(benchmark::kMillisecond);
     }
 
+    for (std::size_t n : N) {
+        // n * 6 * sizeof(float);
+        std::size_t bytes = n * factory::get_size_in_bytes<Snbody, L>();
+        buffer_pointers.emplace_back(new std::byte[bytes]);
+        auto tnbody = factory::buffer_wrapper<Snbody, L>(buffer_pointers.back(), bytes);
+        using wrapper_type = wrapper::wrapper<Snbody, std::span, L>;
+        wrapper_type t_span(tnbody);
+        benchmark::RegisterBenchmark("BM_nbody", BM_nbody<wrapper_type>, t_span)->Arg(n)->Unit(benchmark::kMillisecond);
+    }
+
     benchmark::Initialize(&argc, argv);
     benchmark::RunSpecifiedBenchmarks();
     benchmark::Shutdown();

wrapper/helper.h

@@ -35,6 +35,7 @@ constexpr std::size_t CountMembers() {
     if constexpr (detail::is_aggregate_constructible_from_n<Argument, 2>::value) return 2;
     else if constexpr (detail::is_aggregate_constructible_from_n<Argument,  10>::value) return  10;
     else if constexpr (detail::is_aggregate_constructible_from_n<Argument,  64>::value) return  64;
+    else if constexpr (detail::is_aggregate_constructible_from_n<Argument, 6>::value) return 6;
     else return 100;  // Silence warnings about missing return value
 }
 
@@ -65,7 +66,11 @@ template <
                 m40, m41, m42, m43, m44, m45, m46, m47, m48, m49,
                 m50, m51, m52, m53, m54, m55, m56, m57, m58, m59,
                 m60, m61, m62, m63);
-    } else return void();  // Silence warnings about missing return value
+    } else if constexpr (M == 6) {
+        auto& [m00, m01, m02, m03, m04, m05] = arg;
+        return f(m00, m01, m02, m03, m04, m05);
+    }
+    else return void();  // Silence warnings about missing return value
 }
 
 template <