Merge pull request #50 from krasznaa/AddMatrixType-main-20220211

Add Matrix Type, main branch (2022.02.11.)
acts-project · Feb 14, 2022 · 9e2a527 · 9e2a527
2 parents 9d0d6d0 + 611927d
commit 9e2a527
Show file tree

Hide file tree

Showing 41 changed files with 459 additions and 242 deletions.
diff --git 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-2021 CERN for the benefit of the ACTS project
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -52,18 +52,24 @@ using cmath::theta;
 
 /// @}
 
-/// Function extracting a slice from the matrix used by
-/// @c algebra::array::transform3
+/// Function extracting a slice from a matrix
 template <std::size_t SIZE, std::size_t ROWS, std::size_t COLS,
           typename scalar_t>
 ALGEBRA_HOST_DEVICE inline array::storage_type<scalar_t, SIZE> vector(
-    const array::storage_type<array::storage_type<scalar_t, ROWS>, COLS>& m,
-    std::size_t row, std::size_t col) {
+    const array::matrix_type<scalar_t, ROWS, COLS>& m, std::size_t row,
+    std::size_t col) {
 
   return cmath::vector_getter<std::size_t, array::storage_type, scalar_t,
                               SIZE>()(m, row, col);
 }
 
+/// @name Getter functions on @c algebra::array::matrix_type
+/// @{
+
+using cmath::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp b/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp
@@ -25,49 +25,14 @@
 namespace algebra {
 namespace eigen {
 
-/// Functor used to access elements of Eigen matrices
-struct element_getter {
-  /// Get non-const access to a matrix element
-  template <
-      typename derived_type,
-      std::enable_if_t<std::is_base_of<Eigen::DenseCoeffsBase<
-                                           derived_type, Eigen::WriteAccessors>,
-                                       Eigen::MatrixBase<derived_type> >::value,
-                       bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto& operator()(
-      Eigen::MatrixBase<derived_type>& m, unsigned int row,
-      unsigned int col) const {
-
-    return m(row, col);
-  }
-  /// Get const access to a matrix element
-  template <typename derived_type>
-  ALGEBRA_HOST_DEVICE inline auto operator()(
-      const Eigen::MatrixBase<derived_type>& m, unsigned int row,
-      unsigned int col) const {
-
-    return m(row, col);
-  }
-};  // struct element_getter
-
-/// Functor used to extract a block from Eigen matrices
-struct block_getter {
-  template <std::size_t kROWS, std::size_t kCOLS, typename matrix_type>
-  ALGEBRA_HOST_DEVICE auto operator()(const matrix_type& m, std::size_t row,
-                                      std::size_t col) const {
-
-    return m.template block<kROWS, kCOLS>(row, col);
-  }
-};  // struct block_getter
-
 /// @name cmath based transforms on @c algebra::eigen::storage_type
 /// @{
 
 template <typename T>
 using transform3 = cmath::transform3<
     std::size_t, eigen::storage_type, T,
     typename Eigen::Transform<T, 3, Eigen::Affine>::MatrixType,
-    algebra::eigen::element_getter, algebra::eigen::block_getter>;
+    math::element_getter, math::block_getter>;
 template <typename T>
 using cartesian2 = cmath::cartesian2<transform3<T> >;
 template <typename T>
@@ -101,6 +66,13 @@ ALGEBRA_HOST_DEVICE inline auto vector(const Eigen::MatrixBase<derived_type>& m,
   return m.template block<SIZE, 1>(row, col);
 }
 
+/// @name Getter functions on @c algebra::eigen::matrix_type
+/// @{
+
+using eigen::math::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp
@@ -60,6 +60,13 @@ ALGEBRA_HOST_DEVICE inline auto vector(const Eigen::MatrixBase<derived_type>& m,
   return m.template block<SIZE, 1>(row, col);
 }
 
+/// @name Getter functions on @c algebra::eigen::matrix_type
+/// @{
+
+using eigen::math::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git a/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp b/frontend/smatrix_cmath/include/algebra/smatrix_cmath.hpp
@@ -1,6 +1,6 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2020-2021 CERN for the benefit of the ACTS project
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -18,52 +18,14 @@
 namespace algebra {
 namespace smatrix {
 
-/// Functor used to access elements of Vc matrices
-template <typename scalar_t>
-struct element_getter {
-
-  template <unsigned int ROWS, unsigned int COLS>
-  using matrix_type = ROOT::Math::SMatrix<scalar_t, ROWS, COLS>;
-
-  template <unsigned int ROWS, unsigned int COLS>
-  ALGEBRA_HOST_DEVICE inline scalar_t& operator()(matrix_type<ROWS, COLS>& m,
-                                                  unsigned int row,
-                                                  unsigned int col) const {
-
-    return m(col, row);
-  }
-
-  template <unsigned int ROWS, unsigned int COLS>
-  ALGEBRA_HOST_DEVICE inline scalar_t operator()(
-      const matrix_type<ROWS, COLS>& m, unsigned int row,
-      unsigned int col) const {
-
-    return m(col, row);
-  }
-};  // element_getter
-
-/// Functor used to extract a block from Vc matrices
-template <typename scalar_t>
-struct block_getter {
-
-  template <unsigned int ROWS, unsigned int COLS>
-  using matrix_type = ROOT::Math::SMatrix<scalar_t, ROWS, COLS>;
-
-  template <unsigned int ROWS, unsigned int COLS, class input_matrix_type>
-  ALGEBRA_HOST_DEVICE matrix_type<ROWS, COLS> operator()(
-      const input_matrix_type& m, unsigned int row, unsigned int col) const {
-
-    return m.template Sub<matrix_type<ROWS, COLS> >(row, col);
-  }
-};  // struct block_getter
-
 /// @name cmath based transforms on @c algebra::smatrix::storage_type
 /// @{
 
 template <typename T>
-using transform3 = cmath::transform3<unsigned int, smatrix::storage_type, T,
-                                     ROOT::Math::SMatrix<T, 4, 4>,
-                                     element_getter<T>, block_getter<T> >;
+using transform3 =
+    cmath::transform3<unsigned int, smatrix::storage_type, T,
+                      ROOT::Math::SMatrix<T, 4, 4>, math::element_getter<T>,
+                      math::block_getter<T> >;
 template <typename T>
 using cartesian2 = cmath::cartesian2<transform3<T> >;
 template <typename T>
@@ -99,6 +61,13 @@ ALGEBRA_HOST_DEVICE inline auto vector(
   return m.template SubCol<smatrix::storage_type<scalar_t, SIZE> >(col, row);
 }
 
+/// @name Getter functions on @c algebra::smatrix::matrix_type
+/// @{
+
+using smatrix::math::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git 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-2021 CERN for the benefit of the ACTS project
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -54,6 +54,13 @@ ALGEBRA_HOST_DEVICE inline auto vector(
   return m.template SubCol<smatrix::storage_type<scalar_t, SIZE> >(col, row);
 }
 
+/// @name Getter functions on @c algebra::smatrix::matrix_type
+/// @{
+
+using smatrix::math::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git a/frontend/vc_cmath/include/algebra/vc_cmath.hpp b/frontend/vc_cmath/include/algebra/vc_cmath.hpp
@@ -24,51 +24,6 @@ using algebra::cmath::operator+;
 namespace algebra {
 namespace vc {
 
-/// Functor used to access elements of Vc matrices
-template <typename scalar_t>
-struct element_getter {
-
-  template <std::size_t ROWS, std::size_t COLS>
-  using matrix_type = Vc::array<Vc::array<scalar_t, ROWS>, COLS>;
-
-  template <std::size_t ROWS, std::size_t COLS>
-  ALGEBRA_HOST_DEVICE inline scalar_t& operator()(matrix_type<ROWS, COLS>& m,
-                                                  std::size_t row,
-                                                  std::size_t col) const {
-
-    return m[col][row];
-  }
-
-  template <std::size_t ROWS, std::size_t COLS>
-  ALGEBRA_HOST_DEVICE inline scalar_t operator()(
-      const matrix_type<ROWS, COLS>& m, std::size_t row,
-      std::size_t col) const {
-
-    return m[col][row];
-  }
-};  // element_getter
-
-/// Functor used to extract a block from Vc matrices
-template <typename scalar_t>
-struct block_getter {
-
-  template <std::size_t ROWS, std::size_t COLS>
-  using matrix_type = Vc::array<Vc::array<scalar_t, ROWS>, COLS>;
-
-  template <std::size_t ROWS, std::size_t COLS, class input_matrix_type>
-  ALGEBRA_HOST_DEVICE matrix_type<ROWS, COLS> operator()(
-      const input_matrix_type& m, std::size_t row, std::size_t col) const {
-
-    matrix_type<ROWS, COLS> submatrix{};
-    for (std::size_t icol = col; icol < col + COLS; ++icol) {
-      for (std::size_t irow = row; irow < row + ROWS; ++irow) {
-        submatrix[icol - col][irow - row] = m[icol][irow];
-      }
-    }
-    return submatrix;
-  }
-};  // struct block_getter
-
 /// @name cmath based transforms on @c algebra::vc types
 /// @{
 
@@ -80,8 +35,10 @@ using math::phi;
 template <typename T>
 using transform3 =
     cmath::transform3<std::size_t, vc::storage_type, T,
-                      Vc::array<Vc::array<T, 4>, 4>, element_getter<T>,
-                      block_getter<T>, vc::vector3<T>, vc::point2<T> >;
+                      Vc::array<Vc::array<T, 4>, 4>,
+                      cmath::element_getter<std::size_t, Vc::array, T>,
+                      cmath::block_getter<std::size_t, Vc::array, T>,
+                      vc::vector3<T>, vc::point2<T> >;
 template <typename T>
 using cartesian2 = cmath::cartesian2<transform3<T> >;
 template <typename T>
@@ -117,13 +74,20 @@ using vc::math::theta;
 template <std::size_t SIZE, std::size_t ROWS, std::size_t COLS,
           typename scalar_t>
 ALGEBRA_HOST_DEVICE inline vc::storage_type<scalar_t, SIZE> vector(
-    const Vc::array<Vc::array<scalar_t, ROWS>, COLS>& m, std::size_t row,
+    const vc::matrix_type<scalar_t, ROWS, COLS>& m, std::size_t row,
     std::size_t col) {
 
   return cmath::vector_getter<std::size_t, Vc::array, scalar_t, SIZE,
                               vc::storage_type<scalar_t, SIZE> >()(m, row, col);
 }
 
+/// @name Getter functions on @c algebra::vc::matrix_type
+/// @{
+
+using cmath::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git a/frontend/vc_vc/include/algebra/vc_vc.hpp b/frontend/vc_vc/include/algebra/vc_vc.hpp
@@ -122,6 +122,13 @@ ALGEBRA_HOST_DEVICE inline auto vector(
   }
 }
 
+/// @name Getter functions on @c algebra::vc::matrix_type
+/// @{
+
+using cmath::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp
@@ -1,6 +1,6 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2020-2021 CERN for the benefit of the ACTS project
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -64,6 +64,13 @@ ALGEBRA_HOST_DEVICE inline vecmem::storage_type<scalar_t, SIZE> vector(
                               SIZE>()(m, row, col);
 }
 
+/// @name Getter functions on @c algebra::vecmem::matrix_type
+/// @{
+
+using cmath::element;
+
+/// @}
+
 }  // namespace getter
 
 namespace vector {

diff --git a/math/cmath/include/algebra/math/impl/cmath_getter.hpp b/math/cmath/include/algebra/math/impl/cmath_getter.hpp
@@ -12,6 +12,7 @@
 #include "algebra/qualifiers.hpp"
 
 // System include(s).
+#include <cassert>
 #include <cstddef>
 #include <type_traits>
 
@@ -104,6 +105,8 @@ struct element_getter {
                                                   std::size_t row,
                                                   std::size_t col) const {
 
+    assert(row < ROWS);
+    assert(col < COLS);
     return m[col][row];
   }
 
@@ -113,10 +116,32 @@ struct element_getter {
       const matrix_type<ROWS, COLS> &m, std::size_t row,
       std::size_t col) const {
 
+    assert(row < ROWS);
+    assert(col < COLS);
     return m[col][row];
   }
 };  // struct element_getter
 
+/// Function extracting an element from a matrix (const)
+template <typename size_type, template <typename, size_type> class array_t,
+          typename scalar_t, size_type ROWS, size_type COLS>
+ALGEBRA_HOST_DEVICE inline scalar_t element(
+    const array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
+    std::size_t col) {
+
+  return element_getter<size_type, array_t, scalar_t>()(m, row, col);
+}
+
+/// Function extracting an element from a matrix (non-const)
+template <typename size_type, template <typename, size_type> class array_t,
+          typename scalar_t, size_type ROWS, size_type COLS>
+ALGEBRA_HOST_DEVICE inline scalar_t &element(
+    array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
+    std::size_t col) {
+
+  return element_getter<size_type, array_t, scalar_t>()(m, row, col);
+}
+
 /// "Vector getter", assuming a simple 2D array access
 template <typename size_type, template <typename, size_type> class array_t,
           typename scalar_t, size_type SIZE,
@@ -134,6 +159,8 @@ struct vector_getter {
   ALGEBRA_HOST_DEVICE inline result_type operator()(
       const matrix_type<ROWS, COLS> &m, std::size_t row, std::size_t col) {
 
+    assert(col < COLS);
+    assert(row + SIZE < ROWS);
     result_type subvector{};
     for (std::size_t irow = row; irow < row + SIZE; ++irow) {
       subvector[irow - row] = m[col][irow];