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

Author: niermann999

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

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

math/eigen/include/algebra/math/impl/eigen_transform3.hpp

@@ -65,10 +65,8 @@ struct transform3 {
   /// @name Data objects
   /// @{
 
-  Eigen::Transform<scalar_type, 3, Eigen::Affine> _data =
-      Eigen::Transform<scalar_type, 3, Eigen::Affine>::Identity();
-  Eigen::Transform<scalar_type, 3, Eigen::Affine> _data_inv =
-      Eigen::Transform<scalar_type, 3, Eigen::Affine>::Identity();
+  Eigen::Transform<scalar_type, 3, Eigen::Affine> _data;
+  Eigen::Transform<scalar_type, 3, Eigen::Affine> _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();
     }
   }
 

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;

storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp

@@ -145,8 +145,8 @@ struct block_getter {
 
     ROOT::Math::SVector<scalar_t, SIZE> 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;

tests/common/test_device_basics.hpp

@@ -40,6 +40,7 @@ class test_device_basics : public test_base<T> {
   /// 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<T> {
   /// 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<T> {
   /// 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<T> {
   /// 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<T> {
   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;

tests/common/test_host_basics.hpp

@@ -22,6 +22,8 @@
 #include <cmath>
 #include <cstddef>
 
+using namespace algebra;
+
 /// Test case class, to be specialised for the different plugins - vectors
 template <typename T>
 class test_host_basics_vector : public testing::Test, public test_base<T> {};