From 878388c8f4833276445eb0a3dafb5eda4ff9f644 Mon Sep 17 00:00:00 2001 From: Joana Niermann <53186085+niermann999@users.noreply.github.com> Date: Fri, 6 Dec 2024 15:10:36 +0000 Subject: [PATCH] fix: Eigen Identity matrix and Vc/SMatrix matrix getter and hide arithmetic operators in algebra namespace (#141) The Eigen Identity matrix in the transform class leads to invalid instructions in CUDA. Also fixes small issues in the matrix getters of the Vc and SMatrix plugins. Also moves the arithmetic operators of the cmath plugin to the algebra namespace, since they were picked up in ACTS builds by Eigen types --- frontend/array_cmath/include/algebra/array_cmath.hpp | 4 ++-- frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp | 4 ++-- .../include/algebra/math/impl/eigen_transform3.hpp | 10 ++++++---- .../common/include/algebra/storage/matrix_getter.hpp | 2 +- .../include/algebra/storage/impl/smatrix_getter.hpp | 4 ++-- tests/common/test_device_basics.hpp | 5 +++++ tests/common/test_host_basics.hpp | 2 ++ 7 files changed, 20 insertions(+), 11 deletions(-) diff --git a/frontend/array_cmath/include/algebra/array_cmath.hpp b/frontend/array_cmath/include/algebra/array_cmath.hpp index 581ce4f3..ff4680d9 100644 --- a/frontend/array_cmath/include/algebra/array_cmath.hpp +++ b/frontend/array_cmath/include/algebra/array_cmath.hpp @@ -12,6 +12,8 @@ #include "algebra/math/generic.hpp" #include "algebra/storage/array.hpp" +namespace algebra { + /// @name Operators on @c algebra::array::storage_type /// @{ @@ -21,8 +23,6 @@ using algebra::cmath::operator+; /// @} -namespace algebra { - namespace getter { /// @name Getter functions on @c algebra::array::storage_type diff --git a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp index c3d5a560..914c95f1 100644 --- a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp +++ b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp @@ -12,6 +12,8 @@ #include "algebra/math/generic.hpp" #include "algebra/storage/vecmem.hpp" +namespace algebra { + /// @name Operators on @c algebra::vecmem::storage_type /// @{ @@ -21,8 +23,6 @@ using algebra::cmath::operator+; /// @} -namespace algebra { - namespace getter { /// @name Getter functions on @c algebra::vecmem::matrix_type diff --git a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp index 13c54d85..7aef00f3 100644 --- a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp +++ b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp @@ -65,10 +65,8 @@ struct transform3 { /// @name Data objects /// @{ - Eigen::Transform _data = - Eigen::Transform::Identity(); - Eigen::Transform _data_inv = - Eigen::Transform::Identity(); + Eigen::Transform _data; + Eigen::Transform _data_inv; /// @} @@ -81,6 +79,8 @@ struct transform3 { ALGEBRA_HOST_DEVICE transform3(const vector3 &t, const vector3 &x, const vector3 &y, const vector3 &z, bool get_inverse = true) { + _data.setIdentity(); + auto &matrix = _data.matrix(); matrix.template block<3, 1>(0, 0) = x; matrix.template block<3, 1>(0, 1) = y; @@ -89,6 +89,8 @@ struct transform3 { if (get_inverse) { _data_inv = _data.inverse(); + } else { + _data_inv.setIdentity(); } } diff --git a/storage/common/include/algebra/storage/matrix_getter.hpp b/storage/common/include/algebra/storage/matrix_getter.hpp index 5ac91568..53582880 100644 --- a/storage/common/include/algebra/storage/matrix_getter.hpp +++ b/storage/common/include/algebra/storage/matrix_getter.hpp @@ -241,7 +241,7 @@ ALGEBRA_HOST_DEVICE constexpr void set_block( if constexpr (ROWS == mROW && matrix_t::storage_rows() == input_matrix_t::storage_rows()) { if (row == 0u) { - for (std::size_t j = col; j < mCOL; ++j) { + for (std::size_t j = col; j < col + COLS; ++j) { m[j] = b[j - col]; } return; diff --git a/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp b/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp index 8433779f..e25c1de0 100644 --- a/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp +++ b/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp @@ -145,8 +145,8 @@ struct block_getter { ROOT::Math::SVector ret; - for (std::size_t irow = row; irow < row + SIZE; ++irow) { - ret[irow - row] = m[col][irow]; + for (unsigned int irow = row; irow < row + SIZE; ++irow) { + ret[irow - row] = m(irow, col); } return ret; diff --git a/tests/common/test_device_basics.hpp b/tests/common/test_device_basics.hpp index a4fe602b..561ae346 100644 --- a/tests/common/test_device_basics.hpp +++ b/tests/common/test_device_basics.hpp @@ -40,6 +40,7 @@ class test_device_basics : public test_base { /// Perform various 2D vector operations, and produce a scalar output ALGEBRA_HOST_DEVICE scalar vector_2d_ops(point2 a, point2 b) const { + using namespace algebra; point2 c = a + b; point2 c2 = c * 2.0; @@ -58,6 +59,7 @@ class test_device_basics : public test_base { /// Perform various 3D vector operations, and produce a scalar output ALGEBRA_HOST_DEVICE scalar vector_3d_ops(vector3 a, vector3 b) const { + using namespace algebra; vector3 c = a + b; vector3 c2 = c * 2.0; @@ -78,6 +80,7 @@ class test_device_basics : public test_base { /// Perform some trivial operations on an asymmetrix matrix ALGEBRA_HOST_DEVICE scalar matrix64_ops(const matrix<6, 4>& m) const { + using namespace algebra; matrix<6, 4> m2; for (size_type i = 0; i < 6; ++i) { @@ -139,6 +142,7 @@ class test_device_basics : public test_base { /// Perform some trivial operations on an asymmetrix matrix ALGEBRA_HOST_DEVICE scalar matrix22_ops(const matrix<2, 2>& m22) const { + using namespace algebra; // Test 2 X 2 matrix determinant auto m22_det = algebra::matrix::determinant(m22); @@ -202,6 +206,7 @@ class test_device_basics : public test_base { ALGEBRA_HOST_DEVICE scalar transform3_ops(vector3 t1, vector3 t2, vector3 t3, vector3 a, vector3 b) const { + using namespace algebra; transform3 tr1(t1, t2, t3); transform3 tr2; diff --git a/tests/common/test_host_basics.hpp b/tests/common/test_host_basics.hpp index 1e6dbc3a..cae958ae 100644 --- a/tests/common/test_host_basics.hpp +++ b/tests/common/test_host_basics.hpp @@ -22,6 +22,8 @@ #include #include +using namespace algebra; + /// Test case class, to be specialised for the different plugins - vectors template class test_host_basics_vector : public testing::Test, public test_base {};