From 5b0bd48f43534addfed8294dcfc28248dd0cbc3e Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Sun, 26 Mar 2023 15:34:09 +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 | 10 +- 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 | 5 +- math/vc_soa/CMakeLists.txt | 3 +- .../algebra/math/impl/vc_soa_transform3.hpp | 386 ++++++++++++++++++ math/vc_soa/include/algebra/math/vc_soa.hpp | 1 + tests/common/test_host_basics.hpp | 2 +- tests/vc_soa/vc_soa.cpp | 97 ++++- 18 files changed, 855 insertions(+), 22 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 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 32da28ab..e124ff6d 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 @@ -71,7 +71,7 @@ struct vector_unaryOP_bm : public vector_bm> { 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 and this->m_cfg.do_sleep()) { + if (state.thread_index() == 0 && this->m_cfg.do_sleep()) { std::this_thread::sleep_for(std::chrono::seconds(this->m_cfg.n_sleep())); } @@ -113,7 +113,7 @@ struct vector_binaryOP_bm : public vector_bm> { 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 and this->m_cfg.do_sleep()) { + if (state.thread_index() == 0 && this->m_cfg.do_sleep()) { std::this_thread::sleep_for(std::chrono::seconds(this->m_cfg.n_sleep())); } 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..41af8d04 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; /// @} diff --git a/math/vc_soa/CMakeLists.txt b/math/vc_soa/CMakeLists.txt index d3418d8c..b49bffa0 100644 --- a/math/vc_soa/CMakeLists.txt +++ b/math/vc_soa/CMakeLists.txt @@ -8,6 +8,7 @@ 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 ) 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..a045c4f7 --- /dev/null +++ b/math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp @@ -0,0 +1,386 @@ +/** 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 "algebra/math/cmath.hpp" +#include "algebra/qualifiers.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 + +namespace algebra::vc_soa::math { + +using algebra::storage::operator*; +using algebra::storage::operator/; +using algebra::storage::operator-; +using algebra::storage::operator+; + +namespace internal { + +/// 4x4 matrix type used by @c algebra::vc_soa::math::transform3 that has simd +/// vectors as matrix elements +template class array_t> +struct matrix44 { + + using vector_type = storage::vector<3, Vc::Vector, array_t>; + using value_type = Vc::Vector; + using scalar_type = scalar_t; + + /// Default constructor: Identity with no translation + matrix44() + : x{value_type::One(), value_type::Zero(), value_type::Zero()}, + y{value_type::Zero(), value_type::One(), value_type::Zero()}, + z{value_type::Zero(), value_type::Zero(), value_type::One()}, + t{value_type::Zero(), value_type::Zero(), value_type::Zero()} {} + + /// Identity rotation with translation @param translation + matrix44(const vector_type &v) + : x{value_type::One(), value_type::Zero(), value_type::Zero()}, + y{value_type::Zero(), value_type::One(), value_type::Zero()}, + z{value_type::Zero(), value_type::Zero(), value_type::One()}, + 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 class array_t> +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: + return m.x[0]; + } + } + + /// 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: + return m.x[0]; + } + } + +}; // struct element_getter + +} // namespace internal + +/// Transform wrapper class to ensure standard API within differnt plugins +template