diff --git a/frontend/array_cmath/include/algebra/array_cmath.hpp b/frontend/array_cmath/include/algebra/array_cmath.hpp
index 75d2f4d4..1bc85481 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-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
index 0493b102..e1d53908 100644
--- a/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp
+++ b/frontend/eigen_cmath/include/algebra/eigen_cmath.hpp
@@ -25,41 +25,6 @@
 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
 /// @{
 
@@ -67,7 +32,7 @@ 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
index 1d1d5554..8220c989 100644
--- 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
index 7cc16909..16a7a84e 100644
--- 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
index 2e90be3a..c480b002 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-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
index bd6f1c68..865ce17a 100644
--- 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
index bac282a5..b471b40c 100644
--- 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
index b00bdfe9..5095e1f8 100644
--- a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp
+++ b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp
@@ -1,6 +1,6 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2020-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
index cd9b0463..5200306a 100644
--- 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];
diff --git a/math/eigen/include/algebra/math/impl/eigen_getter.hpp b/math/eigen/include/algebra/math/impl/eigen_getter.hpp
index a01cb1ec..23d89423 100644
--- a/math/eigen/include/algebra/math/impl/eigen_getter.hpp
+++ b/math/eigen/include/algebra/math/impl/eigen_getter.hpp
@@ -93,4 +93,61 @@ ALGEBRA_HOST_DEVICE inline auto eta(
   return algebra::math::atanh(v[2] / v.norm());
 }
 
+/// Functor used to access elements of Eigen matrices
+struct element_getter {
+  /// Get non-const access to a matrix element
+  template <
+      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, std::size_t row,
+      std::size_t 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, std::size_t row,
+      std::size_t col) const {
+
+    return m(row, col);
+  }
+};  // struct element_getter
+
+/// Function extracting an element from a matrix (const)
+template <typename derived_type>
+ALGEBRA_HOST_DEVICE inline auto element(
+    const Eigen::MatrixBase<derived_type> &m, std::size_t row,
+    std::size_t col) {
+
+  return element_getter()(m, row, col);
+}
+
+/// Function extracting an element from a matrix (non-const)
+template <
+    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 &element(Eigen::MatrixBase<derived_type> &m,
+                                         std::size_t row, std::size_t col) {
+
+  return element_getter()(m, row, col);
+}
+
+/// 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
+
 }  // 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 3b885d7f..0a25a719 100644
