From e0c07a4cc1996041511317f3aea5e7909d161c4d Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Fri, 31 Mar 2023 13:54:16 +0200 Subject: [PATCH] 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 | 12 + benchmarks/array/array_transform3.cpp | 51 ++++ .../benchmark/array/data_generator.hpp | 36 ++- .../benchmark/common/benchmark_transform3.hpp | 86 ++++++ .../benchmark/common/benchmark_vector.hpp | 6 +- benchmarks/eigen/eigen_transform3.cpp | 51 ++++ .../benchmark/eigen/data_generator.hpp | 32 ++- .../benchmark/vc_soa/data_generator.hpp | 39 ++- benchmarks/vc_soa/vc_soa_transform3.cpp | 60 ++++ .../include/algebra/eigen_cmath.hpp | 2 + .../include/algebra/eigen_eigen.hpp | 2 + .../include/algebra/fastor_fastor.hpp | 2 + frontend/vc_soa/include/algebra/vc_soa.hpp | 7 +- math/vc_soa/CMakeLists.txt | 5 +- .../algebra/math/impl/vc_soa_transform3.hpp | 271 ++++++++++++++++++ .../algebra/math/impl/vc_soa_vector.hpp | 2 +- math/vc_soa/include/algebra/math/vc_soa.hpp | 5 +- storage/vc_soa/CMakeLists.txt | 5 +- .../algebra/storage/impl/vc_soa_matrix44.hpp | 140 +++++++++ .../vc_soa/include/algebra/storage/vc_soa.hpp | 4 +- tests/common/test_host_basics.hpp | 2 +- tests/vc_soa/vc_soa.cpp | 97 ++++++- 22 files changed, 888 insertions(+), 29 deletions(-) create mode 100644 benchmarks/array/array_transform3.cpp create mode 100644 benchmarks/common/include/benchmark/common/benchmark_transform3.hpp create mode 100644 benchmarks/eigen/eigen_transform3.cpp create mode 100644 benchmarks/vc_soa/vc_soa_transform3.cpp create mode 100644 math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp create mode 100644 storage/vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index 90cf663a..950d8b71 100644 --- a/benchmarks/CMakeLists.txt +++ b/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() diff --git a/benchmarks/array/array_transform3.cpp b/benchmarks/array/array_transform3.cpp new file mode 100644 index 00000000..752b3893 --- /dev/null +++ b/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 + +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(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> v_trf_s{cfg}; + transform3_bm> 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(); +} \ No newline at end of file diff --git a/benchmarks/array/include/benchmark/array/data_generator.hpp b/benchmarks/array/include/benchmark/array/data_generator.hpp index 33708714..e30fc080 100644 --- a/benchmarks/array/include/benchmark/array/data_generator.hpp +++ b/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 #include @@ -16,7 +19,7 @@ namespace algebra { /// Fill an @c std::array based vector with random values template -inline void fill_random(std::vector &collection) { +inline void fill_random_vec(std::vector &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 &collection) { std::generate(collection.begin(), collection.end(), rand_obj); } +/// Fill a @c Vc::Vector based transform3 with random values +template +inline void fill_random_trf(std::vector &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 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 \ No newline at end of file diff --git a/benchmarks/common/include/benchmark/common/benchmark_transform3.hpp b/benchmarks/common/include/benchmark/common/benchmark_transform3.hpp new file mode 100644 index 00000000..fb7298e5 --- /dev/null +++ b/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 +#include +#include +#include +#include + +namespace algebra { + +template +void fill_random_trf(std::vector &); + +/// Benchmark for vector operations +template +struct transform3_bm : public vector_bm { + private: + using base_type = vector_bm; + + public: + /// Prefix for the benchmark name + inline static const std::string bm_name{"transform3"}; + + std::vector 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 diff --git a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp index 2ddd0557..55ec6aee 100644 --- a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp +++ b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp @@ -20,7 +20,7 @@ namespace algebra { template -void fill_random(std::vector &); +void fill_random_vec(std::vector &); /// Benchmark for vector operations template @@ -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 diff --git a/benchmarks/eigen/eigen_transform3.cpp b/benchmarks/eigen/eigen_transform3.cpp new file mode 100644 index 00000000..747a526b --- /dev/null +++ b/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 + +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(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> v_trf_s{cfg}; + transform3_bm> 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(); +} \ No newline at end of file diff --git a/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp index 61a9e504..5b1b3a37 100644 --- a/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp +++ b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp @@ -7,15 +7,18 @@ #pragma once +// Project include(s) +#include "algebra/eigen_eigen.hpp" + // System include(s) #include #include namespace algebra { -/// Fill a @c Eigen3 based vector with random values +/// Fill an @c Eigen3 based vector with random values template -inline void fill_random(std::vector &collection) { +inline void fill_random_vec(std::vector &collection) { auto rand_obj = []() { return vector_t::Random(); }; @@ -23,4 +26,29 @@ inline void fill_random(std::vector &collection) { std::generate(collection.begin(), collection.end(), rand_obj); } +/// Fill a @c Eigen3 based transform3 with random values +template +inline void fill_random_trf(std::vector &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 \ No newline at end of file diff --git a/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp b/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp index 98bc1e85..e0fdd4fe 100644 --- a/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp +++ b/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp @@ -7,6 +7,9 @@ #pragma once +// Project include(s) +#include "algebra/vc_soa.hpp" + // System include(s) #include #include @@ -17,7 +20,7 @@ namespace algebra { /// Fill a @c Vc::SimdArray based vector with random values /*template -inline void fill_random( +inline void fill_random_vec( std::vector, allocator_t>> &collection) { @@ -30,7 +33,7 @@ inline void fill_random( /// Fill a @c Vc::Vector based vector with random values template -inline void fill_random(std::vector &collection) { +inline void fill_random_vec(std::vector &collection) { // Generate a vector of the right type with random values auto rand_obj = []() { using simd_vector_t = typename vector_soa_t::value_type; @@ -45,4 +48,36 @@ inline void fill_random(std::vector &collection) { std::generate(collection.begin(), collection.end(), rand_obj); } +/// Fill a @c Vc::Vector based transform3 with random values +template +inline void fill_random_trf(std::vector &collection) { + // Generate a random, but valid affine transformation + auto rand_obj = []() { + using simd_vector_t = typename transform3_t::value_type; + typename transform3_t::vector3 x_axis, z_axis, t; + x_axis[0] = simd_vector_t::Random(); + x_axis[1] = simd_vector_t::Random(); + x_axis[2] = simd_vector_t::Random(); + x_axis = vector::normalize(x_axis); + + z_axis[0] = simd_vector_t::Random(); + z_axis[1] = simd_vector_t::Random(); + z_axis[2] = simd_vector_t::Random(); + + t[0] = simd_vector_t::Random(); + t[1] = simd_vector_t::Random(); + t[2] = simd_vector_t::Random(); + t = vector::normalize(t); + + // Gram-Schmidt projection + simd_vector_t 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 \ No newline at end of file diff --git a/benchmarks/vc_soa/vc_soa_transform3.cpp b/benchmarks/vc_soa/vc_soa_transform3.cpp new file mode 100644 index 00000000..79191685 --- /dev/null +++ b/benchmarks/vc_soa/vc_soa_transform3.cpp @@ -0,0 +1,60 @@ +/** 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/vc_soa.hpp" +#include "benchmark/common/benchmark_transform3.hpp" +#include "benchmark/vc_soa/data_generator.hpp" + +// Benchmark include +#include + +using namespace algebra; + +/// Run vector benchmarks +int main(int argc, char** argv) { + + constexpr std::size_t n_samples{160000}; + + // + // Prepare benchmarks + // + algebra::benchmark_base::configuration cfg_s{}; + // Reduce the number of samples, since a single SoA struct contains multiple + // vectors + cfg_s.n_samples(n_samples / Vc::float_v::Size) + .n_warmup(static_cast(0.1 * cfg_s.n_samples())); + cfg_s.do_sleep(false); + + // For double precision we need more samples (less vectors per SoA) + algebra::benchmark_base::configuration cfg_d{cfg_s}; + cfg_d.n_samples(n_samples / Vc::double_v::Size) + .n_warmup(static_cast(0.1 * cfg_d.n_samples())); + + transform3_bm> v_trf_s{cfg_s}; + transform3_bm> v_trf_d{cfg_d}; + + std::cout << "Algebra-Plugins 'transform3' benchmark (Vc SoA)\n" + << "-----------------------------------------------\n\n" + << "(single)\n" + << cfg_s << "(double)\n" + << cfg_d; + + // + // 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(); +} \ No newline at end of file diff --git a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp b/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp index 8cc9bf74..0fba8fc3 100644 --- a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp +++ b/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp @@ -5,6 +5,8 @@ * Mozilla Public License Version 2.0 */ +#pragma once + // Project include(s). #include "algebra/math/cmath.hpp" #include "algebra/math/eigen.hpp" diff --git a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp index bcf2e38f..4e502dfa 100644 --- a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp +++ b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp @@ -5,6 +5,8 @@ * Mozilla Public License Version 2.0 */ +#pragma once + // Project include(s). #include "algebra/math/cmath.hpp" #include "algebra/math/eigen.hpp" diff --git a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp index 0bb6005e..463c84f3 100644 --- a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp +++ b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp @@ -5,6 +5,8 @@ * Mozilla Public License Version 2.0 */ +#pragma once + // Project include(s). #include "algebra/math/cmath.hpp" #include "algebra/math/fastor.hpp" diff --git a/frontend/vc_soa/include/algebra/vc_soa.hpp b/frontend/vc_soa/include/algebra/vc_soa.hpp index 08b8052d..e89eb8aa 100644 --- a/frontend/vc_soa/include/algebra/vc_soa.hpp +++ b/frontend/vc_soa/include/algebra/vc_soa.hpp @@ -28,10 +28,11 @@ using algebra::storage::operator+; namespace algebra { namespace vc_soa { -/// @name Vc based transforms on @c algebra::vc::storage_type +/// @name Vc based transforms on @c algebra::vc_soa types /// @{ -//@todo +template +using transform3 = math::transform3; /// @} @@ -77,4 +78,4 @@ using matrix_type = vc_soa::matrix_type; } // namespace matrix -} // namespace algebra \ No newline at end of file +} // namespace algebra diff --git a/math/vc_soa/CMakeLists.txt b/math/vc_soa/CMakeLists.txt index 8494dcca..0a23dc3a 100644 --- a/math/vc_soa/CMakeLists.txt +++ b/math/vc_soa/CMakeLists.txt @@ -8,8 +8,9 @@ algebra_add_library( algebra_vc_soa_math vc_soa_math "include/algebra/math/vc_soa.hpp" "include/algebra/math/impl/vc_soa_getter.hpp" - "include/algebra/math/impl/vc_soa_vector.hpp" ) + "include/algebra/math/impl/vc_soa_vector.hpp" + "include/algebra/math/impl/vc_soa_transform3.hpp" ) target_link_libraries( algebra_vc_soa_math - INTERFACE algebra::common algebra::common_math algebra::common_storage Vc::Vc ) + INTERFACE algebra::common algebra::common_math algebra::common_storage algebra::vc_soa_storage Vc::Vc ) algebra_test_public_headers( algebra_vc_soa_math "algebra/math/vc_soa.hpp" ) diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp new file mode 100644 index 00000000..70c2cb70 --- /dev/null +++ b/math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp @@ -0,0 +1,271 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2023-2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "algebra/math/common.hpp" +#include "algebra/qualifiers.hpp" +#include "algebra/storage/impl/vc_soa_matrix44.hpp" +#include "algebra/storage/vector.hpp" + +// Vc include(s). +#ifdef _MSC_VER +#pragma warning(push, 0) +#endif // MSVC +#include +#ifdef _MSC_VER +#pragma warning(pop) +#endif // MSVC + +// System include(s). +#include +#include +#include + +namespace algebra::vc_soa::math { + +using algebra::storage::operator*; +using algebra::storage::operator/; +using algebra::storage::operator-; +using algebra::storage::operator+; + +/// Transform wrapper class to ensure standard API within differnt plugins +template +struct transform3 { + + /// @name Type definitions for the struct + /// @{ + + /// Scalar type used by the transform + using scalar_type = scalar_t; + /// The type of the matrix elements (in this case: Vc::Vector) + using value_type = Vc::Vector; + + /// 3-element "vector" type (does not observe translations) + using vector3 = storage::vector<3, value_type, std::array>; + /// Point in 3D space (does observe translations) + using point3 = vector3; + /// Point in 2D space + using point2 = storage::vector<2, value_type, std::array>; + + /// 4x4 matrix type + using matrix44 = algebra::vc_soa::matrix44; + + /// Function (object) used for accessing a matrix element + using element_getter = algebra::vc_soa::element_getter; + + /// @} + + /// @name Data objects + /// @{ + + matrix44 _data; + matrix44 _data_inv; + + /// @} + + /// Default constructor: identity + ALGEBRA_HOST_DEVICE + transform3() = default; + + /// Contructor with arguments: t, x, y, z + /// + /// @param t the translation (or origin of the new frame) + /// @param x the x axis of the new frame + /// @param y the y axis of the new frame + /// @param z the z axis of the new frame, normal vector for planes + ALGEBRA_HOST_DEVICE + transform3(const vector3 &t, const vector3 &x, const vector3 &y, + const vector3 &z, [[maybe_unused]] bool get_inverse = true) + : _data{x, y, z, t}, _data_inv{invert(_data)} {} + + /// Contructor with arguments: t, z, x + /// + /// @param t the translation (or origin of the new frame) + /// @param z the z axis of the new frame, normal vector for planes + /// @param x the x axis of the new frame + /// + /// @note y will be constructed by cross product + ALGEBRA_HOST_DEVICE + transform3(const vector3 &t, const vector3 &z, const vector3 &x, + bool get_inverse = true) + : transform3(t, x, cross(z, x), z, get_inverse) {} + + /// Constructor with arguments: translation + /// + /// @param t is the transform + ALGEBRA_HOST_DEVICE + transform3(const vector3 &t) : _data{t}, _data_inv{invert(_data)} {} + + /// Constructor with arguments: matrix + /// + /// @param m is the full 4x4 matrix with simd-vector elements + ALGEBRA_HOST_DEVICE + transform3(const matrix44 &m) : _data{m}, _data_inv{invert(_data)} {} + + /// Defaults + transform3(const transform3 &rhs) = default; + ~transform3() = default; + + /// Equality operator + ALGEBRA_HOST_DEVICE + inline constexpr bool operator==(const transform3 &rhs) const { + return (_data == rhs._data); + } + + /// Matrix access operator + ALGEBRA_HOST_DEVICE + inline const value_type &operator()(std::size_t i, std::size_t j) const { + return element_getter{}(_data, i, j); + } + ALGEBRA_HOST_DEVICE + inline value_type &operator()(std::size_t i, std::size_t j) { + return element_getter{}(_data, i, j); + } + + /// The determinant of a 4x4 matrix + /// + /// @param m is the matrix + /// + /// @return a sacalar determinant - no checking done + ALGEBRA_HOST_DEVICE + inline constexpr value_type determinant(const matrix44 &m) const { + return -m.z[0] * m.y[1] * m.x[2] + m.y[0] * m.z[1] * m.x[2] + + m.z[0] * m.x[1] * m.y[2] - m.x[0] * m.z[1] * m.y[2] - + m.y[0] * m.x[1] * m.z[2] + m.x[0] * m.y[1] * m.z[2]; + } + + /// The inverse of a 4x4 matrix + /// + /// @param m is the matrix + /// + /// @return an inverse matrix + ALGEBRA_HOST_DEVICE + inline constexpr matrix44 invert(const matrix44 &m) const { + matrix44 i; + i.x[0] = -m.z[1] * m.y[2] + m.y[1] * m.z[2]; + i.x[1] = m.z[1] * m.x[2] - m.x[1] * m.z[2]; + i.x[2] = -m.y[1] * m.x[2] + m.x[1] * m.y[2]; + // i.x[3] = 0; + i.y[0] = m.z[0] * m.y[2] - m.y[0] * m.z[2]; + i.y[1] = -m.z[0] * m.x[2] + m.x[0] * m.z[2]; + i.y[2] = m.y[0] * m.x[2] - m.x[0] * m.y[2]; + // i.y[3] = 0; + i.z[0] = -m.z[0] * m.y[1] + m.y[0] * m.z[1]; + i.z[1] = m.z[0] * m.x[1] - m.x[0] * m.z[1]; + i.z[2] = -m.y[0] * m.x[1] + m.x[0] * m.y[1]; + // i.z[3] = 0; + i.t[0] = m.t[0] * m.z[1] * m.y[2] - m.z[0] * m.t[1] * m.y[2] - + m.t[0] * m.y[1] * m.z[2] + m.y[0] * m.t[1] * m.z[2] + + m.z[0] * m.y[1] * m.t[2] - m.y[0] * m.z[1] * m.t[2]; + i.t[1] = m.z[0] * m.t[1] * m.x[2] - m.t[0] * m.z[1] * m.x[2] + + m.t[0] * m.x[1] * m.z[2] - m.x[0] * m.t[1] * m.z[2] - + m.z[0] * m.x[1] * m.t[2] + m.x[0] * m.z[1] * m.t[2]; + i.t[2] = m.t[0] * m.y[1] * m.x[2] - m.y[0] * m.t[1] * m.x[2] - + m.t[0] * m.x[1] * m.y[2] + m.x[0] * m.t[1] * m.y[2] + + m.y[0] * m.x[1] * m.t[2] - m.x[0] * m.y[1] * m.t[2]; + // i.t[3] = 1; + const value_type idet{value_type::One() / determinant(i)}; + + i.x = i.x * idet; + i.y = i.y * idet; + i.z = i.z * idet; + i.t = i.t * idet; + + return i; + } + + /// Rotate a vector into / from a frame + /// + /// @param m is the rotation matrix + /// @param v is the vector to be rotated + ALGEBRA_HOST_DEVICE + inline constexpr auto rotate(const matrix44 &m, const vector3 &v) const { + + return m.x * v[0] + m.y * v[1] + m.z * v[2]; + } + + /// @returns the translation of the transform + ALGEBRA_HOST_DEVICE + inline point3 translation() const { return _data.t; } + + /// @returns the 4x4 matrix of the transform + ALGEBRA_HOST_DEVICE + inline const matrix44 &matrix() const { return _data; } + + /// @returns the 4x4 matrix of the inverse transform + ALGEBRA_HOST_DEVICE + inline const matrix44 &matrix_inverse() const { return _data_inv; } + + /// This method transforms a point from the local 3D cartesian frame + /// to the global 3D cartesian frame + /// + /// @tparam point3_t 3D point + /// + /// @param v is the point to be transformed + /// + /// @return a global point + template < + typename point3_t, + std::enable_if_t, bool> = true> + ALGEBRA_HOST_DEVICE inline point3 point_to_global(const point3_t &v) const { + + return rotate(_data, v) + _data.t; + } + + /// This method transforms a point from the global 3D cartesian frame + /// into the local 3D cartesian frame + /// + /// @tparam point3_t 3D point + /// + /// @param v is the point to be transformed + /// + /// @return a local point + template < + typename point3_t, + std::enable_if_t, bool> = true> + ALGEBRA_HOST_DEVICE inline point3 point_to_local(const point3_t &v) const { + + return rotate(_data_inv, v) + _data_inv.t; + } + + /// This method transforms a vector from the local 3D cartesian frame + /// to the global 3D cartesian frame + /// + /// @tparam vector3_t 3D vector + /// + /// @param v is the vector to be transformed + /// + /// @return a vector in global coordinates + template < + typename vector3_t, + std::enable_if_t, bool> = true> + ALGEBRA_HOST_DEVICE inline vector3 vector_to_global( + const vector3_t &v) const { + + return rotate(_data, v); + } + + /// This method transforms a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame + /// + /// @tparam vector3_t 3D vector + /// + /// @param v is the vector to be transformed + /// + /// @return a vector in local coordinates + template < + typename vector3_t, + std::enable_if_t, bool> = true> + ALGEBRA_HOST_DEVICE inline vector3 vector_to_local(const vector3_t &v) const { + + return rotate(_data_inv, v); + } +}; // struct transform3 + +} // namespace algebra::vc_soa::math diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp index b47a4a98..8413f495 100644 --- a/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp +++ b/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp @@ -87,4 +87,4 @@ normalize(const storage::vector, array_t> &v) { // return Vc::rsqrt(dot(v, v)) * v; } -} // namespace algebra::vc_soa::math \ No newline at end of file +} // namespace algebra::vc_soa::math diff --git a/math/vc_soa/include/algebra/math/vc_soa.hpp b/math/vc_soa/include/algebra/math/vc_soa.hpp index 9c71fd49..fc87b25b 100644 --- a/math/vc_soa/include/algebra/math/vc_soa.hpp +++ b/math/vc_soa/include/algebra/math/vc_soa.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (c) 2020-2022 CERN for the benefit of the ACTS project + * (c) 2023-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -9,4 +9,5 @@ // Project include(s). #include "algebra/math/impl/vc_soa_getter.hpp" -#include "algebra/math/impl/vc_soa_vector.hpp" \ No newline at end of file +#include "algebra/math/impl/vc_soa_transform3.hpp" +#include "algebra/math/impl/vc_soa_vector.hpp" diff --git a/storage/vc_soa/CMakeLists.txt b/storage/vc_soa/CMakeLists.txt index 8eb8b8f8..cb1e761e 100644 --- a/storage/vc_soa/CMakeLists.txt +++ b/storage/vc_soa/CMakeLists.txt @@ -1,12 +1,13 @@ # Algebra plugins library, part of the ACTS project (R&D line) # -# (c) 2021-2022 CERN for the benefit of the ACTS project +# (c) 2024 CERN for the benefit of the ACTS project # # Mozilla Public License Version 2.0 # Set up the library. algebra_add_library( algebra_vc_soa_storage vc_soa_storage - "include/algebra/storage/vc_soa.hpp" ) + "include/algebra/storage/vc_soa.hpp" + "include/algebra/storage/impl/vc_soa_matrix44.hpp" ) target_link_libraries( algebra_vc_soa_storage INTERFACE algebra::common algebra::common_storage Vc::Vc ) algebra_test_public_headers( algebra_vc_soa_storage diff --git a/storage/vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp b/storage/vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp new file mode 100644 index 00000000..86061eb0 --- /dev/null +++ b/storage/vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp @@ -0,0 +1,140 @@ +/** Algebra plugins, part of the ACTS project + * + * (c) 2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ +#pragma once + +// Vc include(s). +#ifdef _MSC_VER +#pragma warning(push, 0) +#endif // MSVC +#include +#ifdef _MSC_VER +#pragma warning(pop) +#endif // MSVC + +// System include(s). +#include +#include +#include + +namespace algebra::vc_soa { + +/// 4x4 matrix type used by @c algebra::vc_soa::math::transform3 that has simd +/// vectors as matrix elements +template +struct matrix44 { + + using vector_type = storage::vector<3, Vc::Vector, std::array>; + using value_type = Vc::Vector; + using scalar_type = scalar_t; + + /// Default constructor: Identity with no translation + matrix44() + : x{1.f, 0.f, 0.f}, + y{0.f, 1.f, 0.f}, + z{0.f, 0.f, 1.f}, + t{0.f, 0.f, 0.f} {} + + /// Identity rotation with translation @param translation + matrix44(const vector_type &v) + : x{1.f, 0.f, 0.f}, y{0.f, 1.f, 0.f}, z{0.f, 0.f, 1.f}, t{v} {} + + /// Construct from given column vectors @param x, @param y, @param z, @param t + matrix44(const vector_type &v_0, const vector_type &v_1, + const vector_type &v_2, const vector_type &v_3) + : x{v_0}, y{v_1}, z{v_2}, t{v_3} {} + + /// Identity rotation with translation from single elemenst @param t_0, + /// @param t_1, @param t_2 + matrix44(const value_type &t_0, const value_type &t_1, const value_type &t_2) + : matrix44{{t_0, t_1, t_2}} {} + + /// Construct from elements (simd vector) of matrix column vectors: + /// - column 0: @param x_0, @param x_1, @param x_2 + /// - column 1: @param y_0, @param y_1, @param y_2 + /// - column 2: @param z_0, @param z_1, @param z_2 + /// - column 3: @param t_0, @param t_1, @param t_2 + matrix44(const value_type &x_0, const value_type &x_1, const value_type &x_2, + const value_type &y_0, const value_type &y_1, const value_type &y_2, + const value_type &z_0, const value_type &z_1, const value_type &z_2, + const value_type &t_0, const value_type &t_1, const value_type &t_2) + : x{{x_0, x_1, x_2}}, + y{{y_0, y_1, y_2}}, + z{{z_0, z_1, z_2}}, + t{{t_0, t_1, t_2}} {} + + /// Equality operator between two matrices + bool operator==(const matrix44 &rhs) const { + return ((x == rhs.x) && (y == rhs.y) && (z == rhs.z) && (t == rhs.t)); + } + + /// Data variables + vector_type x, y, z, t; + +}; // struct matrix44 + +/// Functor used to access elements of matrix44 +template +struct element_getter { + + /// Get const access to a matrix element + ALGEBRA_HOST inline const auto &operator()(const matrix44 &m, + std::size_t row, + std::size_t col) const { + + // Make sure that the indices are valid. + assert(row < 4u); + assert(col < 4u); + + // Return the selected element. + switch (col) { + case 0u: + return m.x[row]; + case 1u: + return m.y[row]; + case 2u: + return m.z[row]; + case 3u: + return m.t[row]; + default: +#ifndef _MSC_VER + __builtin_unreachable(); +#else + return std::numeric_limits::max(); +#endif + } + } + + /// Get const access to a matrix element + ALGEBRA_HOST inline auto &operator()(matrix44 &m, std::size_t row, + std::size_t col) const { + + // Make sure that the indices are valid. + assert(row < 4u); + assert(col < 4u); + + // Return the selected element. + switch (col) { + case 0u: + return m.x[row]; + case 1u: + return m.y[row]; + case 2u: + return m.z[row]; + case 3u: + return m.t[row]; + default: +#ifndef _MSC_VER + __builtin_unreachable(); +#else + return std::numeric_limits::max(); +#endif + } + } + +}; // struct element_getter + +} // namespace algebra::vc_soa diff --git a/storage/vc_soa/include/algebra/storage/vc_soa.hpp b/storage/vc_soa/include/algebra/storage/vc_soa.hpp index ffc44fc3..0e5664ad 100644 --- a/storage/vc_soa/include/algebra/storage/vc_soa.hpp +++ b/storage/vc_soa/include/algebra/storage/vc_soa.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (c) 2023 CERN for the benefit of the ACTS project + * (c) 2023-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -59,4 +59,4 @@ using vector6 = vector_type; template using vector8 = vector_type; -} // namespace algebra::vc_soa \ No newline at end of file +} // namespace algebra::vc_soa diff --git a/tests/common/test_host_basics.hpp b/tests/common/test_host_basics.hpp index f0ca6252..6f093d89 100644 --- a/tests/common/test_host_basics.hpp +++ b/tests/common/test_host_basics.hpp @@ -509,7 +509,7 @@ TYPED_TEST_P(test_host_basics, getter) { // This defines the transform3 test suite TYPED_TEST_P(test_host_basics, transform3) { - // Preparatioon work + // Preparation work typename TypeParam::vector3 z = algebra::vector::normalize(typename TypeParam::vector3{3., 2., 1.}); typename TypeParam::vector3 x = diff --git a/tests/vc_soa/vc_soa.cpp b/tests/vc_soa/vc_soa.cpp index d81c4211..33187f51 100644 --- a/tests/vc_soa/vc_soa.cpp +++ b/tests/vc_soa/vc_soa.cpp @@ -12,13 +12,12 @@ #include #include -#include using namespace algebra; using scalar_t = float; -constexpr float tol{1e-6f}; +constexpr float tol{1e-5f}; /// This test the vector functions on an SoA (Vc::Vector) based vector TEST(test_vc_host, vc_soa_vector) { @@ -27,8 +26,8 @@ TEST(test_vc_host, vc_soa_vector) { // Value type is Vc::Vector using value_t = typename vector3_v::value_type; - vector3_v a{value_t::One(), 2.f * value_t::One(), 3.f * value_t::One()}; - vector3_v b{4.f * value_t::One(), 5.f * value_t::One(), 6.f * value_t::One()}; + vector3_v a{1.f, 2.f, 3.f}; + vector3_v b{4.f, 5.f, 6.f}; EXPECT_TRUE((a[0] == value_t(1.f)).isFull()); EXPECT_TRUE((a[1] == value_t(2.f)).isFull()); @@ -108,10 +107,8 @@ TEST(test_vc_host, vc_soa_vector) { TEST(test_vc_host, vc_soa_getter) { using vector3_v = vc_soa::vector3; - // Value type is Vc::Vector - using value_t = typename vector3_v::value_type; - vector3_v a{value_t::One(), 2.f * value_t::One(), 3.f * value_t::One()}; + vector3_v a{1.f, 2.f, 3.f}; // All results in the vector are the same, so only check the first one @@ -136,4 +133,88 @@ TEST(test_vc_host, vc_soa_getter) { auto v_eta = getter::eta(a); EXPECT_NEAR(v_eta[0], static_cast(std::atanh(1. / std::sqrt(14.) * 3.)), tol); -} \ No newline at end of file +} + +/// This test an SoA (Vc::Vector) based affine transform3 +TEST(test_vc_host, vc_soa_transform3) { + + using vector3 = vc_soa::vector3; + using point3 = vc_soa::point3; + // Value type is Vc::Vector + using value_t = typename vector3::value_type; + using transform3 = vc_soa::transform3; + using transform3 = vc_soa::transform3; + + transform3 idty{}; + + EXPECT_TRUE((idty(0, 0) == value_t::One()).isFull()); + EXPECT_TRUE((idty(1, 0) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(2, 0) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(0, 1) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(1, 1) == value_t::One()).isFull()); + EXPECT_TRUE((idty(2, 1) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(0, 2) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(1, 2) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(2, 2) == value_t::One()).isFull()); + EXPECT_TRUE((idty(0, 3) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(1, 3) == value_t::Zero()).isFull()); + EXPECT_TRUE((idty(2, 3) == value_t::Zero()).isFull()); + + // Preparatioon work + vector3 z = vector::normalize(vector3{3.f, 2.f, 1.f}); + vector3 x = vector::normalize(vector3{2.f, -3.f, 0.f}); + vector3 y = vector::cross(z, x); + point3 t = {2.f, 3.f, 4.f}; + + // Test constructor from t, z, x + transform3 trf1(t, z, x); + ASSERT_TRUE(trf1 == trf1); + transform3 trf2; + trf2 = trf1; + + EXPECT_TRUE((trf2(0, 0) == x[0]).isFull()); + EXPECT_TRUE((trf2(1, 0) == x[1]).isFull()); + EXPECT_TRUE((trf2(2, 0) == x[2]).isFull()); + EXPECT_TRUE((trf2(0, 1) == y[0]).isFull()); + EXPECT_TRUE((trf2(1, 1) == y[1]).isFull()); + EXPECT_TRUE((trf2(2, 1) == y[2]).isFull()); + EXPECT_TRUE((trf2(0, 2) == z[0]).isFull()); + EXPECT_TRUE((trf2(1, 2) == z[1]).isFull()); + EXPECT_TRUE((trf2(2, 2) == z[2]).isFull()); + EXPECT_TRUE((trf2(0, 3) == 2.f * value_t::One()).isFull()); + EXPECT_TRUE((trf2(1, 3) == 3.f * value_t::One()).isFull()); + EXPECT_TRUE((trf2(2, 3) == 4.f * value_t::One()).isFull()); + + // Check that local origin translates into global translation + point3 lzero = {0.f, 0.f, 0.f}; + point3 gzero = trf2.point_to_global(lzero); + EXPECT_TRUE((gzero[0] == t[0]).isFull()); + EXPECT_TRUE((gzero[1] == t[1]).isFull()); + EXPECT_TRUE((gzero[2] == t[2]).isFull()); + + // Check a round trip for point + point3 loc_pt = {3.f, 4.f, 5.f}; + point3 glob_pt = trf2.point_to_global(loc_pt); + point3 loc_pt_r = trf2.point_to_local(glob_pt); + EXPECT_NEAR(loc_pt[0][0], loc_pt_r[0][0], tol); + EXPECT_NEAR(loc_pt[1][0], loc_pt_r[1][0], tol); + EXPECT_NEAR(loc_pt[2][0], loc_pt_r[2][0], tol); + + // Check a point versus vector transform + // vector should not change if transformed by a pure translation + transform3 ttrf(t); + + vector3 glob_vec = {1.f, 1.f, 1.f}; + vector3 loc_vec = ttrf.vector_to_local(glob_vec); + EXPECT_NEAR(glob_vec[0][0], loc_vec[0][0], tol); + EXPECT_NEAR(glob_vec[1][0], loc_vec[1][0], tol); + EXPECT_NEAR(glob_vec[2][0], loc_vec[2][0], tol); + + // Check a round trip for vector + vector3 loc_vecB = {7.f, 8.f, 9.f}; + vector3 glob_vecB = trf2.vector_to_local(loc_vecB); + vector3 loc_vecC = trf2.vector_to_global(glob_vecB); + EXPECT_NEAR(loc_vecB[0][0], loc_vecC[0][0], tol); + EXPECT_NEAR(loc_vecB[1][0], loc_vecC[1][0], tol); + EXPECT_NEAR(loc_vecB[2][0], loc_vecC[2][0], tol); +}