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); +}