diff --git a/README.md b/README.md index 17d15b04..27d4d529 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ for the R&D projects [detray](https://github.com/acts-project/detray) and | Backend | CPU | CUDA | SYCL | | ------------------------------------------------------------------------- | --- | ---- | ---- | -| cmath | ✅ | ✅ | ✅ | +| generic | ✅ | ✅ | ✅ | | [Eigen](https://eigen.tuxfamily.org) | ✅ | ✅ | ✅ | | [SMatrix](https://root.cern.ch/doc/master/group__SMatrixGroup.html) | ✅ | ⚪ | ⚪ | | [VC](https://github.com/VcDevel/Vc) | ✅ | ⚪ | ⚪ | diff --git a/benchmarks/array/array_getter.cpp b/benchmarks/array/array_getter.cpp index 89f7364c..68b1b2d2 100644 --- a/benchmarks/array/array_getter.cpp +++ b/benchmarks/array/array_getter.cpp @@ -27,18 +27,6 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); - using phi_f_t = vector_unaryOP_bm; - using theta_f_t = vector_unaryOP_bm; - using perp_f_t = vector_unaryOP_bm; - using norm_f_t = vector_unaryOP_bm; - using eta_f_t = vector_unaryOP_bm; - - using phi_d_t = vector_unaryOP_bm; - using theta_d_t = vector_unaryOP_bm; - using perp_d_t = vector_unaryOP_bm; - using norm_d_t = vector_unaryOP_bm; - using eta_d_t = vector_unaryOP_bm; - std::cout << "-----------------------------------------------\n" << "Algebra-Plugins 'getter' benchmark (std::array)\n" << "-----------------------------------------------\n\n" @@ -47,7 +35,6 @@ int main(int argc, char** argv) { // // Register all benchmarks // - ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(cfg) ::benchmark::Initialize(&argc, argv); ::benchmark::RunSpecifiedBenchmarks(); diff --git a/benchmarks/array/array_vector.cpp b/benchmarks/array/array_vector.cpp index 5d21f9b0..3b6f388a 100644 --- a/benchmarks/array/array_vector.cpp +++ b/benchmarks/array/array_vector.cpp @@ -27,6 +27,12 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); + using phi_f_t = vector_unaryOP_bm; + using theta_f_t = vector_unaryOP_bm; + using perp_f_t = vector_unaryOP_bm; + using norm_f_t = vector_unaryOP_bm; + using eta_f_t = vector_unaryOP_bm; + using add_f_t = vector_binaryOP_bm; using sub_f_t = vector_binaryOP_bm; using dot_f_t = vector_binaryOP_bm; @@ -34,6 +40,12 @@ int main(int argc, char** argv) { using normlz_f_t = vector_unaryOP_bm; + using phi_d_t = vector_unaryOP_bm; + using theta_d_t = vector_unaryOP_bm; + using perp_d_t = vector_unaryOP_bm; + using norm_d_t = vector_unaryOP_bm; + using eta_d_t = vector_unaryOP_bm; + using add_d_t = vector_binaryOP_bm; using sub_d_t = vector_binaryOP_bm; using dot_d_t = vector_binaryOP_bm; diff --git a/benchmarks/array/include/benchmark/array/data_generator.hpp b/benchmarks/array/include/benchmark/array/data_generator.hpp index cf244088..b4bf5939 100644 --- a/benchmarks/array/include/benchmark/array/data_generator.hpp +++ b/benchmarks/array/include/benchmark/array/data_generator.hpp @@ -53,7 +53,7 @@ inline void fill_random_trf(std::vector &collection) { // Gram-Schmidt projection typename transform3_t::scalar_type coeff = - vector::dot(x_axis, z_axis) / getter::norm(x_axis); + vector::dot(x_axis, z_axis) / vector::norm(x_axis); z_axis = x_axis - coeff * z_axis; return transform3_t{t, x_axis, vector::normalize(z_axis)}; diff --git a/benchmarks/common/include/benchmark/common/benchmark_getter.hpp b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp index bb7bd205..c6c9c9ae 100644 --- a/benchmarks/common/include/benchmark/common/benchmark_getter.hpp +++ b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp @@ -11,36 +11,5 @@ #include "benchmark_vector.hpp" #include "register_benchmark.hpp" -namespace algebra::bench_op { - -// Macro for declaring the predefined materials (with Density effect data) -#define ALGEBRA_PLUGINS_BENCH_GETTER(GETTER_NAME) \ - struct GETTER_NAME { \ - inline static const std::string name{#GETTER_NAME}; \ - template \ - auto operator()(const vector_t &a) const { \ - return algebra::getter::GETTER_NAME(a); \ - } \ - }; - -// Functions to be benchmarked -ALGEBRA_PLUGINS_BENCH_GETTER(phi) -ALGEBRA_PLUGINS_BENCH_GETTER(theta) -ALGEBRA_PLUGINS_BENCH_GETTER(perp) -ALGEBRA_PLUGINS_BENCH_GETTER(norm) -ALGEBRA_PLUGINS_BENCH_GETTER(eta) - -// Macro for registering all getter benchmarks -#define ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(CFG) \ - algebra::register_benchmark(CFG, "_single"); \ - algebra::register_benchmark(CFG, "_double"); \ - algebra::register_benchmark(CFG, "_single"); \ - algebra::register_benchmark(CFG, "_double"); \ - algebra::register_benchmark(CFG, "_single"); \ - algebra::register_benchmark(CFG, "_double"); \ - algebra::register_benchmark(CFG, "_single"); \ - algebra::register_benchmark(CFG, "_double"); \ - algebra::register_benchmark(CFG, "_single"); \ - algebra::register_benchmark(CFG, "_double"); - +namespace algebra::bench_op { /* TODO: Implement */ } // namespace algebra::bench_op diff --git a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp index 93de4c5c..9d672055 100644 --- a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp +++ b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp @@ -152,13 +152,23 @@ struct cross { return algebra::vector::cross(a, b); } }; -struct normalize { - inline static const std::string name{"normalize"}; - template - auto operator()(const vector_t &a) const { - return algebra::vector::normalize(a); - } -}; + +// Macro for declaring vector unary ops +#define ALGEBRA_PLUGINS_BENCH_VECTOR(OP) \ + struct OP { \ + inline static const std::string name{#OP}; \ + template \ + auto operator()(const vector_t &a) const { \ + return algebra::vector::OP(a); \ + } \ + }; + +ALGEBRA_PLUGINS_BENCH_VECTOR(phi) +ALGEBRA_PLUGINS_BENCH_VECTOR(theta) +ALGEBRA_PLUGINS_BENCH_VECTOR(eta) +ALGEBRA_PLUGINS_BENCH_VECTOR(perp) +ALGEBRA_PLUGINS_BENCH_VECTOR(norm) +ALGEBRA_PLUGINS_BENCH_VECTOR(normalize) } // namespace bench_op @@ -173,6 +183,17 @@ struct normalize { algebra::register_benchmark(CFG, "_single"); \ algebra::register_benchmark(CFG, "_double"); \ algebra::register_benchmark(CFG, "_single"); \ - algebra::register_benchmark(CFG, "_double"); + algebra::register_benchmark(CFG, "_double"); \ + \ + algebra::register_benchmark(CFG, "_single"); \ + algebra::register_benchmark(CFG, "_double"); \ + algebra::register_benchmark(CFG, "_single"); \ + algebra::register_benchmark(CFG, "_double"); \ + algebra::register_benchmark(CFG, "_single"); \ + algebra::register_benchmark(CFG, "_double"); \ + algebra::register_benchmark(CFG, "_single"); \ + algebra::register_benchmark(CFG, "_double"); \ + algebra::register_benchmark(CFG, "_single"); \ + algebra::register_benchmark(CFG, "_double"); } // namespace algebra diff --git a/benchmarks/eigen/eigen_getter.cpp b/benchmarks/eigen/eigen_getter.cpp index aa0fd427..4806945d 100644 --- a/benchmarks/eigen/eigen_getter.cpp +++ b/benchmarks/eigen/eigen_getter.cpp @@ -27,18 +27,6 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); - using phi_f_t = vector_unaryOP_bm; - using theta_f_t = vector_unaryOP_bm; - using perp_f_t = vector_unaryOP_bm; - using norm_f_t = vector_unaryOP_bm; - using eta_f_t = vector_unaryOP_bm; - - using phi_d_t = vector_unaryOP_bm; - using theta_d_t = vector_unaryOP_bm; - using perp_d_t = vector_unaryOP_bm; - using norm_d_t = vector_unaryOP_bm; - using eta_d_t = vector_unaryOP_bm; - std::cout << "------------------------------------------\n" << "Algebra-Plugins 'getter' benchmark (Eigen)\n" << "------------------------------------------\n\n" @@ -47,8 +35,6 @@ int main(int argc, char** argv) { // // Register all benchmarks // - ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(cfg) - ::benchmark::Initialize(&argc, argv); ::benchmark::RunSpecifiedBenchmarks(); ::benchmark::Shutdown(); diff --git a/benchmarks/eigen/eigen_vector.cpp b/benchmarks/eigen/eigen_vector.cpp index de791f86..1f162fbe 100644 --- a/benchmarks/eigen/eigen_vector.cpp +++ b/benchmarks/eigen/eigen_vector.cpp @@ -27,6 +27,12 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); + using phi_f_t = vector_unaryOP_bm; + using theta_f_t = vector_unaryOP_bm; + using perp_f_t = vector_unaryOP_bm; + using norm_f_t = vector_unaryOP_bm; + using eta_f_t = vector_unaryOP_bm; + using add_f_t = vector_binaryOP_bm; using sub_f_t = vector_binaryOP_bm; using dot_f_t = vector_binaryOP_bm; @@ -34,6 +40,12 @@ int main(int argc, char** argv) { using normlz_f_t = vector_unaryOP_bm; + using phi_d_t = vector_unaryOP_bm; + using theta_d_t = vector_unaryOP_bm; + using perp_d_t = vector_unaryOP_bm; + using norm_d_t = vector_unaryOP_bm; + using eta_d_t = vector_unaryOP_bm; + using add_d_t = vector_binaryOP_bm; using sub_d_t = vector_binaryOP_bm; using dot_d_t = vector_binaryOP_bm; diff --git a/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp index 21560406..e6963b5c 100644 --- a/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp +++ b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp @@ -41,7 +41,7 @@ inline void fill_random_trf(std::vector &collection) { // Gram-Schmidt projection typename transform3_t::scalar_type coeff = - vector::dot(x_axis, z_axis) / getter::norm(x_axis); + vector::dot(x_axis, z_axis) / vector::norm(x_axis); z_axis = x_axis - coeff * z_axis; return transform3_t{t, x_axis, vector::normalize(z_axis)}; diff --git a/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp b/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp index 2d28a50c..51b3ade3 100644 --- a/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp +++ b/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp @@ -48,7 +48,7 @@ inline void fill_random_trf(std::vector &collection) { // Gram-Schmidt projection typename transform3_t::scalar_type coeff = - vector::dot(x_axis, z_axis) / getter::norm(x_axis); + vector::dot(x_axis, z_axis) / vector::norm(x_axis); z_axis = x_axis - coeff * z_axis; return transform3_t{t, x_axis, vector::normalize(z_axis)}; diff --git a/benchmarks/vc_aos/vc_aos_getter.cpp b/benchmarks/vc_aos/vc_aos_getter.cpp index 8f391020..3b9b0725 100644 --- a/benchmarks/vc_aos/vc_aos_getter.cpp +++ b/benchmarks/vc_aos/vc_aos_getter.cpp @@ -27,18 +27,6 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); - using phi_f_t = vector_unaryOP_bm; - using theta_f_t = vector_unaryOP_bm; - using perp_f_t = vector_unaryOP_bm; - using norm_f_t = vector_unaryOP_bm; - using eta_f_t = vector_unaryOP_bm; - - using phi_d_t = vector_unaryOP_bm; - using theta_d_t = vector_unaryOP_bm; - using perp_d_t = vector_unaryOP_bm; - using norm_d_t = vector_unaryOP_bm; - using eta_d_t = vector_unaryOP_bm; - std::cout << "-------------------------------------------\n" << "Algebra-Plugins 'getter' benchmark (Vc AoS)\n" << "-------------------------------------------\n\n" @@ -47,7 +35,6 @@ int main(int argc, char** argv) { // // Register all benchmarks // - ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(cfg) ::benchmark::Initialize(&argc, argv); ::benchmark::RunSpecifiedBenchmarks(); diff --git a/benchmarks/vc_aos/vc_aos_vector.cpp b/benchmarks/vc_aos/vc_aos_vector.cpp index fbe8261c..6397c1d9 100644 --- a/benchmarks/vc_aos/vc_aos_vector.cpp +++ b/benchmarks/vc_aos/vc_aos_vector.cpp @@ -24,6 +24,12 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); + using phi_f_t = vector_unaryOP_bm; + using theta_f_t = vector_unaryOP_bm; + using perp_f_t = vector_unaryOP_bm; + using norm_f_t = vector_unaryOP_bm; + using eta_f_t = vector_unaryOP_bm; + using add_f_t = vector_binaryOP_bm; using sub_f_t = vector_binaryOP_bm; using dot_f_t = vector_binaryOP_bm; @@ -31,6 +37,12 @@ int main(int argc, char** argv) { using normlz_f_t = vector_unaryOP_bm; + using phi_d_t = vector_unaryOP_bm; + using theta_d_t = vector_unaryOP_bm; + using perp_d_t = vector_unaryOP_bm; + using norm_d_t = vector_unaryOP_bm; + using eta_d_t = vector_unaryOP_bm; + using add_d_t = vector_binaryOP_bm; using sub_d_t = vector_binaryOP_bm; using dot_d_t = vector_binaryOP_bm; @@ -40,7 +52,7 @@ int main(int argc, char** argv) { vector_unaryOP_bm; std::cout << "-------------------------------------------\n" - << "Algebra-Plugins 'vector' benchmark (Vc SoA)\n" + << "Algebra-Plugins 'vector' benchmark (Vc AoS)\n" << "-------------------------------------------\n\n" << cfg; 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 b4d0a591..baaff335 100644 --- a/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp +++ b/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp @@ -55,7 +55,7 @@ inline void fill_random_trf(std::vector &collection) { t = vector::normalize(t); // Gram-Schmidt projection - simd_vector_t coeff = vector::dot(x_axis, z_axis) / getter::norm(x_axis); + simd_vector_t coeff = vector::dot(x_axis, z_axis) / vector::norm(x_axis); z_axis = x_axis - coeff * z_axis; return transform3_t{t, x_axis, vector::normalize(z_axis)}; diff --git a/benchmarks/vc_soa/vc_soa_getter.cpp b/benchmarks/vc_soa/vc_soa_getter.cpp index 28a2ead1..a0c54b94 100644 --- a/benchmarks/vc_soa/vc_soa_getter.cpp +++ b/benchmarks/vc_soa/vc_soa_getter.cpp @@ -36,18 +36,6 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg_d{cfg_s}; cfg_d.n_samples(n_samples / Vc::double_v::Size); - using phi_f_t = vector_unaryOP_bm; - using theta_f_t = vector_unaryOP_bm; - using perp_f_t = vector_unaryOP_bm; - using norm_f_t = vector_unaryOP_bm; - using eta_f_t = vector_unaryOP_bm; - - using phi_d_t = vector_unaryOP_bm; - using theta_d_t = vector_unaryOP_bm; - using perp_d_t = vector_unaryOP_bm; - using norm_d_t = vector_unaryOP_bm; - using eta_d_t = vector_unaryOP_bm; - std::cout << "-------------------------------------------\n" << "Algebra-Plugins 'getter' benchmark (Vc SoA)\n" << "-------------------------------------------\n\n" @@ -58,17 +46,6 @@ int main(int argc, char** argv) { // // Register all benchmarks // - algebra::register_benchmark(cfg_s, "_single"); - algebra::register_benchmark(cfg_d, "_double"); - algebra::register_benchmark(cfg_s, "_single"); - algebra::register_benchmark(cfg_d, "_double"); - algebra::register_benchmark(cfg_s, "_single"); - algebra::register_benchmark(cfg_d, "_double"); - algebra::register_benchmark(cfg_s, "_single"); - algebra::register_benchmark(cfg_d, "_double"); - algebra::register_benchmark(cfg_s, "_single"); - algebra::register_benchmark(cfg_d, "_double"); - ::benchmark::Initialize(&argc, argv); ::benchmark::RunSpecifiedBenchmarks(); ::benchmark::Shutdown(); diff --git a/benchmarks/vc_soa/vc_soa_vector.cpp b/benchmarks/vc_soa/vc_soa_vector.cpp index abdd6413..51b96798 100644 --- a/benchmarks/vc_soa/vc_soa_vector.cpp +++ b/benchmarks/vc_soa/vc_soa_vector.cpp @@ -33,6 +33,12 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg_d{cfg_s}; cfg_d.n_samples(n_samples / Vc::double_v::Size); + using phi_f_t = vector_unaryOP_bm; + using theta_f_t = vector_unaryOP_bm; + using perp_f_t = vector_unaryOP_bm; + using norm_f_t = vector_unaryOP_bm; + using eta_f_t = vector_unaryOP_bm; + using add_f_t = vector_binaryOP_bm; using sub_f_t = vector_binaryOP_bm; using dot_f_t = vector_binaryOP_bm; @@ -40,6 +46,12 @@ int main(int argc, char** argv) { using normlz_f_t = vector_unaryOP_bm; + using phi_d_t = vector_unaryOP_bm; + using theta_d_t = vector_unaryOP_bm; + using perp_d_t = vector_unaryOP_bm; + using norm_d_t = vector_unaryOP_bm; + using eta_d_t = vector_unaryOP_bm; + using add_d_t = vector_binaryOP_bm; using sub_d_t = vector_binaryOP_bm; using dot_d_t = vector_binaryOP_bm; @@ -58,6 +70,17 @@ int main(int argc, char** argv) { // // Register all benchmarks // + algebra::register_benchmark(cfg_s, "_single"); + algebra::register_benchmark(cfg_d, "_double"); + algebra::register_benchmark(cfg_s, "_single"); + algebra::register_benchmark(cfg_d, "_double"); + algebra::register_benchmark(cfg_s, "_single"); + algebra::register_benchmark(cfg_d, "_double"); + algebra::register_benchmark(cfg_s, "_single"); + algebra::register_benchmark(cfg_d, "_double"); + algebra::register_benchmark(cfg_s, "_single"); + algebra::register_benchmark(cfg_d, "_double"); + algebra::register_benchmark(cfg_s, "_single"); algebra::register_benchmark(cfg_d, "_double"); algebra::register_benchmark(cfg_s, "_single"); diff --git a/cmake/algebra-plugins-functions.cmake b/cmake/algebra-plugins-functions.cmake index 21af0427..c19dff0c 100644 --- a/cmake/algebra-plugins-functions.cmake +++ b/cmake/algebra-plugins-functions.cmake @@ -120,6 +120,10 @@ function( algebra_add_benchmark name ) # Create the test executable. set( bench_exe_name "algebra_benchmark_${name}" ) add_executable( ${bench_exe_name} ${ARG_UNPARSED_ARGUMENTS} ) + + target_compile_options( algebra_benchmark_${name} PRIVATE + "-march=native" "-ftree-vectorize") + if( ARG_LINK_LIBRARIES ) target_link_libraries( ${bench_exe_name} PRIVATE ${ARG_LINK_LIBRARIES} ) endif() diff --git a/common/include/algebra/type_traits.hpp b/common/include/algebra/type_traits.hpp index 2ae0519a..5eac0c3e 100644 --- a/common/include/algebra/type_traits.hpp +++ b/common/include/algebra/type_traits.hpp @@ -53,6 +53,15 @@ struct vector {}; template using vector_t = typename vector::type; + +template N, typename T> +using get_vector_t = typename vector::template other_type; + +template +struct matrix {}; + +template ROWS, index_t COLS, typename T> +using get_matrix_t = typename matrix::template other_type; /// @} /// Matrix dimensions @@ -66,10 +75,14 @@ requires std::is_fundamental_v struct dimensions { using size_type = std::size_t; + static constexpr size_type dim{0}; static constexpr size_type rows{1}; static constexpr size_type columns{1}; }; +template +inline constexpr index_t dim{dimensions::sim}; + template inline constexpr index_t rows{dimensions::rows}; @@ -82,10 +95,126 @@ inline constexpr index_t rank{std::min(rows, columns)}; template inline constexpr index_t size{rows * columns}; +template +inline constexpr bool is_matrix{dimensions::dim == 2}; + +template +inline constexpr bool is_vector{dimensions::dim == 1}; + template inline constexpr bool is_square{(rows == columns)}; /// @} +/// Getter types +/// @{ +template +struct element_getter {}; + +template +using element_getter_t = typename element_getter::type; + +template +struct block_getter {}; + +template +using block_getter_t = typename block_getter::type; /// @} } // namespace algebra::trait + +/// Default type trait specializations +/// @{ +#define ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(A) \ + \ + namespace trait { \ + \ + template \ + struct index> { \ + using type = algebra::A::size_type; \ + }; \ + \ + template \ + struct index> { \ + using type = algebra::A::size_type; \ + }; \ + \ + template \ + struct dimensions> { \ + \ + using size_type = index_t>; \ + \ + static constexpr size_type dim{1}; \ + static constexpr size_type rows{N}; \ + static constexpr size_type columns{1}; \ + }; \ + \ + template \ + struct dimensions> { \ + \ + using size_type = index_t>; \ + \ + static constexpr size_type dim{2}; \ + static constexpr size_type rows{ROWS}; \ + static constexpr size_type columns{COLS}; \ + }; \ + \ + template \ + struct value> { \ + using type = T; \ + }; \ + \ + template \ + struct value> { \ + using type = T; \ + }; \ + \ + template \ + struct vector> { \ + \ + template \ + using other_type = A::vector_type; \ + \ + using type = other_type; \ + }; \ + \ + template \ + struct vector> { \ + \ + template \ + using other_type = A::vector_type; \ + \ + using type = other_type; \ + }; \ + \ + template \ + struct matrix> { \ + template \ + using other_type = A::matrix_type; \ + \ + using type = A::matrix_type; \ + }; \ + \ + template \ + struct matrix> { \ + template \ + using other_type = A::matrix_type; \ + \ + using type = other_type; \ + }; \ + \ + template \ + struct element_getter> { \ + using type = A::element_getter; \ + }; \ + \ + template \ + struct element_getter> { \ + using type = A::element_getter; \ + }; \ + \ + template \ + struct block_getter> { \ + using type = A::block_getter; \ + }; \ + \ + } // namespace algebra::trait diff --git a/frontend/CMakeLists.txt b/frontend/CMakeLists.txt index 268ffb92..31c03e5b 100644 --- a/frontend/CMakeLists.txt +++ b/frontend/CMakeLists.txt @@ -1,6 +1,6 @@ # Algebra plugins library, part of the ACTS project (R&D line) # -# (c) 2021-2022 CERN for the benefit of the ACTS project +# (c) 2021-2024 CERN for the benefit of the ACTS project # # Mozilla Public License Version 2.0 @@ -8,18 +8,18 @@ add_subdirectory( array_cmath ) if( ALGEBRA_PLUGINS_INCLUDE_EIGEN ) - add_subdirectory( eigen_cmath ) + add_subdirectory( eigen_generic ) add_subdirectory( eigen_eigen ) endif() if( ALGEBRA_PLUGINS_INCLUDE_SMATRIX ) - add_subdirectory( smatrix_cmath ) + add_subdirectory( smatrix_generic ) add_subdirectory( smatrix_smatrix ) endif() if( ALGEBRA_PLUGINS_INCLUDE_VC ) add_subdirectory( vc_aos ) - add_subdirectory( vc_cmath ) + add_subdirectory( vc_aos_generic ) if( NOT "${CMAKE_CXX_COMPILER_ID}" MATCHES "AppleClang" ) add_subdirectory( vc_soa ) endif() diff --git a/frontend/array_cmath/CMakeLists.txt b/frontend/array_cmath/CMakeLists.txt index 30df0a94..2ee81318 100644 --- a/frontend/array_cmath/CMakeLists.txt +++ b/frontend/array_cmath/CMakeLists.txt @@ -8,6 +8,6 @@ algebra_add_library( algebra_array_cmath array_cmath "include/algebra/array_cmath.hpp" ) target_link_libraries( algebra_array_cmath - INTERFACE algebra::common algebra::array_storage algebra::cmath_math ) + INTERFACE algebra::common algebra::array_storage algebra::cmath_math algebra::generic_math ) algebra_test_public_headers( algebra_array_cmath "algebra/array_cmath.hpp" ) diff --git a/frontend/array_cmath/include/algebra/array_cmath.hpp b/frontend/array_cmath/include/algebra/array_cmath.hpp index 559f200d..177787c3 100644 --- a/frontend/array_cmath/include/algebra/array_cmath.hpp +++ b/frontend/array_cmath/include/algebra/array_cmath.hpp @@ -9,6 +9,7 @@ // Project include(s). #include "algebra/math/cmath.hpp" +#include "algebra/math/generic.hpp" #include "algebra/storage/array.hpp" /// @name Operators on @c algebra::array::storage_type @@ -27,29 +28,10 @@ namespace getter { /// @name Getter functions on @c algebra::array::storage_type /// @{ -using cmath::eta; -using cmath::norm; -using cmath::perp; -using cmath::phi; -using cmath::theta; - -/// @} - -/// Function extracting a slice from a matrix -template -ALGEBRA_HOST_DEVICE inline array::storage_type vector( - const array::matrix_type& m, std::size_t row, - std::size_t col) { - - return cmath::vector_getter()(m, row, col); -} - -/// @name Getter functions on @c algebra::array::matrix_type -/// @{ - -using cmath::element; +using cmath::storage::block; +using cmath::storage::element; +using cmath::storage::set_block; +using cmath::storage::vector; /// @} @@ -60,41 +42,49 @@ namespace vector { /// @name Vector functions on @c algebra::array::storage_type /// @{ -using cmath::cross; +// array specific implementations using cmath::dot; using cmath::normalize; +// generic implementations +using generic::math::cross; +using generic::math::eta; +using generic::math::norm; +using generic::math::perp; +using generic::math::phi; +using generic::math::theta; + /// @} } // namespace vector namespace matrix { -using cmath::block; -using cmath::determinant; +/// @name Matrix functions on @c algebra::array::storage_type +/// @{ + using cmath::identity; -using cmath::inverse; -using cmath::set_block; using cmath::set_identity; using cmath::set_zero; -using cmath::transpose; using cmath::zero; +using generic::math::determinant; +using generic::math::inverse; +using generic::math::transpose; + +/// @} + } // namespace matrix namespace array { -template -using element_getter = - cmath::element_getter; - /// @name cmath based transforms on @c algebra::array /// @{ template -using transform3 = cmath::transform3, - cmath::block_getter>; +using transform3 = + generic::math::transform3; /// @} diff --git a/frontend/eigen_cmath/CMakeLists.txt b/frontend/eigen_cmath/CMakeLists.txt deleted file mode 100644 index d0dc3774..00000000 --- a/frontend/eigen_cmath/CMakeLists.txt +++ /dev/null @@ -1,14 +0,0 @@ -# Algebra plugins library, part of the ACTS project (R&D line) -# -# (c) 2021-2023 CERN for the benefit of the ACTS project -# -# Mozilla Public License Version 2.0 - -# Set up the library. -algebra_add_library( algebra_eigen_cmath eigen_cmath - "include/algebra/eigen_cmath.hpp" ) -target_link_libraries( algebra_eigen_cmath - INTERFACE algebra::common algebra::eigen_storage algebra::cmath_math - algebra::eigen_math Eigen3::Eigen ) -algebra_test_public_headers( algebra_eigen_cmath - "algebra/eigen_cmath.hpp" ) diff --git a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp b/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp deleted file mode 100644 index ff6696ff..00000000 --- a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp +++ /dev/null @@ -1,103 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2020-2022 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/math/eigen.hpp" -#include "algebra/storage/eigen.hpp" - -// Eigen include(s). -#ifdef _MSC_VER -#pragma warning(push, 0) -#endif // MSVC -#include -#ifdef _MSC_VER -#pragma warning(pop) -#endif // MSVC - -// System include(s). -#include - -namespace algebra { - -namespace getter { - -/// @name Getter functions on @c algebra::eigen::storage_type -/// @{ - -using eigen::math::eta; -using eigen::math::norm; -using eigen::math::perp; -using eigen::math::phi; -using eigen::math::theta; - -/// @} - -/// Function extracting a slice from the matrix used by -/// @c algebra::eigen::transform3 -template -ALGEBRA_HOST_DEVICE inline auto vector(const Eigen::MatrixBase& m, - std::size_t row, std::size_t col) { - - return m.template block(static_cast(row), - static_cast(col)); -} - -/// @name Getter functions on @c algebra::eigen::matrix_type -/// @{ - -using eigen::math::element; - -/// @} - -} // namespace getter - -namespace vector { - -/// @name Vector functions on @c algebra::eigen::storage_type -/// @{ - -using eigen::math::cross; -using eigen::math::dot; -using eigen::math::normalize; - -/// @} - -} // namespace vector - -namespace matrix { - -using eigen::math::block; -using eigen::math::determinant; -using eigen::math::identity; -using eigen::math::inverse; -using eigen::math::set_block; -using eigen::math::set_identity; -using eigen::math::set_zero; -using eigen::math::transpose; -using eigen::math::zero; - -} // namespace matrix - -namespace eigen { - -/// @name cmath based transforms on @c algebra::eigen -/// @{ - -template -using transform3 = - cmath::transform3; - -/// @} - -} // namespace eigen - -} // namespace algebra diff --git a/frontend/eigen_eigen/CMakeLists.txt b/frontend/eigen_eigen/CMakeLists.txt index 21b73502..463f30ca 100644 --- a/frontend/eigen_eigen/CMakeLists.txt +++ b/frontend/eigen_eigen/CMakeLists.txt @@ -8,7 +8,7 @@ algebra_add_library( algebra_eigen_eigen eigen_eigen "include/algebra/eigen_eigen.hpp" ) target_link_libraries( algebra_eigen_eigen - INTERFACE algebra::common algebra::eigen_storage algebra::cmath_math + INTERFACE algebra::common algebra::eigen_storage algebra::generic_math algebra::eigen_math Eigen3::Eigen ) algebra_test_public_headers( algebra_eigen_eigen "algebra/eigen_eigen.hpp" ) diff --git a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp index efe599fa..235a2900 100644 --- a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp +++ b/frontend/eigen_eigen/include/algebra/eigen_eigen.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) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -8,8 +8,8 @@ #pragma once // Project include(s). -#include "algebra/math/cmath.hpp" #include "algebra/math/eigen.hpp" +#include "algebra/math/generic.hpp" #include "algebra/storage/eigen.hpp" // Eigen include(s). @@ -25,31 +25,13 @@ namespace algebra { namespace getter { -/// @name Getter functions on @c algebra::eigen::storage_type +/// @name Getter functions on @c algebra::eigen /// @{ -using eigen::math::eta; -using eigen::math::norm; -using eigen::math::perp; -using eigen::math::phi; -using eigen::math::theta; - -/// @} - -/// Function extracting a slice from the matrix used by -/// @c algebra::eigen::transform3 -template -ALGEBRA_HOST_DEVICE inline auto vector(const Eigen::MatrixBase& m, - std::size_t row, std::size_t col) { - - return m.template block(static_cast(row), - static_cast(col)); -} - -/// @name Getter functions on @c algebra::eigen::matrix_type -/// @{ - -using eigen::math::element; +using eigen::storage::block; +using eigen::storage::element; +using eigen::storage::set_block; +using eigen::storage::vector; /// @} @@ -62,24 +44,33 @@ namespace vector { using eigen::math::cross; using eigen::math::dot; +using eigen::math::eta; +using eigen::math::norm; using eigen::math::normalize; +using generic::math::perp; +using generic::math::phi; +using generic::math::theta; + /// @} } // namespace vector namespace matrix { -using eigen::math::block; +/// @name Matrix functions on @c algebra::eigen::storage_type +/// @{ + using eigen::math::determinant; using eigen::math::identity; using eigen::math::inverse; -using eigen::math::set_block; using eigen::math::set_identity; using eigen::math::set_zero; using eigen::math::transpose; using eigen::math::zero; +/// @} + } // namespace matrix namespace eigen { diff --git a/frontend/eigen_generic/CMakeLists.txt b/frontend/eigen_generic/CMakeLists.txt new file mode 100644 index 00000000..4b338dd6 --- /dev/null +++ b/frontend/eigen_generic/CMakeLists.txt @@ -0,0 +1,14 @@ +# Algebra plugins library, part of the ACTS project (R&D line) +# +# (c) 2021-2023 CERN for the benefit of the ACTS project +# +# Mozilla Public License Version 2.0 + +# Set up the library. +algebra_add_library( algebra_eigen_generic eigen_generic + "include/algebra/eigen_generic.hpp" ) +target_link_libraries( algebra_eigen_generic + INTERFACE algebra::common algebra::eigen_storage algebra::generic_math + algebra::eigen_math Eigen3::Eigen ) +algebra_test_public_headers( algebra_eigen_generic + "algebra/eigen_generic.hpp" ) diff --git a/frontend/eigen_generic/include/algebra/eigen_generic.hpp b/frontend/eigen_generic/include/algebra/eigen_generic.hpp new file mode 100644 index 00000000..4ab351c2 --- /dev/null +++ b/frontend/eigen_generic/include/algebra/eigen_generic.hpp @@ -0,0 +1,92 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2020-2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "algebra/math/eigen.hpp" +#include "algebra/math/generic.hpp" +#include "algebra/storage/eigen.hpp" + +// Eigen include(s). +#ifdef _MSC_VER +#pragma warning(push, 0) +#endif // MSVC +#include +#ifdef _MSC_VER +#pragma warning(pop) +#endif // MSVC + +// System include(s). +#include + +namespace algebra { + +namespace getter { + +/// @name Getter functions on @c algebra::eigen::matrix_type +/// @{ + +using eigen::storage::block; +using eigen::storage::element; +using eigen::storage::set_block; +using eigen::storage::vector; + +/// @} + +} // namespace getter + +namespace vector { + +/// @name Vector functions on @c algebra::eigen::storage_type +/// @{ + +using generic::math::cross; +using generic::math::dot; +using generic::math::eta; +using generic::math::norm; +using generic::math::normalize; +using generic::math::perp; +using generic::math::phi; +using generic::math::theta; + +/// @} + +} // namespace vector + +namespace matrix { + +/// @name Matrix functions on @c algebra::eigen::storage_type +/// @{ + +using generic::math::determinant; +using generic::math::identity; +using generic::math::inverse; +using generic::math::set_identity; +using generic::math::set_zero; +using generic::math::transpose; +using generic::math::zero; + +/// @} + +} // namespace matrix + +namespace eigen { + +/// @name generic based transforms on @c algebra::eigen +/// @{ + +template +using transform3 = + generic::math::transform3; + +/// @} + +} // namespace eigen + +} // namespace algebra diff --git a/frontend/fastor_fastor/CMakeLists.txt b/frontend/fastor_fastor/CMakeLists.txt index 2b4c3bdc..91d35878 100644 --- a/frontend/fastor_fastor/CMakeLists.txt +++ b/frontend/fastor_fastor/CMakeLists.txt @@ -8,7 +8,7 @@ algebra_add_library( algebra_fastor_fastor fastor_fastor "include/algebra/fastor_fastor.hpp" ) target_link_libraries( algebra_fastor_fastor - INTERFACE algebra::common algebra::fastor_storage algebra::cmath_math + INTERFACE algebra::common algebra::fastor_storage algebra::generic_math algebra::fastor_math ) algebra_test_public_headers( algebra_fastor_fastor "algebra/fastor_fastor.hpp" ) diff --git a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp index 1dcb4bd7..39b2f21f 100644 --- a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp +++ b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (c) 2022 CERN for the benefit of the ACTS project + * (c) 2022-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -8,8 +8,8 @@ #pragma once // Project include(s). -#include "algebra/math/cmath.hpp" #include "algebra/math/fastor.hpp" +#include "algebra/math/generic.hpp" #include "algebra/storage/fastor.hpp" // Fastor include(s). @@ -25,34 +25,13 @@ namespace algebra { namespace getter { -/// @name Getter functions on @c algebra::fastor::storage_type -/// @{ - -using fastor::math::eta; -using fastor::math::norm; -using fastor::math::perp; -using fastor::math::phi; -using fastor::math::theta; - -/// @} - -/// Function extracting a slice from the matrix used by -/// @c algebra::fastor::transform3 -template -ALGEBRA_HOST_DEVICE inline auto vector( - const Fastor::Tensor& m, std::size_t row, - std::size_t col) { - - return fastor::storage_type( - m(Fastor::seq(static_cast(row), static_cast(row + SIZE)), - static_cast(col))); -} - /// @name Getter functions on @c algebra::fastor::matrix_type /// @{ -using fastor::math::element; +using fastor::storage::block; +using fastor::storage::element; +using fastor::storage::set_block; +using fastor::storage::vector; /// @} @@ -65,7 +44,13 @@ namespace vector { using fastor::math::cross; using fastor::math::dot; +using fastor::math::eta; +using fastor::math::norm; using fastor::math::normalize; +using fastor::math::perp; +using fastor::math::theta; + +using generic::math::phi; /// @} @@ -73,23 +58,31 @@ using fastor::math::normalize; namespace matrix { -using fastor::math::block; +/// @name Matrix functions on @c algebra::fastor::storage_type +/// @{ + using fastor::math::determinant; using fastor::math::identity; using fastor::math::inverse; -using fastor::math::set_block; using fastor::math::set_identity; using fastor::math::set_zero; using fastor::math::transpose; using fastor::math::zero; +/// @} + } // namespace matrix namespace fastor { +/// @name Transform on @c algebra::fastor::storage_type +/// @{ + template using transform3 = math::transform3; +/// @} + } // namespace fastor } // namespace algebra diff --git a/frontend/smatrix_cmath/CMakeLists.txt b/frontend/smatrix_cmath/CMakeLists.txt deleted file mode 100644 index f835e19e..00000000 --- a/frontend/smatrix_cmath/CMakeLists.txt +++ /dev/null @@ -1,17 +0,0 @@ -# Algebra plugins library, part of the ACTS project (R&D line) -# -# (c) 2021-2023 CERN for the benefit of the ACTS project -# -# Mozilla Public License Version 2.0 - -# Find Smatrix. -find_package( ROOT COMPONENTS Smatrix REQUIRED) - -# Set up the library. -algebra_add_library( algebra_smatrix_cmath smatrix_cmath - "include/algebra/smatrix_cmath.hpp" ) -target_link_libraries( algebra_smatrix_cmath - INTERFACE algebra::common algebra::smatrix_storage algebra::cmath_math - algebra::smatrix_math ROOT::Smatrix ) -algebra_test_public_headers( algebra_smatrix_cmath - "algebra/smatrix_cmath.hpp" ) diff --git a/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp b/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp deleted file mode 100644 index 77522811..00000000 --- a/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp +++ /dev/null @@ -1,100 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2020-2022 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/math/smatrix.hpp" -#include "algebra/storage/smatrix.hpp" - -// ROOT/Smatrix include(s). -#include - -namespace algebra { -namespace getter { - -/// @name Getter functions on @c algebra::smatrix::storage_type -/// @{ - -using smatrix::math::eta; -using smatrix::math::norm; -using smatrix::math::perp; -using smatrix::math::phi; -using smatrix::math::theta; - -/// @} - -/// Function extracting a slice from the matrix used by -/// @c algebra::smatrix::transform3 -template -ALGEBRA_HOST_DEVICE inline auto vector( - const ROOT::Math::SMatrix& m, unsigned int row, - unsigned int col) { - - return m.template SubCol>(col, row); -} - -/// @name Getter functions on @c algebra::smatrix::matrix_type -/// @{ - -using smatrix::math::element; - -/// @} - -} // namespace getter - -namespace vector { - -/// @name Vector functions on @c algebra::smatrix::storage_type -/// @{ - -using smatrix::math::cross; -using smatrix::math::dot; -using smatrix::math::normalize; - -/// @} - -} // namespace vector - -namespace matrix { - -using smatrix::math::block; -using smatrix::math::determinant; -using smatrix::math::identity; -using smatrix::math::inverse; -using smatrix::math::set_block; -using smatrix::math::set_identity; -using smatrix::math::set_zero; -using smatrix::math::transpose; -using smatrix::math::zero; - -} // namespace matrix - -namespace smatrix { - -template -using element_getter = smatrix::math::element_getter; - -template -using block_getter = smatrix::math::block_getter; - -/// @name cmath based transforms on @c algebra::smatrix -/// @{ - -template -using transform3 = - cmath::transform3, - block_getter>; - -/// @} - -} // namespace smatrix - -} // namespace algebra diff --git a/frontend/smatrix_generic/CMakeLists.txt b/frontend/smatrix_generic/CMakeLists.txt new file mode 100644 index 00000000..aeb10666 --- /dev/null +++ b/frontend/smatrix_generic/CMakeLists.txt @@ -0,0 +1,17 @@ +# Algebra plugins library, part of the ACTS project (R&D line) +# +# (c) 2021-2023 CERN for the benefit of the ACTS project +# +# Mozilla Public License Version 2.0 + +# Find Smatrix. +find_package( ROOT COMPONENTS Smatrix REQUIRED) + +# Set up the library. +algebra_add_library( algebra_smatrix_generic smatrix_generic + "include/algebra/smatrix_generic.hpp" ) +target_link_libraries( algebra_smatrix_generic + INTERFACE algebra::common algebra::smatrix_storage algebra::generic_math + algebra::smatrix_math ROOT::Smatrix ) +algebra_test_public_headers( algebra_smatrix_generic + "algebra/smatrix_generic.hpp" ) diff --git a/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp b/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp new file mode 100644 index 00000000..9e579407 --- /dev/null +++ b/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp @@ -0,0 +1,83 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2020-2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "algebra/math/generic.hpp" +#include "algebra/math/smatrix.hpp" +#include "algebra/storage/smatrix.hpp" + +// ROOT/Smatrix include(s). +#include + +namespace algebra { + +namespace getter { + +/// @name Getter functions on @c algebra::smatrix::matrix_type +/// @{ + +using smatrix::storage::block; +using smatrix::storage::element; +using smatrix::storage::set_block; +using smatrix::storage::vector; + +/// @} + +} // namespace getter + +namespace vector { + +/// @name Vector functions on @c algebra::smatrix::storage_type +/// @{ + +using generic::math::cross; +using generic::math::dot; +using generic::math::eta; +using generic::math::norm; +using generic::math::normalize; +using generic::math::perp; +using generic::math::phi; +using generic::math::theta; + +/// @} + +} // namespace vector + +namespace matrix { + +/// @name Matrix functions on @c algebra::smatrix::storage_type +/// @{ + +using generic::math::determinant; +using generic::math::identity; +using generic::math::inverse; +using generic::math::set_identity; +using generic::math::set_zero; +using generic::math::transpose; +using generic::math::zero; + +/// @} + +} // namespace matrix + +namespace smatrix { + +/// @name generic based transforms on @c algebra::smatrix +/// @{ + +template +using transform3 = + generic::math::transform3; + +/// @} + +} // namespace smatrix + +} // namespace algebra diff --git a/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp b/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp index e9792aae..da6f14c1 100644 --- a/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp +++ b/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp @@ -15,32 +15,13 @@ namespace algebra { namespace getter { -/// @name Getter functions on @c algebra::smatrix::storage_type -/// @{ - -using smatrix::math::eta; -using smatrix::math::norm; -using smatrix::math::perp; -using smatrix::math::phi; -using smatrix::math::theta; - -/// @} - -/// Function extracting a slice from the matrix used by -/// @c algebra::smatrix::transform3 -template -ALGEBRA_HOST_DEVICE inline auto vector( - const ROOT::Math::SMatrix& m, unsigned int row, - unsigned int col) { - - return m.template SubCol>(col, row); -} - /// @name Getter functions on @c algebra::smatrix::matrix_type /// @{ -using smatrix::math::element; +using smatrix::storage::block; +using smatrix::storage::element; +using smatrix::storage::set_block; +using smatrix::storage::vector; /// @} @@ -53,7 +34,12 @@ namespace vector { using smatrix::math::cross; using smatrix::math::dot; +using smatrix::math::eta; +using smatrix::math::norm; using smatrix::math::normalize; +using smatrix::math::perp; +using smatrix::math::phi; +using smatrix::math::theta; /// @} @@ -61,16 +47,19 @@ using smatrix::math::normalize; namespace matrix { -using smatrix::math::block; +/// @name Matrix functions on @c algebra::smatrix::storage_type +/// @{ + using smatrix::math::determinant; using smatrix::math::identity; using smatrix::math::inverse; -using smatrix::math::set_block; using smatrix::math::set_identity; using smatrix::math::set_zero; using smatrix::math::transpose; using smatrix::math::zero; +/// @} + } // namespace matrix namespace smatrix { diff --git a/frontend/vc_aos/CMakeLists.txt b/frontend/vc_aos/CMakeLists.txt index e00cddff..ceffc48a 100644 --- a/frontend/vc_aos/CMakeLists.txt +++ b/frontend/vc_aos/CMakeLists.txt @@ -8,7 +8,7 @@ algebra_add_library( algebra_vc_aos vc_aos "include/algebra/vc_aos.hpp" ) target_link_libraries( algebra_vc_aos - INTERFACE algebra::common algebra::vc_aos_storage algebra::cmath_math + INTERFACE algebra::common algebra::vc_aos_storage algebra::generic_math algebra::vc_aos_math ) algebra_test_public_headers( algebra_vc_aos "algebra/vc_aos.hpp" ) diff --git a/frontend/vc_aos/include/algebra/vc_aos.hpp b/frontend/vc_aos/include/algebra/vc_aos.hpp index d18c569d..ef52aa6b 100644 --- a/frontend/vc_aos/include/algebra/vc_aos.hpp +++ b/frontend/vc_aos/include/algebra/vc_aos.hpp @@ -8,6 +8,7 @@ #pragma once // Project include(s). +#include "algebra/math/generic.hpp" #include "algebra/math/vc_aos.hpp" #include "algebra/storage/vc_aos.hpp" @@ -16,63 +17,16 @@ #include namespace algebra { -namespace vc_aos { - -/// @name Vc based transforms on @c algebra::vc_aos::storage_type -/// @{ - -template -using transform3 = math::transform3; - -/// @} - -} // namespace vc_aos namespace getter { -/// @name Getter functions on @c algebra::vc_aos types -/// @{ - -using vc_aos::math::eta; -using vc_aos::math::norm; -using vc_aos::math::perp; -using vc_aos::math::phi; -using vc_aos::math::theta; - -/// @} - /// @name Getter functions on @c algebra::vc_aos::matrix_type /// @{ -using vc_aos::math::element; - -/// Function extracting a slice from matrix44 - const -template class array_t, - typename value_t, std::size_t ROW, std::size_t COL> -requires(SIZE <= 4) ALGEBRA_HOST_DEVICE inline const - auto& vector(const storage::matrix& m, - [[maybe_unused]] std::size_t row, - [[maybe_unused]] std::size_t col) { - - assert(row == 0); - assert(col < COL); - - return m[col]; -} - -/// Function extracting a slice from matrix44 - non-const -template class array_t, - typename value_t, std::size_t ROW, std::size_t COL> -requires(SIZE <= 4) ALGEBRA_HOST_DEVICE - inline auto& vector(storage::matrix& m, - [[maybe_unused]] std::size_t row, - [[maybe_unused]] std::size_t col) { - - assert(row == 0); - assert(col < COL); - - return m[col]; -} +using vc_aos::storage::block; +using vc_aos::storage::element; +using vc_aos::storage::set_block; +using vc_aos::storage::vector; /// @} @@ -83,26 +37,51 @@ namespace vector { /// @name Vector functions on @c algebra::vc_aos types /// @{ +// Vc array specific using vc_aos::math::cross; using vc_aos::math::dot; +using vc_aos::math::eta; +using vc_aos::math::norm; using vc_aos::math::normalize; +// No specific vectorized implementation needed +using generic::math::perp; +using generic::math::phi; +using generic::math::theta; + /// @} } // namespace vector namespace matrix { -using vc_aos::math::block; -using vc_aos::math::determinant; +/// @name Matrix functions on @c algebra::vc_aos types +/// @{ + using vc_aos::math::identity; -using vc_aos::math::inverse; -using vc_aos::math::set_block; using vc_aos::math::set_identity; using vc_aos::math::set_zero; using vc_aos::math::transpose; using vc_aos::math::zero; +// Placeholder, until vectorization-friendly version is available +using generic::math::determinant; +using generic::math::inverse; + +/// @} + } // namespace matrix +namespace vc_aos { + +/// @name Vc based transforms on @c algebra::vc_aos::storage_type +/// @{ + +template +using transform3 = math::transform3; + +/// @} + +} // namespace vc_aos + } // namespace algebra diff --git a/frontend/vc_aos_generic/CMakeLists.txt b/frontend/vc_aos_generic/CMakeLists.txt new file mode 100644 index 00000000..8001f597 --- /dev/null +++ b/frontend/vc_aos_generic/CMakeLists.txt @@ -0,0 +1,14 @@ +# Algebra plugins library, part of the ACTS project (R&D line) +# +# (c) 2021-2024 CERN for the benefit of the ACTS project +# +# Mozilla Public License Version 2.0 + +# Set up the library. +algebra_add_library( algebra_vc_aos_generic vc_aos_generic + "include/algebra/vc_aos_generic.hpp" ) +target_link_libraries( algebra_vc_aos_generic + INTERFACE algebra::common algebra::vc_aos_storage algebra::generic_math + algebra::vc_aos_math algebra::cmath_math algebra::array_storage ) +algebra_test_public_headers( algebra_vc_aos_generic + "algebra/vc_aos_generic.hpp" ) diff --git a/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp b/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp new file mode 100644 index 00000000..9a69860f --- /dev/null +++ b/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp @@ -0,0 +1,91 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2020-2024 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/math/generic.hpp" +#include "algebra/math/vc_aos.hpp" +#include "algebra/storage/array.hpp" +#include "algebra/storage/vc_aos.hpp" + +/// @name Operators on @c algebra::vc_aos types +/// @{ + +using algebra::cmath::operator*; +using algebra::cmath::operator-; +using algebra::cmath::operator+; + +/// @} + +namespace algebra { + +namespace getter { + +/// @name Getter functions on @c algebra::vc_aos::matrix_type +/// @{ + +using vc_aos::storage::block; +using vc_aos::storage::element; +using vc_aos::storage::set_block; +using vc_aos::storage::vector; + +/// @} + +} // namespace getter + +namespace vector { + +/// @name Vector functions on @c algebra::vc_aos::storage_type +/// @{ + +using generic::math::cross; +using generic::math::dot; +using generic::math::eta; +using generic::math::norm; +using generic::math::normalize; +using generic::math::perp; +using generic::math::phi; +using generic::math::theta; + +/// @} + +} // namespace vector + +namespace matrix { + +/// @name Matrix functions on @c algebra::vc_aos::storage_type +/// @{ + +using generic::math::determinant; +using generic::math::identity; +using generic::math::inverse; +using generic::math::set_identity; +using generic::math::set_zero; +using generic::math::transpose; +using generic::math::zero; + +/// @} + +} // namespace matrix + +namespace vc_aos { + +/// @name generic based transforms on @c algebra::vc_aos +/// @{ + +template +using transform3 = + generic::math::transform3; + +/// @} + +} // namespace vc_aos + +} // namespace algebra diff --git a/frontend/vc_cmath/CMakeLists.txt b/frontend/vc_cmath/CMakeLists.txt deleted file mode 100644 index a20065a1..00000000 --- a/frontend/vc_cmath/CMakeLists.txt +++ /dev/null @@ -1,14 +0,0 @@ -# Algebra plugins library, part of the ACTS project (R&D line) -# -# (c) 2021-2023 CERN for the benefit of the ACTS project -# -# Mozilla Public License Version 2.0 - -# Set up the library. -algebra_add_library( algebra_vc_cmath vc_cmath - "include/algebra/vc_cmath.hpp" ) -target_link_libraries( algebra_vc_cmath - INTERFACE algebra::common algebra::vc_aos_storage algebra::cmath_math - algebra::vc_aos_math ) -algebra_test_public_headers( algebra_vc_cmath - "algebra/vc_cmath.hpp" ) diff --git a/frontend/vc_cmath/include/algebra/vc_cmath.hpp b/frontend/vc_cmath/include/algebra/vc_cmath.hpp deleted file mode 100644 index 31767d41..00000000 --- a/frontend/vc_cmath/include/algebra/vc_cmath.hpp +++ /dev/null @@ -1,112 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2020-2022 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/math/vc_aos.hpp" -#include "algebra/storage/vc_aos.hpp" - -/// @name Operators on @c algebra::vc_aos types -/// @{ - -using algebra::cmath::operator*; -using algebra::cmath::operator-; -using algebra::cmath::operator+; - -/// @} - -namespace algebra { -namespace getter { - -/// @name Getter functions on @c algebra::vc_aos types -/// @{ - -using cmath::eta; -using cmath::norm; -using cmath::perp; -using cmath::phi; -using cmath::theta; - -using vc_aos::math::eta; -using vc_aos::math::norm; -using vc_aos::math::perp; -using vc_aos::math::phi; -using vc_aos::math::theta; - -/// @| - -/// Function extracting a slice from the matrix used by -/// @c algebra::vc_aos::transform3 -template -ALGEBRA_HOST_DEVICE inline Vc::array vector( - const vc_aos::matrix_type& m, std::size_t row, - std::size_t col) { - - return cmath::vector_getter>()(m, row, col); -} - -/// @name Getter functions on @c algebra::vc_aos::matrix_type -/// @{ - -using cmath::element; - -/// @} - -} // namespace getter - -namespace vector { - -/// @name Vector functions on @c algebra::vc_aos::storage_type -/// @{ - -using cmath::cross; -using cmath::dot; -using cmath::normalize; - -using vc_aos::math::cross; -using vc_aos::math::dot; -using vc_aos::math::normalize; - -/// @} - -} // namespace vector - -namespace matrix { - -using cmath::block; -using cmath::determinant; -using cmath::identity; -using cmath::inverse; -using cmath::set_block; -using cmath::set_identity; -using cmath::set_zero; -using cmath::transpose; -using cmath::zero; - -} // namespace matrix - -namespace vc_aos { - -template -using element_getter = cmath::element_getter; - -/// @name cmath based transforms on @c algebra::vc - -template -using transform3 = - cmath::transform3, cmath::block_getter>; - -/// @} - -} // namespace vc_aos - -} // namespace algebra diff --git a/frontend/vc_soa/CMakeLists.txt b/frontend/vc_soa/CMakeLists.txt index d4980cca..097c4765 100644 --- a/frontend/vc_soa/CMakeLists.txt +++ b/frontend/vc_soa/CMakeLists.txt @@ -8,7 +8,7 @@ algebra_add_library( algebra_vc_soa vc_soa "include/algebra/vc_soa.hpp" ) target_link_libraries( algebra_vc_soa - INTERFACE algebra::common algebra::vc_soa_storage algebra::cmath_math - algebra::vc_soa_math algebra::vc_aos_math ) + INTERFACE algebra::common algebra::vc_soa_storage algebra::vc_soa_math + algebra::vc_aos_math ) algebra_test_public_headers( algebra_vc_soa "algebra/vc_soa.hpp" ) diff --git a/frontend/vc_soa/include/algebra/vc_soa.hpp b/frontend/vc_soa/include/algebra/vc_soa.hpp index 44b62116..a4277882 100644 --- a/frontend/vc_soa/include/algebra/vc_soa.hpp +++ b/frontend/vc_soa/include/algebra/vc_soa.hpp @@ -27,28 +27,16 @@ using algebra::storage::operator+; /// @} namespace algebra { -namespace vc_soa { - -/// @name Vc based transforms on @c algebra::vc_soa types -/// @{ - -template -using transform3 = algebra::vc_aos::math::transform3; - -/// @} - -} // namespace vc_soa namespace getter { /// @name Getter functions on @c algebra::vc_soa types /// @{ -using vc_soa::math::eta; -using vc_soa::math::norm; -using vc_soa::math::perp; -using vc_soa::math::phi; -using vc_soa::math::theta; +using vc_soa::storage::block; +using vc_soa::storage::element; +using vc_soa::storage::set_block; +using vc_soa::storage::vector; /// @} @@ -61,7 +49,12 @@ namespace vector { using vc_soa::math::cross; using vc_soa::math::dot; +using vc_soa::math::eta; +using vc_soa::math::norm; using vc_soa::math::normalize; +using vc_soa::math::perp; +using vc_soa::math::phi; +using vc_soa::math::theta; /// @} @@ -70,11 +63,9 @@ using vc_soa::math::normalize; // Produces clash with matrix typedefs in other plugins namespace matrix { -using vc_soa::math::block; using vc_soa::math::determinant; using vc_soa::math::identity; using vc_soa::math::inverse; -using vc_soa::math::set_block; using vc_soa::math::set_identity; using vc_soa::math::set_zero; using vc_soa::math::transpose; @@ -82,4 +73,17 @@ using vc_soa::math::zero; } // namespace matrix +namespace vc_soa { + +/// @name Vc based transforms on @c algebra::vc_soa types +/// @{ + +template +using transform3 = + algebra::vc_aos::math::transform3; + +/// @} + +} // namespace vc_soa + } // namespace algebra diff --git a/frontend/vecmem_cmath/CMakeLists.txt b/frontend/vecmem_cmath/CMakeLists.txt index 570bd7eb..80f2ad29 100644 --- a/frontend/vecmem_cmath/CMakeLists.txt +++ b/frontend/vecmem_cmath/CMakeLists.txt @@ -8,6 +8,6 @@ algebra_add_library( algebra_vecmem_cmath vecmem_cmath "include/algebra/vecmem_cmath.hpp" ) target_link_libraries( algebra_vecmem_cmath - INTERFACE algebra::common algebra::vecmem_storage algebra::cmath_math ) + INTERFACE algebra::common algebra::vecmem_storage algebra::cmath_math algebra::generic_math ) algebra_test_public_headers( algebra_vecmem_cmath "algebra/vecmem_cmath.hpp" ) diff --git a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp index 585f4140..edf65628 100644 --- a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp +++ b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.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) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -9,6 +9,7 @@ // Project include(s). #include "algebra/math/cmath.hpp" +#include "algebra/math/generic.hpp" #include "algebra/storage/vecmem.hpp" /// @name Operators on @c algebra::vecmem::storage_type @@ -24,33 +25,13 @@ namespace algebra { namespace getter { -/// @name Getter functions on @c algebra::vecmem::storage_type -/// @{ - -using cmath::eta; -using cmath::norm; -using cmath::perp; -using cmath::phi; -using cmath::theta; - -/// @} - -/// Function extracting a slice from the matrix used by -/// @c algebra::vecmem::transform3 -template -ALGEBRA_HOST_DEVICE inline vecmem::storage_type vector( - const vecmem::storage_type, COLS>& m, - std::size_t row, std::size_t col) { - - return cmath::vector_getter()(m, row, col); -} - /// @name Getter functions on @c algebra::vecmem::matrix_type /// @{ -using cmath::element; +using cmath::storage::block; +using cmath::storage::element; +using cmath::storage::set_block; +using cmath::storage::vector; /// @} @@ -61,41 +42,49 @@ namespace vector { /// @name Vector functions on @c algebra::vecmem::storage_type /// @{ -using cmath::cross; +// vecmem::static_array specific implementations using cmath::dot; using cmath::normalize; +// generic implementations +using generic::math::cross; +using generic::math::eta; +using generic::math::norm; +using generic::math::perp; +using generic::math::phi; +using generic::math::theta; + /// @} } // namespace vector namespace matrix { -using cmath::block; -using cmath::determinant; +/// @name Matrix functions on @c algebra::vecmem::storage_type +/// @{ + using cmath::identity; -using cmath::inverse; -using cmath::set_block; using cmath::set_identity; using cmath::set_zero; -using cmath::transpose; using cmath::zero; +using generic::math::determinant; +using generic::math::inverse; +using generic::math::transpose; + +/// @} + } // namespace matrix namespace vecmem { -template -using element_getter = - cmath::element_getter; - /// @name cmath based transforms on @c algebra::vecmem /// @{ template using transform3 = - cmath::transform3, cmath::block_getter>; + generic::math::transform3; /// @} diff --git a/math/CMakeLists.txt b/math/CMakeLists.txt index 42c480ab..b8eb8d20 100644 --- a/math/CMakeLists.txt +++ b/math/CMakeLists.txt @@ -1,12 +1,13 @@ # Algebra plugins library, part of the ACTS project (R&D line) # -# (c) 2021-2023 CERN for the benefit of the ACTS project +# (c) 2021-2024 CERN for the benefit of the ACTS project # # Mozilla Public License Version 2.0 # Set up all enabled libraries. add_subdirectory( common ) add_subdirectory( cmath ) +add_subdirectory( generic ) if( ALGEBRA_PLUGINS_INCLUDE_EIGEN ) add_subdirectory( eigen ) endif() diff --git a/math/cmath/CMakeLists.txt b/math/cmath/CMakeLists.txt index ff99a203..2bcc4057 100644 --- a/math/cmath/CMakeLists.txt +++ b/math/cmath/CMakeLists.txt @@ -1,6 +1,6 @@ # Algebra plugins library, part of the ACTS project (R&D line) # -# (c) 2021-2023 CERN for the benefit of the ACTS project +# (c) 2021-2024 CERN for the benefit of the ACTS project # # Mozilla Public License Version 2.0 @@ -8,21 +8,10 @@ algebra_add_library(algebra_cmath_math cmath_math # impl include "include/algebra/math/cmath.hpp" - "include/algebra/math/impl/cmath_getter.hpp" "include/algebra/math/impl/cmath_matrix.hpp" "include/algebra/math/impl/cmath_operators.hpp" - "include/algebra/math/impl/cmath_transform3.hpp" - "include/algebra/math/impl/cmath_vector.hpp" - # algorithms include - "include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp" - "include/algebra/math/algorithms/matrix/determinant/cofactor.hpp" - "include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp" - "include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp" - "include/algebra/math/algorithms/matrix/inverse/cofactor.hpp" - "include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp" - "include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp" - "include/algebra/math/algorithms/utils/algorithm_finder.hpp") + "include/algebra/math/impl/cmath_vector.hpp") target_link_libraries(algebra_cmath_math - INTERFACE algebra::common algebra::common_math) + INTERFACE algebra::common algebra::common_math algebra::generic_math) algebra_test_public_headers( algebra_cmath_math "algebra/math/cmath.hpp" ) diff --git a/math/cmath/include/algebra/math/cmath.hpp b/math/cmath/include/algebra/math/cmath.hpp index 946de33e..95a2bf3c 100644 --- a/math/cmath/include/algebra/math/cmath.hpp +++ b/math/cmath/include/algebra/math/cmath.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) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -8,17 +8,6 @@ #pragma once // Impl include(s). -#include "algebra/math/impl/cmath_getter.hpp" #include "algebra/math/impl/cmath_matrix.hpp" #include "algebra/math/impl/cmath_operators.hpp" -#include "algebra/math/impl/cmath_transform3.hpp" #include "algebra/math/impl/cmath_vector.hpp" -// Algorithms include(s). -#include "algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp" -#include "algebra/math/algorithms/matrix/determinant/cofactor.hpp" -#include "algebra/math/algorithms/matrix/determinant/hard_coded.hpp" -#include "algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp" -#include "algebra/math/algorithms/matrix/inverse/cofactor.hpp" -#include "algebra/math/algorithms/matrix/inverse/hard_coded.hpp" -#include "algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp" -#include "algebra/math/algorithms/utils/algorithm_finder.hpp" \ No newline at end of file diff --git a/math/cmath/include/algebra/math/impl/cmath_getter.hpp b/math/cmath/include/algebra/math/impl/cmath_getter.hpp deleted file mode 100644 index 03890272..00000000 --- a/math/cmath/include/algebra/math/impl/cmath_getter.hpp +++ /dev/null @@ -1,200 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2020-2022 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" - -// System include(s). -#include -#include -#include - -namespace algebra::cmath { - -/** This method retrieves phi from a vector, vector base with rows >= 2 - * - * @param v the input vector - **/ -template class array_t, - typename scalar_t, size_type N, - std::enable_if_t<(N >= 2) && std::is_scalar_v, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t phi( - const array_t &v) noexcept { - - return algebra::math::atan2(v[1], v[0]); -} - -/** This method retrieves theta from a vector, vector base with rows >= 3 - * - * @param v the input vector - **/ -template class array_t, - typename scalar_t, size_type N, - std::enable_if_t<(N >= 3) && std::is_scalar_v, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t theta( - const array_t &v) noexcept { - - return algebra::math::atan2(algebra::math::sqrt(v[0] * v[0] + v[1] * v[1]), - v[2]); -} - -/** This method retrieves the perpendicular magnitude of a vector with rows >= 2 - * - * @param v the input vector - **/ -template class array_t, - typename scalar_t, size_type N, - std::enable_if_t<(N >= 2) && std::is_scalar_v, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t perp( - const array_t &v) noexcept { - - return algebra::math::sqrt(v[0] * v[0] + v[1] * v[1]); -} - -/** This method retrieves the norm of a vector, no dimension restriction - * - * @param v the input vector - **/ -template class array_t, - typename scalar_t, size_type N, - std::enable_if_t<(N == 2) && std::is_scalar_v, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t norm(const array_t &v) { - - return perp(v); -} - -/** This method retrieves the norm of a vector, no dimension restriction - * - * @param v the input vector - **/ -template class array_t, - typename scalar_t, size_type N, - std::enable_if_t<(N >= 3) && std::is_scalar_v, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t norm(const array_t &v) { - - return algebra::math::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); -} - -/** This method retrieves the pseudo-rapidity from a vector or vector base with - *rows >= 3 - * - * @param v the input vector - **/ -template class array_t, - typename scalar_t, size_type N, - std::enable_if_t<(N >= 3) && std::is_scalar_v, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t eta( - const array_t &v) noexcept { - - return algebra::math::atanh(v[2] / norm(v)); -} - -/// "Element getter", assuming a simple 2D array access -template class array_t, - typename scalar_t> -struct element_getter { - - /// 2D matrix type - template - using matrix_type = array_t, COLS>; - - /// Operator getting a reference to one element of a non-const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t &operator()(matrix_type &m, - std::size_t row, - std::size_t col) const { - - assert(row < ROWS); - assert(col < COLS); - return m[col][row]; - } - - /// Operator getting one value of a const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type &m, std::size_t row, - std::size_t col) const { - - assert(row < ROWS); - assert(col < COLS); - return m[col][row]; - } -}; // struct element_getter - -/// Function extracting an element from a matrix (const) -template class array_t, - typename scalar_t, size_type ROWS, size_type COLS> -ALGEBRA_HOST_DEVICE inline scalar_t element( - const array_t, COLS> &m, std::size_t row, - std::size_t col) { - - return element_getter()(m, row, col); -} - -/// Function extracting an element from a matrix (non-const) -template class array_t, - typename scalar_t, size_type ROWS, size_type COLS> -ALGEBRA_HOST_DEVICE inline scalar_t &element( - array_t, COLS> &m, std::size_t row, - std::size_t col) { - - return element_getter()(m, row, col); -} - -/// "Vector getter", assuming a simple 2D array access -template class array_t, - typename scalar_t, size_type SIZE, - typename result_t = array_t> -struct vector_getter { - - /// Result type - using result_type = result_t; - /// 2D matrix type - template - using matrix_type = array_t, COLS>; - - /// Operator producing a vector out of a const matrix - template - ALGEBRA_HOST_DEVICE inline result_type operator()( - const matrix_type &m, std::size_t row, std::size_t col) { - - assert(col < COLS); - assert(row + SIZE <= ROWS); - result_type subvector{}; - for (std::size_t irow = row; irow < row + SIZE; ++irow) { - subvector[irow - row] = m[col][irow]; - } - return subvector; - } -}; // struct vector_getter - -/// "Block getter", assuming a simple 2D array access -struct block_getter { - - /// Operator producing a sub-matrix from a const matrix - template class array_t> - ALGEBRA_HOST_DEVICE auto operator()( - const array_t, oCOLS> &m, std::size_t row, - std::size_t col) const { - - array_t, COLS> submatrix{}; - - for (std::size_t icol = col; icol < col + COLS; ++icol) { - for (std::size_t irow = row; irow < row + ROWS; ++irow) { - submatrix[icol - col][irow - row] = m[icol][irow]; - } - } - return submatrix; - } -}; // struct block_getter - -} // namespace algebra::cmath diff --git a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp index 93073188..67ea5317 100644 --- a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp +++ b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp @@ -8,50 +8,12 @@ #pragma once // Project include(s). -#include "algebra/math/algorithms/utils/algorithm_finder.hpp" #include "algebra/math/common.hpp" #include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" namespace algebra::cmath { -/// @returns a matrix of dimension @tparam ROW and @tparam COL that contains -/// the submatrix of @param m beginning at row @param row and column -/// @param col -template -ALGEBRA_HOST_DEVICE decltype(auto) block(const input_matrix_type &m, - std::size_t row, std::size_t col) { - - return algebra::cmath::block_getter().template operator()(m, row, - col); -} - -/// Sets a matrix of dimension @tparam ROW and @tparam COL as submatrix of -/// @param m beginning at row @param row and column @param col -template class array_t> -ALGEBRA_HOST_DEVICE void set_block( - input_matrix_type &m, const array_t, COLS> &b, - std::size_t row, std::size_t col) { - for (std::size_t j = 0u; j < COLS; ++j) { - for (std::size_t i = 0u; i < ROWS; ++i) { - m[j + col][i + row] = b[j][i]; - } - } -} - -/// Sets a vector of length @tparam ROW as submatrix of -/// @param m beginning at row @param row and column @param col -template class vector_t, - class input_matrix_type> -ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const vector_t &b, - std::size_t row, std::size_t col) { - for (std::size_t i = 0; i < ROWS; ++i) { - m[col][i + row] = b[i]; - } -} - /// @returns zero matrix of type @tparam matrix_t template requires(std::is_scalar_v) @@ -59,20 +21,6 @@ requires(std::is_scalar_v) return matrix_t{}; } -/// Create zero matrix - cmath transform3 -template -ALGEBRA_HOST_DEVICE inline auto zero() { - matrix_t ret; - - for (std::size_t j = 0; j < algebra::trait::columns; ++j) { - for (std::size_t i = 0; i < algebra::trait::rows; ++i) { - element_getter_t{}(ret, i, j) = 0; - } - } - - return ret; -} - /// @returns identity matrix of type @tparam matrix_t template requires(std::is_scalar_v) @@ -86,22 +34,10 @@ requires(std::is_scalar_v) return ret; } -/// Create identity matrix - cmath transform3 -template -ALGEBRA_HOST_DEVICE inline matrix_t identity() { - auto ret{zero()}; - - for (std::size_t i = 0; i < algebra::trait::rank; ++i) { - element_getter_t{}(ret, i, i) = 1; - } - - return ret; -} - /// Set @param m as zero matrix template class array_t> -ALGEBRA_HOST_DEVICE inline void set_zero( +ALGEBRA_HOST_DEVICE constexpr void set_zero( array_t, COLS> &m) { m = zero, COLS>>(); } @@ -109,50 +45,9 @@ ALGEBRA_HOST_DEVICE inline void set_zero( /// Set @param m as identity matrix template class array_t> -ALGEBRA_HOST_DEVICE inline void set_identity( +ALGEBRA_HOST_DEVICE constexpr void set_identity( array_t, COLS> &m) { m = identity, COLS>>(); } -/// @returns the transpose matrix of @param m -template class array_t> -ALGEBRA_HOST_DEVICE inline array_t, ROWS> transpose( - const array_t, COLS> &m) { - - array_t, ROWS> ret; - - for (std::size_t i = 0; i < ROWS; ++i) { - for (std::size_t j = 0; j < COLS; ++j) { - ret[i][j] = m[j][i]; - } - } - - return ret; -} - -/// @returns the determinant of @param m -template class array_t> -ALGEBRA_HOST_DEVICE inline scalar_t determinant( - const array_t, N> &m) { - - using matrix_t = array_t, N>; - using element_getter_t = element_getter; - - return determinant_t{}(m); -} - -/// @returns the determinant of @param m -template class array_t> -ALGEBRA_HOST_DEVICE inline array_t, N> inverse( - const array_t, N> &m) { - - using matrix_t = array_t, N>; - using element_getter_t = element_getter; - - return inversion_t{}(m); -} - } // namespace algebra::cmath diff --git a/math/cmath/include/algebra/math/impl/cmath_operators.hpp b/math/cmath/include/algebra/math/impl/cmath_operators.hpp index a38581dc..bb4a681f 100644 --- a/math/cmath/include/algebra/math/impl/cmath_operators.hpp +++ b/math/cmath/include/algebra/math/impl/cmath_operators.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (c) 2020-2023 CERN for the benefit of the ACTS project + * (c) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ diff --git a/math/cmath/include/algebra/math/impl/cmath_transform3.hpp b/math/cmath/include/algebra/math/impl/cmath_transform3.hpp deleted file mode 100644 index 950ed8db..00000000 --- a/math/cmath/include/algebra/math/impl/cmath_transform3.hpp +++ /dev/null @@ -1,312 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2020-2024 CERN for the benefit of the ACTS project - * - * Mozilla Public License Version 2.0 - */ - -#pragma once - -// Project include(s). -#include "algebra/math/algorithms/matrix/inverse/hard_coded.hpp" -#include "algebra/math/impl/cmath_getter.hpp" -#include "algebra/math/impl/cmath_matrix.hpp" -#include "algebra/math/impl/cmath_vector.hpp" -#include "algebra/qualifiers.hpp" - -namespace algebra::cmath { - -/** Transform wrapper class to ensure standard API within differnt plugins - **/ -template class matrix_t, - template class array_t, typename element_getter_t, - typename block_getter_t> -struct transform3 { - - /// @name Type definitions for the struct - /// @{ - - /// Size type - using size_type = size_ty; - - /// Scalar type - using scalar_type = scalar_t; - - /// 2D Matrix type - template - using matrix_type = matrix_t; - - /// Array type - template - using array_type = array_t; - - // 4 x 4 Matrix - using matrix44 = matrix_type<4, 4>; - - /// 3-element "vector" type - using vector3 = array_type<3>; - /// Point in 3D space - using point3 = vector3; - /// Point in 2D space - using point2 = array_type<2>; - - /// Function (object) used for accessing a matrix element - using element_getter = element_getter_t; - - /// Function (object) used for accessing a matrix element - using block_getter = block_getter_t; - - /// Matrix inversion algorithm - using matrix_inversion = - cmath::matrix::inverse::hard_coded; - - /// @} - - /// @name Data objects - /// @{ - matrix44 _data{cmath::identity()}; - matrix44 _data_inv{cmath::identity()}; - - /// @} - - /** Default constructor: identity - **/ - constexpr 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, bool get_inverse = true) { - - element_getter{}(_data, 0, 0) = x[0]; - element_getter{}(_data, 1, 0) = x[1]; - element_getter{}(_data, 2, 0) = x[2]; - element_getter{}(_data, 3, 0) = 0.; - element_getter{}(_data, 0, 1) = y[0]; - element_getter{}(_data, 1, 1) = y[1]; - element_getter{}(_data, 2, 1) = y[2]; - element_getter{}(_data, 3, 1) = 0.; - element_getter{}(_data, 0, 2) = z[0]; - element_getter{}(_data, 1, 2) = z[1]; - element_getter{}(_data, 2, 2) = z[2]; - element_getter{}(_data, 3, 2) = 0.; - element_getter{}(_data, 0, 3) = t[0]; - element_getter{}(_data, 1, 3) = t[1]; - element_getter{}(_data, 2, 3) = t[2]; - element_getter{}(_data, 3, 3) = 1.; - - if (get_inverse) { - _data_inv = matrix_inversion{}(_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 - explicit transform3(const vector3 &t) { - - element_getter{}(_data, 0, 0) = 1.; - element_getter{}(_data, 1, 0) = 0.; - element_getter{}(_data, 2, 0) = 0.; - element_getter{}(_data, 3, 0) = 0.; - element_getter{}(_data, 0, 1) = 0.; - element_getter{}(_data, 1, 1) = 1.; - element_getter{}(_data, 2, 1) = 0.; - element_getter{}(_data, 3, 1) = 0.; - element_getter{}(_data, 0, 2) = 0.; - element_getter{}(_data, 1, 2) = 0.; - element_getter{}(_data, 2, 2) = 1.; - element_getter{}(_data, 3, 2) = 0.; - element_getter{}(_data, 0, 3) = t[0]; - element_getter{}(_data, 1, 3) = t[1]; - element_getter{}(_data, 2, 3) = t[2]; - element_getter{}(_data, 3, 3) = 1.; - - _data_inv = matrix_inversion{}(_data); - } - - /** Constructor with arguments: matrix - * - * @param m is the full 4x4 matrix - **/ - ALGEBRA_HOST_DEVICE - explicit transform3(const matrix44 &m) { - - _data = m; - _data_inv = matrix_inversion{}(_data); - } - - /** Constructor with arguments: matrix as array of scalar - * - * @param ma is the full 4x4 matrix 16 array - **/ - ALGEBRA_HOST_DEVICE - explicit transform3(const array_type<16> &ma) { - - element_getter{}(_data, 0, 0) = ma[0]; - element_getter{}(_data, 1, 0) = ma[4]; - element_getter{}(_data, 2, 0) = ma[8]; - element_getter{}(_data, 3, 0) = ma[12]; - element_getter{}(_data, 0, 1) = ma[1]; - element_getter{}(_data, 1, 1) = ma[5]; - element_getter{}(_data, 2, 1) = ma[9]; - element_getter{}(_data, 3, 1) = ma[13]; - element_getter{}(_data, 0, 2) = ma[2]; - element_getter{}(_data, 1, 2) = ma[6]; - element_getter{}(_data, 2, 2) = ma[10]; - element_getter{}(_data, 3, 2) = ma[14]; - element_getter{}(_data, 0, 3) = ma[3]; - element_getter{}(_data, 1, 3) = ma[7]; - element_getter{}(_data, 2, 3) = ma[11]; - element_getter{}(_data, 3, 3) = ma[15]; - - _data_inv = matrix_inversion{}(_data); - } - - /** Equality operator */ - ALGEBRA_HOST_DEVICE - inline bool operator==(const transform3 &rhs) const { - - for (size_type j = 0; j < 4; j++) { - for (size_type i = 0; i < 4; i++) { - if (element_getter{}(_data, i, j) != - element_getter{}(rhs._data, i, j)) { - return false; - } - } - } - - return true; - } - - /** Rotate a vector into / from a frame - * - * @param m is the rotation matrix - * @param v is the vector to be rotated - */ - ALGEBRA_HOST_DEVICE - static inline vector3 rotate(const matrix44 &m, const vector3 &v) { - - vector3 ret{0.f, 0.f, 0.f}; - - ret[0] += element_getter{}(m, 0, 0) * v[0]; - ret[1] += element_getter{}(m, 1, 0) * v[0]; - ret[2] += element_getter{}(m, 2, 0) * v[0]; - - ret[0] += element_getter{}(m, 0, 1) * v[1]; - ret[1] += element_getter{}(m, 1, 1) * v[1]; - ret[2] += element_getter{}(m, 2, 1) * v[1]; - - ret[0] += element_getter{}(m, 0, 2) * v[2]; - ret[1] += element_getter{}(m, 1, 2) * v[2]; - ret[2] += element_getter{}(m, 2, 2) * v[2]; - - return ret; - } - - /** This method retrieves the rotation of a transform */ - ALGEBRA_HOST_DEVICE - auto inline rotation() const { - - return block_getter{}.template operator()<3, 3>(_data, 0, 0); - } - - /** This method retrieves x axis */ - ALGEBRA_HOST_DEVICE - inline point3 x() const { - return {element_getter{}(_data, 0, 0), element_getter{}(_data, 1, 0), - element_getter{}(_data, 2, 0)}; - } - - /** This method retrieves y axis */ - ALGEBRA_HOST_DEVICE - inline point3 y() const { - return {element_getter{}(_data, 0, 1), element_getter{}(_data, 1, 1), - element_getter{}(_data, 2, 1)}; - } - - /** This method retrieves z axis */ - ALGEBRA_HOST_DEVICE - inline point3 z() const { - return {element_getter{}(_data, 0, 2), element_getter{}(_data, 1, 2), - element_getter{}(_data, 2, 2)}; - } - - /** This method retrieves the translation of a transform */ - ALGEBRA_HOST_DEVICE - inline point3 translation() const { - - return {element_getter{}(_data, 0, 3), element_getter{}(_data, 1, 3), - element_getter{}(_data, 2, 3)}; - } - - /** This method retrieves the 4x4 matrix of a transform */ - ALGEBRA_HOST_DEVICE - inline const matrix44 &matrix() const { return _data; } - - /** This method retrieves the 4x4 matrix of an inverse transform */ - ALGEBRA_HOST_DEVICE - inline const matrix44 &matrix_inverse() const { return _data_inv; } - - /** This method transform from a point from the local 3D cartesian frame to - * the global 3D cartesian frame */ - ALGEBRA_HOST_DEVICE inline point3 point_to_global(const point3 &v) const { - - const vector3 rg = rotate(_data, v); - return {rg[0] + element_getter{}(_data, 0, 3), - rg[1] + element_getter{}(_data, 1, 3), - rg[2] + element_getter{}(_data, 2, 3)}; - } - - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ - - ALGEBRA_HOST_DEVICE inline point3 point_to_local(const point3 &v) const { - - const vector3 rg = rotate(_data_inv, v); - return {rg[0] + element_getter{}(_data_inv, 0, 3), - rg[1] + element_getter{}(_data_inv, 1, 3), - rg[2] + element_getter{}(_data_inv, 2, 3)}; - } - - /** This method transform from a vector from the local 3D cartesian frame to - * the global 3D cartesian frame */ - ALGEBRA_HOST_DEVICE inline vector3 vector_to_global(const vector3 &v) const { - - return rotate(_data, v); - } - - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ - ALGEBRA_HOST_DEVICE inline vector3 vector_to_local(const vector3 &v) const { - - return rotate(_data_inv, v); - } - -}; // struct transform3 - -} // namespace algebra::cmath diff --git a/math/cmath/include/algebra/math/impl/cmath_vector.hpp b/math/cmath/include/algebra/math/impl/cmath_vector.hpp index 0ecf3ff7..a4aad3e6 100644 --- a/math/cmath/include/algebra/math/impl/cmath_vector.hpp +++ b/math/cmath/include/algebra/math/impl/cmath_vector.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) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -13,188 +13,126 @@ namespace algebra::cmath { -/** Cross product between two input vectors - 3 Dim - * - * @tparam derived_type_lhs is the first matrix (epresseion) template - * @tparam derived_type_rhs is the second matrix (epresseion) template - * - * @param a the first input vector - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template class array_t, - typename scalar_t, size_type N, - std::enable_if_t= 3 && std::is_scalar_v, bool> = true> -ALGEBRA_HOST_DEVICE inline array_t cross( + typename scalar_t, size_type N> +requires std::is_scalar_v ALGEBRA_HOST_DEVICE inline scalar_t dot( const array_t &a, const array_t &b) { - - return {a[1] * b[2] - b[1] * a[2], a[2] * b[0] - b[2] * a[0], - a[0] * b[1] - b[0] * a[1]}; -} - -/** Cross product between two input vectors - 3 Dim - * - * @tparam derived_type_lhs is the first matrix (epresseion) template - * @tparam derived_type_rhs is the second matrix (epresseion) template - * - * @param a the first input vector - * @param b the second input matrix with single column - * - * @return a vector (expression) representing the cross product - **/ -template class array_t, - typename scalar_t, size_type N, size_type COLS, - std::enable_if_t= 3 && COLS == 1 && std::is_scalar_v, - bool> = true> -ALGEBRA_HOST_DEVICE inline array_t cross( - const array_t &a, - const array_t, COLS> &b) { - - return {a[1] * b[0][2] - b[0][1] * a[2], a[2] * b[0][0] - b[0][2] * a[0], - a[0] * b[0][1] - b[0][0] * a[1]}; -} - -/** Cross product between two input vectors - 3 Dim - * - * @tparam derived_type_lhs is the first matrix (epresseion) template - * @tparam derived_type_rhs is the second matrix (epresseion) template - * - * @param a the first input matrix with single column - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ -template class array_t, - typename scalar_t, size_type N, size_type COLS, - std::enable_if_t= 3 && COLS == 1 && std::is_scalar_v, - bool> = true> -ALGEBRA_HOST_DEVICE inline array_t cross( - const array_t, COLS> &a, - const array_t &b) { - - return {a[0][1] * b[2] - b[1] * a[0][2], a[0][2] * b[0] - b[2] * a[0][0], - a[0][0] * b[1] - b[0] * a[0][1]}; -} - -/** Cross product between two input vectors - 3 Dim - * - * @tparam derived_type_lhs is the first matrix (epresseion) template - * @tparam derived_type_rhs is the second matrix (epresseion) template - * - * @param a the first input matrix with single column - * @param b the second input matrix with single column - * - * @return a vector (expression) representing the cross product - **/ -template class array_t, - typename scalar_t, size_type N, size_type COLS, - std::enable_if_t= 3 && COLS == 1 && std::is_scalar_v, - bool> = true> -ALGEBRA_HOST_DEVICE inline array_t cross( - const array_t, COLS> &a, - const array_t, COLS> &b) { - - return {a[0][1] * b[0][2] - b[0][1] * a[0][2], - a[0][2] * b[0][0] - b[0][2] * a[0][0], - a[0][0] * b[0][1] - b[0][0] * a[0][1]}; -} - -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ -template class array_t, - typename scalar_t, size_type N, - std::enable_if_t, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t dot(const array_t &a, - const array_t &b) { - scalar_t ret = 0; + array_t tmp; + for (size_type i = 0; i < N; i++) { + tmp[i] = a[i] * b[i]; + } + scalar_t ret{0.f}; for (size_type i = 0; i < N; i++) { - ret += a[i] * b[i]; + ret += tmp[i]; } return ret; } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input matrix with single column - * - * @return the scalar dot product value - **/ -template < - typename size_type, template class array_t, - typename scalar_t, size_type N, size_type COLS, - std::enable_if_t, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t dot( - const array_t &a, - const array_t, COLS> &b) { - scalar_t ret = 0; +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input matrix with single column +/// +/// @return the scalar dot product value +template class array_t, + typename scalar_t, size_type N, size_type COLS> +requires(COLS == 1 && std::is_scalar_v) ALGEBRA_HOST_DEVICE + inline scalar_t dot(const array_t &a, + const array_t, COLS> &b) { + array_t tmp; for (size_type i = 0; i < N; i++) { - ret += a[i] * b[0][i]; + tmp[i] = a[i] * b[0][i]; + } + scalar_t ret{0.f}; + for (size_type i = 0; i < N; i++) { + ret += tmp[i]; } return ret; } -/** Dot product between two input vectors - * - * @param a the first input matrix with single column - * @param b the second input vector - * - * @return the scalar dot product value - **/ -template < - typename size_type, template class array_t, - typename scalar_t, size_type N, size_type COLS, - std::enable_if_t, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t dot( - const array_t, COLS> &a, - const array_t &b) { - scalar_t ret = 0; +/// Dot product between two input vectors +/// +/// @param a the first input matrix with single column +/// @param b the second input vector +/// +/// @return the scalar dot product value +template class array_t, + typename scalar_t, size_type N, size_type COLS> +requires(COLS == 1 && std::is_scalar_v) ALGEBRA_HOST_DEVICE + inline scalar_t dot(const array_t, COLS> &a, + const array_t &b) { + array_t tmp; + for (size_type i = 0; i < N; i++) { + tmp[i] = a[0][i] * b[i]; + } + scalar_t ret{0.f}; for (size_type i = 0; i < N; i++) { - ret += a[0][i] * b[i]; + ret += tmp[i]; } return ret; } -/** Dot product between two input vectors - * - * @param a the first input matrix with single column - * @param b the second input matrix with single column - * - * @return the scalar dot product value - **/ -template < - typename size_type, template class array_t, - typename scalar_t, size_type N, size_type COLS, - std::enable_if_t, bool> = true> -ALGEBRA_HOST_DEVICE inline scalar_t dot( - const array_t, COLS> &a, - const array_t, COLS> &b) { - scalar_t ret = 0; +/// Dot product between two input vectors +/// +/// @param a the first input matrix with single column +/// @param b the second input matrix with single column +/// +/// @return the scalar dot product value +template class array_t, + typename scalar_t, size_type N, size_type COLS> +requires(COLS == 1 && std::is_scalar_v) ALGEBRA_HOST_DEVICE + inline scalar_t dot(const array_t, COLS> &a, + const array_t, COLS> &b) { + array_t tmp; for (size_type i = 0; i < N; i++) { - ret += a[0][i] * b[0][i]; + tmp[i] = a[0][i] * b[0][i]; + } + scalar_t ret{0.f}; + for (size_type i = 0; i < N; i++) { + ret += tmp[i]; } return ret; } -/** Get a normalized version of the input vector - * - * @param v the input vector - **/ +/// This method retrieves the norm of a vector, no dimension restriction +/// +/// @param v the input vector +template class array_t, + typename scalar_t, size_type N> +requires(N >= 2 && std::is_scalar_v) ALGEBRA_HOST_DEVICE + inline scalar_t norm(const array_t &v) { + + return algebra::math::sqrt(dot(v, v)); +} + +/// This method retrieves the pseudo-rapidity from a vector or vector base with +/// rows >= 3 +/// +/// @param v the input vector +template class array_t, + typename scalar_t, size_type N> +requires(N >= 3 && std::is_scalar_v) ALGEBRA_HOST_DEVICE + inline scalar_t eta(const array_t &v) noexcept { + + return algebra::math::atanh(v[2] / norm(v)); +} + +/// Get a normalized version of the input vector +/// +/// @param v the input vector template class array_t, - typename scalar_t, size_type N, - std::enable_if_t, bool> = true> -ALGEBRA_HOST_DEVICE inline array_t normalize( - const array_t &v) { + typename scalar_t, size_type N> +requires std::is_scalar_v + ALGEBRA_HOST_DEVICE inline array_t normalize( + const array_t &v) { - const scalar_t oon = - static_cast(1.) / algebra::math::sqrt(dot(v, v)); - return v * oon; + return (static_cast(1.) / norm(v)) * v; } } // namespace algebra::cmath diff --git a/math/common/include/algebra/math/common.hpp b/math/common/include/algebra/math/common.hpp index d1616c24..632413d2 100644 --- a/math/common/include/algebra/math/common.hpp +++ b/math/common/include/algebra/math/common.hpp @@ -30,10 +30,16 @@ namespace math_ns = std; /// Absolute value of arg template -ALGEBRA_HOST_DEVICE inline scalar_t fabs(scalar_t arg) { +ALGEBRA_HOST_DEVICE inline auto fabs(scalar_t arg) { return math_ns::fabs(arg); } +/// Fused multiply add +template +ALGEBRA_HOST_DEVICE inline auto fma(scalar_t x, scalar_t y, scalar_t z) { + return math_ns::fma(x, y, z); +} + /// Arc tangent of y/x template ALGEBRA_HOST_DEVICE inline scalar_t atan2(scalar_t y, scalar_t x) { diff --git a/math/eigen/CMakeLists.txt b/math/eigen/CMakeLists.txt index d74dd643..1b71d08e 100644 --- a/math/eigen/CMakeLists.txt +++ b/math/eigen/CMakeLists.txt @@ -1,13 +1,12 @@ # Algebra plugins library, part of the ACTS project (R&D line) # -# (c) 2021-2023 CERN for the benefit of the ACTS project +# (c) 2021-2024 CERN for the benefit of the ACTS project # # Mozilla Public License Version 2.0 # Set up the library. algebra_add_library(algebra_eigen_math eigen_math "include/algebra/math/eigen.hpp" - "include/algebra/math/impl/eigen_getter.hpp" "include/algebra/math/impl/eigen_matrix.hpp" "include/algebra/math/impl/eigen_transform3.hpp" "include/algebra/math/impl/eigen_vector.hpp") diff --git a/math/eigen/include/algebra/math/eigen.hpp b/math/eigen/include/algebra/math/eigen.hpp index 975cc525..7659f771 100644 --- a/math/eigen/include/algebra/math/eigen.hpp +++ b/math/eigen/include/algebra/math/eigen.hpp @@ -8,7 +8,7 @@ #pragma once // Project include(s). -#include "algebra/math/impl/eigen_getter.hpp" #include "algebra/math/impl/eigen_matrix.hpp" #include "algebra/math/impl/eigen_transform3.hpp" #include "algebra/math/impl/eigen_vector.hpp" +#include "algebra/storage/impl/eigen_getter.hpp" diff --git a/math/eigen/include/algebra/math/impl/eigen_getter.hpp b/math/eigen/include/algebra/math/impl/eigen_getter.hpp deleted file mode 100644 index 975ad0a5..00000000 --- a/math/eigen/include/algebra/math/impl/eigen_getter.hpp +++ /dev/null @@ -1,177 +0,0 @@ -/** Algebra plugins, part of the ACTS project - * - * (c) 2020-2023 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" - -// Eigen include(s). -#ifdef _MSC_VER -#pragma warning(push, 0) -#endif // MSVC -#ifdef __NVCC_DIAG_PRAGMA_SUPPORT__ -#pragma nv_diagnostic push -#pragma nv_diag_suppress 20012 -#endif // __NVCC_DIAG_PRAGMA_SUPPORT__ -#include -#ifdef _MSC_VER -#pragma warning(pop) -#endif // MSVC -#ifdef __NVCC_DIAG_PRAGMA_SUPPORT__ -#pragma nv_diagnostic pop -#endif // __NVCC_DIAG_PRAGMA_SUPPORT__ - -// System include(s). -#include - -namespace algebra::eigen::math { - -/** This method retrieves phi from a vector, vector base with rows >= 2 - * - * @param v the input vector - **/ -template < - typename derived_type, - std::enable_if_t::RowsAtCompileTime >= 2, - bool> = true> -ALGEBRA_HOST_DEVICE inline auto phi( - const Eigen::MatrixBase &v) noexcept { - - return algebra::math::atan2(v[1], v[0]); -} - -/** This method retrieves theta from a vector, vector base with rows >= 3 - * - * @param v the input vector - **/ -template < - typename derived_type, - std::enable_if_t::RowsAtCompileTime >= 3, - bool> = true> -ALGEBRA_HOST_DEVICE inline auto theta( - const Eigen::MatrixBase &v) noexcept { - - return algebra::math::atan2(algebra::math::sqrt(v[0] * v[0] + v[1] * v[1]), - v[2]); -} - -/** This method retrieves the perpendicular magnitude of a vector with rows >= 2 - * - * @param v the input vector - **/ -template < - typename derived_type, - std::enable_if_t::RowsAtCompileTime >= 2, - bool> = true> -ALGEBRA_HOST_DEVICE inline auto perp( - const Eigen::MatrixBase &v) noexcept { - - return algebra::math::sqrt(v[0] * v[0] + v[1] * v[1]); -} - -/** This method retrieves the norm of a vector, no dimension restriction - * - * @param v the input vector - **/ -template -ALGEBRA_HOST_DEVICE inline auto norm(const Eigen::MatrixBase &v) { - - return v.norm(); -} - -/** This method retrieves the pseudo-rapidity from a vector or vector base with - * rows >= 3 - * - * @param v the input vector - **/ -template < - typename derived_type, - std::enable_if_t::RowsAtCompileTime >= 3, - bool> = true> -ALGEBRA_HOST_DEVICE inline auto eta( - const Eigen::MatrixBase &v) noexcept { - - return algebra::math::atanh(v[2] / v.norm()); -} - -/// Functor used to access elements of Eigen matrices -struct element_getter { - /// Get non-const access to a matrix element - template , - Eigen::MatrixBase >::value && - std::is_convertible::value && - std::is_convertible::value, - bool> = true> - ALGEBRA_HOST_DEVICE inline auto &operator()( - Eigen::MatrixBase &m, size_type_1 row, - size_type_2 col) const { - - return m(static_cast(row), static_cast(col)); - } - /// Get const access to a matrix element - template ::value && - std::is_convertible::value, - bool> = true> - ALGEBRA_HOST_DEVICE inline auto operator()( - const Eigen::MatrixBase &m, size_type_1 row, - size_type_2 col) const { - - return m(static_cast(row), static_cast(col)); - } -}; // struct element_getter - -/// Function extracting an element from a matrix (const) -template < - typename derived_type, typename size_type_1, typename size_type_2, - std::enable_if_t::value && - std::is_convertible::value, - bool> = true> -ALGEBRA_HOST_DEVICE inline auto element( - const Eigen::MatrixBase &m, size_type_1 row, - size_type_2 col) { - - return element_getter()(m, row, col); -} - -/// Function extracting an element from a matrix (non-const) -template , - Eigen::MatrixBase >::value && - std::is_convertible::value && - std::is_convertible::value, - bool> = true> -ALGEBRA_HOST_DEVICE inline auto &element(Eigen::MatrixBase &m, - size_type_1 row, size_type_2 col) { - - return element_getter()(m, row, col); -} - -/// Functor used to extract a block from Eigen matrices -struct block_getter { - template ::value && - std::is_convertible::value, - bool> = true> - ALGEBRA_HOST_DEVICE auto operator()(const matrix_type &m, size_type_1 row, - size_type_2 col) const { - - return m.template block(row, col); - } -}; // struct block_getter - -} // namespace algebra::eigen::math \ No newline at end of file diff --git a/math/eigen/include/algebra/math/impl/eigen_matrix.hpp b/math/eigen/include/algebra/math/impl/eigen_matrix.hpp index 30ee8d67..35455ffe 100644 --- a/math/eigen/include/algebra/math/impl/eigen_matrix.hpp +++ b/math/eigen/include/algebra/math/impl/eigen_matrix.hpp @@ -29,70 +29,32 @@ namespace algebra::eigen::math { -/// Operator getting a block of a const matrix -template -requires(std::is_convertible_v - &&std::is_convertible_v) - ALGEBRA_HOST_DEVICE auto block(const Eigen::MatrixBase &m, - size_type_1 row, size_type_2 col) { - return m.template block(static_cast(row), - static_cast(col)); -} - -/// Operator getting a block of a const matrix -template -requires(std::is_convertible_v - &&std::is_convertible_v) - ALGEBRA_HOST_DEVICE auto block(Eigen::MatrixBase &m, - size_type_1 row, size_type_2 col) { - return m.template block(static_cast(row), - static_cast(col)); -} - -/// Operator setting a block -template -requires(std::is_convertible_v - &&std::is_convertible_v) - ALGEBRA_HOST_DEVICE - void set_block(Eigen::MatrixBase &m, - const Eigen::MatrixBase &b, size_type_1 row, - size_type_2 col) { - using block_t = Eigen::MatrixBase; - constexpr auto R{block_t::RowsAtCompileTime}; - constexpr auto C{block_t::ColsAtCompileTime}; - m.template block(static_cast(row), - static_cast(col)) = b; -} - -// Create zero matrix +/// Create zero matrix template ALGEBRA_HOST_DEVICE inline matrix_t zero() { return matrix_t::Zero(); } -// Create identity matrix +/// Create identity matrix template ALGEBRA_HOST_DEVICE inline matrix_t identity() { return matrix_t::Identity(); } -// Set input matrix as zero matrix +/// Set input matrix as zero matrix template ALGEBRA_HOST_DEVICE inline void set_zero(Eigen::MatrixBase &m) { m.setZero(); } -// Set input matrix as identity matrix +/// Set input matrix as identity matrix template ALGEBRA_HOST_DEVICE inline void set_identity( Eigen::MatrixBase &m) { m.setIdentity(); } -// Create transpose matrix +/// Create transpose matrix template ALGEBRA_HOST_DEVICE inline matrix_type< typename Eigen::MatrixBase::value_type, diff --git a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp index 02b3378f..0c155944 100644 --- a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp +++ b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp @@ -34,28 +34,19 @@ namespace algebra::eigen::math { -/** Transform wrapper class to ensure standard API within differnt plugins */ +/// Transform wrapper class to ensure standard API within differnt plugins template struct transform3 { /// @name Type definitions for the struct /// @{ - /// Size type - using size_type = int; + /// Scalar type used by the transform using scalar_type = scalar_t; - /// 2D Matrix type - template - using matrix_type = Eigen::Matrix; - - /// 4x4 matrix type - using matrix44 = - typename Eigen::Transform::MatrixType; - /// Array type used by the transform template - using array_type = Eigen::Matrix; + using array_type = Eigen::Matrix; /// 3-element "vector" type using vector3 = array_type<3>; @@ -64,6 +55,10 @@ struct transform3 { /// Point in 2D space using point2 = array_type<2>; + /// 4x4 matrix type + using matrix44 = + typename Eigen::Transform::MatrixType; + /// @} /// @name Data objects @@ -76,14 +71,12 @@ struct transform3 { /// @} - /** 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 - * - **/ + /// 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, bool get_inverse = true) { @@ -98,22 +91,19 @@ struct transform3 { } } - /** 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 - * - **/ + /// 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 ALGEBRA_HOST_DEVICE transform3(const vector3 &t, const vector3 &z, const vector3 &x, bool get_inverse = true) : transform3(t, x, z.cross(x), z, get_inverse) {} - /** Constructor with arguments: translation - * - * @param t is the transform - **/ + /// Constructor with arguments: translation + /// + /// @param t is the transform ALGEBRA_HOST_DEVICE explicit transform3(const vector3 &t) { _data.setIdentity(); @@ -124,10 +114,9 @@ struct transform3 { _data_inv = _data.inverse(); } - /** Constructor with arguments: matrix - * - * @param m is the full 4x4 matrix - **/ + /// Constructor with arguments: matrix + /// + /// @param m is the full 4x4 matrix ALGEBRA_HOST_DEVICE explicit transform3(const matrix44 &m) { @@ -136,10 +125,9 @@ struct transform3 { _data_inv = _data.inverse(); } - /** Constructor with arguments: matrix as std::aray of scalar - * - * @param ma is the full 4x4 matrix as a 16 array - **/ + /// Constructor with arguments: matrix as std::aray of scalar + /// + /// @param ma is the full 4x4 matrix as a 16 array ALGEBRA_HOST_DEVICE explicit transform3(const array_type<16> &ma) { @@ -149,29 +137,28 @@ struct transform3 { _data_inv = _data.inverse(); } - /** Default constructor: set contents to identity matrices */ + /// Default constructor: set contents to identity matrices ALGEBRA_HOST_DEVICE transform3() { _data.setIdentity(); _data_inv.setIdentity(); } - /** Default contructors */ + /// Default contructors transform3(const transform3 &rhs) = default; ~transform3() = default; - /** Equality operator */ + /// Equality operator ALGEBRA_HOST_DEVICE inline bool operator==(const transform3 &rhs) const { return (_data.isApprox(rhs._data)); } - /** Rotate a vector into / from a frame - * - * @param m is the rotation matrix - * @param v is the vector to be rotated - */ + /// Rotate a vector into / from a frame + /// + /// @param m is the rotation matrix + /// @param v is the vector to be rotated template < typename derived_type, std::enable_if_t::RowsAtCompileTime == 3, @@ -185,42 +172,42 @@ struct transform3 { return m.matrix().template block<3, 3>(0, 0) * v; } - /** This method retrieves the rotation of a transform **/ + /// This method retrieves the rotation of a transform ALGEBRA_HOST_DEVICE inline auto rotation() const { return _data.matrix().template block<3, 3>(0, 0); } - /** This method retrieves x axis */ + /// This method retrieves x axis ALGEBRA_HOST_DEVICE inline point3 x() const { return _data.matrix().template block<3, 1>(0, 0); } - /** This method retrieves y axis */ + /// This method retrieves y axis ALGEBRA_HOST_DEVICE inline point3 y() const { return _data.matrix().template block<3, 1>(0, 1); } - /** This method retrieves z axis */ + /// This method retrieves z axis ALGEBRA_HOST_DEVICE inline point3 z() const { return _data.matrix().template block<3, 1>(0, 2); } - /** This method retrieves the translation of a transform **/ + /// This method retrieves the translation of a transform ALGEBRA_HOST_DEVICE inline auto translation() const { return _data.matrix().template block<3, 1>(0, 3); } - /** This method retrieves the 4x4 matrix of a transform */ + /// This method retrieves the 4x4 matrix of a transform ALGEBRA_HOST_DEVICE inline const matrix44 &matrix() const { return _data.matrix(); } - /** This method retrieves the 4x4 matrix of an inverse transform */ + /// This method retrieves the 4x4 matrix of an inverse transform ALGEBRA_HOST_DEVICE inline const matrix44 &matrix_inverse() const { return _data_inv.matrix(); } - /** This method transform from a point from the local 3D cartesian frame to - * the global 3D cartesian frame */ + /// This method transform from a point from the local 3D cartesian frame to + /// the global 3D cartesian frame template < typename derived_type, std::enable_if_t::RowsAtCompileTime == 3, @@ -233,8 +220,8 @@ struct transform3 { return (_data * v); } - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame template < typename derived_type, std::enable_if_t::RowsAtCompileTime == 3, @@ -247,8 +234,8 @@ struct transform3 { return (_data_inv * v); } - /** This method transform from a vector from the local 3D cartesian frame to - * the global 3D cartesian frame */ + /// This method transform from a vector from the local 3D cartesian frame to + /// the global 3D cartesian frame template < typename derived_type, std::enable_if_t::RowsAtCompileTime == 3, @@ -261,8 +248,8 @@ struct transform3 { return (_data.linear() * v); } - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame template < typename derived_type, std::enable_if_t::RowsAtCompileTime == 3, diff --git a/math/eigen/include/algebra/math/impl/eigen_vector.hpp b/math/eigen/include/algebra/math/impl/eigen_vector.hpp index 003c004d..3486cf81 100644 --- a/math/eigen/include/algebra/math/impl/eigen_vector.hpp +++ b/math/eigen/include/algebra/math/impl/eigen_vector.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (c) 2020-2023 CERN for the benefit of the ACTS project + * (c) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -28,28 +28,47 @@ namespace algebra::eigen::math { -/** Get a normalized version of the input vector - * - * @tparam derived_type is the matrix template - * - * @param v the input vector - **/ +/// This method retrieves the norm of a vector, no dimension restriction +/// +/// @param v the input vector +template +ALGEBRA_HOST_DEVICE inline auto norm(const Eigen::MatrixBase &v) { + + return v.norm(); +} + +/// This method retrieves the pseudo-rapidity from a vector or vector base with +/// rows >= 3 +/// +/// @param v the input vector +template +requires(Eigen::MatrixBase::RowsAtCompileTime >= + 3) ALGEBRA_HOST_DEVICE + inline auto eta(const Eigen::MatrixBase &v) noexcept { + + return algebra::math::atanh(v[2] / v.norm()); +} + +/// Get a normalized version of the input vector +/// +/// @tparam derived_type is the matrix template +/// +/// @param v the input vector template ALGEBRA_HOST_DEVICE inline auto normalize( const Eigen::MatrixBase &v) { return v.normalized(); } -/** Dot product between two input vectors - * - * @tparam derived_type_lhs is the first matrix (expression) template - * @tparam derived_type_rhs is the second matrix (expression) template - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @tparam derived_type_lhs is the first matrix (expression) template +/// @tparam derived_type_rhs is the second matrix (expression) template +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST_DEVICE inline auto dot( const Eigen::MatrixBase &a, @@ -57,16 +76,15 @@ ALGEBRA_HOST_DEVICE inline auto dot( return a.dot(b); } -/** Cross product between two input vectors - * - * @tparam derived_type_lhs is the first matrix (expression) template - * @tparam derived_type_rhs is the second matrix (expression) template - * - * @param a the first input vector - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between two input vectors +/// +/// @tparam derived_type_lhs is the first matrix (expression) template +/// @tparam derived_type_rhs is the second matrix (expression) template +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST_DEVICE inline auto cross( const Eigen::MatrixBase &a, diff --git a/math/fastor/CMakeLists.txt b/math/fastor/CMakeLists.txt index 643f0ef3..5852ab9b 100644 --- a/math/fastor/CMakeLists.txt +++ b/math/fastor/CMakeLists.txt @@ -1,13 +1,12 @@ # Algebra plugins library, part of the ACTS project (R&D line) # -# (c) 2022-2023 CERN for the benefit of the ACTS project +# (c) 2022-2024 CERN for the benefit of the ACTS project # # Mozilla Public License Version 2.0 # Set up the library. algebra_add_library(algebra_fastor_math fastor_math "include/algebra/math/fastor.hpp" - "include/algebra/math/impl/fastor_getter.hpp" "include/algebra/math/impl/fastor_matrix.hpp" "include/algebra/math/impl/fastor_transform3.hpp" "include/algebra/math/impl/fastor_vector.hpp") diff --git a/math/fastor/include/algebra/math/fastor.hpp b/math/fastor/include/algebra/math/fastor.hpp index ad454a68..0e3e7d5d 100644 --- a/math/fastor/include/algebra/math/fastor.hpp +++ b/math/fastor/include/algebra/math/fastor.hpp @@ -8,7 +8,6 @@ #pragma once // Project include(s). -#include "algebra/math/impl/fastor_getter.hpp" #include "algebra/math/impl/fastor_matrix.hpp" #include "algebra/math/impl/fastor_transform3.hpp" #include "algebra/math/impl/fastor_vector.hpp" diff --git a/math/fastor/include/algebra/math/impl/fastor_getter.hpp b/math/fastor/include/algebra/math/impl/fastor_getter.hpp deleted file mode 100644 index 38d42b5a..00000000 --- a/math/fastor/include/algebra/math/impl/fastor_getter.hpp +++ /dev/null @@ -1,144 +0,0 @@ -/** Algebra plugins, part of the ACTS project - * - * (c) 2022 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" - -// Fastor include(s). -#ifdef _MSC_VER -#pragma warning(disable : 4244 4701 4702) -#endif // MSVC -#include -#ifdef _MSC_VER -#pragma warning(default : 4244 4701 4702) -#endif // MSVC - -// System include(s). -#include // for the std::size_t type -#include - -namespace algebra::fastor::math { - -/** This method retrieves phi from a vector, vector base with rows >= 2 - * - * @param v the input vector - **/ -template = 2, bool> = true> -ALGEBRA_HOST inline scalar_t phi( - const Fastor::Tensor &v) noexcept { - - return algebra::math::atan2(v[1], v[0]); -} - -/** This method retrieves theta from a vector, vector base with rows >= 3 - * - * @param v the input vector - **/ -template = 3, bool> = true> -ALGEBRA_HOST inline scalar_t theta( - const Fastor::Tensor &v) noexcept { - - return algebra::math::atan2(Fastor::norm(v(Fastor::fseq<0, 2>())), v[2]); -} - -/** This method retrieves the perpenticular magnitude of a vector with rows >= 2 - * - * @param v the input vector - **/ -template = 2, bool> = true> -ALGEBRA_HOST inline scalar_t perp( - const Fastor::Tensor &v) noexcept { - - return algebra::math::sqrt( - Fastor::inner(v(Fastor::fseq<0, 2>()), v(Fastor::fseq<0, 2>()))); -} - -/** This method retrieves the norm of a vector, no dimension restriction - * - * @param v the input vector - **/ -template -ALGEBRA_HOST inline scalar_t norm(const Fastor::Tensor &v) { - - return Fastor::norm(v); -} - -/** This method retrieves the pseudo-rapidity from a vector or vector base with - *rows >= 3 - * - * @param v the input vector - **/ -template = 3, bool> = true> -ALGEBRA_HOST inline scalar_t eta( - const Fastor::Tensor &v) noexcept { - - return algebra::math::atanh(v[2] / Fastor::norm(v)); -} - -/// Functor used to access elements of Fastor matrices -template -struct element_getter { - - template - using matrix_type = Fastor::Tensor; - - template - ALGEBRA_HOST_DEVICE inline scalar_t &operator()(matrix_type &m, - std::size_t row, - std::size_t col) const { - assert(row < ROWS); - assert(col < COLS); - return m(row, col); - } - - template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type &m, std::size_t row, - std::size_t col) const { - - assert(row < ROWS); - assert(col < COLS); - return m(row, col); - } -}; // element_getter - -/// Function extracting an element from a matrix (const) -template -ALGEBRA_HOST_DEVICE inline scalar_t element( - const Fastor::Tensor &m, std::size_t row, - std::size_t col) { - - return element_getter()(m, row, col); -} - -/// Function extracting an element from a matrix (non-const) -template -ALGEBRA_HOST_DEVICE inline scalar_t &element( - Fastor::Tensor &m, std::size_t row, std::size_t col) { - - return element_getter()(m, row, col); -} - -/// Functor used to extract a block from Fastor matrices -template -struct block_getter { - - template - using matrix_type = Fastor::Tensor; - - template - ALGEBRA_HOST_DEVICE matrix_type operator()( - const input_matrix_type &m, std::size_t row, std::size_t col) const { - - return m(Fastor::seq(row, row + ROWS), Fastor::seq(col, col + COLS)); - } -}; // struct block_getter - -} // namespace algebra::fastor::math diff --git a/math/fastor/include/algebra/math/impl/fastor_matrix.hpp b/math/fastor/include/algebra/math/impl/fastor_matrix.hpp index b3d421af..bfeff7fb 100644 --- a/math/fastor/include/algebra/math/impl/fastor_matrix.hpp +++ b/math/fastor/include/algebra/math/impl/fastor_matrix.hpp @@ -20,55 +20,15 @@ #pragma warning(default : 4244 4701 4702) #endif // MSVC -// System include(s). -#include // for the std::size_t type - namespace algebra::fastor::math { -/// Operator getting a block of a const matrix -template -ALGEBRA_HOST_DEVICE - matrix_type, ROWS, COLS> - block(const input_matrix_type &m, std::size_t row, std::size_t col) { - // In `Fastor::seq`, the last element is not included. - // `Fastor::seq` takes `int`s as input, but `row`, `col`, `ROWS`, and `COLS` - // have type `std::size_t`. - return m(Fastor::seq(static_cast(row), static_cast(row + ROWS)), - Fastor::seq(static_cast(col), static_cast(col + COLS))); -} - -/// Operator setting a block with a matrix -template -ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const matrix_type &b, - std::size_t row, std::size_t col) { - // In `Fastor::seq`, the last element is not included. - // `Fastor::seq` takes `int`s as input, but `ROWS` and `COLS` have type - // `std::size_t`, which is `std::size_t`. - m(Fastor::seq(static_cast(row), static_cast(row + ROWS)), - Fastor::seq(static_cast(col), static_cast(col + COLS))) = b; -} - -/// Operator setting a block with a vector -template -ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const Fastor::Tensor &b, - std::size_t row, std::size_t col) { - // In `Fastor::seq`, the last element is not included. - // `Fastor::seq` takes `int`s as input, but `ROWS` and `COLS` have type - // `std::size_t`, which is `std::size_t`. - m(Fastor::seq(static_cast(row), static_cast(row + ROWS)), - static_cast(col)) = b; -} - -// Create zero matrix +/// Create zero matrix template ALGEBRA_HOST_DEVICE inline matrix_t zero() { return matrix_t(0); } -// Create identity matrix +/// Create identity matrix template ALGEBRA_HOST_DEVICE inline matrix_t identity() { using scalar_t = algebra::trait::value_t; @@ -88,13 +48,13 @@ ALGEBRA_HOST_DEVICE inline matrix_t identity() { } } -// Set input matrix as zero matrix +/// Set input matrix as zero matrix template ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type &m) { m.zeros(); } -// Set input matrix as identity matrix +/// Set input matrix as identity matrix template ALGEBRA_HOST_DEVICE inline void set_identity( matrix_type &m) { @@ -102,7 +62,7 @@ ALGEBRA_HOST_DEVICE inline void set_identity( m = identity>(); } -// Create transpose matrix +/// Create transpose matrix template ALGEBRA_HOST_DEVICE inline matrix_type transpose( const matrix_type &m) { diff --git a/math/fastor/include/algebra/math/impl/fastor_transform3.hpp b/math/fastor/include/algebra/math/impl/fastor_transform3.hpp index 459a1480..b5fedecf 100644 --- a/math/fastor/include/algebra/math/impl/fastor_transform3.hpp +++ b/math/fastor/include/algebra/math/impl/fastor_transform3.hpp @@ -1,6 +1,6 @@ /** Algebra plugins, part of the ACTS project * - * (c) 2022 CERN for the benefit of the ACTS project + * (c) 2022-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -24,19 +24,18 @@ namespace algebra::fastor::math { -/** Transform wrapper class to ensure standard API within different plugins - * - **/ +/// Transform wrapper class to ensure standard API towards the Fastor library template struct transform3 { /// @name Type definitions for the struct /// @{ + /// Scalar type used by the transform + using scalar_type = scalar_t; + /// Array type used by the transform template using array_type = Fastor::Tensor; - /// Scalar type used by the transform - using scalar_type = scalar_t; /// 3-element "vector" type using vector3 = array_type<3>; @@ -48,16 +47,6 @@ struct transform3 { /// 4x4 matrix type using matrix44 = Fastor::Tensor; - /// Function (object) used for accessing a matrix element - using element_getter = algebra::fastor::math::element_getter; - - /// Size type - using size_type = std::size_t; - - /// 2D Matrix type - template - using matrix_type = Fastor::Tensor; - /// @} /// @name Data objects @@ -68,14 +57,12 @@ struct transform3 { /// @} - /** 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 - * - **/ + /// 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, bool get_inverse = true) { @@ -97,22 +84,19 @@ struct transform3 { } } - /** 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 - * - **/ + /// 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 ALGEBRA_HOST transform3(const vector3 &t, const vector3 &z, const vector3 &x, bool get_inverse = true) : transform3(t, x, Fastor::cross(z, x), z, get_inverse) {} - /** Constructor with arguments: translation - * - * @param t is the translation - **/ + /// Constructor with arguments: translation + /// + /// @param t is the translation ALGEBRA_HOST explicit transform3(const vector3 &t) { @@ -127,10 +111,9 @@ struct transform3 { _data_inv = Fastor::inverse(_data); } - /** Constructor with arguments: matrix - * - * @param m is the full 4x4 matrix - **/ + /// Constructor with arguments: matrix + /// + /// @param m is the full 4x4 matrix ALGEBRA_HOST explicit transform3(const matrix44 &m) { _data = m; @@ -138,11 +121,10 @@ struct transform3 { _data_inv = Fastor::inverse(_data); } - /** Constructor with arguments: matrix as Fastor::Tensor of - * scalars - * - * @param ma is the full 4x4 matrix as a 16-element array - **/ + /// Constructor with arguments: matrix as Fastor::Tensor of + /// scalars + /// + /// @param ma is the full 4x4 matrix as a 16-element array ALGEBRA_HOST explicit transform3(const array_type<16> &ma) { _data = ma; @@ -150,19 +132,19 @@ struct transform3 { _data_inv = Fastor::inverse(_data); } - /** Default contructors */ + /// Default contructors transform3() = default; transform3(const transform3 &rhs) = default; ~transform3() = default; - /** Equality operator */ + /// Equality operator ALGEBRA_HOST inline bool operator==(const transform3 &rhs) const { return Fastor::isequal(_data, rhs._data); } - /** This method retrieves the rotation of a transform */ + /// This method retrieves the rotation of a transform ALGEBRA_HOST inline auto rotation() const { @@ -170,32 +152,32 @@ struct transform3 { _data(Fastor::fseq<0, 3>(), Fastor::fseq<0, 3>())); } - /** This method retrieves x axis */ + /// This method retrieves x axis ALGEBRA_HOST_DEVICE inline point3 x() const { return _data(Fastor::fseq<0, 3>(), 0); } - /** This method retrieves y axis */ + /// This method retrieves y axis ALGEBRA_HOST_DEVICE inline point3 y() const { return _data(Fastor::fseq<0, 3>(), 1); } - /** This method retrieves z axis */ + /// This method retrieves z axis ALGEBRA_HOST_DEVICE inline point3 z() const { return _data(Fastor::fseq<0, 3>(), 2); } - /** This method retrieves the translation of a transform */ + /// This method retrieves the translation of a transform ALGEBRA_HOST inline vector3 translation() const { return _data(Fastor::fseq<0, 3>(), 3); } - /** This method retrieves the 4x4 matrix of a transform */ + /// This method retrieves the 4x4 matrix of a transform ALGEBRA_HOST inline matrix44 matrix() const { return _data; } - /** This method retrieves the 4x4 matrix of an inverse transform */ + /// This method retrieves the 4x4 matrix of an inverse transform ALGEBRA_HOST inline matrix44 matrix_inverse() const { return _data_inv; } - /** This method transform from a point from the local 3D cartesian frame to - * the global 3D cartesian frame */ + /// This method transform from a point from the local 3D cartesian frame to + /// the global 3D cartesian frame ALGEBRA_HOST inline const point3 point_to_global(const point3 &v) const { @@ -206,8 +188,8 @@ struct transform3 { Fastor::matmul(_data, vector_4)(Fastor::fseq<0, 3>())); } - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame ALGEBRA_HOST inline const point3 point_to_local(const point3 &v) const { @@ -218,8 +200,8 @@ struct transform3 { Fastor::matmul(_data_inv, vector_4)(Fastor::fseq<0, 3>())); } - /** This method transform from a vector from the local 3D cartesian frame to - * the global 3D cartesian frame */ + /// This method transform from a vector from the local 3D cartesian frame to + /// the global 3D cartesian frame ALGEBRA_HOST inline const point3 vector_to_global(const vector3 &v) const { @@ -230,8 +212,8 @@ struct transform3 { Fastor::matmul(_data, vector_4)(Fastor::fseq<0, 3>())); } - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame ALGEBRA_HOST inline const point3 vector_to_local(const vector3 &v) const { diff --git a/math/fastor/include/algebra/math/impl/fastor_vector.hpp b/math/fastor/include/algebra/math/impl/fastor_vector.hpp index 15e9a5d5..08c2d603 100644 --- a/math/fastor/include/algebra/math/impl/fastor_vector.hpp +++ b/math/fastor/include/algebra/math/impl/fastor_vector.hpp @@ -1,6 +1,6 @@ /** Algebra plugins, part of the ACTS project * - * (c) 2022 CERN for the benefit of the ACTS project + * (c) 2022-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -21,10 +21,50 @@ namespace algebra::fastor::math { -/** Get a normalized version of the input vector - * - * @param v the input vector - **/ +/// This method retrieves theta from a vector, vector base with rows >= 3 +/// +/// @param v the input vector +template +requires(N >= 3) ALGEBRA_HOST inline scalar_t + theta(const Fastor::Tensor &v) noexcept { + + return algebra::math::atan2(Fastor::norm(v(Fastor::fseq<0, 2>())), v[2]); +} + +/// This method retrieves the perpenticular magnitude of a vector with rows >= 2 +/// +/// @param v the input vector +template +requires(N >= 2) ALGEBRA_HOST inline scalar_t + perp(const Fastor::Tensor &v) noexcept { + + return algebra::math::sqrt( + Fastor::inner(v(Fastor::fseq<0, 2>()), v(Fastor::fseq<0, 2>()))); +} + +/// This method retrieves the norm of a vector, no dimension restriction +/// +/// @param v the input vector +template +ALGEBRA_HOST inline scalar_t norm(const Fastor::Tensor &v) { + + return Fastor::norm(v); +} + +/// This method retrieves the pseudo-rapidity from a vector or vector base with +/// rows >= 3 +/// +/// @param v the input vector +template +requires(N >= 3) ALGEBRA_HOST inline scalar_t + eta(const Fastor::Tensor &v) noexcept { + + return algebra::math::atanh(v[2] / Fastor::norm(v)); +} + +/// Get a normalized version of the input vector +/// +/// @param v the input vector template ALGEBRA_HOST inline Fastor::Tensor normalize( const Fastor::Tensor &v) { @@ -32,26 +72,24 @@ ALGEBRA_HOST inline Fastor::Tensor normalize( return (static_cast(1.0) / Fastor::norm(v)) * v; } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST_DEVICE inline scalar_t dot(const Fastor::Tensor &a, const Fastor::Tensor &b) { return Fastor::inner(a, b); } -/** Dot product between Tensor and Tensor - * - * @param a the first input vector - * @param b the second input Tensor - * - * @return the scalar dot product value - **/ +/// Dot product between Tensor and Tensor +/// +/// @param a the first input vector +/// @param b the second input Tensor +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor &a, const Fastor::Tensor &b) { @@ -63,13 +101,12 @@ ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor &a, Fastor::Tensor(b(Fastor::fseq<0, N>(), 0))); } -/** Dot product between Tensor and Tensor - * - * @param a the second input Tensor - * @param b the first input vector - * - * @return the scalar dot product value - **/ +/// Dot product between Tensor and Tensor +/// +/// @param a the second input Tensor +/// @param b the first input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor &a, const Fastor::Tensor &b) { @@ -78,13 +115,12 @@ ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor &a, b); } -/** Dot product between two Tensor - * - * @param a the second input Tensor - * @param b the first input Tensor - * - * @return the scalar dot product value - **/ +/// Dot product between two Tensor +/// +/// @param a the second input Tensor +/// @param b the first input Tensor +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor &a, const Fastor::Tensor &b) { @@ -93,13 +129,12 @@ ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor &a, Fastor::Tensor(b(Fastor::fseq<0, 3>(), 0))); } -/** Cross product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST_DEVICE inline Fastor::Tensor cross( const Fastor::Tensor &a, @@ -107,13 +142,12 @@ ALGEBRA_HOST_DEVICE inline Fastor::Tensor cross( return Fastor::cross(a, b); } -/** Cross product between Tensor and Tensor - * - * @param a the first input vector - * @param b the second input Tensor - * - * @return a vector representing the cross product - **/ +/// Cross product between Tensor and Tensor +/// +/// @param a the first input vector +/// @param b the second input Tensor +/// +/// @return a vector representing the cross product template ALGEBRA_HOST inline Fastor::Tensor cross( const Fastor::Tensor &a, @@ -126,13 +160,12 @@ ALGEBRA_HOST inline Fastor::Tensor cross( Fastor::Tensor(b(Fastor::fseq<0, 3>(), 0))); } -/** Cross product between Tensor and Tensor - * - * @param a the second input Tensor - * @param b the first input vector - * - * @return a vector representing the cross product - **/ +/// Cross product between Tensor and Tensor +/// +/// @param a the second input Tensor +/// @param b the first input vector +/// +/// @return a vector representing the cross product template ALGEBRA_HOST inline Fastor::Tensor cross( const Fastor::Tensor &a, @@ -142,13 +175,12 @@ ALGEBRA_HOST inline Fastor::Tensor cross( b); } -/** Cross product between two Tensor - * - * @param a the second input Tensor - * @param b the first input Tensor - * - * @return a vector representing the cross product - **/ +/// Cross product between two Tensor +/// +/// @param a the second input Tensor +/// @param b the first input Tensor +/// +/// @return a vector representing the cross product template ALGEBRA_HOST inline Fastor::Tensor cross( const Fastor::Tensor &a, diff --git a/math/generic/CMakeLists.txt b/math/generic/CMakeLists.txt new file mode 100644 index 00000000..098949a8 --- /dev/null +++ b/math/generic/CMakeLists.txt @@ -0,0 +1,26 @@ +# Algebra plugins library, part of the ACTS project (R&D line) +# +# (c) 2021-2024 CERN for the benefit of the ACTS project +# +# Mozilla Public License Version 2.0 + +# Set up the library. +algebra_add_library(algebra_generic_math generic_math + # impl include + "include/algebra/math/generic.hpp" + "include/algebra/math/impl/generic_matrix.hpp" + "include/algebra/math/impl/generic_transform3.hpp" + "include/algebra/math/impl/generic_vector.hpp" + # algorithms include + "include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp" + "include/algebra/math/algorithms/matrix/determinant/cofactor.hpp" + "include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp" + "include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp" + "include/algebra/math/algorithms/matrix/inverse/cofactor.hpp" + "include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp" + "include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp" + "include/algebra/math/algorithms/utils/algorithm_finder.hpp") +target_link_libraries(algebra_generic_math + INTERFACE algebra::common algebra::common_math) +algebra_test_public_headers( algebra_generic_math + "algebra/math/generic.hpp" ) diff --git a/math/common/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp b/math/generic/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp similarity index 94% rename from math/common/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp rename to math/generic/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp index 3c91a7d5..9d929c42 100644 --- a/math/common/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp +++ b/math/generic/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp @@ -14,7 +14,7 @@ // System include(s). #include -namespace algebra::cmath::matrix::decomposition { +namespace algebra::generic::matrix::decomposition { /// "Partial Pivot LU Decomposition", assuming a N X N matrix template @@ -73,7 +73,7 @@ struct partial_pivot_lud { max_idx = i; for (size_type k = i; k < N; k++) { - abs_val = math::fabs(element_getter()(lu, k, i)); + abs_val = algebra::math::fabs(element_getter()(lu, k, i)); if (abs_val > max_val) { @@ -119,4 +119,4 @@ struct partial_pivot_lud { } }; -} // namespace algebra::cmath::matrix::decomposition +} // namespace algebra::generic::matrix::decomposition diff --git a/math/common/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp b/math/generic/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp similarity index 96% rename from math/common/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp rename to math/generic/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp index f0b83788..461c5338 100644 --- a/math/common/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp +++ b/math/generic/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp @@ -14,7 +14,7 @@ // System include(s) #include -namespace algebra::cmath::matrix::determinant { +namespace algebra::generic::matrix::determinant { /// "Determinant getter", assuming a N X N matrix template @@ -100,4 +100,4 @@ struct cofactor { }; }; -} // namespace algebra::cmath::matrix::determinant +} // namespace algebra::generic::matrix::determinant diff --git a/math/common/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp b/math/generic/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp similarity index 97% rename from math/common/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp rename to math/generic/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp index ef754b86..14262203 100644 --- a/math/common/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp +++ b/math/generic/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp @@ -14,7 +14,7 @@ // System include(s) #include -namespace algebra::cmath::matrix::determinant { +namespace algebra::generic::matrix::determinant { /// "Determinant getter", assuming a N X N matrix template @@ -91,4 +91,4 @@ struct hard_coded { } }; -} // namespace algebra::cmath::matrix::determinant +} // namespace algebra::generic::matrix::determinant diff --git a/math/common/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp b/math/generic/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp similarity index 88% rename from math/common/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp rename to math/generic/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp index 472fca4d..a9a17f33 100644 --- a/math/common/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp +++ b/math/generic/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp @@ -12,7 +12,7 @@ #include "algebra/qualifiers.hpp" #include "algebra/type_traits.hpp" -namespace algebra::cmath::matrix::determinant { +namespace algebra::generic::matrix::determinant { /// "Partial Pivot LU Decomposition", assuming a N X N matrix template @@ -25,7 +25,7 @@ struct partial_pivot_lud { using element_getter = element_getter_t; using decomposition_t = - typename algebra::cmath::matrix::decomposition::partial_pivot_lud< + typename algebra::generic::matrix::decomposition::partial_pivot_lud< matrix_t, element_getter_t>; ALGEBRA_HOST_DEVICE inline scalar_type operator()(const matrix_t& m) const { @@ -49,4 +49,4 @@ struct partial_pivot_lud { } }; -} // namespace algebra::cmath::matrix::determinant +} // namespace algebra::generic::matrix::determinant diff --git a/math/common/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp b/math/generic/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp similarity index 97% rename from math/common/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp rename to math/generic/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp index a68d701b..40eeed50 100644 --- a/math/common/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp +++ b/math/generic/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp @@ -15,7 +15,7 @@ // System include(s) #include -namespace algebra::cmath::matrix { +namespace algebra::generic::matrix { namespace adjoint { @@ -133,4 +133,4 @@ struct cofactor { } // namespace inverse -} // namespace algebra::cmath::matrix +} // namespace algebra::generic::matrix diff --git a/math/common/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp b/math/generic/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp similarity index 99% rename from math/common/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp rename to math/generic/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp index 8736dd9c..ba817a1c 100644 --- a/math/common/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp +++ b/math/generic/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp @@ -15,7 +15,7 @@ // System include(s) #include -namespace algebra::cmath::matrix::inverse { +namespace algebra::generic::matrix::inverse { /// "inverse getter", assuming a N X N matrix template @@ -272,4 +272,4 @@ struct hard_coded { } }; -} // namespace algebra::cmath::matrix::inverse +} // namespace algebra::generic::matrix::inverse diff --git a/math/common/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp b/math/generic/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp similarity index 92% rename from math/common/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp rename to math/generic/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp index 6bd64ea9..aa29c4f1 100644 --- a/math/common/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp +++ b/math/generic/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp @@ -12,7 +12,7 @@ #include "algebra/qualifiers.hpp" #include "algebra/type_traits.hpp" -namespace algebra::cmath::matrix::inverse { +namespace algebra::generic::matrix::inverse { /// "Partial Pivot LU Decomposition", assuming a N X N matrix template @@ -25,7 +25,7 @@ struct partial_pivot_lud { using element_getter = element_getter_t; using decomposition_t = - typename algebra::cmath::matrix::decomposition::partial_pivot_lud< + typename algebra::generic::matrix::decomposition::partial_pivot_lud< matrix_t, element_getter_t>; ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t& m) const { @@ -70,4 +70,4 @@ struct partial_pivot_lud { } }; -} // namespace algebra::cmath::matrix::inverse +} // namespace algebra::generic::matrix::inverse diff --git a/math/common/include/algebra/math/algorithms/utils/algorithm_finder.hpp b/math/generic/include/algebra/math/algorithms/utils/algorithm_finder.hpp similarity index 76% rename from math/common/include/algebra/math/algorithms/utils/algorithm_finder.hpp rename to math/generic/include/algebra/math/algorithms/utils/algorithm_finder.hpp index 2a70c9a9..918b7d6c 100644 --- a/math/common/include/algebra/math/algorithms/utils/algorithm_finder.hpp +++ b/math/generic/include/algebra/math/algorithms/utils/algorithm_finder.hpp @@ -13,7 +13,7 @@ #include "algebra/math/algorithms/matrix/inverse/hard_coded.hpp" #include "algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp" -namespace algebra::cmath { +namespace algebra::generic { /// Get the type of determinant algorithm acording to matrix dimension /// @{ @@ -32,9 +32,11 @@ struct determinant_selector<4, Args...> { using type = matrix::determinant::hard_coded; }; -template +/// @tparam M matrix type +template using determinant_t = - typename determinant_selector::type; + typename determinant_selector, M, + algebra::trait::element_getter_t>::type; /// @} /// Get the type of inversion algorithm acording to matrix dimension @@ -54,9 +56,11 @@ struct inversion_selector<4, Args...> { using type = matrix::inverse::hard_coded; }; -template +/// @tparam M matrix type +template using inversion_t = - typename inversion_selector::type; -/// @} + typename inversion_selector, M, + algebra::trait::element_getter_t>::type; +/// @}generic -} // namespace algebra::cmath +} // namespace algebra::generic diff --git a/math/generic/include/algebra/math/generic.hpp b/math/generic/include/algebra/math/generic.hpp new file mode 100644 index 00000000..e3d7596c --- /dev/null +++ b/math/generic/include/algebra/math/generic.hpp @@ -0,0 +1,23 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2020-2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Impl include(s). +#include "algebra/math/impl/generic_matrix.hpp" +#include "algebra/math/impl/generic_transform3.hpp" +#include "algebra/math/impl/generic_vector.hpp" + +// Algorithms include(s). +#include "algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp" +#include "algebra/math/algorithms/matrix/determinant/cofactor.hpp" +#include "algebra/math/algorithms/matrix/determinant/hard_coded.hpp" +#include "algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp" +#include "algebra/math/algorithms/matrix/inverse/cofactor.hpp" +#include "algebra/math/algorithms/matrix/inverse/hard_coded.hpp" +#include "algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp" +#include "algebra/math/algorithms/utils/algorithm_finder.hpp" diff --git a/math/generic/include/algebra/math/impl/generic_matrix.hpp b/math/generic/include/algebra/math/impl/generic_matrix.hpp new file mode 100644 index 00000000..41c49864 --- /dev/null +++ b/math/generic/include/algebra/math/impl/generic_matrix.hpp @@ -0,0 +1,100 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2022-2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "algebra/math/algorithms/utils/algorithm_finder.hpp" +#include "algebra/math/common.hpp" +#include "algebra/qualifiers.hpp" + +namespace algebra::generic::math { + +/// Create zero matrix - generic transform3 +template +ALGEBRA_HOST_DEVICE inline matrix_t zero() { + + using index_t = algebra::trait::index_t; + using element_getter_t = algebra::trait::element_getter_t; + + matrix_t ret; + + for (index_t j = 0; j < algebra::trait::columns; ++j) { + for (index_t i = 0; i < algebra::trait::rows; ++i) { + element_getter_t{}(ret, i, j) = 0; + } + } + + return ret; +} + +/// Create identity matrix - generic transform3 +template +ALGEBRA_HOST_DEVICE inline matrix_t identity() { + + using index_t = algebra::trait::index_t; + using element_getter_t = algebra::trait::element_getter_t; + + auto ret{zero()}; + + for (index_t i = 0; i < algebra::trait::rank; ++i) { + element_getter_t{}(ret, i, i) = 1; + } + + return ret; +} + +/// Set @param m as zero matrix +template +ALGEBRA_HOST_DEVICE inline void set_zero(matrix_t &m) { + m = zero(); +} + +/// Set @param m as identity matrix +template +ALGEBRA_HOST_DEVICE inline void set_identity(matrix_t &m) { + m = identity(); +} + +/// @returns the transpose matrix of @param m +template +ALGEBRA_HOST_DEVICE inline auto transpose(const matrix_t &m) { + + using index_t = algebra::trait::index_t; + using value_t = algebra::trait::value_t; + using element_getter_t = algebra::trait::element_getter_t; + + constexpr index_t rows{algebra::trait::rows}; + constexpr index_t columns{algebra::trait::columns}; + + algebra::trait::get_matrix_t ret; + + for (index_t i = 0; i < rows; ++i) { + for (index_t j = 0; j < columns; ++j) { + element_getter_t{}(ret, j, i) = element_getter_t{}(m, i, j); + } + } + + return ret; +} + +/// @returns the determinant of @param m +template +ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t determinant( + const matrix_t &m) { + + return determinant_t{}(m); +} + +/// @returns the determinant of @param m +template +ALGEBRA_HOST_DEVICE inline matrix_t inverse(const matrix_t &m) { + + return inversion_t{}(m); +} + +} // namespace algebra::generic::math diff --git a/math/generic/include/algebra/math/impl/generic_transform3.hpp b/math/generic/include/algebra/math/impl/generic_transform3.hpp new file mode 100644 index 00000000..ee740e3b --- /dev/null +++ b/math/generic/include/algebra/math/impl/generic_transform3.hpp @@ -0,0 +1,299 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2020-2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "algebra/math/algorithms/matrix/inverse/hard_coded.hpp" +#include "algebra/math/impl/generic_matrix.hpp" +#include "algebra/math/impl/generic_vector.hpp" +#include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" + +namespace algebra::generic::math { + +/// Transform wrapper class to ensure standard API within differnt plugins +template class matrix_t, + template class array_t> +struct transform3 { + + /// @name Type definitions for the struct + /// @{ + + /// Scalar type + using scalar_type = scalar_t; + + /// Array type used by the transform + template + using array_type = array_t; + + // 4 x 4 Matrix (keep four rows for better alignment) + using matrix44 = matrix_t; + + /// 3-element "vector" type + using vector3 = algebra::trait::get_vector_t; + /// Point in 3D space + using point3 = vector3; + /// Point in 2D space + using point2 = algebra::trait::get_vector_t; + + /// Function (object) used for accessing a matrix element + using element_getter = algebra::trait::element_getter_t; + + /// Function (object) used for accessing a matrix element + using block_getter = algebra::trait::block_getter_t; + + /// Matrix inversion algorithm + using matrix_inversion = + generic::matrix::inverse::hard_coded; + + /// @} + + /// @name Data objects + /// @{ + matrix44 _data{generic::math::identity()}; + matrix44 _data_inv{generic::math::identity()}; + + /// @} + + /// Default constructor: identity + constexpr 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, bool get_inverse = true) { + + element_getter{}(_data, 0, 0) = element_getter{}(x, 0); + element_getter{}(_data, 1, 0) = element_getter{}(x, 1); + element_getter{}(_data, 2, 0) = element_getter{}(x, 2); + element_getter{}(_data, 3, 0) = 0.f; + element_getter{}(_data, 0, 1) = element_getter{}(y, 0); + element_getter{}(_data, 1, 1) = element_getter{}(y, 1); + element_getter{}(_data, 2, 1) = element_getter{}(y, 2); + element_getter{}(_data, 3, 1) = 0.f; + element_getter{}(_data, 0, 2) = element_getter{}(z, 0); + element_getter{}(_data, 1, 2) = element_getter{}(z, 1); + element_getter{}(_data, 2, 2) = element_getter{}(z, 2); + element_getter{}(_data, 3, 2) = 0.f; + element_getter{}(_data, 0, 3) = element_getter{}(t, 0); + element_getter{}(_data, 1, 3) = element_getter{}(t, 1); + element_getter{}(_data, 2, 3) = element_getter{}(t, 2); + element_getter{}(_data, 3, 3) = 1.f; + + if (get_inverse) { + _data_inv = matrix_inversion{}(_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 + explicit transform3(const vector3 &t) { + + element_getter{}(_data, 0, 0) = 1.f; + element_getter{}(_data, 1, 0) = 0.f; + element_getter{}(_data, 2, 0) = 0.f; + element_getter{}(_data, 3, 0) = 0.f; + element_getter{}(_data, 0, 1) = 0.f; + element_getter{}(_data, 1, 1) = 1.f; + element_getter{}(_data, 2, 1) = 0.f; + element_getter{}(_data, 3, 1) = 0.f; + element_getter{}(_data, 0, 2) = 0.f; + element_getter{}(_data, 1, 2) = 0.f; + element_getter{}(_data, 2, 2) = 1.f; + element_getter{}(_data, 3, 2) = 0.f; + element_getter{}(_data, 0, 3) = element_getter{}(t, 0); + element_getter{}(_data, 1, 3) = element_getter{}(t, 1); + element_getter{}(_data, 2, 3) = element_getter{}(t, 2); + element_getter{}(_data, 3, 3) = 1.f; + + _data_inv = matrix_inversion{}(_data); + } + + /// Constructor with arguments: matrix + /// + /// @param m is the full 4x4 matrix + ALGEBRA_HOST_DEVICE + explicit transform3(const matrix44 &m) : _data{m} { + _data_inv = matrix_inversion{}(_data); + } + + /// Constructor with arguments: matrix as array of scalar + /// + /// @param ma is the full 4x4 matrix 16 array + ALGEBRA_HOST_DEVICE + explicit transform3(const array_type<16> &ma) { + + element_getter{}(_data, 0, 0) = ma[0]; + element_getter{}(_data, 1, 0) = ma[4]; + element_getter{}(_data, 2, 0) = ma[8]; + element_getter{}(_data, 3, 0) = 0.f; + element_getter{}(_data, 0, 1) = ma[1]; + element_getter{}(_data, 1, 1) = ma[5]; + element_getter{}(_data, 2, 1) = ma[9]; + element_getter{}(_data, 3, 1) = 0.f; + element_getter{}(_data, 0, 2) = ma[2]; + element_getter{}(_data, 1, 2) = ma[6]; + element_getter{}(_data, 2, 2) = ma[10]; + element_getter{}(_data, 3, 2) = 0.f; + element_getter{}(_data, 0, 3) = ma[3]; + element_getter{}(_data, 1, 3) = ma[7]; + element_getter{}(_data, 2, 3) = ma[11]; + element_getter{}(_data, 3, 3) = 1.f; + + _data_inv = matrix_inversion{}(_data); + } + + /// Equality operator + ALGEBRA_HOST_DEVICE + inline bool operator==(const transform3 &rhs) const { + + for (index_t j = 0; j < 4; j++) { + // Check only the rows that can differ + for (index_t i = 0; i < 3; i++) { + if (element_getter{}(_data, i, j) != + element_getter{}(rhs._data, i, j)) { + return false; + } + } + } + + return true; + } + + /// Rotate a vector into / from a frame + /// + /// @param m is the rotation matrix + /// @param v is the vector to be rotated + ALGEBRA_HOST_DEVICE + static inline vector3 rotate(const matrix44 &m, const vector3 &v) { + + vector3 ret{0.f, 0.f, 0.f}; + + element_getter{}(ret, 0) += + element_getter{}(m, 0, 0) * element_getter{}(v, 0); + element_getter{}(ret, 1) += + element_getter{}(m, 1, 0) * element_getter{}(v, 0); + element_getter{}(ret, 2) += + element_getter{}(m, 2, 0) * element_getter{}(v, 0); + + element_getter{}(ret, 0) += + element_getter{}(m, 0, 1) * element_getter{}(v, 1); + element_getter{}(ret, 1) += + element_getter{}(m, 1, 1) * element_getter{}(v, 1); + element_getter{}(ret, 2) += + element_getter{}(m, 2, 1) * element_getter{}(v, 1); + + element_getter{}(ret, 0) += + element_getter{}(m, 0, 2) * element_getter{}(v, 2); + element_getter{}(ret, 1) += + element_getter{}(m, 1, 2) * element_getter{}(v, 2); + element_getter{}(ret, 2) += + element_getter{}(m, 2, 2) * element_getter{}(v, 2); + + return ret; + } + + /// This method retrieves the rotation of a transform + ALGEBRA_HOST_DEVICE + auto inline rotation() const { + return block_getter{}.template operator()<3, 3>(_data, 0, 0); + } + + /// This method retrieves x axis + ALGEBRA_HOST_DEVICE + inline point3 x() const { + return {element_getter{}(_data, 0, 0), element_getter{}(_data, 1, 0), + element_getter{}(_data, 2, 0)}; + } + + /// This method retrieves y axis + ALGEBRA_HOST_DEVICE + inline point3 y() const { + return {element_getter{}(_data, 0, 1), element_getter{}(_data, 1, 1), + element_getter{}(_data, 2, 1)}; + } + + /// This method retrieves z axis + ALGEBRA_HOST_DEVICE + inline point3 z() const { + return {element_getter{}(_data, 0, 2), element_getter{}(_data, 1, 2), + element_getter{}(_data, 2, 2)}; + } + + /// This method retrieves the translation of a transform + ALGEBRA_HOST_DEVICE + inline point3 translation() const { + return {element_getter{}(_data, 0, 3), element_getter{}(_data, 1, 3), + element_getter{}(_data, 2, 3)}; + } + + /// This method retrieves the 4x4 matrix of a transform + ALGEBRA_HOST_DEVICE + inline const matrix44 &matrix() const { return _data; } + + /// This method retrieves the 4x4 matrix of an inverse transform + ALGEBRA_HOST_DEVICE + inline const matrix44 &matrix_inverse() const { return _data_inv; } + + /// This method transform from a point from the local 3D cartesian frame to + /// the global 3D cartesian frame + ALGEBRA_HOST_DEVICE inline point3 point_to_global(const point3 &v) const { + + const vector3 rg = rotate(_data, v); + + return {element_getter{}(rg, 0) + element_getter{}(_data, 0, 3), + element_getter{}(rg, 1) + element_getter{}(_data, 1, 3), + element_getter{}(rg, 2) + element_getter{}(_data, 2, 3)}; + } + + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame + ALGEBRA_HOST_DEVICE inline point3 point_to_local(const point3 &v) const { + + const vector3 rg = rotate(_data_inv, v); + + return {element_getter{}(rg, 0) + element_getter{}(_data_inv, 0, 3), + element_getter{}(rg, 1) + element_getter{}(_data_inv, 1, 3), + element_getter{}(rg, 2) + element_getter{}(_data_inv, 2, 3)}; + } + + /// This method transform from a vector from the local 3D cartesian frame to + /// the global 3D cartesian frame + ALGEBRA_HOST_DEVICE inline vector3 vector_to_global(const vector3 &v) const { + return rotate(_data, v); + } + + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame + ALGEBRA_HOST_DEVICE inline vector3 vector_to_local(const vector3 &v) const { + return rotate(_data_inv, v); + } + +}; // struct transform3 + +} // namespace algebra::generic::math diff --git a/math/generic/include/algebra/math/impl/generic_vector.hpp b/math/generic/include/algebra/math/impl/generic_vector.hpp new file mode 100644 index 00000000..51256f95 --- /dev/null +++ b/math/generic/include/algebra/math/impl/generic_vector.hpp @@ -0,0 +1,153 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2020-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/type_traits.hpp" + +namespace algebra::generic::math { + +/// This method retrieves phi from a vector with rows >= 2 +/// +/// @param v the input vector +template +requires algebra::trait::is_vector + ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t phi( + const vector_t &v) noexcept { + + using element_getter_t = algebra::trait::element_getter_t; + + return algebra::math::atan2(element_getter_t{}(v, 1), + element_getter_t{}(v, 0)); +} + +/// This method retrieves the perpendicular magnitude of a vector with rows >= 2 +/// +/// @param v the input vector +template +requires algebra::trait::is_vector + ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t perp( + const vector_t &v) noexcept { + + using element_getter_t = algebra::trait::element_getter_t; + + return algebra::math::sqrt( + algebra::math::fma(element_getter_t{}(v, 0), element_getter_t{}(v, 0), + element_getter_t{}(v, 1) * element_getter_t{}(v, 1))); +} + +/// This method retrieves theta from a vector with rows >= 3 +/// +/// @param v the input vector +template +requires algebra::trait::is_vector + ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t theta( + const vector_t &v) noexcept { + + using element_getter_t = algebra::trait::element_getter_t; + + return algebra::math::atan2(perp(v), element_getter_t{}(v, 2)); +} + +/// Cross product between two input vectors - 3 Dim +/// +/// @tparam vector1_t first vector or column matrix type +/// @tparam vector2_t second vector or column matrix type +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return a vector representing the cross product +template +requires(algebra::trait::columns == 1 && + algebra::trait::columns == 1) ALGEBRA_HOST_DEVICE + inline algebra::trait::vector_t cross(const vector1_t &a, + const vector2_t &b) { + + using element_getter_t = algebra::trait::element_getter_t; + + return { + algebra::math::fma(element_getter_t{}(a, 1), element_getter_t{}(b, 2), + -element_getter_t{}(b, 1) * element_getter_t{}(a, 2)), + algebra::math::fma(element_getter_t{}(a, 2), element_getter_t{}(b, 0), + -element_getter_t{}(b, 2) * element_getter_t{}(a, 0)), + algebra::math::fma(element_getter_t{}(a, 0), element_getter_t{}(b, 1), + -element_getter_t{}(b, 0) * element_getter_t{}(a, 1))}; +} + +/// Dot product between two input vectors +/// +/// @tparam vector1_t first vector or column matrix type +/// @tparam vector2_t second vector or column matrix type +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value +template +requires(algebra::trait::columns == 1 && + algebra::trait::columns == 1) ALGEBRA_HOST_DEVICE + inline algebra::trait::scalar_t dot(const vector1_t &a, + const vector2_t &b) { + + using scalar_t = algebra::trait::scalar_t; + using index_t = algebra::trait::index_t; + using element_getter_t = algebra::trait::element_getter_t; + + scalar_t ret{element_getter_t{}(a, 0) * element_getter_t{}(b, 0)}; + + for (index_t i = 1; i < algebra::trait::size; i++) { + ret = std::fma(element_getter_t{}(a, i), element_getter_t{}(b, i), ret); + } + + return ret; +} + +/// This method retrieves the norm of a vector with rows >= 3 +/// +/// @param v the input vector +template +requires algebra::trait::is_vector + ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t norm( + const vector_t &v) { + + return algebra::math::sqrt(dot(v, v)); +} + +/// This method retrieves the pseudo-rapidity from a vector or vector base with +/// rows >= 3 +/// +/// @param v the input vector +template +requires algebra::trait::is_vector + ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t eta( + const vector_t &v) noexcept { + + using element_getter_t = algebra::trait::element_getter_t; + + return algebra::math::atanh(element_getter_t{}(v, 2) / norm(v)); +} + +/// Get a normalized version of the input vector +/// +/// @tparam vector_t vector or column matrix type +/// +/// @param v the input vector +/// +/// @returns the normalized vector +template +ALGEBRA_HOST_DEVICE inline vector_t normalize(const vector_t &v) { + + using scalar_t = algebra::trait::scalar_t; + + return v * static_cast(1.) / norm(v); +} + +} // namespace algebra::generic::math diff --git a/math/smatrix/CMakeLists.txt b/math/smatrix/CMakeLists.txt index 0c47132d..d69f6527 100644 --- a/math/smatrix/CMakeLists.txt +++ b/math/smatrix/CMakeLists.txt @@ -1,6 +1,6 @@ # Algebra plugins library, part of the ACTS project (R&D line) # -# (c) 2021-2023 CERN for the benefit of the ACTS project +# (c) 2021-2024 CERN for the benefit of the ACTS project # # Mozilla Public License Version 2.0 @@ -11,7 +11,6 @@ find_package(ROOT COMPONENTS Core MathCore Smatrix REQUIRED) algebra_add_library(algebra_smatrix_math smatrix_math "include/algebra/math/smatrix.hpp" "include/algebra/math/impl/smatrix_errorcheck.hpp" - "include/algebra/math/impl/smatrix_getter.hpp" "include/algebra/math/impl/smatrix_matrix.hpp" "include/algebra/math/impl/smatrix_transform3.hpp" "include/algebra/math/impl/smatrix_vector.hpp") diff --git a/math/smatrix/include/algebra/math/impl/smatrix_getter.hpp b/math/smatrix/include/algebra/math/impl/smatrix_getter.hpp deleted file mode 100644 index 32d4c895..00000000 --- a/math/smatrix/include/algebra/math/impl/smatrix_getter.hpp +++ /dev/null @@ -1,184 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2020-2022 CERN for the benefit of the ACTS project - * - * Mozilla Public License Version 2.0 - */ - -#pragma once - -// Project include(s). -#include "algebra/qualifiers.hpp" - -// ROOT/Smatrix include(s). -#include -#include -#include -#include -#include - -// System include(s). -#include - -namespace algebra::smatrix::math { - -/** This method retrieves phi from a vector, vector base with rows >= 2 - * - * @param v the input vector - **/ -template = 2, bool> = true> -ALGEBRA_HOST inline scalar_t phi( - const ROOT::Math::SVector &v) noexcept { - - return static_cast(TMath::ATan2(v[1], v[0])); -} - -template = 2, bool> = true> -ALGEBRA_HOST inline scalar_t phi( - const ROOT::Math::VecExpr &v) noexcept { - - return static_cast(TMath::ATan2(v.apply(1), v.apply(0))); -} - -/** This method retrieves theta from a vector, vector base with rows >= 3 - * - * @param v the input vector - **/ -template = 3, bool> = true> -ALGEBRA_HOST inline scalar_t theta( - const ROOT::Math::SVector &v) noexcept { - - return static_cast( - TMath::ATan2(TMath::Sqrt(v[0] * v[0] + v[1] * v[1]), v[2])); -} - -template = 3, bool> = true> -ALGEBRA_HOST inline scalar_t theta( - const ROOT::Math::VecExpr &v) noexcept { - - return static_cast(TMath::ATan2( - TMath::Sqrt(v.apply(0) * v.apply(0) + v.apply(1) * v.apply(1)), - v.apply(2))); -} - -/** This method retrieves the norm of a vector, no dimension restriction - * - * @param v the input vector - **/ -template -ALGEBRA_HOST inline scalar_t norm(const ROOT::Math::SVector &v) { - - return static_cast(TMath::Sqrt(ROOT::Math::Dot(v, v))); -} - -template -ALGEBRA_HOST inline scalar_t norm( - const ROOT::Math::VecExpr &v) { - - return static_cast(TMath::Sqrt(ROOT::Math::Dot(v, v))); -} - -/** This method retrieves the pseudo-rapidity from a vector or vector base with - * rows >= 3 - * - * @param v the input vector - **/ -template = 3, bool> = true> -ALGEBRA_HOST inline scalar_t eta( - const ROOT::Math::SVector &v) noexcept { - - return static_cast(TMath::ATanH(v[2] / norm(v))); -} - -template = 3, bool> = true> -ALGEBRA_HOST inline scalar_t eta( - const ROOT::Math::VecExpr &v) noexcept { - - return static_cast(TMath::ATanH(v.apply(2) / norm(v))); -} - -/** This method retrieves the perpendicular magnitude of a vector with rows >= 2 - * - * @param v the input vector - **/ -template = 2, bool> = true> -ALGEBRA_HOST inline scalar_t perp( - const ROOT::Math::SVector &v) noexcept { - - return static_cast(TMath::Sqrt(v[0] * v[0] + v[1] * v[1])); -} - -template = 2, bool> = true> -ALGEBRA_HOST inline scalar_t perp( - const ROOT::Math::VecExpr &v) noexcept { - - return static_cast( - TMath::Sqrt(v.apply(0) * v.apply(0) + v.apply(1) * v.apply(1))); -} - -/// Functor used to access elements of SMatrix matrices -template -struct element_getter { - - template - using matrix_type = ROOT::Math::SMatrix; - - template - ALGEBRA_HOST_DEVICE inline scalar_t &operator()(matrix_type &m, - unsigned int row, - unsigned int col) const { - - assert(row < ROWS); - assert(col < COLS); - return m(row, col); - } - - template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type &m, unsigned int row, - unsigned int col) const { - - assert(row < ROWS); - assert(col < COLS); - return m(row, col); - } -}; // element_getter - -/// Function extracting an element from a matrix (const) -template -ALGEBRA_HOST_DEVICE inline scalar_t element( - const ROOT::Math::SMatrix &m, unsigned int row, - unsigned int col) { - - return element_getter()(m, row, col); -} - -/// Function extracting an element from a matrix (non-const) -template -ALGEBRA_HOST_DEVICE inline scalar_t &element( - ROOT::Math::SMatrix &m, unsigned int row, - unsigned int col) { - - return element_getter()(m, row, col); -} - -/// Functor used to extract a block from SMatrix matrices -template -struct block_getter { - - template - using matrix_type = ROOT::Math::SMatrix; - - template - ALGEBRA_HOST_DEVICE matrix_type operator()( - const input_matrix_type &m, unsigned int row, unsigned int col) const { - - return m.template Sub >(row, col); - } -}; // struct block_getter - -} // namespace algebra::smatrix::math diff --git a/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp b/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp index a5c9f0f8..20fc20bf 100644 --- a/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp +++ b/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp @@ -16,45 +16,13 @@ namespace algebra::smatrix::math { -/// Operator getting a block of a const matrix -template -ALGEBRA_HOST_DEVICE - matrix_type, ROWS, COLS> - block(const input_matrix_type &m, unsigned int row, unsigned int col) { - using scalar_t = algebra::trait::value_t; - return m.template Sub >(row, col); -} - -/// Operator setting a block with a matrix -template -ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const matrix_type &b, - unsigned int row, unsigned int col) { - for (unsigned int i = 0; i < ROWS; ++i) { - for (unsigned int j = 0; j < COLS; ++j) { - m(i + row, j + col) = b(i, j); - } - } -} - -/// Operator setting a block with a vector -template -ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const storage_type &b, - unsigned int row, unsigned int col) { - for (unsigned int i = 0; i < ROWS; ++i) { - m(i + row, col) = b[i]; - } -} - -// Create zero matrix +/// Create zero matrix template ALGEBRA_HOST_DEVICE inline matrix_t zero() { return matrix_t(); } -// Create identity matrix +/// Create identity matrix template ALGEBRA_HOST_DEVICE inline matrix_t identity() { return matrix_t(ROOT::Math::SMatrixIdentity()); diff --git a/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp b/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp index 9733aa2f..3bf1a7aa 100644 --- a/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp +++ b/math/smatrix/include/algebra/math/impl/smatrix_transform3.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) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -17,20 +17,19 @@ namespace algebra::smatrix::math { -/** Transform wrapper class to ensure standard API within differnt plugins - * - **/ +/// Transform wrapper class to ensure standard API towards the ROOT::SMatrix lib template struct transform3 { /// @name Type definitions for the struct /// @{ + /// Scalar type used by the transform + using scalar_type = scalar_t; + /// Array type used by the transform template using array_type = ROOT::Math::SVector; - /// Scalar type used by the transform - using scalar_type = scalar_t; /// 3-element "vector" type using vector3 = array_type<3>; @@ -42,13 +41,6 @@ struct transform3 { /// 4x4 matrix type using matrix44 = ROOT::Math::SMatrix; - /// Size type - using size_type = unsigned int; - - /// 2D Matrix type - template - using matrix_type = ROOT::Math::SMatrix; - /// @} /// @name Data objects @@ -59,14 +51,12 @@ struct transform3 { /// @} - /** 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 - * - **/ + /// 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, bool get_inverse = true) { @@ -90,22 +80,19 @@ struct transform3 { } } - /** 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 - * - **/ + /// 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 ALGEBRA_HOST transform3(const vector3 &t, const vector3 &z, const vector3 &x, bool get_inverse = true) : transform3(t, x, ROOT::Math::Cross(z, x), z, get_inverse) {} - /** Constructor with arguments: translation - * - * @param t is the translation - **/ + /// Constructor with arguments: translation + /// + /// @param t is the translation ALGEBRA_HOST explicit transform3(const vector3 &t) { @@ -118,10 +105,9 @@ struct transform3 { SMATRIX_CHECK(ifail); } - /** Constructor with arguments: matrix - * - * @param m is the full 4x4 matrix - **/ + /// Constructor with arguments: matrix + /// + /// @param m is the full 4x4 matrix ALGEBRA_HOST explicit transform3(const matrix44 &m) { _data = m; @@ -131,11 +117,10 @@ struct transform3 { SMATRIX_CHECK(ifail); } - /** Constructor with arguments: matrix as ROOT::Math::SVector of - * scalars - * - * @param ma is the full 4x4 matrix as a 16-element array - **/ + /// Constructor with arguments: matrix as ROOT::Math::SVector of + /// scalars + /// + /// @param ma is the full 4x4 matrix as a 16-element array ALGEBRA_HOST explicit transform3(const array_type<16> &ma) { @@ -162,54 +147,54 @@ struct transform3 { // from ROOT in this place... } - /** Default contructors */ + /// Default contructors transform3() = default; transform3(const transform3 &rhs) = default; ~transform3() = default; - /** Equality operator */ + /// Equality operator ALGEBRA_HOST inline bool operator==(const transform3 &rhs) const { return _data == rhs._data; } - /** This method retrieves the rotation of a transform */ + /// This method retrieves the rotation of a transform ALGEBRA_HOST inline auto rotation() const { return (_data.template Sub >(0, 0)); } - /** This method retrieves x axis */ + /// This method retrieves x axis ALGEBRA_HOST_DEVICE inline point3 x() const { return (_data.template SubCol(0, 0)); } - /** This method retrieves y axis */ + /// This method retrieves y axis ALGEBRA_HOST_DEVICE inline point3 y() const { return (_data.template SubCol(1, 0)); } - /** This method retrieves z axis */ + /// This method retrieves z axis ALGEBRA_HOST_DEVICE inline point3 z() const { return (_data.template SubCol(2, 0)); } - /** This method retrieves the translation of a transform */ + /// This method retrieves the translation of a transform ALGEBRA_HOST inline vector3 translation() const { return (_data.template SubCol(3, 0)); } - /** This method retrieves the 4x4 matrix of a transform */ + /// This method retrieves the 4x4 matrix of a transform ALGEBRA_HOST inline matrix44 matrix() const { return _data; } - /** This method retrieves the 4x4 matrix of an inverse transform */ + /// This method retrieves the 4x4 matrix of an inverse transform ALGEBRA_HOST inline matrix44 matrix_inverse() const { return _data_inv; } - /** This method transform from a point from the local 3D cartesian frame to - * the global 3D cartesian frame */ + /// This method transform from a point from the local 3D cartesian frame to + /// the global 3D cartesian frame ALGEBRA_HOST inline const point3 point_to_global(const point3 &v) const { @@ -220,8 +205,8 @@ struct transform3 { .template Sub(0); } - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame ALGEBRA_HOST inline const point3 point_to_local(const point3 &v) const { @@ -232,8 +217,8 @@ struct transform3 { .template Sub(0); } - /** This method transform from a vector from the local 3D cartesian frame to - * the global 3D cartesian frame */ + /// This method transform from a vector from the local 3D cartesian frame to + /// the global 3D cartesian frame ALGEBRA_HOST inline const point3 vector_to_global(const vector3 &v) const { @@ -243,8 +228,8 @@ struct transform3 { .template Sub(0); } - /** This method transform from a vector from the global 3D cartesian frame - * into the local 3D cartesian frame */ + /// This method transform from a vector from the global 3D cartesian frame + /// into the local 3D cartesian frame ALGEBRA_HOST inline const point3 vector_to_local(const vector3 &v) const { diff --git a/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp b/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp index ca7cd7a3..bb5fca1d 100644 --- a/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp +++ b/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (c) 2020-2021 CERN for the benefit of the ACTS project + * (c) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -18,20 +18,108 @@ namespace algebra::smatrix::math { -/** Get a normalized version of the input vector - * - * @param v the input vector - **/ +/// This method retrieves phi from a vector, vector base with rows >= 2 +/// +/// @param v the input vector +template +requires(N >= 2) ALGEBRA_HOST inline scalar_t + phi(const ROOT::Math::SVector &v) noexcept { + + return static_cast(TMath::ATan2(v[1], v[0])); +} + +template +requires(N >= 2) ALGEBRA_HOST inline scalar_t + phi(const ROOT::Math::VecExpr &v) noexcept { + + return static_cast(TMath::ATan2(v.apply(1), v.apply(0))); +} + +/// This method retrieves theta from a vector, vector base with rows >= 3 +/// +/// @param v the input vector +template +requires(N >= 3) ALGEBRA_HOST inline scalar_t + theta(const ROOT::Math::SVector &v) noexcept { + + return static_cast( + TMath::ATan2(TMath::Sqrt(v[0] * v[0] + v[1] * v[1]), v[2])); +} + +template +requires(N >= 3) ALGEBRA_HOST inline scalar_t + theta(const ROOT::Math::VecExpr &v) noexcept { + + return static_cast(TMath::ATan2( + TMath::Sqrt(v.apply(0) * v.apply(0) + v.apply(1) * v.apply(1)), + v.apply(2))); +} + +/// This method retrieves the norm of a vector, no dimension restriction +/// +/// @param v the input vector +template +ALGEBRA_HOST inline scalar_t norm(const ROOT::Math::SVector &v) { + + return static_cast(TMath::Sqrt(ROOT::Math::Dot(v, v))); +} + +template +ALGEBRA_HOST inline scalar_t norm( + const ROOT::Math::VecExpr &v) { + + return static_cast(TMath::Sqrt(ROOT::Math::Dot(v, v))); +} + +/// This method retrieves the pseudo-rapidity from a vector or vector base with +/// rows >= 3 +/// +/// @param v the input vector +template +requires(N >= 3) ALGEBRA_HOST inline scalar_t + eta(const ROOT::Math::SVector &v) noexcept { + + return static_cast(TMath::ATanH(v[2] / norm(v))); +} + +template +requires(N >= 3) ALGEBRA_HOST inline scalar_t + eta(const ROOT::Math::VecExpr &v) noexcept { + + return static_cast(TMath::ATanH(v.apply(2) / norm(v))); +} + +/// This method retrieves the perpendicular magnitude of a vector with rows >= 2 +/// +/// @param v the input vector +template +requires(N >= 2) ALGEBRA_HOST inline scalar_t + perp(const ROOT::Math::SVector &v) noexcept { + + return static_cast(TMath::Sqrt(v[0] * v[0] + v[1] * v[1])); +} + +template +requires(N >= 2) ALGEBRA_HOST inline scalar_t + perp(const ROOT::Math::VecExpr &v) noexcept { + + return static_cast( + TMath::Sqrt(v.apply(0) * v.apply(0) + v.apply(1) * v.apply(1))); +} + +/// Get a normalized version of the input vector +/// +/// @param v the input vector template ALGEBRA_HOST inline ROOT::Math::SVector normalize( const ROOT::Math::SVector &v) { return ROOT::Math::Unit(v); } -/** Get a normalized version of the input vector - * - * @param v the input vector - **/ + +/// Get a normalized version of the input vector +/// +/// @param v the input vector template ALGEBRA_HOST inline ROOT::Math::SVector normalize( const ROOT::Math::VecExpr &v) { @@ -39,39 +127,38 @@ ALGEBRA_HOST inline ROOT::Math::SVector normalize( return ROOT::Math::Unit(v); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector &a, const ROOT::Math::SVector &b) { return ROOT::Math::Dot(a, b); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ + +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector &a, const ROOT::Math::VecExpr &b) { return ROOT::Math::Dot(a, b); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ + +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr &a, const ROOT::Math::SVector &b) { @@ -79,13 +166,12 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr &a, return ROOT::Math::Dot(a, b); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr &a, const ROOT::Math::VecExpr &b) { @@ -93,13 +179,12 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr &a, return ROOT::Math::Dot(a, b); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix &a, const ROOT::Math::VecExpr &b) { @@ -107,26 +192,24 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix &a, return ROOT::Math::Dot(a.Col(0), b); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr &a, const ROOT::Math::SMatrix &b) { return dot(b, a); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix &a, const ROOT::Math::SVector &b) { @@ -134,26 +217,24 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix &a, return ROOT::Math::Dot(a.Col(0), b); } -/** Dot product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return the scalar dot product value - **/ +/// Dot product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector &a, const ROOT::Math::SMatrix &b) { return dot(b, a); } -/** Cross product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST inline ROOT::Math::SVector cross( const ROOT::Math::SVector &a, @@ -162,13 +243,12 @@ ALGEBRA_HOST inline ROOT::Math::SVector cross( return ROOT::Math::Cross(a, b); } -/** Cross product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST inline ROOT::Math::SVector cross( const ROOT::Math::SVector &a, @@ -177,13 +257,12 @@ ALGEBRA_HOST inline ROOT::Math::SVector cross( return ROOT::Math::Cross(a, b); } -/** Cross product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST inline ROOT::Math::SVector cross( const ROOT::Math::VecExpr &a, @@ -192,13 +271,12 @@ ALGEBRA_HOST inline ROOT::Math::SVector cross( return ROOT::Math::Cross(a, b); } -/** Cross product between two input vectors - * - * @param a the first input vector - * @param b the second input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between two input vectors +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST inline ROOT::Math::SVector cross( const ROOT::Math::VecExpr &a, @@ -207,13 +285,12 @@ ALGEBRA_HOST inline ROOT::Math::SVector cross( return ROOT::Math::Cross(a, b); } -/** Cross product between vector3 and matrix<3,1> - * - * @param a the first input vector - * @param b the second input matrix<3,1> - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between vector3 and matrix<3,1> +/// +/// @param a the first input vector +/// @param b the second input matrix<3,1> +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST inline ROOT::Math::SVector cross( const ROOT::Math::SVector &a, @@ -222,13 +299,12 @@ ALGEBRA_HOST inline ROOT::Math::SVector cross( return ROOT::Math::Cross(a, b.Col(0)); } -/** Cross product between matrix<3,1> and vector3 - * - * @param a the second input matrix<3,1> - * @param b the first input vector - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between matrix<3,1> and vector3 +/// +/// @param a the second input matrix<3,1> +/// @param b the first input vector +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST inline ROOT::Math::SVector cross( const ROOT::Math::SMatrix &a, @@ -237,13 +313,12 @@ ALGEBRA_HOST inline ROOT::Math::SVector cross( return ROOT::Math::Cross(a.Col(0), b); } -/** Cross product between two matrix<3,1> - * - * @param a the second input matrix<3,1> - * @param b the first input matrix<3,1> - * - * @return a vector (expression) representing the cross product - **/ +/// Cross product between two matrix<3,1> +/// +/// @param a the second input matrix<3,1> +/// @param b the first input matrix<3,1> +/// +/// @return a vector (expression) representing the cross product template ALGEBRA_HOST inline ROOT::Math::SVector cross( const ROOT::Math::SMatrix &a, diff --git a/math/smatrix/include/algebra/math/smatrix.hpp b/math/smatrix/include/algebra/math/smatrix.hpp index 418bb931..e55507f4 100644 --- a/math/smatrix/include/algebra/math/smatrix.hpp +++ b/math/smatrix/include/algebra/math/smatrix.hpp @@ -8,7 +8,6 @@ #pragma once // Project include(s). -#include "algebra/math/impl/smatrix_getter.hpp" #include "algebra/math/impl/smatrix_matrix.hpp" #include "algebra/math/impl/smatrix_transform3.hpp" #include "algebra/math/impl/smatrix_vector.hpp" diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_getter.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_getter.hpp deleted file mode 100644 index 15744041..00000000 --- a/math/vc_aos/include/algebra/math/impl/vc_aos_getter.hpp +++ /dev/null @@ -1,104 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2020-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/math/impl/vc_aos_vector.hpp" -#include "algebra/qualifiers.hpp" -#include "algebra/storage/matrix.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 - -namespace algebra::vc_aos::math { - -/// This method retrieves phi from a vector, vector base with rows >= 2 -/// -/// @param v the input vector -template -requires(Vc::is_simd_vector::value || - algebra::detail::is_storage_vector_v) ALGEBRA_HOST_DEVICE - inline auto phi(const vector_t &v) noexcept { - - return algebra::math::atan2(v[1], v[0]); -} - -/// This method retrieves the perpendicular magnitude of a vector with rows >= 2 -/// -/// @param v the input vector -template -requires(Vc::is_simd_vector::value || - algebra::detail::is_storage_vector_v) ALGEBRA_HOST_DEVICE - inline auto perp(const vector_t &v) noexcept { - - return algebra::math::sqrt(v[0] * v[0] + v[1] * v[1]); -} - -/// This method retrieves theta from a vector, vector base with rows >= 3 -/// -/// @param v the input vector -template -requires(Vc::is_simd_vector::value || - algebra::detail::is_storage_vector_v) ALGEBRA_HOST_DEVICE - inline auto theta(const vector_t &v) noexcept { - - return algebra::math::atan2(perp(v), v[2]); -} - -/// This method retrieves the norm of a vector, no dimension restriction -/// -/// @param v the input vector -template -requires(Vc::is_simd_vector::value || - algebra::detail::is_storage_vector_v) ALGEBRA_HOST_DEVICE - inline auto norm(const vector_t &v) { - - return algebra::math::sqrt(dot(v, v)); -} - -/// This method retrieves the pseudo-rapidity from a vector or vector base with -/// rows >= 3 -/// -/// @param v the input vector -template -requires(Vc::is_simd_vector::value || - algebra::detail::is_storage_vector_v) ALGEBRA_HOST_DEVICE - inline auto eta(const vector_t &v) noexcept { - - return algebra::math::atanh(v[2] / norm(v)); -} - -/// Function extracting an element from a matrix (const) -template -requires(Vc::is_simd_vector::value || - algebra::detail::is_storage_vector_v) - ALGEBRA_HOST_DEVICE inline decltype(auto) - element(const matrix_t &m, std::size_t row, std::size_t col) { - - return m[col][row]; -} - -/// Function extracting an element from a matrix (non-const) -template -requires(Vc::is_simd_vector::value || - algebra::detail::is_storage_vector_v) - ALGEBRA_HOST_DEVICE inline decltype(auto) - element(matrix_t &m, std::size_t row, std::size_t col) { - - return m[col][row]; -} - -} // namespace algebra::vc_aos::math diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp index 3abae506..1df40a6e 100644 --- a/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp +++ b/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp @@ -13,9 +13,7 @@ namespace algebra::vc_aos::math { -using storage::block; using storage::identity; -using storage::set_block; using storage::set_identity; using storage::set_zero; using storage::transpose; diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp index 32680afc..f31126e5 100644 --- a/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp +++ b/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp @@ -11,6 +11,7 @@ #include "algebra/math/common.hpp" #include "algebra/qualifiers.hpp" #include "algebra/storage/matrix.hpp" +#include "algebra/storage/matrix_getter.hpp" #include "algebra/storage/vector.hpp" // Vc include(s). @@ -48,8 +49,7 @@ struct transform3 { public: /// @name Type definitions for the struct /// @{ - /// Index type - using size_type = std::size_t; + /// Scalar type used by the transform using scalar_type = scalar_t; /// The type of the matrix elements (scalar for AoS, Vc::Vector for SoA) diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp index 2fa41738..112e9e3b 100644 --- a/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp +++ b/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp @@ -47,6 +47,18 @@ requires( return (a * b).sum(); } +/// This method retrieves the norm of a vector, no dimension restriction +/// +/// @param v the input vector +template ::value || + algebra::detail::is_storage_vector_v), + bool> = true> +ALGEBRA_HOST_DEVICE inline auto norm(const vector_t &v) { + + return algebra::math::sqrt(dot(v, v)); +} + /// Get a normalized version of the input vector /// /// @tparam vector_type generic input vector type @@ -57,7 +69,20 @@ requires(Vc::is_simd_vector::value || algebra::detail::is_storage_vector_v) ALGEBRA_HOST_DEVICE inline auto normalize(const vector_type &v) { - return v / algebra::math::sqrt(dot(v, v)); + return v / norm(v); +} + +/// This method retrieves the pseudo-rapidity from a vector or vector base with +/// rows >= 3 +/// +/// @param v the input vector +template ::value || + algebra::detail::is_storage_vector_v), + bool> = true> +ALGEBRA_HOST_DEVICE inline auto eta(const vector_t &v) noexcept { + + return algebra::math::atanh(v[2] / norm(v)); } /// Cross product between two input vectors - 3 Dim @@ -78,8 +103,9 @@ requires( inline auto cross(const vector_type1 &a, const vector_type2 &b) -> decltype(a * b - a * b) { - return {a[1] * b[2] - b[1] * a[2], a[2] * b[0] - b[2] * a[0], - a[0] * b[1] - b[0] * a[1], 0.f}; + return {algebra::math::fma(a[1], b[2], -b[1] * a[2]), + algebra::math::fma(a[2], b[0], -b[2] * a[0]), + algebra::math::fma(a[0], b[1], -b[0] * a[1]), 0.f}; } /// Elementwise sum diff --git a/math/vc_aos/include/algebra/math/vc_aos.hpp b/math/vc_aos/include/algebra/math/vc_aos.hpp index ec2cfa23..1421eba5 100644 --- a/math/vc_aos/include/algebra/math/vc_aos.hpp +++ b/math/vc_aos/include/algebra/math/vc_aos.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) 2020-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -8,7 +8,6 @@ #pragma once // Project include(s). -#include "algebra/math/impl/vc_aos_getter.hpp" #include "algebra/math/impl/vc_aos_matrix.hpp" #include "algebra/math/impl/vc_aos_transform3.hpp" #include "algebra/math/impl/vc_aos_vector.hpp" diff --git a/math/vc_soa/CMakeLists.txt b/math/vc_soa/CMakeLists.txt index 9e277e98..d66838b8 100644 --- a/math/vc_soa/CMakeLists.txt +++ b/math/vc_soa/CMakeLists.txt @@ -7,7 +7,6 @@ # Set up the library. 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") target_link_libraries( algebra_vc_soa_math INTERFACE algebra::common algebra::common_math algebra::common_storage algebra::vc_soa_storage Vc::Vc ) diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_getter.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_getter.hpp deleted file mode 100644 index 67842e05..00000000 --- a/math/vc_soa/include/algebra/math/impl/vc_soa_getter.hpp +++ /dev/null @@ -1,123 +0,0 @@ -/** 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/impl/vc_soa_vector.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 { - -/// This method retrieves phi from a vector, vector base with rows >= 2 -/// -/// @tparam N dimension of the vector -/// @tparam scalar_t scalar type -/// @tparam array_t array type that holds the vector elements -/// -/// @param v the input vector -template class array_t, - std::enable_if_t= 2, bool> = true> -ALGEBRA_HOST_DEVICE inline auto phi( - const storage::vector, array_t> &v) { - - return Vc::atan2(v[1], v[0]); -} - -/// This method retrieves the perpenticular magnitude of a vector with rows >= 2 -/// -/// @tparam N dimension of the vector -/// @tparam scalar_t scalar type -/// @tparam array_t array type that holds the vector elements -/// -/// @param v the input vector -template class array_t, - std::enable_if_t= 2, bool> = true> -ALGEBRA_HOST_DEVICE inline auto perp( - const storage::vector, array_t> &v) { - - auto tmp = Vc::Vector::Zero(); - tmp = Vc::fma(v[0], v[0], tmp); - tmp = Vc::fma(v[1], v[1], tmp); - - return Vc::sqrt(tmp); -} - -/// This method retrieves theta from a vector, vector base with rows >= 3 -/// -/// @tparam N dimension of the vector -/// @tparam scalar_t scalar type -/// @tparam array_t array type that holds the vector elements -/// -/// @param v the input vector -template class array_t, - std::enable_if_t= 3, bool> = true> -ALGEBRA_HOST_DEVICE inline auto theta( - const storage::vector, array_t> &v) { - - return Vc::atan2(perp(v), v[2]); -} - -/// This method retrieves the norm of a vector, no dimension restriction -/// -/// @tparam N dimension of the vector -/// @tparam scalar_t scalar type -/// @tparam array_t array type that holds the vector elements -/// -/// @param v the input vector -template class array_t> -ALGEBRA_HOST_DEVICE inline auto norm( - const storage::vector, array_t> &v) { - - return Vc::sqrt(dot(v, v)); -} - -/// This method retrieves the pseudo-rapidity from a vector or vector base with -/// rows >= 3 -/// -/// @tparam N dimension of the vector -/// @tparam scalar_t scalar type -/// @tparam array_t array type that holds the vector elements -/// -/// @param v the input vector -template class array_t, - std::enable_if_t= 3, bool> = true> -ALGEBRA_HOST_DEVICE inline auto eta( - const storage::vector, array_t> &v) { - - // atanh does not exist in Vc - auto atanh_func = [](scalar_t e) { return std::atanh(e); }; - - return (v[2] / norm(v)).apply(atanh_func); - - // Faster, but less accurate - // return (Vc::reciprocal(norm(v)) * v[2]).apply(atanh_func); - - // Even faster, but even less accurate - // return (Vc::rsqrt(dot(v, v)) * v[2]).apply(atanh_func); -} - -} // namespace algebra::vc_soa::math diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp index e3383e6b..5e30e421 100644 --- a/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp +++ b/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp @@ -14,9 +14,7 @@ namespace algebra::vc_soa::math { -using storage::block; using storage::identity; -using storage::set_block; using storage::set_identity; using storage::set_zero; using storage::transpose; 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 afb4d4c8..07fdf65c 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 @@ -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 */ @@ -22,6 +22,51 @@ namespace algebra::vc_soa::math { +/// This method retrieves phi from a vector, vector base with rows >= 2 +/// +/// @tparam N dimension of the vector +/// @tparam scalar_t scalar type +/// @tparam array_t array type that holds the vector elements +/// +/// @param v the input vector +template class array_t> +requires(N >= 2) ALGEBRA_HOST_DEVICE inline auto phi( + const storage::vector, array_t> &v) { + + return Vc::atan2(v[1], v[0]); +} + +/// This method retrieves the perpenticular magnitude of a vector with rows >= 2 +/// +/// @tparam N dimension of the vector +/// @tparam scalar_t scalar type +/// @tparam array_t array type that holds the vector elements +/// +/// @param v the input vector +template class array_t> +requires(N >= 2) ALGEBRA_HOST_DEVICE inline auto perp( + const storage::vector, array_t> &v) { + + return Vc::sqrt(Vc::fma(v[0], v[0], v[1] * v[1])); +} + +/// This method retrieves theta from a vector, vector base with rows >= 3 +/// +/// @tparam N dimension of the vector +/// @tparam scalar_t scalar type +/// @tparam array_t array type that holds the vector elements +/// +/// @param v the input vector +template class array_t> +requires(N >= 3) ALGEBRA_HOST_DEVICE inline auto theta( + const storage::vector, array_t> &v) { + + return Vc::atan2(perp(v), v[2]); +} + /// Cross product between two input vectors - 3 Dim /// /// @tparam N dimension of the vector @@ -33,11 +78,11 @@ namespace algebra::vc_soa::math { /// /// @return a vector (expression) representing the cross product template class array_t, - std::enable_if_t = true> -ALGEBRA_HOST_DEVICE inline storage::vector, array_t> -cross(const storage::vector, array_t> &a, - const storage::vector, array_t> &b) { + template class array_t> +requires(N == 3) ALGEBRA_HOST_DEVICE + inline storage::vector, array_t> cross( + const storage::vector, array_t> &a, + const storage::vector, array_t> &b) { return {Vc::fma(a[1], b[2], -b[1] * a[2]), Vc::fma(a[2], b[0], -b[2] * a[0]), Vc::fma(a[0], b[1], -b[0] * a[1])}; @@ -59,13 +104,30 @@ ALGEBRA_HOST_DEVICE inline Vc::Vector dot( const storage::vector, array_t> &a, const storage::vector, array_t> &b) { - auto ret = Vc::Vector::Zero(); - for (unsigned int i{0u}; i < N; i++) { + auto ret = a[0] * b[0]; + + for (unsigned int i{1u}; i < N; i++) { ret = Vc::fma(a[i], b[i], ret); } + return ret; } +/// This method retrieves the norm of a vector, no dimension restriction +/// +/// @tparam N dimension of the vector +/// @tparam scalar_t scalar type +/// @tparam array_t array type that holds the vector elements +/// +/// @param v the input vector +template class array_t> +ALGEBRA_HOST_DEVICE inline auto norm( + const storage::vector, array_t> &v) { + + return Vc::sqrt(dot(v, v)); +} + /// Get a normalized version of the input vector /// /// @tparam N dimension of the vector @@ -78,15 +140,40 @@ template , array_t> normalize(const storage::vector, array_t> &v) { - return (Vc::Vector::One() / Vc::sqrt(dot(v, v))) * v; + return (Vc::Vector::One() / norm(v)) * v; // Less accurate, but faster - // return Vc::reciprocal(Vc::sqrt(dot(v, v))) * v; + // return Vc::reciprocal(norm(v)) * v; // Even less accurate, but even faster // return Vc::rsqrt(dot(v, v)) * v; } +/// This method retrieves the pseudo-rapidity from a vector or vector base with +/// rows >= 3 +/// +/// @tparam N dimension of the vector +/// @tparam scalar_t scalar type +/// @tparam array_t array type that holds the vector elements +/// +/// @param v the input vector +template class array_t> +requires(N >= 3) ALGEBRA_HOST_DEVICE inline auto eta( + const storage::vector, array_t> &v) { + + // atanh does not exist in Vc + auto atanh_func = [](scalar_t e) { return std::atanh(e); }; + + return (v[2] / norm(v)).apply(atanh_func); + + // Faster, but less accurate + // return (Vc::reciprocal(norm(v)) * v[2]).apply(atanh_func); + + // Even faster, but even less accurate + // return (Vc::rsqrt(dot(v, v)) * v[2]).apply(atanh_func); +} + /// Elementwise sum /// /// @tparam N dimension of the vector diff --git a/math/vc_soa/include/algebra/math/vc_soa.hpp b/math/vc_soa/include/algebra/math/vc_soa.hpp index b469df8f..0a4b1854 100644 --- a/math/vc_soa/include/algebra/math/vc_soa.hpp +++ b/math/vc_soa/include/algebra/math/vc_soa.hpp @@ -8,6 +8,5 @@ #pragma once // Project include(s). -#include "algebra/math/impl/vc_soa_getter.hpp" #include "algebra/math/impl/vc_soa_matrix.hpp" #include "algebra/math/impl/vc_soa_vector.hpp" diff --git a/storage/CMakeLists.txt b/storage/CMakeLists.txt index 019ff865..ce257fd0 100644 --- a/storage/CMakeLists.txt +++ b/storage/CMakeLists.txt @@ -6,6 +6,7 @@ # Set up all enabled libraries. add_subdirectory( array ) +add_subdirectory( cmath ) add_subdirectory( common ) if( ALGEBRA_PLUGINS_INCLUDE_EIGEN ) add_subdirectory( eigen ) diff --git a/storage/array/CMakeLists.txt b/storage/array/CMakeLists.txt index 3a33bf36..8cc612bf 100644 --- a/storage/array/CMakeLists.txt +++ b/storage/array/CMakeLists.txt @@ -8,6 +8,6 @@ algebra_add_library( algebra_array_storage array_storage "include/algebra/storage/array.hpp" ) target_link_libraries( algebra_array_storage - INTERFACE algebra::common ) + INTERFACE algebra::common algebra::common_math algebra::cmath_storage ) algebra_test_public_headers( algebra_array_storage "algebra/storage/array.hpp" ) diff --git a/storage/array/include/algebra/storage/array.hpp b/storage/array/include/algebra/storage/array.hpp index 667977d1..ed79c674 100644 --- a/storage/array/include/algebra/storage/array.hpp +++ b/storage/array/include/algebra/storage/array.hpp @@ -8,6 +8,7 @@ #pragma once // Project include(s) +#include "algebra/storage/impl/cmath_getter.hpp" #include "algebra/type_traits.hpp" // System include(s). @@ -43,37 +44,13 @@ using vector2 = storage_type; template using point2 = vector2; -} // namespace array - -namespace trait { - -/// Type trait specializations -/// @{ -template -struct index, COLS>> { - using type = algebra::array::size_type; -}; - -template -struct dimensions, COLS>> { +/// Element Getter +using element_getter = cmath::storage::element_getter; +/// Block Getter +using block_getter = cmath::storage::block_getter; - using size_type = index_t, COLS>>; - - static constexpr size_type rows{ROWS}; - static constexpr size_type columns{COLS}; -}; - -template -struct value, COLS>> { - using type = T; -}; - -template -struct vector, COLS>> { - using type = std::array; -}; -/// @} +} // namespace array -} // namespace trait +ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(array) } // namespace algebra diff --git a/storage/cmath/CMakeLists.txt b/storage/cmath/CMakeLists.txt new file mode 100644 index 00000000..97c69c80 --- /dev/null +++ b/storage/cmath/CMakeLists.txt @@ -0,0 +1,13 @@ +# Algebra plugins library, part of the ACTS project (R&D line) +# +# (c) 2024 CERN for the benefit of the ACTS project +# +# Mozilla Public License Version 2.0 + +# Set up the library. +algebra_add_library( algebra_cmath_storage cmath_storage + "include/algebra/storage/impl/cmath_getter.hpp" ) +target_link_libraries( algebra_cmath_storage + INTERFACE algebra::common algebra::common_math ) +algebra_test_public_headers( algebra_cmath_storage + "algebra/storage/impl/cmath_getter.hpp" ) diff --git a/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp b/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp new file mode 100644 index 00000000..4ab385a5 --- /dev/null +++ b/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp @@ -0,0 +1,234 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2020-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" + +// System include(s). +#include +#include +#include + +namespace algebra::cmath::storage { + +/// "Element getter", assuming a simple 2D array access +struct element_getter { + + /// Operator getting a reference to one element of a non-const matrix + template class array_t> + ALGEBRA_HOST_DEVICE inline scalar_t &operator()( + array_t, COLS> &m, std::size_t row, + std::size_t col) const { + + assert(row < ROWS); + assert(col < COLS); + return m[col][row]; + } + + /// Operator getting one value of a const matrix + template class array_t> + ALGEBRA_HOST_DEVICE inline scalar_t operator()( + const array_t, COLS> &m, std::size_t row, + std::size_t col) const { + + assert(row < ROWS); + assert(col < COLS); + return m[col][row]; + } + + /// Operator getting a reference to one element of a non-const matrix + template class array_t> + ALGEBRA_HOST_DEVICE inline scalar_t &operator()( + array_t, 1> &m, std::size_t row) const { + + assert(row < N); + return m[0][row]; + } + + /// Operator getting a reference to one element of a const matrix + template class array_t> + ALGEBRA_HOST_DEVICE inline scalar_t operator()( + const array_t, 1> &m, std::size_t row) const { + + assert(row < N); + return m[0][row]; + } + + /// Operator getting a reference to one element of a non-const matrix + template class array_t> + ALGEBRA_HOST_DEVICE inline scalar_t &operator()(array_t &m, + std::size_t row) const { + + assert(row < N); + return m[row]; + } + + /// Operator getting a reference to one element of a const matrix + template class array_t> + ALGEBRA_HOST_DEVICE inline scalar_t operator()(const array_t &m, + std::size_t row) const { + + assert(row < N); + return m[row]; + } + +}; // struct element_getter + +/// Function extracting an element from a matrix (const) +template