feat: add a Vc SoA implementation for transform3 (#97)

This PR is based on the previous development for an SoA Vc-based algebra plugin and adds the transform3 implementation, including a test and benchmarks. Like the current vc_vc plugin, it uses the vector3 type as column vectors in the 4x4 matrix type that is used by the transform3. Elements that are known to be equal to zero or one are optimized away in the inversion and determinant calculations.

benchmarks/CMakeLists.txt

@@ -33,6 +33,10 @@ algebra_add_benchmark( array_vector
    "array/array_vector.cpp"
    LINK_LIBRARIES benchmark::benchmark algebra::bench_common
                    algebra::bench_array algebra::array_cmath )
+algebra_add_benchmark( array_transform3
+   "array/array_transform3.cpp"
+   LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                   algebra::bench_array algebra::array_cmath )
 
 if( ALGEBRA_PLUGINS_INCLUDE_EIGEN )
    add_library( algebra_bench_eigen INTERFACE )
@@ -50,6 +54,10 @@ if( ALGEBRA_PLUGINS_INCLUDE_EIGEN )
       "eigen/eigen_vector.cpp"
       LINK_LIBRARIES benchmark::benchmark algebra::bench_common
                      algebra::bench_eigen algebra::eigen_eigen )
+   algebra_add_benchmark( eigen_transform3
+      "eigen/eigen_transform3.cpp"
+      LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                     algebra::bench_eigen algebra::eigen_eigen )
 endif()
 
 if( ALGEBRA_PLUGINS_INCLUDE_VC )
@@ -70,5 +78,9 @@ if( ALGEBRA_PLUGINS_INCLUDE_VC )
          "vc_soa/vc_soa_vector.cpp"
          LINK_LIBRARIES benchmark::benchmark algebra::bench_common
                         algebra::bench_vc_soa algebra::vc_soa )
+      algebra_add_benchmark( vc_soa_transform3
+         "vc_soa/vc_soa_transform3.cpp"
+         LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                        algebra::bench_vc_soa algebra::vc_soa )
    endif()
 endif()

benchmarks/array/array_transform3.cpp

@@ -0,0 +1,51 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/array_cmath.hpp"
+#include "benchmark/array/data_generator.hpp"
+#include "benchmark/common/benchmark_transform3.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+  constexpr std::size_t n_warmup{static_cast<std::size_t>(0.1 * n_samples)};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg{};
+  cfg.n_samples(n_samples).n_warmup(n_warmup);
+  cfg.do_sleep(false);
+
+  transform3_bm<array::transform3<float>> v_trf_s{cfg};
+  transform3_bm<array::transform3<double>> v_trf_d{cfg};
+
+  std::cout << "Algebra-Plugins 'transform3' benchmark (std::array)\n"
+            << "---------------------------------------------------\n\n"
+            << cfg;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_trf_s.name() + "_single").c_str(), v_trf_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_trf_d.name() + "_double").c_str(), v_trf_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}

benchmarks/array/include/benchmark/array/data_generator.hpp

@@ -7,6 +7,9 @@
 
 #pragma once
 
+// Project include(s)
+#include "algebra/array_cmath.hpp"
+
 // System include(s)
 #include <algorithm>
 #include <random>