--- a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp
+++ b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp
@@ -25,36 +25,6 @@
 
 namespace algebra::eigen::math {
 
-namespace internal {
-
-/// Functor used to access elements of Eigen matrices
-template <typename scalar_t>
-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 scalar_t &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 scalar_t operator()(
-      const Eigen::MatrixBase<derived_type> &m, unsigned int row,
-      unsigned int col) const {
-
-    return m(row, col);
-  }
-};  // struct element_getter
-
-}  // namespace internal
-
 /** Transform wrapper class to ensure standard API within differnt plugins */
 template <typename scalar_t>
 struct transform3 {
@@ -80,7 +50,7 @@ struct transform3 {
       typename Eigen::Transform<scalar_type, 3, Eigen::Affine>::MatrixType;
 
   /// Function (object) used for accessing a matrix element
-  using element_getter = internal::element_getter<scalar_type>;
+  using element_getter = algebra::eigen::math::element_getter;
 
   /// @}
 
diff --git a/math/smatrix/include/algebra/math/impl/smatrix_getter.hpp b/math/smatrix/include/algebra/math/impl/smatrix_getter.hpp
index 53997b73..951969a8 100644
--- a/math/smatrix/include/algebra/math/impl/smatrix_getter.hpp
+++ b/math/smatrix/include/algebra/math/impl/smatrix_getter.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
  */
@@ -13,9 +13,13 @@
 // ROOT/Smatrix include(s).
 #include <Math/Expression.h>
 #include <Math/Functions.h>
+#include <Math/SMatrix.h>
 #include <Math/SVector.h>
 #include <TMath.h>
 
+// System include(s).
+#include <cassert>
+
 namespace algebra::smatrix::math {
 
 /** This method retrieves phi from a vector, vector base with rows >= 2
@@ -114,4 +118,65 @@ ALGEBRA_HOST inline scalar_t perp(
   return TMath::Sqrt(v.apply(0) * v.apply(0) + v.apply(1) * v.apply(1));
 }
 
+/// 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 {
+
+    assert(row < ROWS);
+    assert(col < COLS);
+    return m(row, col);
+  }
+
+  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 {
+
+    assert(row < ROWS);
+    assert(col < COLS);
+    return m(row, col);
+  }
+};  // element_getter
+
+/// Function extracting an element from a matrix (const)
+template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
+ALGEBRA_HOST_DEVICE inline scalar_t element(
+    const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, std::size_t row,
+    std::size_t col) {
+
+  return element_getter<scalar_t>()(m, row, col);
+}
+
+/// Function extracting an element from a matrix (non-const)
+template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
+ALGEBRA_HOST_DEVICE inline scalar_t &element(
+    ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, std::size_t row,
+    std::size_t col) {
+
+  return element_getter<scalar_t>()(m, row, col);
+}
+
+/// Functor used to extract a block from SMatrix 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
+
 }  // 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 32ef072a..f16df305 100644
--- a/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp
+++ b/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp
@@ -1,6 +1,6 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2020-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
  */
@@ -17,25 +17,6 @@
 
 namespace algebra::smatrix::math {
 
-namespace internal {
-
-/// Functor used to access elements of SMatrix matrices
-template <typename scalar_t>
-struct element_getter {
-  /// The matrix type used by this struct
-  template <unsigned int ROWS, unsigned int COLS>
-  using matrix_type = ROOT::Math::SMatrix<scalar_t, ROWS, COLS>;
-  /// Get const access to a matrix element
-  template <unsigned int ROWS, unsigned int COLS>
-  ALGEBRA_HOST inline scalar_t operator()(const matrix_type<ROWS, COLS> &m,
-                                          unsigned int row,
-                                          unsigned int col) const {
-    return m(row, col);
-  }
-};  // element_getter
-
-}  // namespace internal
-
 /** Transform wrapper class to ensure standard API within differnt plugins
  *
  **/
@@ -62,7 +43,7 @@ struct transform3 {
   using matrix44 = ROOT::Math::SMatrix<scalar_type, 4, 4>;
 
   /// Function (object) used for accessing a matrix element
-  using element_getter = internal::element_getter<scalar_type>;
+  using element_getter = algebra::smatrix::math::element_getter<scalar_t>;
 
   /// @}
 
diff --git a/storage/array/include/algebra/storage/array.hpp b/storage/array/include/algebra/storage/array.hpp
index 50d66c86..e2a18180 100644
--- a/storage/array/include/algebra/storage/array.hpp
+++ b/storage/array/include/algebra/storage/array.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
  */
@@ -16,6 +16,9 @@ namespace algebra::array {
 /// Array type used in the Array storage model
 template <typename T, std::size_t N>
 using storage_type = std::array<T, N>;
+/// Matrix type used in the Array storage model
+template <typename T, std::size_t ROWS, std::size_t COLS>
+using matrix_type = storage_type<storage_type<T, ROWS>, COLS>;
 
 /// 3-element "vector" type, using @c std::array
 template <typename T>
diff --git a/storage/eigen/include/algebra/storage/eigen.hpp b/storage/eigen/include/algebra/storage/eigen.hpp
index b9922d8e..37e1d35f 100644
--- a/storage/eigen/include/algebra/storage/eigen.hpp
+++ b/storage/eigen/include/algebra/storage/eigen.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,6 +18,9 @@ namespace algebra::eigen {
 /// Array type used in the Eigen storage model
 template <typename T, std::size_t N>
 using storage_type = array<T, N>;
+/// Matrix type used in the Eigen storage model
+template <typename T, std::size_t ROWS, std::size_t COLS>
+using matrix_type = Eigen::Matrix<T, ROWS, COLS>;
 
 /// 3-element "vector" type, using @c algebra::eigen::array
 template <typename T>
diff --git a/storage/smatrix/include/algebra/storage/smatrix.hpp b/storage/smatrix/include/algebra/storage/smatrix.hpp
index 6e7b4e43..7089d0bd 100644
--- a/storage/smatrix/include/algebra/storage/smatrix.hpp
+++ b/storage/smatrix/include/algebra/storage/smatrix.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // ROOT/Smatrix include(s).
+#include <Math/SMatrix.h>
 #include <Math/SVector.h>
 
 // System include(s).
@@ -15,9 +16,12 @@
 
 namespace algebra::smatrix {
 
-/// Array type used in the Smatrix storage model
+/// Array type used in the SMatrix storage model
 template <typename T, unsigned int N>
 using storage_type = ROOT::Math::SVector<T, N>;
+/// Matrix type used in the SMatrix storage model
+template <typename T, unsigned int ROWS, unsigned int COLS>
+using matrix_type = ROOT::Math::SMatrix<T, ROWS, COLS>;
 
 /// 3-element "vector" type, using @c ROOT::Math::SVector
 template <typename T>
diff --git a/storage/vc/include/algebra/storage/vc.hpp b/storage/vc/include/algebra/storage/vc.hpp
index 26f3e25c..5601e0bf 100644
--- a/storage/vc/include/algebra/storage/vc.hpp
+++ b/storage/vc/include/algebra/storage/vc.hpp
@@ -28,6 +28,9 @@ namespace algebra::vc {
 /// Array type used in the Vc storage model
 template <typename T, std::size_t N>
 using storage_type = Vc::SimdArray<T, N>;
+/// Matrix type used in the Vc storage model
+template <typename T, std::size_t ROWS, std::size_t COLS>
+using matrix_type = Vc::array<Vc::array<T, ROWS>, COLS>;
 
 /// 3-element "vector" type, using @c algebra::vc::array4
 template <typename T>
diff --git a/storage/vecmem/include/algebra/storage/vecmem.hpp b/storage/vecmem/include/algebra/storage/vecmem.hpp
index 8fc287f2..4a571920 100644
--- a/storage/vecmem/include/algebra/storage/vecmem.hpp
+++ b/storage/vecmem/include/algebra/storage/vecmem.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,6 +18,9 @@ namespace algebra::vecmem {
 /// Array type used in the VecMem storage model
 template <typename T, std::size_t N>
 using storage_type = ::vecmem::static_array<T, N>;
+/// Matrix type used in the VecMem storage model
+template <typename T, std::size_t ROWS, std::size_t COLS>
+using matrix_type = storage_type<storage_type<T, ROWS>, COLS>;
 
 /// 3-element "vector" type, using @c vecmem::static_array
 template <typename T>
diff --git a/tests/accelerator/common/test_basics_base.hpp b/tests/accelerator/common/test_basics_base.hpp
index 53641f2a..8306d415 100644
--- a/tests/accelerator/common/test_basics_base.hpp
+++ b/tests/accelerator/common/test_basics_base.hpp
@@ -51,6 +51,10 @@ class test_basics_base : public testing::Test, public test_base<T> {
     m_v2 = std::make_unique<vecmem::vector<typename T::vector3> >(s_arraySize,
                                                                   &m_resource);
 
+    m_m1 =
+        std::make_unique<vecmem::vector<typename T::template matrix<6, 4> > >(
+            s_arraySize, &m_resource);
+
     m_output_host = std::make_unique<vecmem::vector<typename T::scalar> >(
         s_arraySize, &m_resource);
     m_output_device = std::make_unique<vecmem::vector<typename T::scalar> >(
@@ -81,6 +85,13 @@ class test_basics_base : public testing::Test, public test_base<T> {
                      static_cast<typename T::scalar>(i * 2.3),
                      static_cast<typename T::scalar>((i + 2) * 3.4)};
 
+      for (typename T::size_type j = 0; j < 6; ++j) {
+        for (typename T::size_type k = 0; k < 4; ++k) {
+          algebra::getter::element(m_m1->at(i), j, k) =
+              static_cast<typename T::scalar>(j * 20.3 + k * 10.5);
+        }
+      }
+
       m_output_host->at(i) = 0;
       m_output_device->at(i) = 0;
     }
@@ -122,6 +133,7 @@ class test_basics_base : public testing::Test, public test_base<T> {
   std::unique_ptr<vecmem::vector<typename T::vector3> > m_t1, m_t2, m_t3;
   std::unique_ptr<vecmem::vector<typename T::point2> > m_p1, m_p2;
   std::unique_ptr<vecmem::vector<typename T::vector3> > m_v1, m_v2;
+  std::unique_ptr<vecmem::vector<typename T::template matrix<6, 4> > > m_m1;
 
   /// @}
 
diff --git a/tests/accelerator/common/test_basics_functors.hpp b/tests/accelerator/common/test_basics_functors.hpp
index 62462715..96d3518a 100644
--- a/tests/accelerator/common/test_basics_functors.hpp
+++ b/tests/accelerator/common/test_basics_functors.hpp
@@ -66,6 +66,26 @@ class vector_3d_ops_functor : public functor_base<T> {
   }
 };
 
+/// Functor running @c test_device_basics::matrix64_ops
+template <typename T>
+class matrix64_ops_functor : public functor_base<T> {
+
+ public:
+  ALGEBRA_HOST_DEVICE void operator()(
+      std::size_t i,
+      vecmem::data::vector_view<const typename T::template matrix<6, 4> > m,
+      vecmem::data::vector_view<typename T::scalar> output) const {
+
+    // Create the VecMem vector(s).
+    vecmem::device_vector<const typename T::template matrix<6, 4> > vec_m(m);
+    vecmem::device_vector<typename T::scalar> vec_output(output);
+
+    // Perform the operation.
+    auto ii = static_cast<typename decltype(vec_output)::size_type>(i);
+    vec_output[ii] = this->m_tester.matrix64_ops(vec_m[ii]);
+  }
+};
+
 /// Functor running @c test_device_basics::transform3_ops
 template <typename T>
 class transform3_ops_functor : public functor_base<T> {
diff --git a/tests/accelerator/cuda/CMakeLists.txt b/tests/accelerator/cuda/CMakeLists.txt
index 39512c6a..e9f9be18 100644
--- a/tests/accelerator/cuda/CMakeLists.txt
+++ b/tests/accelerator/cuda/CMakeLists.txt
@@ -1,6 +1,6 @@
 # Algebra plugins library, part of the ACTS project (R&D line)
 #
-# (c) 2021 CERN for the benefit of the ACTS project
+# (c) 2021-2022 CERN for the benefit of the ACTS project
 #
 # Mozilla Public License Version 2.0
 
@@ -28,6 +28,15 @@ target_compile_options( algebra_tests_cuda_common INTERFACE
    $<$<COMPILE_LANGUAGE:CUDA>:-Xcudafe --diag_suppress=177>
    $<$<COMPILE_LANGUAGE:CUDA>:--ftz=false --prec-div=true>
    $<$<COMPILE_LANGUAGE:CUDA>:--prec-sqrt=true --fmad=false> )
+if( "${CUDAToolkit_VERSION}" VERSION_GREATER_EQUAL "11.5" )
+   # Replace all "#pragma diag_suppress" calls in the Eigen code with
+   # "#pragma nv_diag_suppress" calls. To be removed once the Eigen
+   # code starts behaving better. Note that this is only needed for
+   # Windows, as these warnings are suppressed on Linux by declaring
+   # Eigen as a "system include".
+   target_compile_definitions( algebra_tests_cuda_common
+      INTERFACE diag_suppress=nv_diag_suppress )
+endif()
 add_library( algebra::tests_cuda_common ALIAS algebra_tests_cuda_common )
 
 # Set up all of the (available) CUDA tests.
diff --git a/tests/accelerator/cuda/array_cmath.cu b/tests/accelerator/cuda/array_cmath.cu
index 9de77a7a..c06141da 100644
--- a/tests/accelerator/cuda/array_cmath.cu
+++ b/tests/accelerator/cuda/array_cmath.cu
@@ -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
  */
@@ -38,12 +38,14 @@ typedef testing::Types<
         float, algebra::array::point2<float>, algebra::array::point3<float>,
         algebra::array::vector2<float>, algebra::array::vector3<float>,
         algebra::array::transform3<float>, algebra::array::cartesian2<float>,
-        algebra::array::polar2<float>, algebra::array::cylindrical2<float> >,
+        algebra::array::polar2<float>, algebra::array::cylindrical2<float>,
+        std::size_t, algebra::array::matrix_type>,
     test_types<
         double, algebra::array::point2<double>, algebra::array::point3<double>,
         algebra::array::vector2<double>, algebra::array::vector3<double>,
         algebra::array::transform3<double>, algebra::array::cartesian2<double>,
-        algebra::array::polar2<double>, algebra::array::cylindrical2<double> > >
+        algebra::array::polar2<double>, algebra::array::cylindrical2<double>,
+        std::size_t, algebra::array::matrix_type> >
     array_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_cuda_basics,
                                array_cmath_types, test_specialisation_name);
diff --git a/tests/accelerator/cuda/common/test_cuda_basics.cuh b/tests/accelerator/cuda/common/test_cuda_basics.cuh
index d969b17b..556e54c6 100644
--- a/tests/accelerator/cuda/common/test_cuda_basics.cuh
+++ b/tests/accelerator/cuda/common/test_cuda_basics.cuh
@@ -76,6 +76,21 @@ TYPED_TEST_P(test_cuda_basics, vector_3d_ops) {
   this->compareOutputs();
 }
 
+/// Test for handling matrices
+TYPED_TEST_P(test_cuda_basics, matrix64_ops) {
+
+  // Run the test on the host, and on the/a device.
+  execute_host_test<matrix64_ops_functor<TypeParam> >(
+      this->m_m1->size(), vecmem::get_data(*(this->m_m1)),
+      vecmem::get_data(*(this->m_output_host)));
+  execute_cuda_test<matrix64_ops_functor<TypeParam> >(
+      this->m_m1->size(), vecmem::get_data(*(this->m_m1)),
+      vecmem::get_data(*(this->m_output_device)));
+
+  // Compare the outputs.
+  this->compareOutputs();
+}
+
 /// Test for some operations with @c transform3
 TYPED_TEST_P(test_cuda_basics, transform3) {
 
@@ -153,4 +168,5 @@ TYPED_TEST_P(test_cuda_basics, polar2) {
 }
 
 REGISTER_TYPED_TEST_SUITE_P(test_cuda_basics, vector_2d_ops, vector_3d_ops,
-                            transform3, cartesian2, cylindrical2, polar2);
+                            matrix64_ops, transform3, cartesian2, cylindrical2,
+                            polar2);
diff --git a/tests/accelerator/cuda/eigen_cmath.cu b/tests/accelerator/cuda/eigen_cmath.cu
index 48ae6a7f..83d1e19b 100644
--- a/tests/accelerator/cuda/eigen_cmath.cu
+++ b/tests/accelerator/cuda/eigen_cmath.cu
@@ -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
  */
@@ -38,12 +38,14 @@ typedef testing::Types<
         float, algebra::eigen::point2<float>, algebra::eigen::point3<float>,
         algebra::eigen::vector2<float>, algebra::eigen::vector3<float>,
         algebra::eigen::transform3<float>, algebra::eigen::cartesian2<float>,
-        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float> >,
+        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float>,
+        std::size_t, algebra::eigen::matrix_type>,
     test_types<
         double, algebra::eigen::point2<double>, algebra::eigen::point3<double>,
         algebra::eigen::vector2<double>, algebra::eigen::vector3<double>,
         algebra::eigen::transform3<double>, algebra::eigen::cartesian2<double>,
-        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double> > >
+        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double>,
+        std::size_t, algebra::eigen::matrix_type> >
     eigen_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_cuda_basics,
                                eigen_cmath_types, test_specialisation_name);
diff --git a/tests/accelerator/cuda/eigen_eigen.cu b/tests/accelerator/cuda/eigen_eigen.cu
index 1ae9eaf6..93caacc7 100644
--- a/tests/accelerator/cuda/eigen_eigen.cu
+++ b/tests/accelerator/cuda/eigen_eigen.cu
@@ -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
  */
@@ -38,12 +38,14 @@ typedef testing::Types<
         float, algebra::eigen::point2<float>, algebra::eigen::point3<float>,
         algebra::eigen::vector2<float>, algebra::eigen::vector3<float>,
         algebra::eigen::transform3<float>, algebra::eigen::cartesian2<float>,
-        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float> >,
+        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float>,
+        std::size_t, algebra::eigen::matrix_type>,
     test_types<
         double, algebra::eigen::point2<double>, algebra::eigen::point3<double>,
         algebra::eigen::vector2<double>, algebra::eigen::vector3<double>,
         algebra::eigen::transform3<double>, algebra::eigen::cartesian2<double>,
-        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double> > >
+        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double>,
+        std::size_t, algebra::eigen::matrix_type> >
     eigen_eigen_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_cuda_basics,
                                eigen_eigen_types, test_specialisation_name);
diff --git a/tests/accelerator/cuda/vecmem_cmath.cu b/tests/accelerator/cuda/vecmem_cmath.cu
index 73f80bc2..7179f1cf 100644
--- a/tests/accelerator/cuda/vecmem_cmath.cu
+++ b/tests/accelerator/cuda/vecmem_cmath.cu
@@ -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
  */
@@ -38,13 +38,15 @@ typedef testing::Types<
         float, algebra::vecmem::point2<float>, algebra::vecmem::point3<float>,
         algebra::vecmem::vector2<float>, algebra::vecmem::vector3<float>,
         algebra::vecmem::transform3<float>, algebra::vecmem::cartesian2<float>,
-        algebra::vecmem::polar2<float>, algebra::vecmem::cylindrical2<float> >,
+        algebra::vecmem::polar2<float>, algebra::vecmem::cylindrical2<float>,
+        std::size_t, algebra::vecmem::matrix_type>,
     test_types<
         double, algebra::vecmem::point2<double>,
         algebra::vecmem::point3<double>, algebra::vecmem::vector2<double>,
         algebra::vecmem::vector3<double>, algebra::vecmem::transform3<double>,
         algebra::vecmem::cartesian2<double>, algebra::vecmem::polar2<double>,
-        algebra::vecmem::cylindrical2<double> > >
+        algebra::vecmem::cylindrical2<double>, std::size_t,
+        algebra::vecmem::matrix_type> >
     vecmem_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_cuda_basics,
                                vecmem_cmath_types, test_specialisation_name);
diff --git a/tests/accelerator/sycl/array_cmath.sycl b/tests/accelerator/sycl/array_cmath.sycl
index 21ecaddd..112c2f40 100644
--- a/tests/accelerator/sycl/array_cmath.sycl
+++ b/tests/accelerator/sycl/array_cmath.sycl
@@ -38,12 +38,14 @@ typedef testing::Types<
         float, algebra::array::point2<float>, algebra::array::point3<float>,
         algebra::array::vector2<float>, algebra::array::vector3<float>,
         algebra::array::transform3<float>, algebra::array::cartesian2<float>,
-        algebra::array::polar2<float>, algebra::array::cylindrical2<float> >,
+        algebra::array::polar2<float>, algebra::array::cylindrical2<float>,
+        std::size_t, algebra::array::matrix_type>,
     test_types<
         double, algebra::array::point2<double>, algebra::array::point3<double>,
         algebra::array::vector2<double>, algebra::array::vector3<double>,
         algebra::array::transform3<double>, algebra::array::cartesian2<double>,
-        algebra::array::polar2<double>, algebra::array::cylindrical2<double> > >
+        algebra::array::polar2<double>, algebra::array::cylindrical2<double>,
+        std::size_t, algebra::array::matrix_type> >
     array_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_sycl_basics,
                                array_cmath_types, test_specialisation_name);
diff --git a/tests/accelerator/sycl/common/test_sycl_basics.hpp b/tests/accelerator/sycl/common/test_sycl_basics.hpp
index a7ffb718..9388f039 100644
--- a/tests/accelerator/sycl/common/test_sycl_basics.hpp
+++ b/tests/accelerator/sycl/common/test_sycl_basics.hpp
@@ -81,6 +81,21 @@ TYPED_TEST_P(test_sycl_basics, vector_3d_ops) {
   this->compareOutputs();
 }
 
+/// Test for handling matrices
+TYPED_TEST_P(test_sycl_basics, matrix64_ops) {
+
+  // Run the test on the host, and on the/a device.
+  execute_host_test<matrix64_ops_functor<TypeParam> >(
+      this->m_m1->size(), vecmem::get_data(*(this->m_m1)),
+      vecmem::get_data(*(this->m_output_host)));
+  execute_sycl_test<matrix64_ops_functor<TypeParam> >(
+      this->m_queue, this->m_m1->size(), vecmem::get_data(*(this->m_m1)),
+      vecmem::get_data(*(this->m_output_device)));
+
+  // Compare the outputs.
+  this->compareOutputs();
+}
+
 /// Test for some operations with @c transform3
 TYPED_TEST_P(test_sycl_basics, transform3) {
 
@@ -158,4 +173,5 @@ TYPED_TEST_P(test_sycl_basics, polar2) {
 }
 
 REGISTER_TYPED_TEST_SUITE_P(test_sycl_basics, vector_2d_ops, vector_3d_ops,
-                            transform3, cartesian2, cylindrical2, polar2);
+                            matrix64_ops, transform3, cartesian2, cylindrical2,
+                            polar2);
diff --git a/tests/accelerator/sycl/vecmem_cmath.sycl b/tests/accelerator/sycl/vecmem_cmath.sycl
index 1a173c4b..a211aa71 100644
--- a/tests/accelerator/sycl/vecmem_cmath.sycl
+++ b/tests/accelerator/sycl/vecmem_cmath.sycl
@@ -38,13 +38,15 @@ typedef testing::Types<
         float, algebra::vecmem::point2<float>, algebra::vecmem::point3<float>,
         algebra::vecmem::vector2<float>, algebra::vecmem::vector3<float>,
         algebra::vecmem::transform3<float>, algebra::vecmem::cartesian2<float>,
-        algebra::vecmem::polar2<float>, algebra::vecmem::cylindrical2<float> >,
+        algebra::vecmem::polar2<float>, algebra::vecmem::cylindrical2<float>,
+        std::size_t, algebra::vecmem::matrix_type>,
     test_types<
         double, algebra::vecmem::point2<double>,
         algebra::vecmem::point3<double>, algebra::vecmem::vector2<double>,
         algebra::vecmem::vector3<double>, algebra::vecmem::transform3<double>,
         algebra::vecmem::cartesian2<double>, algebra::vecmem::polar2<double>,
-        algebra::vecmem::cylindrical2<double> > >
+        algebra::vecmem::cylindrical2<double>, std::size_t,
+        algebra::vecmem::matrix_type> >
     vecmem_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_sycl_basics,
                                vecmem_cmath_types, test_specialisation_name);
diff --git a/tests/array/array_cmath.cpp b/tests/array/array_cmath.cpp
index b83e0870..4d43d256 100644
--- a/tests/array/array_cmath.cpp
+++ b/tests/array/array_cmath.cpp
@@ -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
  */
@@ -38,12 +38,14 @@ typedef testing::Types<
         float, algebra::array::point2<float>, algebra::array::point3<float>,
         algebra::array::vector2<float>, algebra::array::vector3<float>,
         algebra::array::transform3<float>, algebra::array::cartesian2<float>,
-        algebra::array::polar2<float>, algebra::array::cylindrical2<float> >,
+        algebra::array::polar2<float>, algebra::array::cylindrical2<float>,
+        std::size_t, algebra::array::matrix_type>,
     test_types<
         double, algebra::array::point2<double>, algebra::array::point3<double>,
         algebra::array::vector2<double>, algebra::array::vector3<double>,
         algebra::array::transform3<double>, algebra::array::cartesian2<double>,
-        algebra::array::polar2<double>, algebra::array::cylindrical2<double> > >
+        algebra::array::polar2<double>, algebra::array::cylindrical2<double>,
+        std::size_t, algebra::array::matrix_type> >
     array_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics,
                                array_cmath_types, test_specialisation_name);
diff --git a/tests/common/test_base.hpp b/tests/common/test_base.hpp
index e1357200..c88d2bea 100644
--- a/tests/common/test_base.hpp
+++ b/tests/common/test_base.hpp
@@ -20,10 +20,12 @@ class test_base {};
 /// Test base class, using a @c test_types type argument
 template <typename scalar_t, typename point2_t, typename point3_t,
           typename vector2_t, typename vector3_t, typename transform3_t,
-          typename cartesian2_t, typename polar2_t, typename cylindrical2_t>
+          typename cartesian2_t, typename polar2_t, typename cylindrical2_t,
+          typename size_ty,
+          template <typename, size_ty, size_ty> class matrix_t>
 class test_base<
     test_types<scalar_t, point2_t, point3_t, vector2_t, vector3_t, transform3_t,
-               cartesian2_t, polar2_t, cylindrical2_t> > {
+               cartesian2_t, polar2_t, cylindrical2_t, size_ty, matrix_t> > {
 
  public:
   /// @name Type definitions
@@ -38,6 +40,9 @@ class test_base<
   using cartesian2 = cartesian2_t;
   using polar2 = polar2_t;
   using cylindrical2 = cylindrical2_t;
+  using size_type = size_ty;
+  template <size_type ROWS, size_type COLS>
+  using matrix = matrix_t<scalar, ROWS, COLS>;
 
   /// @}
 
diff --git a/tests/common/test_device_basics.hpp b/tests/common/test_device_basics.hpp
index d8d05be1..1d90d9e7 100644
--- a/tests/common/test_device_basics.hpp
+++ b/tests/common/test_device_basics.hpp
@@ -34,6 +34,9 @@ class test_device_basics : public test_base<T> {
   using cartesian2 = typename test_base<T>::cartesian2;
   using polar2 = typename test_base<T>::polar2;
   using cylindrical2 = typename test_base<T>::cylindrical2;
+  using size_type = typename test_base<T>::size_type;
+  template <size_type ROWS, size_type COLS>
+  using matrix = typename test_base<T>::template matrix<ROWS, COLS>;
 
   /// @}
 
@@ -75,6 +78,27 @@ class test_device_basics : public test_base<T> {
     return (phi + perp + norm1 + dot + norm3);
   }
 
+  /// Perform some trivial operations on an asymmetrix matrix
+  ALGEBRA_HOST_DEVICE
+  scalar matrix64_ops(const matrix<6, 4>& m) const {
+
+    matrix<6, 4> m2;
+    for (size_type i = 0; i < 6; ++i) {
+      for (size_type j = 0; j < 4; ++j) {
+        algebra::getter::element(m2, i, j) = algebra::getter::element(m, i, j);
+      }
+    }
+
+    scalar result = 0.;
+    for (size_type i = 0; i < 6; ++i) {
+      for (size_type j = 0; j < 4; ++j) {
+        result += 0.6f * algebra::getter::element(m, i, j) +
+                  0.7f * algebra::getter::element(m2, i, j);
+      }
+    }
+    return result;
+  }
+
   /// Perform various operations using the @c transform3 type
   ALGEBRA_HOST_DEVICE
   scalar transform3_ops(vector3 t1, vector3 t2, vector3 t3, vector3 a,
diff --git a/tests/common/test_host_basics.hpp b/tests/common/test_host_basics.hpp
index da493cae..9ac9776b 100644
--- a/tests/common/test_host_basics.hpp
+++ b/tests/common/test_host_basics.hpp
@@ -107,6 +107,35 @@ TYPED_TEST_P(test_host_basics, vector3) {
   ASSERT_NEAR(norm, std::sqrt(3.), this->m_epsilon);
 }
 
+// Test generic access to a 6x4 matrix
+TYPED_TEST_P(test_host_basics, matrix64) {
+
+  // Create the matrix.
+  static constexpr typename TypeParam::size_type ROWS = 6;
+  static constexpr typename TypeParam::size_type COLS = 4;
+  typename TypeParam::template matrix<ROWS, COLS> m;
+
+  // Fill it.
+  for (typename TypeParam::size_type i = 0; i < ROWS; ++i) {
+    for (typename TypeParam::size_type j = 0; j < COLS; ++j) {
+      algebra::getter::element(m, i, j) =
+          static_cast<typename TypeParam::scalar>(0.5 * i + j);
+    }
+  }
+
+  // Check its content.
+  const typename TypeParam::template matrix<ROWS, COLS>& m_const_ref = m;
+  for (typename TypeParam::size_type i = 0; i < ROWS; ++i) {
+    for (typename TypeParam::size_type j = 0; j < COLS; ++j) {
+      const typename TypeParam::scalar ref =
+          static_cast<typename TypeParam::scalar>(0.5 * i + j);
+      ASSERT_NEAR(algebra::getter::element(m, i, j), ref, this->m_epsilon);
+      ASSERT_NEAR(algebra::getter::element(m_const_ref, i, j), ref,
+                  this->m_epsilon);
+    }
+  }
+}
+
 // This defines the vector operation test suite
 TYPED_TEST_P(test_host_basics, getter) {
 
@@ -293,6 +322,6 @@ TYPED_TEST_P(test_host_basics, local_transformations) {
   ASSERT_NEAR(polfrom2[1], polfrom3[1], this->m_epsilon);
 }
 
-REGISTER_TYPED_TEST_SUITE_P(test_host_basics, local_vectors, vector3, getter,
-                            transform3, global_transformations,
+REGISTER_TYPED_TEST_SUITE_P(test_host_basics, local_vectors, vector3, matrix64,
+                            getter, transform3, global_transformations,
                             local_transformations);
diff --git a/tests/common/test_types.hpp b/tests/common/test_types.hpp
index 66fe9f0b..2e51599f 100644
--- a/tests/common/test_types.hpp
+++ b/tests/common/test_types.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
  */
@@ -10,7 +10,9 @@
 /// Simple struct holding the types that describe a given plugin
 template <typename scalar_t, typename point2_t, typename point3_t,
           typename vector2_t, typename vector3_t, typename transform3_t,
-          typename cartesian2_t, typename polar2_t, typename cylindrical2_t>
+          typename cartesian2_t, typename polar2_t, typename cylindrical2_t,
+          typename size_ty,
+          template <typename, size_ty, size_ty> class matrix_t>
 struct test_types {
 
   using scalar = scalar_t;
@@ -22,5 +24,8 @@ struct test_types {
   using cartesian2 = cartesian2_t;
   using polar2 = polar2_t;
   using cylindrical2 = cylindrical2_t;
+  using size_type = size_ty;
+  template <size_type ROWS, size_type COLS>
+  using matrix = matrix_t<scalar, ROWS, COLS>;
 
 };  // struct test_types
diff --git a/tests/eigen/eigen_cmath.cpp b/tests/eigen/eigen_cmath.cpp
index 67bc1534..bf8bd4aa 100644
--- a/tests/eigen/eigen_cmath.cpp
+++ b/tests/eigen/eigen_cmath.cpp
@@ -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
  */
@@ -38,12 +38,14 @@ typedef testing::Types<
         float, algebra::eigen::point2<float>, algebra::eigen::point3<float>,
         algebra::eigen::vector2<float>, algebra::eigen::vector3<float>,
         algebra::eigen::transform3<float>, algebra::eigen::cartesian2<float>,
-        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float> >,
+        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float>,
+        std::size_t, algebra::eigen::matrix_type>,
     test_types<
         double, algebra::eigen::point2<double>, algebra::eigen::point3<double>,
         algebra::eigen::vector2<double>, algebra::eigen::vector3<double>,
         algebra::eigen::transform3<double>, algebra::eigen::cartesian2<double>,
-        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double> > >
+        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double>,
+        std::size_t, algebra::eigen::matrix_type> >
     eigen_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics,
                                eigen_cmath_types, test_specialisation_name);
diff --git a/tests/eigen/eigen_eigen.cpp b/tests/eigen/eigen_eigen.cpp
index 74a2997f..ddee1a09 100644
--- a/tests/eigen/eigen_eigen.cpp
+++ b/tests/eigen/eigen_eigen.cpp
@@ -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
  */
@@ -38,12 +38,14 @@ typedef testing::Types<
         float, algebra::eigen::point2<float>, algebra::eigen::point3<float>,
         algebra::eigen::vector2<float>, algebra::eigen::vector3<float>,
         algebra::eigen::transform3<float>, algebra::eigen::cartesian2<float>,
-        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float> >,
+        algebra::eigen::polar2<float>, algebra::eigen::cylindrical2<float>,
+        std::size_t, algebra::eigen::matrix_type>,
     test_types<
         double, algebra::eigen::point2<double>, algebra::eigen::point3<double>,
         algebra::eigen::vector2<double>, algebra::eigen::vector3<double>,
         algebra::eigen::transform3<double>, algebra::eigen::cartesian2<double>,
-        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double> > >
+        algebra::eigen::polar2<double>, algebra::eigen::cylindrical2<double>,
+        std::size_t, algebra::eigen::matrix_type> >
     eigen_eigen_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics,
                                eigen_eigen_types, test_specialisation_name);
diff --git a/tests/smatrix/smatrix_cmath.cpp b/tests/smatrix/smatrix_cmath.cpp
index 325d7fa6..70ea8934 100644
--- a/tests/smatrix/smatrix_cmath.cpp
+++ b/tests/smatrix/smatrix_cmath.cpp
@@ -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
  */
@@ -39,13 +39,15 @@ typedef testing::Types<
         algebra::smatrix::vector2<float>, algebra::smatrix::vector3<float>,
         algebra::smatrix::transform3<float>,
         algebra::smatrix::cartesian2<float>, algebra::smatrix::polar2<float>,
-        algebra::smatrix::cylindrical2<float> >,
+        algebra::smatrix::cylindrical2<float>, unsigned int,
+        algebra::smatrix::matrix_type>,
     test_types<
         double, algebra::smatrix::point2<double>,
         algebra::smatrix::point3<double>, algebra::smatrix::vector2<double>,
         algebra::smatrix::vector3<double>, algebra::smatrix::transform3<double>,
         algebra::smatrix::cartesian2<double>, algebra::smatrix::polar2<double>,
-        algebra::smatrix::cylindrical2<double> > >
+        algebra::smatrix::cylindrical2<double>, unsigned int,
+        algebra::smatrix::matrix_type> >
     smatrix_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics,
                                smatrix_cmath_types, test_specialisation_name);
diff --git a/tests/smatrix/smatrix_smatrix.cpp b/tests/smatrix/smatrix_smatrix.cpp
index 6161609c..d48fa65f 100644
--- a/tests/smatrix/smatrix_smatrix.cpp
+++ b/tests/smatrix/smatrix_smatrix.cpp
@@ -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
  */
@@ -39,13 +39,15 @@ typedef testing::Types<
         algebra::smatrix::vector2<float>, algebra::smatrix::vector3<float>,
         algebra::smatrix::transform3<float>,
         algebra::smatrix::cartesian2<float>, algebra::smatrix::polar2<float>,
-        algebra::smatrix::cylindrical2<float> >,
+        algebra::smatrix::cylindrical2<float>, unsigned int,
+        algebra::smatrix::matrix_type>,
     test_types<
         double, algebra::smatrix::point2<double>,
         algebra::smatrix::point3<double>, algebra::smatrix::vector2<double>,
         algebra::smatrix::vector3<double>, algebra::smatrix::transform3<double>,
         algebra::smatrix::cartesian2<double>, algebra::smatrix::polar2<double>,
-        algebra::smatrix::cylindrical2<double> > >
+        algebra::smatrix::cylindrical2<double>, unsigned int,
+        algebra::smatrix::matrix_type> >
     smatrix_smatrix_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics,
                                smatrix_smatrix_types, test_specialisation_name);
diff --git a/tests/vc/vc_cmath.cpp b/tests/vc/vc_cmath.cpp
index aa1e69b1..7490eeb6 100644
--- a/tests/vc/vc_cmath.cpp
+++ b/tests/vc/vc_cmath.cpp
@@ -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
  */
@@ -37,12 +37,13 @@ typedef testing::Types<
     test_types<float, algebra::vc::point2<float>, algebra::vc::point3<float>,
                algebra::vc::vector2<float>, algebra::vc::vector3<float>,
                algebra::vc::transform3<float>, algebra::vc::cartesian2<float>,
-               algebra::vc::polar2<float>, algebra::vc::cylindrical2<float> >,
+               algebra::vc::polar2<float>, algebra::vc::cylindrical2<float>,
+               std::size_t, algebra::vc::matrix_type>,
     test_types<double, algebra::vc::point2<double>, algebra::vc::point3<double>,
                algebra::vc::vector2<double>, algebra::vc::vector3<double>,
                algebra::vc::transform3<double>, algebra::vc::cartesian2<double>,
-               algebra::vc::polar2<double>,
-               algebra::vc::cylindrical2<double> > >
+               algebra::vc::polar2<double>, algebra::vc::cylindrical2<double>,
+               std::size_t, algebra::vc::matrix_type> >
     vc_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics,
                                vc_cmath_types, test_specialisation_name);
diff --git a/tests/vc/vc_vc.cpp b/tests/vc/vc_vc.cpp
index 11e2ff2c..50ad2443 100644
--- a/tests/vc/vc_vc.cpp
+++ b/tests/vc/vc_vc.cpp
@@ -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
  */
@@ -37,12 +37,13 @@ typedef testing::Types<
     test_types<float, algebra::vc::point2<float>, algebra::vc::point3<float>,
                algebra::vc::vector2<float>, algebra::vc::vector3<float>,
                algebra::vc::transform3<float>, algebra::vc::cartesian2<float>,
-               algebra::vc::polar2<float>, algebra::vc::cylindrical2<float> >,
+               algebra::vc::polar2<float>, algebra::vc::cylindrical2<float>,
+               std::size_t, algebra::vc::matrix_type>,
     test_types<double, algebra::vc::point2<double>, algebra::vc::point3<double>,
                algebra::vc::vector2<double>, algebra::vc::vector3<double>,
                algebra::vc::transform3<double>, algebra::vc::cartesian2<double>,
-               algebra::vc::polar2<double>,
-               algebra::vc::cylindrical2<double> > >
+               algebra::vc::polar2<double>, algebra::vc::cylindrical2<double>,
+               std::size_t, algebra::vc::matrix_type> >
     vc_vc_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics, vc_vc_types,
                                test_specialisation_name);
diff --git a/tests/vecmem/vecmem_cmath.cpp b/tests/vecmem/vecmem_cmath.cpp
index c1bb8572..00e781c3 100644
--- a/tests/vecmem/vecmem_cmath.cpp
+++ b/tests/vecmem/vecmem_cmath.cpp
@@ -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
  */
@@ -38,13 +38,15 @@ typedef testing::Types<
         float, algebra::vecmem::point2<float>, algebra::vecmem::point3<float>,
         algebra::vecmem::vector2<float>, algebra::vecmem::vector3<float>,
         algebra::vecmem::transform3<float>, algebra::vecmem::cartesian2<float>,
-        algebra::vecmem::polar2<float>, algebra::vecmem::cylindrical2<float> >,
+        algebra::vecmem::polar2<float>, algebra::vecmem::cylindrical2<float>,
+        std::size_t, algebra::vecmem::matrix_type>,
     test_types<
         double, algebra::vecmem::point2<double>,
         algebra::vecmem::point3<double>, algebra::vecmem::vector2<double>,
         algebra::vecmem::vector3<double>, algebra::vecmem::transform3<double>,
         algebra::vecmem::cartesian2<double>, algebra::vecmem::polar2<double>,
-        algebra::vecmem::cylindrical2<double> > >
+        algebra::vecmem::cylindrical2<double>, std::size_t,
+        algebra::vecmem::matrix_type> >
     vecmem_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_host_basics,
                                vecmem_cmath_types, test_specialisation_name);