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