@@ -16,7 +19,7 @@ namespace algebra {
 
 /// Fill an @c std::array based vector with random values
 template <typename vector_t>
-inline void fill_random(std::vector<vector_t> &collection) {
+inline void fill_random_vec(std::vector<vector_t> &collection) {
 
   // Generate a vector of the right type with random values
   std::random_device rd;
@@ -29,4 +32,35 @@ inline void fill_random(std::vector<vector_t> &collection) {
   std::generate(collection.begin(), collection.end(), rand_obj);
 }
 
+/// Fill a @c Vc::Vector based transform3 with random values
+template <typename transform3_t>
+inline void fill_random_trf(std::vector<transform3_t> &collection) {
+
+  using vector_t = typename transform3_t::vector3;
+
+  // Generate a random, but valid affine transformation
+  std::random_device rd;
+  std::mt19937 mt(rd());
+  std::uniform_real_distribution<typename transform3_t::scalar_type> dist(0.f,
+                                                                          1.f);
+
+  auto rand_obj = [&]() {
+    vector_t x_axis, z_axis, t;
+
+    x_axis = vector::normalize(vector_t{dist(mt), dist(mt), dist(mt)});
+    z_axis = {dist(mt), dist(mt), dist(mt)};
+    t = vector::normalize(vector_t{dist(mt), dist(mt), dist(mt)});
+
+    // Gram-Schmidt projection
+    typename transform3_t::scalar_type coeff =
+        vector::dot(x_axis, z_axis) / getter::norm(x_axis);
+    z_axis = x_axis - coeff * z_axis;
+
+    return transform3_t{t, x_axis, vector::normalize(z_axis)};
+  };
+
+  collection.resize(collection.capacity());
+  std::generate(collection.begin(), collection.end(), rand_obj);
+}
+
 }  // namespace algebra

benchmarks/common/include/benchmark/common/benchmark_transform3.hpp

@@ -0,0 +1,86 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "benchmark_vector.hpp"
+
+// System include(s)
+#include <chrono>
+#include <iostream>
+#include <string>
+#include <thread>
+#include <vector>
+
+namespace algebra {
+
+template <typename transform3_t>
+void fill_random_trf(std::vector<transform3_t> &);
+
+/// Benchmark for vector operations
+template <typename transform3_t>
+struct transform3_bm : public vector_bm<typename transform3_t::vector3> {
+ private:
+  using base_type = vector_bm<typename transform3_t::vector3>;
+
+ public:
+  /// Prefix for the benchmark name
+  inline static const std::string bm_name{"transform3"};
+
+  std::vector<transform3_t> trfs;
+
+  /// No default construction: Cannot prepare data
+  transform3_bm() = delete;
+  std::string name() const override { return base_type::name + "_" + bm_name; }
+
+  /// Construct from an externally provided configuration @param cfg
+  transform3_bm(benchmark_base::configuration cfg) : base_type{cfg} {
+
+    const std::size_t n_data{this->m_cfg.n_samples() + this->m_cfg.n_warmup()};
+
+    trfs.reserve(n_data);
+
+    fill_random_trf(trfs);
+  }
+
+  /// Clear state
+  virtual ~transform3_bm() { trfs.clear(); }
+
+  /// Benchmark case
+  void operator()(::benchmark::State &state) override {
+
+    const std::size_t n_samples{this->m_cfg.n_samples()};
+    const std::size_t n_warmup{this->m_cfg.n_warmup()};
+
+    // Spin down before benchmark (Thread zero is counting the clock)
+    if (state.thread_index() == 0 && this->m_cfg.do_sleep()) {
+      std::this_thread::sleep_for(std::chrono::seconds(this->m_cfg.n_sleep()));
+    }
+
+    // Run the benchmark
+    for (auto _ : state) {
+      // Warm-up
+      state.PauseTiming();
+      if (this->m_cfg.do_warmup()) {
+        for (std::size_t i{0u}; i < n_warmup; ++i) {
+          ::benchmark::DoNotOptimize(
+              this->trfs[i].vector_to_global(this->a[i]));
+          benchmark::ClobberMemory();
+        }
+      }
+      state.ResumeTiming();
+
+      for (std::size_t i{n_warmup}; i < n_samples + n_warmup; ++i) {
+        ::benchmark::DoNotOptimize(this->trfs[i].vector_to_global(this->a[i]));
+        benchmark::ClobberMemory();
+      }
+    }
+  }
+};
+
+}  // namespace algebra

benchmarks/common/include/benchmark/common/benchmark_vector.hpp

@@ -20,7 +20,7 @@
 namespace algebra {
 
 template <typename vector_t>
-void fill_random(std::vector<vector_t> &);
+void fill_random_vec(std::vector<vector_t> &);
 
 /// Benchmark for vector operations
 template <typename vector_t>
@@ -42,8 +42,8 @@ struct vector_bm : public benchmark_base {
     a.reserve(n_data);
     b.reserve(n_data);
 
-    fill_random(a);
-    fill_random(b);
+    fill_random_vec(a);
+    fill_random_vec(b);
   }
 
   /// Clear state

benchmarks/eigen/eigen_transform3.cpp

@@ -0,0 +1,51 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/eigen_eigen.hpp"
+#include "benchmark/common/benchmark_transform3.hpp"
+#include "benchmark/eigen/data_generator.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+  constexpr std::size_t n_warmup{static_cast<std::size_t>(0.1 * n_samples)};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg{};
+  cfg.n_samples(n_samples).n_warmup(n_warmup);
+  cfg.do_sleep(false);
+
+  transform3_bm<eigen::transform3<float>> v_trf_s{cfg};
+  transform3_bm<eigen::transform3<double>> v_trf_d{cfg};
+
+  std::cout << "Algebra-Plugins 'transform3' benchmark (Eigen3)\n"
+            << "-----------------------------------------------\n\n"
+            << cfg;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_trf_s.name() + "_single").c_str(), v_trf_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_trf_d.name() + "_double").c_str(), v_trf_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}

benchmarks/eigen/include/benchmark/eigen/data_generator.hpp

@@ -7,20 +7,48 @@
 
 #pragma once
 
+// Project include(s)
+#include "algebra/eigen_eigen.hpp"
+
 // System include(s)
 #include <algorithm>
 #include <random>
 
 namespace algebra {
 
-/// Fill a @c Eigen3 based vector with random values
+/// Fill an @c Eigen3 based vector with random values
 template <typename vector_t>
-inline void fill_random(std::vector<vector_t> &collection) {
+inline void fill_random_vec(std::vector<vector_t> &collection) {
 
   auto rand_obj = []() { return vector_t::Random(); };
 
   collection.resize(collection.capacity());
   std::generate(collection.begin(), collection.end(), rand_obj);
 }
 
+/// Fill a @c Eigen3 based transform3 with random values
+template <typename transform3_t>
+inline void fill_random_trf(std::vector<transform3_t> &collection) {
+
+  using vector_t = typename transform3_t::vector3;
+
+  auto rand_obj = []() {
+    vector_t x_axis, z_axis, t;
+
+    x_axis = vector::normalize(vector_t::Random());
+    z_axis = vector_t::Random();
+    t = vector::normalize(vector_t::Random());
+
+    // Gram-Schmidt projection
+    typename transform3_t::scalar_type coeff =
+        vector::dot(x_axis, z_axis) / getter::norm(x_axis);
+    z_axis = x_axis - coeff * z_axis;
+
+    return transform3_t{t, x_axis, vector::normalize(z_axis)};
+  };
+
+  collection.resize(collection.capacity());
+  std::generate(collection.begin(), collection.end(), rand_obj);
+}
+
 }  // namespace algebra