From c7c9e73bf9c1132691ebe2785cbb32bbc0646741 Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Fri, 12 Apr 2024 13:27:07 +0200 Subject: [PATCH] Cleanup matrix actors --- benchmarks/array/array_matrix.cpp | 13 + .../benchmark/common/benchmark_matrix.hpp | 55 +-- benchmarks/eigen/eigen_matrix.cpp | 13 + benchmarks/vc_aos/vc_aos_matrix.cpp | 13 + benchmarks/vc_soa/vc_soa_matrix.cpp | 20 ++ common/include/algebra/type_traits.hpp | 91 +++++ extern/vecmem/CMakeLists.txt | 2 +- .../include/algebra/array_cmath.hpp | 108 +----- .../include/algebra/eigen_cmath.hpp | 98 +----- .../include/algebra/eigen_eigen.hpp | 13 +- .../include/algebra/fastor_fastor.hpp | 13 +- .../include/algebra/smatrix_cmath.hpp | 104 ++---- .../include/algebra/smatrix_smatrix.hpp | 15 +- frontend/vc_aos/include/algebra/vc_aos.hpp | 15 +- .../vc_cmath/include/algebra/vc_cmath.hpp | 106 +----- frontend/vc_soa/include/algebra/vc_soa.hpp | 24 +- .../include/algebra/vecmem_cmath.hpp | 107 +----- .../matrix/determinant/partial_pivot_lud.hpp | 56 ---- .../algorithms/utils/algorithm_finder.hpp | 57 ---- .../algebra/math/impl/cmath_getter.hpp | 18 +- .../algebra/math/impl/cmath_matrix.hpp | 312 +++++++----------- .../algebra/math/impl/cmath_transform3.hpp | 246 ++++++-------- .../decomposition/partial_pivot_lud.hpp | 54 ++- .../matrix/determinant/cofactor.hpp | 34 +- .../matrix/determinant/hard_coded.hpp | 28 +- .../matrix/determinant/partial_pivot_lud.hpp | 52 +++ .../algorithms/matrix/inverse/cofactor.hpp | 61 ++-- .../algorithms/matrix/inverse/hard_coded.hpp | 38 +-- .../matrix/inverse/partial_pivot_lud.hpp | 34 +- .../algorithms/utils/algorithm_finder.hpp | 62 ++++ math/common/include/algebra/math/common.hpp | 13 + .../algebra/math/impl/eigen_matrix.hpp | 215 +++++------- .../algebra/math/impl/eigen_transform3.hpp | 41 +-- .../algebra/math/impl/fastor_getter.hpp | 1 - .../algebra/math/impl/fastor_matrix.hpp | 228 ++++++------- .../algebra/math/impl/fastor_transform3.hpp | 15 +- math/smatrix/CMakeLists.txt | 3 +- .../algebra/math/impl/smatrix_matrix.hpp | 212 +++++------- .../algebra/math/impl/smatrix_transform3.hpp | 18 +- .../algebra/math/impl/vc_aos_matrix.hpp | 144 ++------ .../algebra/math/impl/vc_aos_transform3.hpp | 15 +- .../algebra/math/impl/vc_aos_vector.hpp | 61 ++-- .../algebra/math/impl/vc_soa_matrix.hpp | 41 +++ .../algebra/math/impl/vc_soa_vector.hpp | 23 ++ math/vc_soa/include/algebra/math/vc_soa.hpp | 1 + .../array/include/algebra/storage/array.hpp | 44 ++- .../common/include/algebra/storage/matrix.hpp | 14 +- .../common/include/algebra/storage/vector.hpp | 4 +- .../eigen/include/algebra/storage/eigen.hpp | 70 +++- .../fastor/include/algebra/storage/fastor.hpp | 68 +++- .../include/algebra/storage/smatrix.hpp | 44 ++- .../vc_aos/include/algebra/storage/vc_aos.hpp | 42 ++- .../vc_soa/include/algebra/storage/vc_soa.hpp | 46 ++- .../vecmem/include/algebra/storage/vecmem.hpp | 46 ++- tests/accelerator/cuda/CMakeLists.txt | 6 + tests/accelerator/cuda/array_cmath.cu | 10 +- tests/accelerator/cuda/eigen_cmath.cu | 10 +- tests/accelerator/cuda/eigen_eigen.cu | 18 +- tests/accelerator/cuda/vecmem_cmath.cu | 10 +- tests/accelerator/sycl/CMakeLists.txt | 4 + tests/accelerator/sycl/array_cmath.sycl | 10 +- tests/accelerator/sycl/eigen_cmath.sycl | 10 +- tests/accelerator/sycl/eigen_eigen.sycl | 18 +- tests/accelerator/sycl/vecmem_cmath.sycl | 10 +- tests/array/array_cmath.cpp | 10 +- tests/common/test_base.hpp | 11 +- tests/common/test_device_basics.hpp | 29 +- tests/common/test_host_basics.hpp | 137 +++++--- tests/common/test_types.hpp | 9 +- tests/eigen/eigen_cmath.cpp | 10 +- tests/eigen/eigen_eigen.cpp | 18 +- tests/fastor/fastor_fastor.cpp | 13 +- tests/smatrix/smatrix_cmath.cpp | 10 +- tests/smatrix/smatrix_smatrix.cpp | 13 +- tests/vc_aos/vc_aos.cpp | 14 +- tests/vc_aos/vc_cmath.cpp | 4 +- tests/vc_soa/vc_soa.cpp | 46 +++ tests/vecmem/vecmem_cmath.cpp | 10 +- 78 files changed, 1824 insertions(+), 1887 deletions(-) create mode 100644 common/include/algebra/type_traits.hpp delete mode 100644 math/cmath/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp delete mode 100644 math/cmath/include/algebra/math/algorithms/utils/algorithm_finder.hpp rename math/{cmath => common}/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp (61%) rename math/{cmath => common}/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp (72%) rename math/{cmath => common}/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp (83%) create mode 100644 math/common/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp rename math/{cmath => common}/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp (60%) rename math/{cmath => common}/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp (91%) rename math/{cmath => common}/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp (63%) create mode 100644 math/common/include/algebra/math/algorithms/utils/algorithm_finder.hpp create mode 100644 math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp diff --git a/benchmarks/array/array_matrix.cpp b/benchmarks/array/array_matrix.cpp index ac96ba1a..eb284b17 100644 --- a/benchmarks/array/array_matrix.cpp +++ b/benchmarks/array/array_matrix.cpp @@ -27,6 +27,19 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); + using mat44_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_add_f_t = matrix_binaryOP_bm, bench_op::add>; using mat44_add_d_t = diff --git a/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp b/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp index 3da3333b..e17ce107 100644 --- a/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp +++ b/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp @@ -182,29 +182,44 @@ struct mul { } }; +struct transpose { + inline static const std::string name{"transpose"}; + template + auto operator()(const matrix_t& a) const { + return algebra::matrix::transpose(a); + } +}; + } // namespace bench_op // Macro for registering all vector benchmarks -#define ALGEBRA_PLUGINS_REGISTER_MATRIX_BENCH(CFG) \ - algebra::register_benchmark(CFG, "_4x4_add_single"); \ - algebra::register_benchmark(CFG, "_4x4_add_double"); \ - algebra::register_benchmark(CFG, "_6x6_add_single"); \ - algebra::register_benchmark(CFG, "_6x6_add_double"); \ - algebra::register_benchmark(CFG, "_8x8_add_single"); \ - algebra::register_benchmark(CFG, "_8x8_add_double"); \ - \ - algebra::register_benchmark(CFG, "_4x4_mul_single"); \ - algebra::register_benchmark(CFG, "_4x4_mul_double"); \ - algebra::register_benchmark(CFG, "_6x6_mul_single"); \ - algebra::register_benchmark(CFG, "_6x6_mul_double"); \ - algebra::register_benchmark(CFG, "_8x8_mul_single"); \ - algebra::register_benchmark(CFG, "_8x8_mul_double"); \ - \ - algebra::register_benchmark(CFG, "_4x4_vec_single"); \ - algebra::register_benchmark(CFG, "_4x4_vec_double"); \ - algebra::register_benchmark(CFG, "_6x6_vec_single"); \ - algebra::register_benchmark(CFG, "_6x6_vec_double"); \ - algebra::register_benchmark(CFG, "_8x8_vec_single"); \ +#define ALGEBRA_PLUGINS_REGISTER_MATRIX_BENCH(CFG) \ + algebra::register_benchmark(cfg, "_4x4_transpose_single"); \ + algebra::register_benchmark(cfg, "_4x4_transpose_double"); \ + algebra::register_benchmark(cfg, "_6x6_transpose_single"); \ + algebra::register_benchmark(cfg, "_6x6_transpose_double"); \ + algebra::register_benchmark(cfg, "_8x8_transpose_single"); \ + algebra::register_benchmark(cfg, "_8x8_transpose_double"); \ + \ + algebra::register_benchmark(CFG, "_4x4_add_single"); \ + algebra::register_benchmark(CFG, "_4x4_add_double"); \ + algebra::register_benchmark(CFG, "_6x6_add_single"); \ + algebra::register_benchmark(CFG, "_6x6_add_double"); \ + algebra::register_benchmark(CFG, "_8x8_add_single"); \ + algebra::register_benchmark(CFG, "_8x8_add_double"); \ + \ + algebra::register_benchmark(CFG, "_4x4_mul_single"); \ + algebra::register_benchmark(CFG, "_4x4_mul_double"); \ + algebra::register_benchmark(CFG, "_6x6_mul_single"); \ + algebra::register_benchmark(CFG, "_6x6_mul_double"); \ + algebra::register_benchmark(CFG, "_8x8_mul_single"); \ + algebra::register_benchmark(CFG, "_8x8_mul_double"); \ + \ + algebra::register_benchmark(CFG, "_4x4_vec_single"); \ + algebra::register_benchmark(CFG, "_4x4_vec_double"); \ + algebra::register_benchmark(CFG, "_6x6_vec_single"); \ + algebra::register_benchmark(CFG, "_6x6_vec_double"); \ + algebra::register_benchmark(CFG, "_8x8_vec_single"); \ algebra::register_benchmark(CFG, "_8x8_vec_double"); } // namespace algebra diff --git a/benchmarks/eigen/eigen_matrix.cpp b/benchmarks/eigen/eigen_matrix.cpp index b7d30ada..09126733 100644 --- a/benchmarks/eigen/eigen_matrix.cpp +++ b/benchmarks/eigen/eigen_matrix.cpp @@ -27,6 +27,19 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); + using mat44_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_add_f_t = matrix_binaryOP_bm, bench_op::add>; using mat44_add_d_t = diff --git a/benchmarks/vc_aos/vc_aos_matrix.cpp b/benchmarks/vc_aos/vc_aos_matrix.cpp index 95be41d7..a3b7a1d0 100644 --- a/benchmarks/vc_aos/vc_aos_matrix.cpp +++ b/benchmarks/vc_aos/vc_aos_matrix.cpp @@ -27,6 +27,19 @@ int main(int argc, char** argv) { algebra::benchmark_base::configuration cfg{}; cfg.n_samples(100000); + using mat44_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_add_f_t = matrix_binaryOP_bm, bench_op::add>; using mat44_add_d_t = diff --git a/benchmarks/vc_soa/vc_soa_matrix.cpp b/benchmarks/vc_soa/vc_soa_matrix.cpp index 7ff39abb..362d1cb9 100644 --- a/benchmarks/vc_soa/vc_soa_matrix.cpp +++ b/benchmarks/vc_soa/vc_soa_matrix.cpp @@ -36,6 +36,19 @@ 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 mat44_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat66_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_f_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat88_transp_d_t = + matrix_unaryOP_bm, bench_op::transpose>; + using mat44_add_f_t = matrix_binaryOP_bm, bench_op::add>; using mat44_add_d_t = @@ -85,6 +98,13 @@ int main(int argc, char** argv) { // // Register all benchmarks // + algebra::register_benchmark(cfg_s, "_4x4_transpose_single"); + algebra::register_benchmark(cfg_d, "_4x4_transpose_double"); + algebra::register_benchmark(cfg_s, "_6x6_transpose_single"); + algebra::register_benchmark(cfg_d, "_6x6_transpose_double"); + algebra::register_benchmark(cfg_s, "_8x8_transpose_single"); + algebra::register_benchmark(cfg_d, "_8x8_transpose_double"); + algebra::register_benchmark(cfg_s, "_4x4_add_single"); algebra::register_benchmark(cfg_d, "_4x4_add_double"); algebra::register_benchmark(cfg_s, "_6x6_add_single"); diff --git a/common/include/algebra/type_traits.hpp b/common/include/algebra/type_traits.hpp new file mode 100644 index 00000000..2ae0519a --- /dev/null +++ b/common/include/algebra/type_traits.hpp @@ -0,0 +1,91 @@ +/** Algebra plugins library, part of the ACTS project + * + * (c) 2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// System include(s). +#include +#include +#include + +namespace algebra::trait { + +/// Matrix traits +/// @{ + +/// Type of the matrix indices +/// @{ +template +struct index {}; + +template +using index_t = typename index::type; +/// @} + +/// Value type that is used with the matrix (e.g. float or double precision) +/// @{ +template +struct value {}; + +template +using value_t = typename value::type; +/// @} + +/// Scalar type that is used with the matrix (can be multiple values in SoA) +/// @{ +template +struct scalar { + using type = value_t; +}; + +template +using scalar_t = typename scalar::type; +/// @} + +/// Vector type that is compatible with the matrix +/// @{ +template +struct vector {}; + +template +using vector_t = typename vector::type; +/// @} + +/// Matrix dimensions +/// @{ +template +struct dimensions {}; + +/// Specilization for scalar types +template +requires std::is_fundamental_v struct dimensions { + + using size_type = std::size_t; + + static constexpr size_type rows{1}; + static constexpr size_type columns{1}; +}; + +template +inline constexpr index_t rows{dimensions::rows}; + +template +inline constexpr index_t columns{dimensions::columns}; + +template +inline constexpr index_t rank{std::min(rows, columns)}; + +template +inline constexpr index_t size{rows * columns}; + +template +inline constexpr bool is_square{(rows == columns)}; +/// @} + +/// @} + +} // namespace algebra::trait diff --git a/extern/vecmem/CMakeLists.txt b/extern/vecmem/CMakeLists.txt index c235f696..673a8ec0 100644 --- a/extern/vecmem/CMakeLists.txt +++ b/extern/vecmem/CMakeLists.txt @@ -13,7 +13,7 @@ message( STATUS "Building VecMem as part of the Algebra Plugins project" ) # Declare where to get VecMem from. set( ALGEBRA_PLUGINS_VECMEM_SOURCE - "URL;https://github.com/acts-project/vecmem/archive/refs/tags/v1.5.0.tar.gz;URL_MD5;3cc5a3bb14b93f611513535173a6be28" + "URL;https://github.com/acts-project/vecmem/archive/refs/tags/v1.8.0.tar.gz;URL_MD5;afddf52d9568964f25062e1c887246b7" CACHE STRING "Source for VecMem, when built as part of this project" ) mark_as_advanced( ALGEBRA_PLUGINS_VECMEM_SOURCE ) diff --git a/frontend/array_cmath/include/algebra/array_cmath.hpp b/frontend/array_cmath/include/algebra/array_cmath.hpp index cb031fb6..559f200d 100644 --- a/frontend/array_cmath/include/algebra/array_cmath.hpp +++ b/frontend/array_cmath/include/algebra/array_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 */ @@ -70,103 +70,31 @@ using cmath::normalize; namespace matrix { -using size_type = array::size_type; - -template -using array_type = array::storage_type; - -template -using matrix_type = array::matrix_type; - -template -using element_getter = cmath::element_getter; - -template -using block_getter = cmath::block_getter; - -// matrix actor -template -using actor = - cmath::matrix::actor, block_getter>; - -namespace determinant { - -// determinant aggregation -template -using actor = - cmath::matrix::determinant::actor; - -// determinant::cofactor -template -using cofactor = - cmath::matrix::determinant::cofactor, Ds...>; - -// determinant::partial_pivot_lud -template -using partial_pivot_lud = cmath::matrix::determinant::partial_pivot_lud< - size_type, matrix_type, scalar_t, element_getter, Ds...>; - -// determinant::hard_coded -template -using hard_coded = - cmath::matrix::determinant::hard_coded, Ds...>; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace determinant - -namespace inverse { - -// inverion aggregation -template -using actor = - cmath::matrix::inverse::actor; - -// inverse::cofactor -template -using cofactor = - cmath::matrix::inverse::cofactor, Ds...>; - -// inverse::partial_pivot_lud -template -using partial_pivot_lud = - cmath::matrix::inverse::partial_pivot_lud, Ds...>; - -// inverse::hard_coded -template -using hard_coded = - cmath::matrix::inverse::hard_coded, Ds...>; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace inverse +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 array { -/// @name cmath based transforms on @c algebra::matrix::actor +template +using element_getter = + cmath::element_getter; + +/// @name cmath based transforms on @c algebra::array /// @{ template -using transform3_actor = matrix::actor, - matrix::inverse::preset0>; -template -using transform3 = cmath::transform3>; +using transform3 = cmath::transform3, + cmath::block_getter>; /// @} diff --git a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp b/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp index 0fba8fc3..ff6696ff 100644 --- a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp +++ b/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp @@ -73,98 +73,28 @@ using eigen::math::normalize; namespace matrix { -using size_type = eigen::size_type; -template -using array_type = eigen::storage_type; -template -using matrix_type = eigen::matrix_type; -using element_getter = eigen::math::element_getter; -using block_getter = eigen::math::block_getter; - -// matrix actor -template -using actor = cmath::matrix::actor; - -namespace determinant { - -// determinant aggregation -template -using actor = - cmath::matrix::determinant::actor; - -// determinant::cofactor -template -using cofactor = - cmath::matrix::determinant::cofactor; - -// determinant::partial_pivot_lud -template -using partial_pivot_lud = cmath::matrix::determinant::partial_pivot_lud< - size_type, matrix_type, scalar_t, element_getter, Ds...>; - -// determinant::hard_coded -template -using hard_coded = - cmath::matrix::determinant::hard_coded; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace determinant - -namespace inverse { - -// inverion aggregation -template -using actor = - cmath::matrix::inverse::actor; - -// inverse::cofactor -template -using cofactor = - cmath::matrix::inverse::cofactor; - -// inverse::partial_pivot_lud -template -using partial_pivot_lud = - cmath::matrix::inverse::partial_pivot_lud; - -// inverse::hard_coded -template -using hard_coded = - cmath::matrix::inverse::hard_coded; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace inverse +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::matrix::actor +/// @name cmath based transforms on @c algebra::eigen /// @{ template -using transform3_actor = - algebra::matrix::actor, - algebra::matrix::inverse::preset0>; - -template -using transform3 = cmath::transform3>; +using transform3 = + cmath::transform3; /// @} diff --git a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp index 4e502dfa..efe599fa 100644 --- a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp +++ b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp @@ -70,15 +70,22 @@ using eigen::math::normalize; namespace matrix { -template -using actor = eigen::matrix::actor; +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 { template -using transform3 = math::transform3>; +using transform3 = math::transform3; } // namespace eigen diff --git a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp index 463c84f3..1dcb4bd7 100644 --- a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp +++ b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp @@ -73,15 +73,22 @@ using fastor::math::normalize; namespace matrix { -template -using actor = fastor::matrix::actor; +using fastor::math::block; +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 { template -using transform3 = math::transform3>; +using transform3 = math::transform3; } // namespace fastor diff --git a/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp b/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp index cf4a3be5..77522811 100644 --- a/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp +++ b/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp @@ -64,100 +64,34 @@ using smatrix::math::normalize; namespace matrix { -using size_type = smatrix::size_type; -template -using array_type = smatrix::storage_type; -template -using matrix_type = smatrix::matrix_type; -template -using element_getter = smatrix::math::element_getter; -template -using block_getter = smatrix::math::block_getter; - -// matrix actor -template -using actor = - cmath::matrix::actor, block_getter>; - -namespace determinant { - -// determinant aggregation -template -using actor = - cmath::matrix::determinant::actor; - -// determinant::cofactor -template -using cofactor = - cmath::matrix::determinant::cofactor, Ds...>; - -// determinant::partial_pivot_lud -template -using partial_pivot_lud = cmath::matrix::determinant::partial_pivot_lud< - size_type, matrix_type, scalar_t, element_getter, Ds...>; - -// determinant::hard_coded -template -using hard_coded = - cmath::matrix::determinant::hard_coded, Ds...>; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace determinant - -namespace inverse { - -// inverion aggregation -template -using actor = - cmath::matrix::inverse::actor; - -// inverse::cofactor -template -using cofactor = - cmath::matrix::inverse::cofactor, Ds...>; +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; -// inverse::partial_pivot_lud -template -using partial_pivot_lud = - cmath::matrix::inverse::partial_pivot_lud, Ds...>; +} // namespace matrix -// inverse::hard_coded -template -using hard_coded = - cmath::matrix::inverse::hard_coded, Ds...>; +namespace smatrix { -// preset(s) as standard option(s) for user's convenience template -using preset0 = - actor, hard_coded>; - -} // namespace inverse - -} // namespace matrix +using element_getter = smatrix::math::element_getter; -namespace smatrix { +template +using block_getter = smatrix::math::block_getter; -/// @name cmath based transforms on @c algebra::matrix::actor +/// @name cmath based transforms on @c algebra::smatrix /// @{ template -using transform3_actor = - algebra::matrix::actor, - algebra::matrix::inverse::preset0>; -template -using transform3 = cmath::transform3>; +using transform3 = + cmath::transform3, + block_getter>; /// @} diff --git a/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp b/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp index 4c4107e1..e9792aae 100644 --- a/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp +++ b/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.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 */ @@ -61,8 +61,15 @@ using smatrix::math::normalize; namespace matrix { -template -using actor = smatrix::matrix::actor; +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 @@ -72,7 +79,7 @@ namespace smatrix { /// @{ template -using transform3 = math::transform3>; +using transform3 = math::transform3; /// @} diff --git a/frontend/vc_aos/include/algebra/vc_aos.hpp b/frontend/vc_aos/include/algebra/vc_aos.hpp index cbf39786..d18c569d 100644 --- a/frontend/vc_aos/include/algebra/vc_aos.hpp +++ b/frontend/vc_aos/include/algebra/vc_aos.hpp @@ -93,12 +93,15 @@ using vc_aos::math::normalize; namespace matrix { -template -using actor = vc_aos::matrix::actor; - -using storage::identity; -using storage::transpose; -using storage::zero; +using vc_aos::math::block; +using vc_aos::math::determinant; +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; } // namespace matrix diff --git a/frontend/vc_cmath/include/algebra/vc_cmath.hpp b/frontend/vc_cmath/include/algebra/vc_cmath.hpp index 6fecd4da..31767d41 100644 --- a/frontend/vc_cmath/include/algebra/vc_cmath.hpp +++ b/frontend/vc_cmath/include/algebra/vc_cmath.hpp @@ -81,105 +81,29 @@ using vc_aos::math::normalize; namespace matrix { -using size_type = vc_aos::size_type; - -template -using array_type = Vc::array; - -template -using matrix_type = vc_aos::matrix_type; - -template -using element_getter = cmath::element_getter; - -template -using block_getter = cmath::block_getter; - -// matrix actor -template -using actor = - cmath::matrix::actor, block_getter>; - -namespace determinant { - -// determinant aggregation -template -using actor = - cmath::matrix::determinant::actor; - -// determinant::cofactor -template -using cofactor = - cmath::matrix::determinant::cofactor, Ds...>; - -// determinant::partial_pivot_lud -template -using partial_pivot_lud = cmath::matrix::determinant::partial_pivot_lud< - size_type, matrix_type, scalar_t, element_getter, Ds...>; - -// determinant::hard_coded -template -using hard_coded = - cmath::matrix::determinant::hard_coded, Ds...>; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace determinant - -namespace inverse { - -// inverion aggregation -template -using actor = - cmath::matrix::inverse::actor; - -// inverse::cofactor -template -using cofactor = - cmath::matrix::inverse::cofactor, Ds...>; - -// inverse::partial_pivot_lud -template -using partial_pivot_lud = - cmath::matrix::inverse::partial_pivot_lud, Ds...>; - -// inverse::hard_coded -template -using hard_coded = - cmath::matrix::inverse::hard_coded, Ds...>; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace inverse +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 { -/// @name cmath based transforms on @c algebra::matrix::actor -/// @{ +template +using element_getter = cmath::element_getter; -template -using transform3_actor = - algebra::matrix::actor, - algebra::matrix::inverse::preset0>; +/// @name cmath based transforms on @c algebra::vc template -using transform3 = cmath::transform3>; +using transform3 = + cmath::transform3, cmath::block_getter>; /// @} diff --git a/frontend/vc_soa/include/algebra/vc_soa.hpp b/frontend/vc_soa/include/algebra/vc_soa.hpp index 89a15ab6..44b62116 100644 --- a/frontend/vc_soa/include/algebra/vc_soa.hpp +++ b/frontend/vc_soa/include/algebra/vc_soa.hpp @@ -68,16 +68,18 @@ using vc_soa::math::normalize; } // namespace vector // Produces clash with matrix typedefs in other plugins -/*namespace matrix { - -using size_type = vc_soa::size_type; - -template -using array_type = vc_soa::storage_type; - -template -using matrix_type = vc_soa::matrix_type; - -}*/ // namespace matrix +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; +using vc_soa::math::zero; + +} // namespace matrix } // namespace algebra diff --git a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp index 06b2f7f0..585f4140 100644 --- a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp +++ b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp @@ -71,104 +71,31 @@ using cmath::normalize; namespace matrix { -using size_type = std::size_t; - -template -using array_type = vecmem::storage_type; - -template -using matrix_type = vecmem::matrix_type; - -template -using element_getter = cmath::element_getter; - -template -using block_getter = cmath::block_getter; - -// matrix actor -template -using actor = - cmath::matrix::actor, block_getter>; - -namespace determinant { - -// determinant aggregation -template -using actor = - cmath::matrix::determinant::actor; - -// determinant::cofactor -template -using cofactor = - cmath::matrix::determinant::cofactor, Ds...>; - -// determinant::partial_pivot_lud -template -using partial_pivot_lud = cmath::matrix::determinant::partial_pivot_lud< - size_type, matrix_type, scalar_t, element_getter, Ds...>; - -// determinant::hard_coded -template -using hard_coded = - cmath::matrix::determinant::hard_coded, Ds...>; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace determinant - -namespace inverse { - -// inverion aggregation -template -using actor = - cmath::matrix::inverse::actor; - -// inverse::cofactor -template -using cofactor = - cmath::matrix::inverse::cofactor, Ds...>; - -// inverse::partial_pivot_lud -template -using partial_pivot_lud = - cmath::matrix::inverse::partial_pivot_lud, Ds...>; - -// inverse::hard_coded -template -using hard_coded = - cmath::matrix::inverse::hard_coded, Ds...>; - -// preset(s) as standard option(s) for user's convenience -template -using preset0 = - actor, hard_coded>; - -} // namespace inverse +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 vecmem { -/// @name cmath based transforms on @c algebra::matrix::actor -/// @{ +template +using element_getter = + cmath::element_getter; -template -using transform3_actor = matrix::actor, - matrix::inverse::preset0>; +/// @name cmath based transforms on @c algebra::vecmem +/// @{ template -using transform3 = cmath::transform3>; +using transform3 = + cmath::transform3, cmath::block_getter>; /// @} diff --git a/math/cmath/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp b/math/cmath/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp deleted file mode 100644 index b7b5bd3f..00000000 --- a/math/cmath/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp +++ /dev/null @@ -1,56 +0,0 @@ -/** Algebra plugins library, 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/algorithms/matrix/decomposition/partial_pivot_lud.hpp" -#include "algebra/qualifiers.hpp" - -namespace algebra::cmath::matrix::determinant { - -/// "Partial Pivot LU Decomposition", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t, size_type... Ds> -struct partial_pivot_lud { - - using _dims = std::integer_sequence; - - /// Function (object) used for accessing a matrix element - using element_getter = element_getter_t; - - /// 2D matrix type - template - using matrix_type = matrix_t; - - using decomposition_t = - typename algebra::cmath::matrix::decomposition::partial_pivot_lud< - size_type, matrix_t, scalar_t, element_getter_t>; - - template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type& m) const { - - const typename decomposition_t::template lud decomp_res = - decomposition_t()(m); - - // Get the LU decomposition matrix equal to (L - I) + U - const auto& lu = decomp_res.lu; - const size_type n_pivot = static_cast(decomp_res.n_pivot); - - scalar_t det = element_getter_t()(lu, 0, 0); - - for (size_type i = 1; i < N; i++) { - det *= element_getter_t()(lu, i, i); - } - - return (n_pivot - N) % 2 == 0 ? det : -det; - } -}; - -} // namespace algebra::cmath::matrix::determinant \ No newline at end of file diff --git a/math/cmath/include/algebra/math/algorithms/utils/algorithm_finder.hpp b/math/cmath/include/algebra/math/algorithms/utils/algorithm_finder.hpp deleted file mode 100644 index 3d2c9d47..00000000 --- a/math/cmath/include/algebra/math/algorithms/utils/algorithm_finder.hpp +++ /dev/null @@ -1,57 +0,0 @@ -/** Algebra plugins library, part of the ACTS project - * - * (c) 2022 CERN for the benefit of the ACTS project - * - * Mozilla Public License Version 2.0 - */ - -#pragma once - -// System include -#include - -namespace algebra::cmath { - -template -constexpr bool is_in([[maybe_unused]] size_type i, - std::integer_sequence) { - return ((i == Is) || ...); -} - -template -struct find_algorithm_or_default; - -template -struct find_algorithm_or_default { - - using algorithm_type = A; - static constexpr bool found = false; -}; - -template -struct find_algorithm_or_default { - - private: - using next = find_algorithm_or_default; - - public: - using algorithm_type = - typename std::conditional_t; - - static constexpr bool found = is_in(N, typename A::_dims{}) || next::found; -}; - -template -struct find_algorithm { - private: - using helper = - find_algorithm_or_default; - - public: - using algorithm_type = typename helper::algorithm_type; - static constexpr bool found = helper::found; -}; - -} // namespace algebra::cmath diff --git a/math/cmath/include/algebra/math/impl/cmath_getter.hpp b/math/cmath/include/algebra/math/impl/cmath_getter.hpp index ef0c31b5..03890272 100644 --- a/math/cmath/include/algebra/math/impl/cmath_getter.hpp +++ b/math/cmath/include/algebra/math/impl/cmath_getter.hpp @@ -176,20 +176,18 @@ struct vector_getter { }; // struct vector_getter /// "Block getter", assuming a simple 2D array access -template class array_t, - typename scalar_t> struct block_getter { - /// 2D matrix type - template - using matrix_type = array_t, COLS>; - /// Operator producing a sub-matrix from a const matrix - template - ALGEBRA_HOST_DEVICE matrix_type operator()( - const input_matrix_type &m, std::size_t row, std::size_t col) const { + 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{}; - matrix_type 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]; diff --git a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp index c1cb3b42..93073188 100644 --- a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp +++ b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (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 */ @@ -9,216 +9,150 @@ // Project include(s). #include "algebra/math/algorithms/utils/algorithm_finder.hpp" +#include "algebra/math/common.hpp" #include "algebra/qualifiers.hpp" -namespace algebra::cmath::matrix { - -/// "Matrix actor", assuming a simple 2D matrix -template class array_t, - template class matrix_t, - typename scalar_t, class determinant_actor_t, class inverse_actor_t, - class element_getter_t, class block_getter_t> -struct actor { - - /// Size type - using size_ty = size_type; - - /// Scalar type - using scalar_type = scalar_t; - - /// Function (object) used for accessing a matrix element - using element_getter = element_getter_t; - - /// Function (object) used for accessing a matrix block - using block_getter = block_getter_t; - - /// 2D matrix type - template - using matrix_type = matrix_t; - - /// vector type - template - using array_type = array_t; - - /// Operator getting a reference to one element of a non-const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t &element(matrix_type &m, - size_type row, - size_type col) const { - return element_getter()(m, row, col); - } - - /// Operator getting one value of a const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t element(const matrix_type &m, - size_type row, - size_type col) const { - return element_getter()(m, row, col); - } - - /// Operator getting a block of a const matrix - template - ALGEBRA_HOST_DEVICE matrix_type block(const input_matrix_type &m, - size_type row, - size_type col) const { - return block_getter().template operator()(m, row, col); - } - - /// Operator setting a block with a vector matrix - template - ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const matrix_type &b, - size_type row, size_type col) const { - for (size_type j = 0; j < COLS; ++j) { - for (size_type i = 0; i < ROWS; ++i) { - element_getter()(m, i + row, j + col) = element_getter()(b, i, j); - } - } - } - - /// Operator setting a block with a vector - template class vector_t, - class input_matrix_type> - ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const vector_t &b, - size_type row, size_type col) const { - for (size_type i = 0; i < ROWS; ++i) { - element_getter()(m, i + row, col) = b[i]; - } - } - - // Create zero matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type zero() const { - matrix_type ret; - - for (size_type j = 0; j < COLS; ++j) { - for (size_type i = 0; i < ROWS; ++i) { - element_getter()(ret, i, j) = 0; - } - } - - return ret; - } - - // Create identity matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type identity() const { - matrix_type ret; - - for (size_type j = 0; j < COLS; ++j) { - for (size_type i = 0; i < ROWS; ++i) { - if (i == j) { - element_getter()(ret, i, j) = 1; - } else { - element_getter()(ret, i, j) = 0; - } - } +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]; } - - return ret; } - - // Set input matrix as zero matrix - template - ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type &m) const { - - for (size_type j = 0; j < COLS; ++j) { - for (size_type i = 0; i < ROWS; ++i) { - element_getter()(m, i, j) = 0; - } - } +} + +/// 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]; } - - // Set input matrix as identity matrix - template - ALGEBRA_HOST_DEVICE inline void set_identity( - matrix_type &m) const { - - for (size_type j = 0; j < COLS; ++j) { - for (size_type i = 0; i < ROWS; ++i) { - if (i == j) { - element_getter()(m, i, j) = 1; - } else { - element_getter()(m, i, j) = 0; - } - } +} + +/// @returns zero matrix of type @tparam matrix_t +template +requires(std::is_scalar_v) + ALGEBRA_HOST_DEVICE inline matrix_t zero() { + 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; } } - // Create transpose matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type transpose( - const matrix_type &m) const { + return ret; +} - matrix_type ret; +/// @returns identity matrix of type @tparam matrix_t +template +requires(std::is_scalar_v) + ALGEBRA_HOST_DEVICE inline auto identity() { + auto ret{zero()}; - for (size_type i = 0; i < ROWS; ++i) { - for (size_type j = 0; j < COLS; ++j) { - element_getter()(ret, j, i) = element_getter()(m, i, j); - } - } - - return ret; + for (std::size_t i = 0; i < algebra::trait::rank; ++i) { + ret[i][i] = 1; } - // Get determinant - template - ALGEBRA_HOST_DEVICE inline scalar_t determinant( - const matrix_type &m) const { + return ret; +} - return determinant_actor_t()(m); - } +/// Create identity matrix - cmath transform3 +template +ALGEBRA_HOST_DEVICE inline matrix_t identity() { + auto ret{zero()}; - // Create inverse matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type inverse( - const matrix_type &m) const { - - return inverse_actor_t()(m); + for (std::size_t i = 0; i < algebra::trait::rank; ++i) { + element_getter_t{}(ret, i, i) = 1; } -}; - -namespace determinant { -template class matrix_t, - typename scalar_t, class... As> -struct actor { - - /// 2D matrix type - template - using matrix_type = matrix_t; - - template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type &m) const { - - return typename find_algorithm::algorithm_type()(m); + return ret; +} + +/// Set @param m as zero matrix +template class array_t> +ALGEBRA_HOST_DEVICE inline void set_zero( + array_t, COLS> &m) { + m = zero, COLS>>(); +} + +/// Set @param m as identity matrix +template class array_t> +ALGEBRA_HOST_DEVICE inline 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]; + } } -}; -} // namespace determinant + return ret; +} -namespace inverse { +/// @returns the determinant of @param m +template class array_t> +ALGEBRA_HOST_DEVICE inline scalar_t determinant( + const array_t, N> &m) { -template class matrix_t, - typename scalar_t, class... As> -struct actor { + using matrix_t = array_t, N>; + using element_getter_t = element_getter; - /// 2D matrix type - template - using matrix_type = matrix_t; + return determinant_t{}(m); +} - template - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type &m) const { +/// @returns the determinant of @param m +template class array_t> +ALGEBRA_HOST_DEVICE inline array_t, N> inverse( + const array_t, N> &m) { - return typename find_algorithm::algorithm_type()(m); - } -}; + using matrix_t = array_t, N>; + using element_getter_t = element_getter; -} // namespace inverse + return inversion_t{}(m); +} -} // namespace algebra::cmath::matrix +} // namespace algebra::cmath diff --git a/math/cmath/include/algebra/math/impl/cmath_transform3.hpp b/math/cmath/include/algebra/math/impl/cmath_transform3.hpp index 262daf95..950ed8db 100644 --- a/math/cmath/include/algebra/math/impl/cmath_transform3.hpp +++ b/math/cmath/include/algebra/math/impl/cmath_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 */ @@ -8,7 +8,9 @@ #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" @@ -16,28 +18,28 @@ namespace algebra::cmath { /** Transform wrapper class to ensure standard API within differnt plugins **/ -template +template class matrix_t, + template class array_t, typename element_getter_t, + typename block_getter_t> struct transform3 { /// @name Type definitions for the struct /// @{ - /// Matrix actor - using matrix_actor = matrix_actor_t; - /// Size type - using size_type = typename matrix_actor_t::size_ty; + using size_type = size_ty; /// Scalar type - using scalar_type = typename matrix_actor_t::scalar_type; + using scalar_type = scalar_t; /// 2D Matrix type template - using matrix_type = typename matrix_actor::template matrix_type; + using matrix_type = matrix_t; /// Array type template - using array_type = typename matrix_actor::template array_type; + using array_type = array_t; // 4 x 4 Matrix using matrix44 = matrix_type<4, 4>; @@ -50,21 +52,28 @@ struct transform3 { using point2 = array_type<2>; /// Function (object) used for accessing a matrix element - using element_getter = typename matrix_actor::element_getter; + using element_getter = element_getter_t; + + /// Function (object) used for accessing a matrix element + using block_getter = block_getter_t; - /// Function (object) used for accessing a matrix block - using block_getter = typename matrix_actor::block_getter; + /// Matrix inversion algorithm + using matrix_inversion = + cmath::matrix::inverse::hard_coded; /// @} /// @name Data objects /// @{ - - matrix44 _data = matrix_actor().template identity<4, 4>(); - matrix44 _data_inv = matrix_actor().template identity<4, 4>(); + 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) @@ -77,25 +86,25 @@ struct transform3 { transform3(const vector3 &t, const vector3 &x, const vector3 &y, const vector3 &z, bool get_inverse = true) { - matrix_actor().element(_data, 0, 0) = x[0]; - matrix_actor().element(_data, 1, 0) = x[1]; - matrix_actor().element(_data, 2, 0) = x[2]; - matrix_actor().element(_data, 3, 0) = 0.; - matrix_actor().element(_data, 0, 1) = y[0]; - matrix_actor().element(_data, 1, 1) = y[1]; - matrix_actor().element(_data, 2, 1) = y[2]; - matrix_actor().element(_data, 3, 1) = 0.; - matrix_actor().element(_data, 0, 2) = z[0]; - matrix_actor().element(_data, 1, 2) = z[1]; - matrix_actor().element(_data, 2, 2) = z[2]; - matrix_actor().element(_data, 3, 2) = 0.; - matrix_actor().element(_data, 0, 3) = t[0]; - matrix_actor().element(_data, 1, 3) = t[1]; - matrix_actor().element(_data, 2, 3) = t[2]; - matrix_actor().element(_data, 3, 3) = 1.; + 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_actor().inverse(_data); + _data_inv = matrix_inversion{}(_data); } } @@ -118,26 +127,26 @@ struct transform3 { * @param t is the transform **/ ALGEBRA_HOST_DEVICE - transform3(const vector3 &t) { - - matrix_actor().element(_data, 0, 0) = 1.; - matrix_actor().element(_data, 1, 0) = 0.; - matrix_actor().element(_data, 2, 0) = 0.; - matrix_actor().element(_data, 3, 0) = 0.; - matrix_actor().element(_data, 0, 1) = 0.; - matrix_actor().element(_data, 1, 1) = 1.; - matrix_actor().element(_data, 2, 1) = 0.; - matrix_actor().element(_data, 3, 1) = 0.; - matrix_actor().element(_data, 0, 2) = 0.; - matrix_actor().element(_data, 1, 2) = 0.; - matrix_actor().element(_data, 2, 2) = 1.; - matrix_actor().element(_data, 3, 2) = 0.; - matrix_actor().element(_data, 0, 3) = t[0]; - matrix_actor().element(_data, 1, 3) = t[1]; - matrix_actor().element(_data, 2, 3) = t[2]; - matrix_actor().element(_data, 3, 3) = 1.; - - _data_inv = matrix_actor().inverse(_data); + 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 @@ -145,10 +154,10 @@ struct transform3 { * @param m is the full 4x4 matrix **/ ALGEBRA_HOST_DEVICE - transform3(const matrix44 &m) { + explicit transform3(const matrix44 &m) { _data = m; - _data_inv = matrix_actor().inverse(_data); + _data_inv = matrix_inversion{}(_data); } /** Constructor with arguments: matrix as array of scalar @@ -156,35 +165,26 @@ struct transform3 { * @param ma is the full 4x4 matrix 16 array **/ ALGEBRA_HOST_DEVICE - transform3(const array_type<16> &ma) { - - matrix_actor().element(_data, 0, 0) = ma[0]; - matrix_actor().element(_data, 1, 0) = ma[4]; - matrix_actor().element(_data, 2, 0) = ma[8]; - matrix_actor().element(_data, 3, 0) = ma[12]; - matrix_actor().element(_data, 0, 1) = ma[1]; - matrix_actor().element(_data, 1, 1) = ma[5]; - matrix_actor().element(_data, 2, 1) = ma[9]; - matrix_actor().element(_data, 3, 1) = ma[13]; - matrix_actor().element(_data, 0, 2) = ma[2]; - matrix_actor().element(_data, 1, 2) = ma[6]; - matrix_actor().element(_data, 2, 2) = ma[10]; - matrix_actor().element(_data, 3, 2) = ma[14]; - matrix_actor().element(_data, 0, 3) = ma[3]; - matrix_actor().element(_data, 1, 3) = ma[7]; - matrix_actor().element(_data, 2, 3) = ma[11]; - matrix_actor().element(_data, 3, 3) = ma[15]; - - _data_inv = matrix_actor().inverse(_data); - } - - /** Constructor with arguments: identity - * - **/ - ALGEBRA_HOST_DEVICE - transform3() { - matrix_actor().set_identity(_data); - _data_inv = _data; + 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 */ @@ -193,8 +193,8 @@ struct transform3 { for (size_type j = 0; j < 4; j++) { for (size_type i = 0; i < 4; i++) { - if (matrix_actor().element(_data, i, j) != - matrix_actor().element(rhs._data, i, j)) { + if (element_getter{}(_data, i, j) != + element_getter{}(rhs._data, i, j)) { return false; } } @@ -203,30 +203,6 @@ struct transform3 { return true; } - /** The determinant of a 4x4 matrix - * - * @param m is the matrix - * - * @return a sacalar determinant - no checking done - */ - ALGEBRA_HOST_DEVICE - static inline scalar_type determinant(const matrix44 &m) { - - return matrix_actor().determinant(m); - } - - /** The inverse of a 4x4 matrix - * - * @param m is the matrix - * - * @return an inverse matrix - */ - ALGEBRA_HOST_DEVICE - static inline matrix44 invert(const matrix44 &m) { - - return matrix_actor().inverse(m); - } - /** Rotate a vector into / from a frame * * @param m is the rotation matrix @@ -237,17 +213,17 @@ struct transform3 { vector3 ret{0.f, 0.f, 0.f}; - ret[0] += matrix_actor().element(m, 0, 0) * v[0]; - ret[1] += matrix_actor().element(m, 1, 0) * v[0]; - ret[2] += matrix_actor().element(m, 2, 0) * v[0]; + 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] += matrix_actor().element(m, 0, 1) * v[1]; - ret[1] += matrix_actor().element(m, 1, 1) * v[1]; - ret[2] += matrix_actor().element(m, 2, 1) * v[1]; + 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] += matrix_actor().element(m, 0, 2) * v[2]; - ret[1] += matrix_actor().element(m, 1, 2) * v[2]; - ret[2] += matrix_actor().element(m, 2, 2) * v[2]; + 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; } @@ -256,40 +232,36 @@ struct transform3 { ALGEBRA_HOST_DEVICE auto inline rotation() const { - return matrix_actor().template block<3, 3>(_data, 0, 0); + return block_getter{}.template operator()<3, 3>(_data, 0, 0); } /** This method retrieves x axis */ ALGEBRA_HOST_DEVICE inline point3 x() const { - return {matrix_actor().element(_data, 0, 0), - matrix_actor().element(_data, 1, 0), - matrix_actor().element(_data, 2, 0)}; + 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 {matrix_actor().element(_data, 0, 1), - matrix_actor().element(_data, 1, 1), - matrix_actor().element(_data, 2, 1)}; + 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 {matrix_actor().element(_data, 0, 2), - matrix_actor().element(_data, 1, 2), - matrix_actor().element(_data, 2, 2)}; + 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 {matrix_actor().element(_data, 0, 3), - matrix_actor().element(_data, 1, 3), - matrix_actor().element(_data, 2, 3)}; + 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 */ @@ -305,9 +277,9 @@ struct transform3 { ALGEBRA_HOST_DEVICE inline point3 point_to_global(const point3 &v) const { const vector3 rg = rotate(_data, v); - return {rg[0] + matrix_actor().element(_data, 0, 3), - rg[1] + matrix_actor().element(_data, 1, 3), - rg[2] + matrix_actor().element(_data, 2, 3)}; + 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 @@ -316,9 +288,9 @@ struct transform3 { ALGEBRA_HOST_DEVICE inline point3 point_to_local(const point3 &v) const { const vector3 rg = rotate(_data_inv, v); - return {rg[0] + matrix_actor().element(_data_inv, 0, 3), - rg[1] + matrix_actor().element(_data_inv, 1, 3), - rg[2] + matrix_actor().element(_data_inv, 2, 3)}; + 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 diff --git a/math/cmath/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp b/math/common/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp similarity index 61% rename from math/cmath/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp rename to math/common/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp index b7e7bda5..3c91a7d5 100644 --- a/math/cmath/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp +++ b/math/common/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp @@ -17,57 +17,55 @@ namespace algebra::cmath::matrix::decomposition { /// "Partial Pivot LU Decomposition", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t, size_type... Ds> +template struct partial_pivot_lud { - using _dims = std::integer_sequence; + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; + using vector_type = algebra::trait::vector_t; /// Function (object) used for accessing a matrix element using element_getter = element_getter_t; - /// 2D matrix type - template - using matrix_type = matrix_t; - template struct lud { // LU decomposition matrix, equal to (L - I) + U, where the diagonal // components of L is always 1 - matrix_type lu; + matrix_t lu; // Permutation vector - matrix_type<1, N> P; + vector_type P; // Number of pivots int n_pivot = 0; }; - template - ALGEBRA_HOST_DEVICE inline lud operator()( - const matrix_type& m) const { + ALGEBRA_HOST_DEVICE inline lud> operator()( + const matrix_t& m) const { + + constexpr size_type N{algebra::trait::rank}; + // LU decomposition matrix - matrix_type lu = m; + matrix_t lu = m; // Permutation - matrix_type<1, N> P; + vector_type P; // Max index and value size_type max_idx; - scalar_t max_val; - scalar_t abs_val; + scalar_type max_val; + scalar_type abs_val; // Number of pivoting int n_pivot = N; // Rows for swapping - matrix_type<1, N> row_0; - matrix_type<1, N> row_1; + vector_type row_0; + vector_type row_1; // Unit permutation matrix, P[N] initialized with N for (size_type i = 0; i < N; i++) { - element_getter_t()(P, 0, i) = static_cast(i); + P[i] = static_cast(i); } for (size_type i = 0; i < N; i++) { @@ -86,19 +84,19 @@ struct partial_pivot_lud { if (max_idx != i) { // Pivoting P - auto j = element_getter_t()(P, 0, i); + auto j = P[i]; - element_getter_t()(P, 0, i) = element_getter_t()(P, 0, max_idx); - element_getter_t()(P, 0, max_idx) = j; + P[i] = P[max_idx]; + P[max_idx] = j; // Pivoting rows of A for (size_type q = 0; q < N; q++) { - element_getter_t()(row_0, 0, q) = element_getter_t()(lu, i, q); - element_getter_t()(row_1, 0, q) = element_getter_t()(lu, max_idx, q); + row_0[q] = element_getter_t()(lu, i, q); + row_1[q] = element_getter_t()(lu, max_idx, q); } for (size_type q = 0; q < N; q++) { - element_getter_t()(lu, i, q) = element_getter_t()(row_1, 0, q); - element_getter_t()(lu, max_idx, q) = element_getter_t()(row_0, 0, q); + element_getter_t()(lu, i, q) = row_1[q]; + element_getter_t()(lu, max_idx, q) = row_0[q]; } // counting pivots starting from N (for determinant) @@ -121,4 +119,4 @@ struct partial_pivot_lud { } }; -} // namespace algebra::cmath::matrix::decomposition \ No newline at end of file +} // namespace algebra::cmath::matrix::decomposition diff --git a/math/cmath/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp b/math/common/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp similarity index 72% rename from math/cmath/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp rename to math/common/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp index 283cf887..f0b83788 100644 --- a/math/cmath/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp +++ b/math/common/include/algebra/math/algorithms/matrix/determinant/cofactor.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 */ @@ -9,6 +9,7 @@ // Project include(s). #include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" // System include(s) #include @@ -16,24 +17,17 @@ namespace algebra::cmath::matrix::determinant { /// "Determinant getter", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t, size_type... Ds> +template struct cofactor { - using _dims = std::integer_sequence; + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; /// Function (object) used for accessing a matrix element using element_getter = element_getter_t; - /// 2D matrix type - template - using matrix_type = matrix_t; - - template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type &m) const { - return determinant_getter_helper()(m); + ALGEBRA_HOST_DEVICE inline scalar_type operator()(const matrix_t &m) const { + return determinant_getter_helper>()(m); } template @@ -42,7 +36,7 @@ struct cofactor { template struct determinant_getter_helper> { template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( + ALGEBRA_HOST_DEVICE inline scalar_type operator()( const input_matrix_type &m) const { return element_getter()(m, 0, 0); } @@ -52,13 +46,13 @@ struct cofactor { struct determinant_getter_helper> { template - ALGEBRA_HOST_DEVICE inline scalar_t operator()( + ALGEBRA_HOST_DEVICE inline scalar_type operator()( const input_matrix_type &m) const { - scalar_t D = 0; + scalar_type D = 0; // To store cofactors - matrix_type temp; + matrix_t temp; // To store sign multiplier int sign = 1; @@ -79,11 +73,11 @@ struct cofactor { template ALGEBRA_HOST_DEVICE inline void get_cofactor(const input_matrix_type &m, - matrix_type &temp, - size_type p, + matrix_t &temp, size_type p, size_type q) const { - size_type i = 0, j = 0; + size_type i = 0; + size_type j = 0; // Looping for each element of the matrix for (size_type row = 0; row < N; row++) { diff --git a/math/cmath/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp b/math/common/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp similarity index 83% rename from math/cmath/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp rename to math/common/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp index 8bfa5cf1..ef754b86 100644 --- a/math/cmath/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp +++ b/math/common/include/algebra/math/algorithms/matrix/determinant/hard_coded.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 */ @@ -9,6 +9,7 @@ // Project include(s). #include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" // System include(s) #include @@ -16,33 +17,28 @@ namespace algebra::cmath::matrix::determinant { /// "Determinant getter", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t, size_type... Ds> +template struct hard_coded { - using _dims = std::integer_sequence; + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; /// Function (object) used for accessing a matrix element using element_getter = element_getter_t; - /// 2D matrix type - template - using matrix_type = matrix_t; - // 2 X 2 matrix determinant - template = true> - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type &m) const { + template + requires(algebra::trait::rank == 2) ALGEBRA_HOST_DEVICE inline scalar_type + operator()(const matrix_t &m) const { return element_getter()(m, 0, 0) * element_getter()(m, 1, 1) - element_getter()(m, 0, 1) * element_getter()(m, 1, 0); } // 4 X 4 matrix determinant - template = true> - ALGEBRA_HOST_DEVICE inline scalar_t operator()( - const matrix_type &m) const { + template + requires(algebra::trait::rank == 4) ALGEBRA_HOST_DEVICE inline scalar_type + operator()(const matrix_t &m) const { return element_getter()(m, 0, 3) * element_getter()(m, 1, 2) * element_getter()(m, 2, 1) * element_getter()(m, 3, 0) - @@ -95,4 +91,4 @@ struct hard_coded { } }; -} // namespace algebra::cmath::matrix::determinant \ No newline at end of file +} // namespace algebra::cmath::matrix::determinant diff --git a/math/common/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp b/math/common/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp new file mode 100644 index 00000000..472fca4d --- /dev/null +++ b/math/common/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp @@ -0,0 +1,52 @@ +/** 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/matrix/decomposition/partial_pivot_lud.hpp" +#include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" + +namespace algebra::cmath::matrix::determinant { + +/// "Partial Pivot LU Decomposition", assuming a N X N matrix +template +struct partial_pivot_lud { + + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; + + /// Function (object) used for accessing a matrix element + using element_getter = element_getter_t; + + using decomposition_t = + typename algebra::cmath::matrix::decomposition::partial_pivot_lud< + matrix_t, element_getter_t>; + + ALGEBRA_HOST_DEVICE inline scalar_type operator()(const matrix_t& m) const { + + constexpr size_type N{algebra::trait::rank}; + + const typename decomposition_t::template lud> + decomp_res = decomposition_t()(m); + + // Get the LU decomposition matrix equal to (L - I) + U + const auto& lu = decomp_res.lu; + const auto n_pivot = static_cast(decomp_res.n_pivot); + + scalar_type det = element_getter()(lu, 0, 0); + + for (size_type i = 1; i < N; i++) { + det *= element_getter()(lu, i, i); + } + + return (n_pivot - N) % 2 == 0 ? det : -det; + } +}; + +} // namespace algebra::cmath::matrix::determinant diff --git a/math/cmath/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp b/math/common/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp similarity index 60% rename from math/cmath/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp rename to math/common/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp index b9d67cd4..a68d701b 100644 --- a/math/cmath/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp +++ b/math/common/include/algebra/math/algorithms/matrix/inverse/cofactor.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 */ @@ -10,6 +10,7 @@ // Project include(s). #include "algebra/math/algorithms/matrix/determinant/cofactor.hpp" #include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" // System include(s) #include @@ -19,22 +20,17 @@ namespace algebra::cmath::matrix { namespace adjoint { /// "Adjoint getter", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t> +template struct cofactor { + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; + /// Function (object) used for accessing a matrix element using element_getter = element_getter_t; - /// 2D matrix type - template - using matrix_type = matrix_t; - - template - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type &m) const { - return adjoint_getter_helper()(m); + ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t &m) const { + return adjoint_getter_helper>()(m); } template @@ -42,10 +38,9 @@ struct cofactor { template struct adjoint_getter_helper> { - - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type & /*m*/) const { - matrix_type ret; + ALGEBRA_HOST_DEVICE inline matrix_t operator()( + const matrix_t & /*m*/) const { + matrix_t ret; element_getter()(ret, 0, 0) = 1; return ret; } @@ -55,18 +50,17 @@ struct cofactor { struct adjoint_getter_helper> { using determinant_getter = - determinant::cofactor; + determinant::cofactor; - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type &m) const { + ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t &m) const { - matrix_type adj; + matrix_t adj; // temp is used to store cofactors of m int sign = 1; // To store cofactors - matrix_type temp; + matrix_t temp; for (size_type i = 0; i < N; i++) { for (size_type j = 0; j < N; j++) { @@ -97,34 +91,27 @@ struct cofactor { namespace inverse { /// "inverse getter", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t, size_type... Ds> +template struct cofactor { - using _dims = std::integer_sequence; + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; /// Function (object) used for accessing a matrix element using element_getter = element_getter_t; - /// 2D matrix type - template - using matrix_type = matrix_t; + using determinant_getter = determinant::cofactor; - using determinant_getter = - determinant::cofactor; + using adjoint_getter = adjoint::cofactor; - using adjoint_getter = - adjoint::cofactor; + ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t &m) const { - template - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type &m) const { + constexpr size_type N{algebra::trait::rank}; - matrix_type ret; + matrix_t ret; // Find determinant of A - scalar_t det = determinant_getter()(m); + scalar_type det = determinant_getter()(m); // TODO: handle singular matrix error // if (det == 0) { diff --git a/math/cmath/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp b/math/common/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp similarity index 91% rename from math/cmath/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp rename to math/common/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp index 82818c9c..8736dd9c 100644 --- a/math/cmath/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp +++ b/math/common/include/algebra/math/algorithms/matrix/inverse/hard_coded.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 */ @@ -10,6 +10,7 @@ // Project include(s). #include "algebra/math/algorithms/matrix/determinant/hard_coded.hpp" #include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" // System include(s) #include @@ -17,31 +18,26 @@ namespace algebra::cmath::matrix::inverse { /// "inverse getter", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t, size_type... Ds> +template struct hard_coded { - using _dims = std::integer_sequence; + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; /// Function (object) used for accessing a matrix element using element_getter = element_getter_t; - /// 2D matrix type - template - using matrix_type = matrix_t; - using determinant_getter = - determinant::hard_coded; + determinant::hard_coded; // 2 X 2 matrix inverse - template = true> - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type &m) const { + template + requires(algebra::trait::rank == 2) ALGEBRA_HOST_DEVICE inline matrix_t + operator()(const matrix_t &m) const { - matrix_type ret; + matrix_t ret; - scalar_t det = determinant_getter()(m); + scalar_type det = determinant_getter()(m); element_getter()(ret, 0, 0) = element_getter()(m, 1, 1) / det; element_getter()(ret, 0, 1) = -1 * element_getter()(m, 0, 1) / det; @@ -52,11 +48,11 @@ struct hard_coded { } // 4 X 4 matrix inverse - template = true> - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type &m) const { + template + requires(algebra::trait::rank == 4) ALGEBRA_HOST_DEVICE inline matrix_t + operator()(const matrix_t &m) const { - matrix_type ret; + matrix_t ret; element_getter()(ret, 0, 0) = element_getter()(m, 1, 2) * element_getter()(m, 2, 3) * element_getter()(m, 3, 1) - @@ -266,7 +262,7 @@ struct hard_coded { element_getter()(m, 0, 0) * element_getter()(m, 1, 1) * element_getter()(m, 2, 2); - scalar_t idet = static_cast(1.) / determinant_getter()(ret); + scalar_type idet = static_cast(1.) / determinant_getter()(ret); for (unsigned int c = 0; c < 4; ++c) { for (unsigned int r = 0; r < 4; ++r) { element_getter()(ret, c, r) *= idet; @@ -276,4 +272,4 @@ struct hard_coded { } }; -} // namespace algebra::cmath::matrix::inverse \ No newline at end of file +} // namespace algebra::cmath::matrix::inverse diff --git a/math/cmath/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp b/math/common/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp similarity index 63% rename from math/cmath/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp rename to math/common/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp index 327b69cb..6bd64ea9 100644 --- a/math/cmath/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp +++ b/math/common/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (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 */ @@ -10,31 +10,28 @@ // Project include(s). #include "algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp" #include "algebra/qualifiers.hpp" +#include "algebra/type_traits.hpp" namespace algebra::cmath::matrix::inverse { /// "Partial Pivot LU Decomposition", assuming a N X N matrix -template class matrix_t, - typename scalar_t, class element_getter_t, size_type... Ds> +template struct partial_pivot_lud { - using _dims = std::integer_sequence; + using scalar_type = algebra::trait::value_t; + using size_type = algebra::trait::index_t; /// Function (object) used for accessing a matrix element using element_getter = element_getter_t; - /// 2D matrix type - template - using matrix_type = matrix_t; - using decomposition_t = typename algebra::cmath::matrix::decomposition::partial_pivot_lud< - size_type, matrix_t, scalar_t, element_getter_t>; + matrix_t, element_getter_t>; + + ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t& m) const { + + constexpr size_type N{algebra::trait::rank}; - template - ALGEBRA_HOST_DEVICE inline matrix_type operator()( - const matrix_type& m) const { const typename decomposition_t::template lud decomp_res = decomposition_t()(m); @@ -45,15 +42,14 @@ struct partial_pivot_lud { const auto& P = decomp_res.P; // Inverse matrix - matrix_type inv; + matrix_t inv; // Calculate inv(A) = inv(U) * inv(L) * P; for (size_type j = 0; j < N; j++) { for (size_type i = 0; i < N; i++) { - element_getter_t()(inv, i, j) = - static_cast(element_getter_t()(P, 0, i)) == j - ? static_cast(1.0) - : static_cast(0.0); + element_getter_t()(inv, i, j) = static_cast(P[i]) == j + ? static_cast(1.0) + : static_cast(0.0); for (size_type k = 0; k < i; k++) { element_getter_t()(inv, i, j) -= @@ -74,4 +70,4 @@ struct partial_pivot_lud { } }; -} // namespace algebra::cmath::matrix::inverse \ No newline at end of file +} // namespace algebra::cmath::matrix::inverse diff --git a/math/common/include/algebra/math/algorithms/utils/algorithm_finder.hpp b/math/common/include/algebra/math/algorithms/utils/algorithm_finder.hpp new file mode 100644 index 00000000..2a70c9a9 --- /dev/null +++ b/math/common/include/algebra/math/algorithms/utils/algorithm_finder.hpp @@ -0,0 +1,62 @@ +/** 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/matrix/determinant/hard_coded.hpp" +#include "algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp" +#include "algebra/math/algorithms/matrix/inverse/hard_coded.hpp" +#include "algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp" + +namespace algebra::cmath { + +/// Get the type of determinant algorithm acording to matrix dimension +/// @{ +template +struct determinant_selector { + using type = matrix::determinant::partial_pivot_lud; +}; + +template +struct determinant_selector<2, Args...> { + using type = matrix::determinant::hard_coded; +}; + +template +struct determinant_selector<4, Args...> { + using type = matrix::determinant::hard_coded; +}; + +template +using determinant_t = + typename determinant_selector::type; +/// @} + +/// Get the type of inversion algorithm acording to matrix dimension +/// @{ +template +struct inversion_selector { + using type = matrix::inverse::partial_pivot_lud; +}; + +template +struct inversion_selector<2, Args...> { + using type = matrix::inverse::hard_coded; +}; + +template +struct inversion_selector<4, Args...> { + using type = matrix::inverse::hard_coded; +}; + +template +using inversion_t = + typename inversion_selector::type; +/// @} + +} // namespace algebra::cmath diff --git a/math/common/include/algebra/math/common.hpp b/math/common/include/algebra/math/common.hpp index aa1635ed..d1616c24 100644 --- a/math/common/include/algebra/math/common.hpp +++ b/math/common/include/algebra/math/common.hpp @@ -16,6 +16,7 @@ #endif // System include(s). +#include #include namespace algebra::math { @@ -51,4 +52,16 @@ ALGEBRA_HOST_DEVICE inline scalar_t atanh(scalar_t arg) { return math_ns::atanh(arg); } +/// Minimum of two values +template +ALGEBRA_HOST_DEVICE inline scalar_t min(scalar_t a, scalar_t b) { + return math_ns::min(a, b); +} + +/// Maximum of two values +template +ALGEBRA_HOST_DEVICE inline scalar_t max(scalar_t a, scalar_t b) { + return math_ns::max(a, b); +} + } // namespace algebra::math diff --git a/math/eigen/include/algebra/math/impl/eigen_matrix.hpp b/math/eigen/include/algebra/math/impl/eigen_matrix.hpp index c9f6770a..30ee8d67 100644 --- a/math/eigen/include/algebra/math/impl/eigen_matrix.hpp +++ b/math/eigen/include/algebra/math/impl/eigen_matrix.hpp @@ -1,6 +1,6 @@ /** Algebra plugins library, part of the ACTS project * - * (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 */ @@ -27,126 +27,93 @@ #pragma nv_diagnostic pop #endif // __NVCC_DIAG_PRAGMA_SUPPORT__ -namespace algebra::eigen::matrix { - -/// "Matrix actor", assuming an Eigen matrix -template -struct actor { - - /// Size type - using size_ty = int; - - /// Scalar type - using scalar_type = scalar_t; - - /// 2D matrix type - template - using matrix_type = algebra::eigen::matrix_type; - - /// Array type - template - using array_type = storage_type; - - /// 3-element "vector" type - using vector3 = array_type<3>; - - /// Operator getting a reference to one element of a non-const matrix - template ::value && - std::is_convertible::value, - bool> = true> - ALGEBRA_HOST_DEVICE inline scalar_t &element(matrix_type &m, - size_type_1 row, - size_type_2 col) const { - return m(static_cast(row), static_cast(col)); - } - - /// Operator getting one value of a const matrix - template ::value && - std::is_convertible::value, - bool> = true> - ALGEBRA_HOST_DEVICE inline scalar_t element(const matrix_type &m, - size_type_1 row, - size_type_2 col) const { - return m(static_cast(row), static_cast(col)); - } - - /// Operator getting a block of a const matrix - template ::value && - std::is_convertible::value, - bool> = true> - ALGEBRA_HOST_DEVICE matrix_type block(const input_matrix_type &m, - size_type_1 row, - size_type_2 col) const { - return m.template block(static_cast(row), - static_cast(col)); - } - - /// Operator setting a block - template ::value && - std::is_convertible::value, - bool> = true> - ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const matrix_type &b, - size_type_1 row, size_type_2 col) const { - m.template block(static_cast(row), - static_cast(col)) = b; - } - - // Create zero matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type zero() const { - return matrix_type::Zero(); - } - - // Create identity matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type identity() const { - return matrix_type::Identity(); - } - - // Set input matrix as zero matrix - template - ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type &m) const { - m.setZero(); - } - - // Set input matrix as identity matrix - template - ALGEBRA_HOST_DEVICE inline void set_identity( - matrix_type &m) const { - m.setIdentity(); - } - - // Create transpose matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type transpose( - const matrix_type &m) const { - return m.transpose(); - } - - // Get determinant - template - ALGEBRA_HOST_DEVICE inline scalar_t determinant( - const matrix_type &m) const { - return m.determinant(); - } - - // Create inverse matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type inverse( - const matrix_type &m) const { - return m.inverse(); - } -}; - -} // namespace algebra::eigen::matrix +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 +template +ALGEBRA_HOST_DEVICE inline matrix_t zero() { + return matrix_t::Zero(); +} + +// Create identity matrix +template +ALGEBRA_HOST_DEVICE inline matrix_t identity() { + return matrix_t::Identity(); +} + +// 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 +template +ALGEBRA_HOST_DEVICE inline void set_identity( + Eigen::MatrixBase &m) { + m.setIdentity(); +} + +// Create transpose matrix +template +ALGEBRA_HOST_DEVICE inline matrix_type< + typename Eigen::MatrixBase::value_type, + Eigen::MatrixBase::ColsAtCompileTime, + Eigen::MatrixBase::RowsAtCompileTime> +transpose(const Eigen::MatrixBase &m) { + return m.transpose(); +} + +/// @returns the determinant of @param m +template +ALGEBRA_HOST_DEVICE inline typename Eigen::MatrixBase::value_type +determinant(const Eigen::MatrixBase &m) { + return m.determinant(); +} + +/// @returns the inverse of @param m +template +ALGEBRA_HOST_DEVICE inline auto inverse( + const Eigen::MatrixBase &m) { + return m.inverse(); +} + +} // namespace algebra::eigen::math diff --git a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp index fc064fc9..02b3378f 100644 --- a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp +++ b/math/eigen/include/algebra/math/impl/eigen_transform3.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 */ @@ -35,17 +35,27 @@ namespace algebra::eigen::math { /** Transform wrapper class to ensure standard API within differnt plugins */ -template +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; - /// Scalar type used by the transform - using scalar_type = scalar_t; /// 3-element "vector" type using vector3 = array_type<3>; @@ -54,23 +64,6 @@ struct transform3 { /// Point in 2D space using point2 = array_type<2>; - /// 4x4 matrix type - using matrix44 = - typename Eigen::Transform::MatrixType; - - /// Function (object) used for accessing a matrix element - using element_getter = algebra::eigen::math::element_getter; - - /// Size type - using size_type = typename matrix_actor_t::size_ty; - - /// Matrix actor - using matrix_actor = matrix_actor_t; - - /// 2D Matrix type - template - using matrix_type = typename matrix_actor::template matrix_type; - /// @} /// @name Data objects @@ -122,7 +115,7 @@ struct transform3 { * @param t is the transform **/ ALGEBRA_HOST_DEVICE - transform3(const vector3 &t) { + explicit transform3(const vector3 &t) { _data.setIdentity(); auto &matrix = _data.matrix(); @@ -136,7 +129,7 @@ struct transform3 { * @param m is the full 4x4 matrix **/ ALGEBRA_HOST_DEVICE - transform3(const matrix44 &m) { + explicit transform3(const matrix44 &m) { _data.matrix() = m; @@ -148,7 +141,7 @@ struct transform3 { * @param ma is the full 4x4 matrix as a 16 array **/ ALGEBRA_HOST_DEVICE - transform3(const array_type<16> &ma) { + explicit transform3(const array_type<16> &ma) { _data.matrix() << ma[0], ma[1], ma[2], ma[3], ma[4], ma[5], ma[6], ma[7], ma[8], ma[9], ma[10], ma[11], ma[12], ma[13], ma[14], ma[15]; diff --git a/math/fastor/include/algebra/math/impl/fastor_getter.hpp b/math/fastor/include/algebra/math/impl/fastor_getter.hpp index a1854cbc..38d42b5a 100644 --- a/math/fastor/include/algebra/math/impl/fastor_getter.hpp +++ b/math/fastor/include/algebra/math/impl/fastor_getter.hpp @@ -93,7 +93,6 @@ struct element_getter { 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); diff --git a/math/fastor/include/algebra/math/impl/fastor_matrix.hpp b/math/fastor/include/algebra/math/impl/fastor_matrix.hpp index efc5bfa6..b3d421af 100644 --- a/math/fastor/include/algebra/math/impl/fastor_matrix.hpp +++ b/math/fastor/include/algebra/math/impl/fastor_matrix.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 */ @@ -23,147 +23,107 @@ // System include(s). #include // for the std::size_t type -namespace algebra::fastor::matrix { - -/// "Matrix actor", assuming a Fastor matrix -template -struct actor { - - /// Size type - using size_ty = std::size_t; - - /// Scalar type - using scalar_type = scalar_t; - - /// 2D matrix type - template - using matrix_type = algebra::fastor::matrix_type; - - /// Vector type - template - using vector_type = Fastor::Tensor; - - /// Array type - template - using array_type = storage_type; - - /// 3-element "vector" type - using vector3 = array_type<3>; - - /// Operator getting a reference to one element of a non-const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t &element(matrix_type &m, - size_ty row, size_ty col) const { - return m(row, col); - } - - /// Operator getting one value of a const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t element(const matrix_type &m, - size_ty row, size_ty col) const { - return m(row, col); - } - - /// Operator getting a block of a const matrix - template - ALGEBRA_HOST_DEVICE matrix_type block(const input_matrix_type &m, - size_ty row, - size_ty col) const { - // In `Fastor::seq`, the last element is not included. - // `Fastor::seq` takes `int`s as input, but `row`, `col`, `ROWS`, and `COLS` - // have type `size_ty`, which is `std::size_t`. - return m(Fastor::seq(static_cast(row), static_cast(row + ROWS)), - Fastor::seq(static_cast(col), static_cast(col + COLS))); +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 +template +ALGEBRA_HOST_DEVICE inline matrix_t zero() { + return matrix_t(0); +} + +// Create identity matrix +template +ALGEBRA_HOST_DEVICE inline matrix_t identity() { + using scalar_t = algebra::trait::value_t; + constexpr auto rows{algebra::trait::rows}; + constexpr auto cols{algebra::trait::columns}; + + if constexpr (rows >= cols) { + matrix_type identity_matrix; + identity_matrix.eye2(); + return matrix_t( + identity_matrix(Fastor::fseq<0, rows>(), Fastor::fseq<0, cols>())); + } else { + matrix_type identity_matrix; + identity_matrix.eye2(); + return matrix_t( + identity_matrix(Fastor::fseq<0, rows>(), Fastor::fseq<0, cols>())); } +} - /// Operator setting a block with a matrix - template - ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const matrix_type &b, - size_ty row, size_ty col) const { - // In `Fastor::seq`, the last element is not included. - // `Fastor::seq` takes `int`s as input, but `ROWS` and `COLS` have type - // `size_ty`, 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; - } +// Set input matrix as zero matrix +template +ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type &m) { + m.zeros(); +} - /// Operator setting a block with a vector - template - ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m, - const vector_type &b, size_ty row, - size_ty col) const { - // In `Fastor::seq`, the last element is not included. - // `Fastor::seq` takes `int`s as input, but `ROWS` and `COLS` have type - // `size_ty`, which is `std::size_t`. - m(Fastor::seq(static_cast(row), static_cast(row + ROWS)), - static_cast(col)) = b; - } +// Set input matrix as identity matrix +template +ALGEBRA_HOST_DEVICE inline void set_identity( + matrix_type &m) { - // Create zero matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type zero() const { - return matrix_type(0); - } + m = identity>(); +} - // Create identity matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type identity() const { - // There are 2 identity tensor methods in Fastor, eye() and eye2(). The - // former is for arbitrary order tensors, whereas the latter is specifically - // for second order tensors. As such, I chose to use eye2() here because it - // does less and hence would be faster. - - // eye2() only works for square matrices. The idea is to take the largest - // dimension of the matrix, make an identity matrix of that dimension, and - // then return the appropriately sized submatrix of it. - if constexpr (ROWS >= COLS) { - matrix_type identity_matrix; - identity_matrix.eye2(); - return matrix_type( - identity_matrix(Fastor::fseq<0, ROWS>(), Fastor::fseq<0, COLS>())); - } else { - matrix_type identity_matrix; - identity_matrix.eye2(); - return matrix_type( - identity_matrix(Fastor::fseq<0, ROWS>(), Fastor::fseq<0, COLS>())); - } - } +// Create transpose matrix +template +ALGEBRA_HOST_DEVICE inline matrix_type transpose( + const matrix_type &m) { - // Set input matrix as zero matrix - template - ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type &m) const { - m.zeros(); - } + return Fastor::transpose(m); +} - // Set input matrix as identity matrix - template - ALGEBRA_HOST_DEVICE inline void set_identity( - matrix_type &m) const { +/// @returns the determinant of @param m +template +ALGEBRA_HOST_DEVICE inline scalar_t determinant( + const matrix_type &m) { - m = identity(); - } + return Fastor::determinant(m); +} - // Create transpose matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type transpose( - const matrix_type &m) const { - return Fastor::transpose(m); - } +/// @returns the inverse of @param m +template +ALGEBRA_HOST_DEVICE inline matrix_type inverse( + const matrix_type &m) { - // Get determinant - template - ALGEBRA_HOST_DEVICE inline scalar_t determinant( - const matrix_type &m) const { - return Fastor::determinant(m); - } - - // Create inverse matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type inverse( - const matrix_type &m) const { - return Fastor::inverse(m); - } -}; + return Fastor::inverse(m); +} -} // namespace algebra::fastor::matrix +} // namespace algebra::fastor::math diff --git a/math/fastor/include/algebra/math/impl/fastor_transform3.hpp b/math/fastor/include/algebra/math/impl/fastor_transform3.hpp index 224a5485..459a1480 100644 --- a/math/fastor/include/algebra/math/impl/fastor_transform3.hpp +++ b/math/fastor/include/algebra/math/impl/fastor_transform3.hpp @@ -27,7 +27,7 @@ namespace algebra::fastor::math { /** Transform wrapper class to ensure standard API within different plugins * **/ -template +template struct transform3 { /// @name Type definitions for the struct /// @{ @@ -52,14 +52,11 @@ struct transform3 { using element_getter = algebra::fastor::math::element_getter; /// Size type - using size_type = typename matrix_actor_t::size_ty; - - /// Matrix actor - using matrix_actor = matrix_actor_t; + using size_type = std::size_t; /// 2D Matrix type template - using matrix_type = typename matrix_actor::template matrix_type; + using matrix_type = Fastor::Tensor; /// @} @@ -117,7 +114,7 @@ struct transform3 { * @param t is the translation **/ ALGEBRA_HOST - transform3(const vector3 &t) { + explicit transform3(const vector3 &t) { // The matrix needs to be initialized to the identity matrix first. In this // case, the `transform3` requires `_data` to look just like an identity @@ -135,7 +132,7 @@ struct transform3 { * @param m is the full 4x4 matrix **/ ALGEBRA_HOST - transform3(const matrix44 &m) { + explicit transform3(const matrix44 &m) { _data = m; _data_inv = Fastor::inverse(_data); @@ -147,7 +144,7 @@ struct transform3 { * @param ma is the full 4x4 matrix as a 16-element array **/ ALGEBRA_HOST - transform3(const array_type<16> &ma) { + explicit transform3(const array_type<16> &ma) { _data = ma; _data_inv = Fastor::inverse(_data); diff --git a/math/smatrix/CMakeLists.txt b/math/smatrix/CMakeLists.txt index ebd2ec89..0c47132d 100644 --- a/math/smatrix/CMakeLists.txt +++ b/math/smatrix/CMakeLists.txt @@ -16,6 +16,7 @@ algebra_add_library(algebra_smatrix_math smatrix_math "include/algebra/math/impl/smatrix_transform3.hpp" "include/algebra/math/impl/smatrix_vector.hpp") target_link_libraries(algebra_smatrix_math - INTERFACE algebra::common ROOT::Core ROOT::MathCore ROOT::Smatrix) + INTERFACE algebra::common algebra::smatrix_storage ROOT::Core ROOT::MathCore + ROOT::Smatrix) algebra_test_public_headers( algebra_smatrix_math "algebra/math/smatrix.hpp" ) diff --git a/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp b/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp index 34a68e45..a5c9f0f8 100644 --- a/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp +++ b/math/smatrix/include/algebra/math/impl/smatrix_matrix.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,143 +9,109 @@ // Project include(s). #include "algebra/qualifiers.hpp" +#include "algebra/storage/smatrix.hpp" // ROOT/Smatrix include(s). #include -namespace algebra::smatrix::matrix { - -/// "Matrix actor", assuming an SMatrix matrix -template -struct actor { - - /// Size type - using size_ty = unsigned int; - - /// Scalar_type - using scalar_type = scalar_t; - - /// 2D matrix type - template - using matrix_type = ROOT::Math::SMatrix; - - template - using vector_type = ROOT::Math::SVector; - - template - using array_type = vector_type; - - using vector3 = array_type<3>; - - /// Operator getting a reference to one element of a non-const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t &element(matrix_type &m, - unsigned int row, - unsigned int col) const { - return m(row, col); - } - - /// Operator getting one value of a const matrix - template - ALGEBRA_HOST_DEVICE inline scalar_t element(const matrix_type &m, - unsigned int row, - unsigned int col) const { - return m(row, col); - } - - /// Operator getting a block of a const matrix - template - ALGEBRA_HOST_DEVICE matrix_type block(const input_matrix_type &m, - unsigned int row, - unsigned int col) const { - 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) const { - for (unsigned int i = 0; i < ROWS; ++i) { - for (unsigned int j = 0; j < COLS; ++j) { - m(i + row, j + col) = b(i, j); - } +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 vector_type &b, - unsigned int row, unsigned int col) const { - for (unsigned int i = 0; i < ROWS; ++i) { - m(i + row, col) = b[i]; - } - } - - // Create zero matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type zero() const { - return matrix_type(); +} + +/// 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 identity matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type identity() const { - return matrix_type(ROOT::Math::SMatrixIdentity()); +} + +// Create zero matrix +template +ALGEBRA_HOST_DEVICE inline matrix_t zero() { + return matrix_t(); +} + +// Create identity matrix +template +ALGEBRA_HOST_DEVICE inline matrix_t identity() { + return matrix_t(ROOT::Math::SMatrixIdentity()); +} + +// Set input matrix as zero matrix +template +ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type &m) { + + for (unsigned int i = 0; i < ROWS; ++i) { + for (unsigned int j = 0; j < COLS; ++j) { + m(i, j) = 0; + } } - - // Set input matrix as zero matrix - template - ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type &m) const { - - for (unsigned int i = 0; i < ROWS; ++i) { - for (unsigned int j = 0; j < COLS; ++j) { +} + +// Set input matrix as identity matrix +template +ALGEBRA_HOST_DEVICE inline void set_identity( + matrix_type &m) { + + for (unsigned int i = 0; i < ROWS; ++i) { + for (unsigned int j = 0; j < COLS; ++j) { + if (i == j) { + m(i, j) = 1; + } else { m(i, j) = 0; } } } +} - // Set input matrix as identity matrix - template - ALGEBRA_HOST_DEVICE inline void set_identity( - matrix_type &m) const { - - for (unsigned int i = 0; i < ROWS; ++i) { - for (unsigned int j = 0; j < COLS; ++j) { - if (i == j) { - m(i, j) = 1; - } else { - m(i, j) = 0; - } - } - } - } +// Create transpose matrix +template +ALGEBRA_HOST_DEVICE inline matrix_type transpose( + const matrix_type &m) { + return ROOT::Math::Transpose(m); +} - // Create transpose matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type transpose( - const matrix_type &m) const { - return ROOT::Math::Transpose(m); - } +/// @returns the determinant of @param m +template +ALGEBRA_HOST_DEVICE inline scalar_t determinant( + const matrix_type &m) { - // Get determinant - template - ALGEBRA_HOST_DEVICE inline scalar_t determinant( - const matrix_type &m) const { - scalar_t det; - [[maybe_unused]] bool success = m.Det2(det); + scalar_t det; + [[maybe_unused]] bool success = m.Det2(det); - return det; - } + return det; +} - // Create inverse matrix - template - ALGEBRA_HOST_DEVICE inline matrix_type inverse( - const matrix_type &m) const { - int ifail = 0; - return m.Inverse(ifail); - } -}; +/// @returns the inverse of @param m +template +ALGEBRA_HOST_DEVICE inline matrix_type inverse( + const matrix_type &m) { + + int ifail = 0; + return m.Inverse(ifail); +} -} // namespace algebra::smatrix::matrix +} // namespace algebra::smatrix::math diff --git a/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp b/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp index 6db9bd79..9733aa2f 100644 --- a/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp +++ b/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp @@ -20,7 +20,7 @@ namespace algebra::smatrix::math { /** Transform wrapper class to ensure standard API within differnt plugins * **/ -template +template struct transform3 { /// @name Type definitions for the struct @@ -42,18 +42,12 @@ struct transform3 { /// 4x4 matrix type using matrix44 = ROOT::Math::SMatrix; - /// Function (object) used for accessing a matrix element - using element_getter = algebra::smatrix::math::element_getter; - /// Size type - using size_type = typename matrix_actor_t::size_ty; - - /// Matrix actor - using matrix_actor = matrix_actor_t; + using size_type = unsigned int; /// 2D Matrix type template - using matrix_type = typename matrix_actor::template matrix_type; + using matrix_type = ROOT::Math::SMatrix; /// @} @@ -113,7 +107,7 @@ struct transform3 { * @param t is the translation **/ ALGEBRA_HOST - transform3(const vector3 &t) { + explicit transform3(const vector3 &t) { _data(0, 3) = t[0]; _data(1, 3) = t[1]; @@ -129,7 +123,7 @@ struct transform3 { * @param m is the full 4x4 matrix **/ ALGEBRA_HOST - transform3(const matrix44 &m) { + explicit transform3(const matrix44 &m) { _data = m; int ifail = 0; @@ -143,7 +137,7 @@ struct transform3 { * @param ma is the full 4x4 matrix as a 16-element array **/ ALGEBRA_HOST - transform3(const array_type<16> &ma) { + explicit transform3(const array_type<16> &ma) { _data(0, 0) = ma[0]; _data(1, 0) = ma[4]; 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 17e617d8..3abae506 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 @@ -11,117 +11,33 @@ #include "algebra/qualifiers.hpp" #include "algebra/storage/matrix.hpp" -namespace algebra::vc_aos { - -namespace matrix { - -/// Explicitly vectorized matrix implementation -template