From 44d9055d7d2071212e4dc2875bb78a94ce7cb17f Mon Sep 17 00:00:00 2001
From: Joana Niermann <joana.niermann@cern.ch>
Date: Wed, 24 Apr 2024 16:35:49 +0200
Subject: [PATCH] Make AoS and SoA share the matrix44 and transform3 tpyes

---
 frontend/vc_aos/include/algebra/vc_aos.hpp    |  42 ++-
 frontend/vc_soa/CMakeLists.txt                |   2 +-
 frontend/vc_soa/include/algebra/vc_soa.hpp    |  10 +-
 .../algebra/math/impl/cmath_getter.hpp        |  18 +-
 math/vc_aos/CMakeLists.txt                    |   3 +-
 .../algebra/math/impl/vc_aos_transform3.hpp   | 278 +++++------------
 math/vc_soa/CMakeLists.txt                    |   5 +-
 .../algebra/math/impl/vc_soa_transform3.hpp   | 283 ------------------
 math/vc_soa/include/algebra/math/vc_soa.hpp   |   1 -
 storage/common/CMakeLists.txt                 |   3 +-
 .../include/algebra/storage/matrix44.hpp}     |  76 +++--
 .../common/include/algebra/storage/vector.hpp |   2 +
 .../vc_aos/include/algebra/storage/vc_aos.hpp |   1 +
 storage/vc_soa/CMakeLists.txt                 |   3 +-
 14 files changed, 189 insertions(+), 538 deletions(-)
 delete mode 100644 math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp
 rename storage/{vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp => common/include/algebra/storage/matrix44.hpp} (54%)

diff --git a/frontend/vc_aos/include/algebra/vc_aos.hpp b/frontend/vc_aos/include/algebra/vc_aos.hpp
index 4245f9d1..a32c3926 100644
--- a/frontend/vc_aos/include/algebra/vc_aos.hpp
+++ b/frontend/vc_aos/include/algebra/vc_aos.hpp
@@ -23,7 +23,7 @@ namespace vc_aos {
 /// @{
 
 template <typename T>
-using transform3 = math::transform3<storage_type, T, vector3<T>, point2<T>>;
+using transform3 = math::transform3<storage_type, T>;
 
 /// @}
 
@@ -42,11 +42,17 @@ using vc_aos::math::theta;
 
 /// @}
 
-/// Function extracting a slice from the matrix used by
-/// @c algebra::vc_aos::transform3<float>
-template <std::size_t SIZE, std::enable_if_t<SIZE <= 4, bool> = true>
-ALGEBRA_HOST_DEVICE inline vc_aos::vector3<float> vector(
-    const vc_aos::transform3<float>::matrix44& m,
+/// @name Getter functions on @c algebra::vc_aos::matrix_type
+/// @{
+
+using cmath::element;
+
+/// Function extracting a slice from matrix44 - const
+template <std::size_t SIZE, template <typename, std::size_t> class array_t,
+          typename value_t, std::size_t N,
+          std::enable_if_t<SIZE <= 4, bool> = true>
+ALGEBRA_HOST_DEVICE inline const auto& vector(
+    const storage::matrix44<array_t, value_t, N>& m,
     std::size_t
 #ifndef NDEBUG
         row
@@ -66,15 +72,20 @@ ALGEBRA_HOST_DEVICE inline vc_aos::vector3<float> vector(
     case 3:
       return m.t;
     default:
+#ifndef _MSC_VER
+      __builtin_unreachable();
+#else
       return m.x;
+#endif
   }
 }
 
-/// Function extracting a slice from the matrix used by
-/// @c algebra::vc_aos::transform3<double>
-template <std::size_t SIZE, std::enable_if_t<SIZE <= 4, bool> = true>
-ALGEBRA_HOST_DEVICE inline vc_aos::vector3<double> vector(
-    const vc_aos::transform3<double>::matrix44& m,
+/// Function extracting a slice from matrix44 - non-const
+template <std::size_t SIZE, template <typename, std::size_t> class array_t,
+          typename value_t, std::size_t N,
+          std::enable_if_t<SIZE <= 4, bool> = true>
+ALGEBRA_HOST_DEVICE inline auto& vector(
+    storage::matrix44<array_t, value_t, N>& m,
     std::size_t
 #ifndef NDEBUG
         row
@@ -94,15 +105,14 @@ ALGEBRA_HOST_DEVICE inline vc_aos::vector3<double> vector(
     case 3:
       return m.t;
     default:
+#ifndef _MSC_VER
+      __builtin_unreachable();
+#else
       return m.x;
+#endif
   }
 }
 
-/// @name Getter functions on @c algebra::vc_aos::matrix_type
-/// @{
-
-using cmath::element;
-
 /// @}
 
 }  // namespace getter
diff --git a/frontend/vc_soa/CMakeLists.txt b/frontend/vc_soa/CMakeLists.txt
index b6344424..d4980cca 100644
--- a/frontend/vc_soa/CMakeLists.txt
+++ b/frontend/vc_soa/CMakeLists.txt
@@ -9,6 +9,6 @@ algebra_add_library( algebra_vc_soa vc_soa
    "include/algebra/vc_soa.hpp" )
 target_link_libraries( algebra_vc_soa
    INTERFACE algebra::common algebra::vc_soa_storage algebra::cmath_math
-             algebra::vc_soa_math )
+             algebra::vc_soa_math algebra::vc_aos_math )
 algebra_test_public_headers( algebra_vc_soa
    "algebra/vc_soa.hpp" )
diff --git a/frontend/vc_soa/include/algebra/vc_soa.hpp b/frontend/vc_soa/include/algebra/vc_soa.hpp
index e89eb8aa..89a15ab6 100644
--- a/frontend/vc_soa/include/algebra/vc_soa.hpp
+++ b/frontend/vc_soa/include/algebra/vc_soa.hpp
@@ -1,6 +1,6 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2023 CERN for the benefit of the ACTS project
+ * (c) 2023-2024 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/math/impl/vc_aos_transform3.hpp"
 #include "algebra/math/vc_soa.hpp"
 #include "algebra/storage/vc_soa.hpp"
 
@@ -32,7 +33,7 @@ namespace vc_soa {
 /// @{
 
 template <typename T>
-using transform3 = math::transform3<T>;
+using transform3 = algebra::vc_aos::math::transform3<storage_type, T>;
 
 /// @}
 
@@ -66,7 +67,8 @@ using vc_soa::math::normalize;
 
 }  // namespace vector
 
-namespace matrix {
+// Produces clash with matrix typedefs in other plugins
+/*namespace matrix {
 
 using size_type = vc_soa::size_type;
 
@@ -76,6 +78,6 @@ using array_type = vc_soa::storage_type<T, N>;
 template <typename T, size_type ROWS, size_type COLS>
 using matrix_type = vc_soa::matrix_type<T, ROWS, COLS>;
 
-}  // namespace matrix
+}*/  // namespace matrix
 
 }  // namespace algebra
diff --git a/math/cmath/include/algebra/math/impl/cmath_getter.hpp b/math/cmath/include/algebra/math/impl/cmath_getter.hpp
index 6d7b2199..ef0c31b5 100644
--- a/math/cmath/include/algebra/math/impl/cmath_getter.hpp
+++ b/math/cmath/include/algebra/math/impl/cmath_getter.hpp
@@ -23,7 +23,8 @@ namespace algebra::cmath {
  * @param v the input vector
  **/
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N, std::enable_if_t<N >= 2, bool> = true>
+          typename scalar_t, size_type N,
+          std::enable_if_t<(N >= 2) && std::is_scalar_v<scalar_t>, bool> = true>
 ALGEBRA_HOST_DEVICE inline scalar_t phi(
     const array_t<scalar_t, N> &v) noexcept {
 
@@ -35,7 +36,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t phi(
  * @param v the input vector
  **/
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N, std::enable_if_t<N >= 3, bool> = true>
+          typename scalar_t, size_type N,
+          std::enable_if_t<(N >= 3) && std::is_scalar_v<scalar_t>, bool> = true>
 ALGEBRA_HOST_DEVICE inline scalar_t theta(
     const array_t<scalar_t, N> &v) noexcept {
 
@@ -48,7 +50,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t theta(
  * @param v the input vector
  **/
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N, std::enable_if_t<N >= 2, bool> = true>
+          typename scalar_t, size_type N,
+          std::enable_if_t<(N >= 2) && std::is_scalar_v<scalar_t>, bool> = true>
 ALGEBRA_HOST_DEVICE inline scalar_t perp(
     const array_t<scalar_t, N> &v) noexcept {
 
@@ -60,7 +63,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t perp(
  * @param v the input vector
  **/
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N, std::enable_if_t<N == 2, bool> = true>
+          typename scalar_t, size_type N,
+          std::enable_if_t<(N == 2) && std::is_scalar_v<scalar_t>, bool> = true>
 ALGEBRA_HOST_DEVICE inline scalar_t norm(const array_t<scalar_t, N> &v) {
 
   return perp(v);
@@ -71,7 +75,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t norm(const array_t<scalar_t, N> &v) {
  * @param v the input vector
  **/
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N, std::enable_if_t<N >= 3, bool> = true>
+          typename scalar_t, size_type N,
+          std::enable_if_t<(N >= 3) && std::is_scalar_v<scalar_t>, bool> = true>
 ALGEBRA_HOST_DEVICE inline scalar_t norm(const array_t<scalar_t, N> &v) {
 
   return algebra::math::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
@@ -83,7 +88,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t norm(const array_t<scalar_t, N> &v) {
  * @param v the input vector
  **/
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N, std::enable_if_t<N >= 3, bool> = true>
+          typename scalar_t, size_type N,
+          std::enable_if_t<(N >= 3) && std::is_scalar_v<scalar_t>, bool> = true>
 ALGEBRA_HOST_DEVICE inline scalar_t eta(
     const array_t<scalar_t, N> &v) noexcept {
 
diff --git a/math/vc_aos/CMakeLists.txt b/math/vc_aos/CMakeLists.txt
index 64a723c7..8e6a63a9 100644
--- a/math/vc_aos/CMakeLists.txt
+++ b/math/vc_aos/CMakeLists.txt
@@ -11,7 +11,6 @@ algebra_add_library( algebra_vc_aos_math vc_aos_math
    "include/algebra/math/impl/vc_aos_transform3.hpp"
    "include/algebra/math/impl/vc_aos_vector.hpp" )
 target_link_libraries( algebra_vc_aos_math
-   INTERFACE Vc::Vc algebra::common algebra::common_math algebra::vc_storage
-             algebra::cmath_math )
+   INTERFACE Vc::Vc algebra::common algebra::common_math algebra::vc_storage )
 algebra_test_public_headers( algebra_vc_aos_math
    "algebra/math/vc_aos.hpp" )
diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp
index d40f8ff8..be3eea06 100644
--- a/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp
+++ b/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp
@@ -8,9 +8,10 @@
 #pragma once
 
 // Project include(s).
-#include "algebra/math/cmath.hpp"
-#include "algebra/math/impl/vc_aos_vector.hpp"
+#include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
+#include "algebra/storage/matrix44.hpp"
+#include "algebra/storage/vector.hpp"
 
 // Vc include(s).
 #ifdef _MSC_VER
@@ -26,95 +27,46 @@
 
 namespace algebra::vc_aos::math {
 
-namespace internal {
-
-/// 4x4 matrix type used by @c algebra::vc_aos::math::transform3
-template <typename scalar_t>
-struct matrix44 {
-
-  /// Equality operator between two matrices
-  bool operator==(const matrix44 &rhs) const {
-    return ((x == rhs.x).isFull() && (y == rhs.y).isFull() &&
-            (z == rhs.z).isFull() && (t == rhs.t).isFull());
-  }
-
-  /// Data variables
-  Vc::SimdArray<scalar_t, 4> x, y, z, t;
-
-};  // struct matrix44
-
-/// Functor used to access elements of Vc matrices
-template <template <std::size_t> class array_t, typename scalar_t>
-struct element_getter {
-
-  /// Get const access to a matrix element
-  ALGEBRA_HOST inline scalar_t operator()(const array_t<16> &m,
-                                          unsigned int row,
-                                          unsigned int col) const {
-
-    // Make sure that the indices are valid.
-    assert(row < 4);
-    assert(col < 4);
-
-    // Return the selected element.
-    return m[col * 4 + row];
-  }
-
-  /// Get const access to a matrix element
-  ALGEBRA_HOST inline scalar_t operator()(const matrix44<scalar_t> &m,
-                                          unsigned int row,
-                                          unsigned int col) const {
-
-    // Make sure that the indices are valid.
-    assert(row < 4);
-    assert(col < 4);
-
-    // Return the selected element.
-    switch (row) {
-      case 0:
-        return m.x[col];
-      case 1:
-        return m.y[col];
-      case 2:
-        return m.z[col];
-      case 3:
-        return m.t[col];
-      default:
-        return 0;
-    }
-  }
-
-};  // struct element_getter
-
-}  // namespace internal
+using algebra::storage::operator*;
+using algebra::storage::operator/;
+using algebra::storage::operator-;
+using algebra::storage::operator+;
 
 /// Transform wrapper class to ensure standard API within differnt plugins
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          typename vector3_t = array_t<scalar_t, 3>,
-          typename point2_t = array_t<scalar_t, 2>>
+template <template <typename, std::size_t> class array_t, typename scalar_t>
 struct transform3 {
 
   /// @name Type definitions for the struct
   /// @{
 
-  /// Array type used by the transform
-  template <std::size_t N>
-  using array_type = array_t<scalar_t, N>;
   /// Scalar type used by the transform
   using scalar_type = scalar_t;
+  /// The type of the matrix elements (scalar for AoS, Vc::Vector for SoA)
+  using value_type =
+      std::conditional_t<Vc::is_simd_vector<array_t<scalar_t, 4>>::value,
+                         scalar_t, Vc::Vector<scalar_t>>;
 
-  /// 3-element "vector" type
-  using vector3 = vector3_t;
-  /// Point in 3D space
+  // AoS stores four elements per vector for alignment
+  using storage_dim =
+      std::conditional_t<Vc::is_simd_vector<value_type>::value,
+                         std::integral_constant<std::size_t, 3u>,
+                         std::integral_constant<std::size_t, 4u>>;
+
+  template <std::size_t N>
+  using array_type = array_t<value_type, N>;
+
+  /// 3-element "vector" type (does not observe translations)
+  using vector3 = storage::vector<storage_dim::value, value_type, array_t>;
+  /// Point in 3D space (does observe translations)
   using point3 = vector3;
   /// Point in 2D space
-  using point2 = point2_t;
+  using point2 = storage::vector<2, value_type, array_t>;
 
   /// 4x4 matrix type
-  using matrix44 = internal::matrix44<scalar_type>;
+  using matrix44 = storage::matrix44<array_t, value_type, storage_dim::value>;
 
   /// Function (object) used for accessing a matrix element
-  using element_getter = internal::element_getter<array_type, scalar_type>;
+  using element_getter = storage::element_getter;
 
   /// @}
 
@@ -125,6 +77,9 @@ struct transform3 {
   matrix44 _data_inv;
 
   /// @}
+  /// Default constructor: identity
+  ALGEBRA_HOST_DEVICE
+  transform3() = default;
 
   /// Contructor with arguments: t, x, y, z
   ///
@@ -134,20 +89,8 @@ struct transform3 {
   /// @param z the z axis of the new frame, normal vector for planes
   ALGEBRA_HOST_DEVICE
   transform3(const vector3 &t, const vector3 &x, const vector3 &y,
-             const vector3 &z, bool get_inverse = true) {
-    _data.x = x;
-    _data.y = y;
-    _data.z = z;
-    _data.t = t;
-    _data.x[3] = 0.f;
-    _data.y[3] = 0.f;
-    _data.z[3] = 0.f;
-    _data.t[3] = 1.f;
-
-    if (get_inverse) {
-      _data_inv = invert(_data);
-    }
-  }
+             const vector3 &z, [[maybe_unused]] bool get_inverse = true)
+      : _data{x, y, z, t}, _data_inv{invert(_data)} {}
 
   /// Contructor with arguments: t, z, x
   ///
@@ -159,32 +102,22 @@ struct transform3 {
   ALGEBRA_HOST_DEVICE
   transform3(const vector3 &t, const vector3 &z, const vector3 &x,
              bool get_inverse = true)
-      : transform3(t, x, cross(z, x), z, get_inverse) {}
+      : transform3(t, x,
+                   vector3(z[1] * x[2] - x[1] * z[2], z[2] * x[0] - x[2] * z[0],
+                           z[0] * x[1] - x[0] * z[1]),
+                   z, get_inverse) {}
 
   /// Constructor with arguments: translation
   ///
   /// @param t is the transform
   ALGEBRA_HOST_DEVICE
-  transform3(const vector3 &t) {
-
-    _data.x = {1.f, 0.f, 0.f, 0.f};
-    _data.y = {0.f, 1.f, 0.f, 0.f};
-    _data.z = {0.f, 0.f, 1.f, 0.f};
-    _data.t = t;
-    _data.t[3] = 1.f;
-
-    _data_inv = invert(_data);
-  }
+  transform3(const vector3 &t) : _data{t}, _data_inv{invert(_data)} {}
 
   /// Constructor with arguments: matrix
   ///
-  /// @param m is the full 4x4 matrix
+  /// @param m is the full 4x4 matrix with simd-vector elements
   ALGEBRA_HOST_DEVICE
-  transform3(const matrix44 &m) {
-
-    _data = m;
-    _data_inv = invert(_data);
-  }
+  transform3(const matrix44 &m) : _data{m}, _data_inv{invert(_data)} {}
 
   /// Constructor with arguments: matrix as std::aray of scalar
   ///
@@ -192,67 +125,43 @@ struct transform3 {
   ALGEBRA_HOST_DEVICE
   transform3(const array_type<16> &ma) {
 
-    _data.x = {ma[0], ma[4], ma[8], ma[12]};
-    _data.y = {ma[1], ma[5], ma[9], ma[13]};
-    _data.z = {ma[2], ma[6], ma[10], ma[14]};
-    _data.t = {ma[3], ma[7], ma[11], ma[15]};
-
+    _data.x = typename matrix44::vector_type{ma[0], ma[4], ma[8], ma[12]};
+    _data.y = typename matrix44::vector_type{ma[1], ma[5], ma[9], ma[13]};
+    _data.z = typename matrix44::vector_type{ma[2], ma[6], ma[10], ma[14]};
+    _data.t = typename matrix44::vector_type{ma[3], ma[7], ma[11], ma[15]};
     _data_inv = invert(_data);
   }
 
-  /// Constructor with arguments: identity
-  ALGEBRA_HOST_DEVICE
-  transform3() {
-
-    _data.x = {1.f, 0.f, 0.f, 0.f};
-    _data.y = {0.f, 1.f, 0.f, 0.f};
-    _data.z = {0.f, 0.f, 1.f, 0.f};
-    _data.t = {0.f, 0.f, 0.f, 1.f};
-
-    _data_inv = _data;
-  }
-
-  /// Default contructors
+  /// Defaults
   transform3(const transform3 &rhs) = default;
   ~transform3() = default;
 
   /// Equality operator
   ALGEBRA_HOST_DEVICE
-  inline bool operator==(const transform3 &rhs) const {
+  inline constexpr bool operator==(const transform3 &rhs) const {
     return (_data == rhs._data);
   }
 
+  /// Matrix access operator
+  ALGEBRA_HOST_DEVICE
+  inline const value_type &operator()(std::size_t i, std::size_t j) const {
+    return element_getter{}(_data, i, j);
+  }
+  ALGEBRA_HOST_DEVICE
+  inline value_type &operator()(std::size_t i, std::size_t j) {
+    return element_getter{}(_data, i, j);
+  }
+
   /// 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 m.t[0] * m.z[1] * m.y[2] * m.x[3] -
-           m.z[0] * m.t[1] * m.y[2] * m.x[3] -
-           m.t[0] * m.y[1] * m.z[2] * m.x[3] +
-           m.y[0] * m.t[1] * m.z[2] * m.x[3] +
-           m.z[0] * m.y[1] * m.t[2] * m.x[3] -
-           m.y[0] * m.z[1] * m.t[2] * m.x[3] -
-           m.t[0] * m.z[1] * m.x[2] * m.y[3] +
-           m.z[0] * m.t[1] * m.x[2] * m.y[3] +
-           m.t[0] * m.x[1] * m.z[2] * m.y[3] -
-           m.x[0] * m.t[1] * m.z[2] * m.y[3] -
-           m.z[0] * m.x[1] * m.t[2] * m.y[3] +
-           m.x[0] * m.z[1] * m.t[2] * m.y[3] +
-           m.t[0] * m.y[1] * m.x[2] * m.z[3] -
-           m.y[0] * m.t[1] * m.x[2] * m.z[3] -
-           m.t[0] * m.x[1] * m.y[2] * m.z[3] +
-           m.x[0] * m.t[1] * m.y[2] * m.z[3] +
-           m.y[0] * m.x[1] * m.t[2] * m.z[3] -
-           m.x[0] * m.y[1] * m.t[2] * m.z[3] -
-           m.z[0] * m.y[1] * m.x[2] * m.t[3] +
-           m.y[0] * m.z[1] * m.x[2] * m.t[3] +
-           m.z[0] * m.x[1] * m.y[2] * m.t[3] -
-           m.x[0] * m.z[1] * m.y[2] * m.t[3] -
-           m.y[0] * m.x[1] * m.z[2] * m.t[3] +
-           m.x[0] * m.y[1] * m.z[2] * m.t[3];
+  inline constexpr value_type determinant(const matrix44 &m) const {
+    return -m.z[0] * m.y[1] * m.x[2] + m.y[0] * m.z[1] * m.x[2] +
+           m.z[0] * m.x[1] * m.y[2] - m.x[0] * m.z[1] * m.y[2] -
+           m.y[0] * m.x[1] * m.z[2] + m.x[0] * m.y[1] * m.z[2];
   }
 
   /// The inverse of a 4x4 matrix
@@ -261,44 +170,20 @@ struct transform3 {
   ///
   /// @return an inverse matrix
   ALGEBRA_HOST_DEVICE
-  static inline matrix44 invert(const matrix44 &m) {
+  inline constexpr matrix44 invert(const matrix44 &m) const {
     matrix44 i;
-    i.x[0] = m.z[1] * m.t[2] * m.y[3] - m.t[1] * m.z[2] * m.y[3] +
-             m.t[1] * m.y[2] * m.z[3] - m.y[1] * m.t[2] * m.z[3] -
-             m.z[1] * m.y[2] * m.t[3] + m.y[1] * m.z[2] * m.t[3];
-    i.x[1] = m.t[1] * m.z[2] * m.x[3] - m.z[1] * m.t[2] * m.x[3] -
-             m.t[1] * m.x[2] * m.z[3] + m.x[1] * m.t[2] * m.z[3] +
-             m.z[1] * m.x[2] * m.t[3] - m.x[1] * m.z[2] * m.t[3];
-    i.x[2] = m.y[1] * m.t[2] * m.x[3] - m.t[1] * m.y[2] * m.x[3] +
-             m.t[1] * m.x[2] * m.y[3] - m.x[1] * m.t[2] * m.y[3] -
-             m.y[1] * m.x[2] * m.t[3] + m.x[1] * m.y[2] * m.t[3];
-    i.x[3] = m.z[1] * m.y[2] * m.x[3] - m.y[1] * m.z[2] * m.x[3] -
-             m.z[1] * m.x[2] * m.y[3] + m.x[1] * m.z[2] * m.y[3] +
-             m.y[1] * m.x[2] * m.z[3] - m.x[1] * m.y[2] * m.z[3];
-    i.y[0] = m.t[0] * m.z[2] * m.y[3] - m.z[0] * m.t[2] * m.y[3] -
-             m.t[0] * m.y[2] * m.z[3] + m.y[0] * m.t[2] * m.z[3] +
-             m.z[0] * m.y[2] * m.t[3] - m.y[0] * m.z[2] * m.t[3];
-    i.y[1] = m.z[0] * m.t[2] * m.x[3] - m.t[0] * m.z[2] * m.x[3] +
-             m.t[0] * m.x[2] * m.z[3] - m.x[0] * m.t[2] * m.z[3] -
-             m.z[0] * m.x[2] * m.t[3] + m.x[0] * m.z[2] * m.t[3];
-    i.y[2] = m.t[0] * m.y[2] * m.x[3] - m.y[0] * m.t[2] * m.x[3] -
-             m.t[0] * m.x[2] * m.y[3] + m.x[0] * m.t[2] * m.y[3] +
-             m.y[0] * m.x[2] * m.t[3] - m.x[0] * m.y[2] * m.t[3];
-    i.y[3] = m.y[0] * m.z[2] * m.x[3] - m.z[0] * m.y[2] * m.x[3] +
-             m.z[0] * m.x[2] * m.y[3] - m.x[0] * m.z[2] * m.y[3] -
-             m.y[0] * m.x[2] * m.z[3] + m.x[0] * m.y[2] * m.z[3];
-    i.z[0] = m.z[0] * m.t[1] * m.y[3] - m.t[0] * m.z[1] * m.y[3] +
-             m.t[0] * m.y[1] * m.z[3] - m.y[0] * m.t[1] * m.z[3] -
-             m.z[0] * m.y[1] * m.t[3] + m.y[0] * m.z[1] * m.t[3];
-    i.z[1] = m.t[0] * m.z[1] * m.x[3] - m.z[0] * m.t[1] * m.x[3] -
-             m.t[0] * m.x[1] * m.z[3] + m.x[0] * m.t[1] * m.z[3] +
-             m.z[0] * m.x[1] * m.t[3] - m.x[0] * m.z[1] * m.t[3];
-    i.z[2] = m.y[0] * m.t[1] * m.x[3] - m.t[0] * m.y[1] * m.x[3] +
-             m.t[0] * m.x[1] * m.y[3] - m.x[0] * m.t[1] * m.y[3] -
-             m.y[0] * m.x[1] * m.t[3] + m.x[0] * m.y[1] * m.t[3];
-    i.z[3] = m.z[0] * m.y[1] * m.x[3] - m.y[0] * m.z[1] * m.x[3] -
-             m.z[0] * m.x[1] * m.y[3] + m.x[0] * m.z[1] * m.y[3] +
-             m.y[0] * m.x[1] * m.z[3] - m.x[0] * m.y[1] * m.z[3];
+    i.x[0] = -m.z[1] * m.y[2] + m.y[1] * m.z[2];
+    i.x[1] = m.z[1] * m.x[2] - m.x[1] * m.z[2];
+    i.x[2] = -m.y[1] * m.x[2] + m.x[1] * m.y[2];
+    // i.x[3] = 0;
+    i.y[0] = m.z[0] * m.y[2] - m.y[0] * m.z[2];
+    i.y[1] = -m.z[0] * m.x[2] + m.x[0] * m.z[2];
+    i.y[2] = m.y[0] * m.x[2] - m.x[0] * m.y[2];
+    // i.y[3] = 0;
+    i.z[0] = -m.z[0] * m.y[1] + m.y[0] * m.z[1];
+    i.z[1] = m.z[0] * m.x[1] - m.x[0] * m.z[1];
+    i.z[2] = -m.y[0] * m.x[1] + m.x[0] * m.y[1];
+    // i.z[3] = 0;
     i.t[0] = m.t[0] * m.z[1] * m.y[2] - m.z[0] * m.t[1] * m.y[2] -
              m.t[0] * m.y[1] * m.z[2] + m.y[0] * m.t[1] * m.z[2] +
              m.z[0] * m.y[1] * m.t[2] - m.y[0] * m.z[1] * m.t[2];
@@ -308,15 +193,13 @@ struct transform3 {
     i.t[2] = m.t[0] * m.y[1] * m.x[2] - m.y[0] * m.t[1] * m.x[2] -
              m.t[0] * m.x[1] * m.y[2] + m.x[0] * m.t[1] * m.y[2] +
              m.y[0] * m.x[1] * m.t[2] - m.x[0] * m.y[1] * m.t[2];
-    i.t[3] = m.y[0] * m.z[1] * m.x[2] - m.z[0] * m.y[1] * m.x[2] +
-             m.z[0] * m.x[1] * m.y[2] - m.x[0] * m.z[1] * m.y[2] -
-             m.y[0] * m.x[1] * m.z[2] + m.x[0] * m.y[1] * m.z[2];
-    scalar_type idet = static_cast<scalar_type>(1.) / determinant(i);
+    // i.t[3] = 1;
+    const value_type idet{value_type(1.f) / determinant(i)};
 
-    i.x *= idet;
-    i.y *= idet;
-    i.z *= idet;
-    i.t *= idet;
+    i.x = i.x * idet;
+    i.y = i.y * idet;
+    i.z = i.z * idet;
+    i.t = i.t * idet;
 
     return i;
   }
@@ -325,9 +208,8 @@ struct transform3 {
   ///
   /// @param m is the rotation matrix
   /// @param v is the vector to be rotated
-  template <typename vector3_type>
-  ALGEBRA_HOST_DEVICE static inline auto rotate(const matrix44 &m,
-                                                const vector3_type &v) {
+  ALGEBRA_HOST_DEVICE
+  inline constexpr auto rotate(const matrix44 &m, const vector3 &v) const {
 
     return m.x * v[0] + m.y * v[1] + m.z * v[2];
   }
diff --git a/math/vc_soa/CMakeLists.txt b/math/vc_soa/CMakeLists.txt
index 0a23dc3a..9e277e98 100644
--- a/math/vc_soa/CMakeLists.txt
+++ b/math/vc_soa/CMakeLists.txt
@@ -1,6 +1,6 @@
 # Algebra plugins library, part of the ACTS project (R&D line)
 #
-# (c) 2021-2022 CERN for the benefit of the ACTS project
+# (c) 2021-2024 CERN for the benefit of the ACTS project
 #
 # Mozilla Public License Version 2.0
 
@@ -8,8 +8,7 @@
 algebra_add_library( algebra_vc_soa_math vc_soa_math
    "include/algebra/math/vc_soa.hpp"
    "include/algebra/math/impl/vc_soa_getter.hpp"
-   "include/algebra/math/impl/vc_soa_vector.hpp"
-   "include/algebra/math/impl/vc_soa_transform3.hpp" )
+   "include/algebra/math/impl/vc_soa_vector.hpp")
 target_link_libraries( algebra_vc_soa_math
    INTERFACE algebra::common algebra::common_math algebra::common_storage  algebra::vc_soa_storage Vc::Vc )
 algebra_test_public_headers( algebra_vc_soa_math
diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp
deleted file mode 100644
index aee22c02..00000000
--- a/math/vc_soa/include/algebra/math/impl/vc_soa_transform3.hpp
+++ /dev/null
@@ -1,283 +0,0 @@
-/** Algebra plugins library, part of the ACTS project
- *
- * (c) 2023-2024 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-// Project include(s).
-#include "algebra/math/common.hpp"
-#include "algebra/qualifiers.hpp"
-#include "algebra/storage/impl/vc_soa_matrix44.hpp"
-#include "algebra/storage/vector.hpp"
-
-// Vc include(s).
-#ifdef _MSC_VER
-#pragma warning(push, 0)
-#endif  // MSVC
-#include <Vc/Vc>
-#ifdef _MSC_VER
-#pragma warning(pop)
-#endif  // MSVC
-
-// System include(s).
-#include <array>
-#include <cassert>
-#include <type_traits>
-
-namespace algebra::vc_soa::math {
-
-using algebra::storage::operator*;
-using algebra::storage::operator/;
-using algebra::storage::operator-;
-using algebra::storage::operator+;
-
-/// Transform wrapper class to ensure standard API within differnt plugins
-template <typename scalar_t>
-struct transform3 {
-
-  /// @name Type definitions for the struct
-  /// @{
-
-  /// Scalar type used by the transform
-  using scalar_type = scalar_t;
-  /// The type of the matrix elements (in this case: Vc::Vector)
-  using value_type = Vc::Vector<scalar_t>;
-
-  /// 3-element "vector" type (does not observe translations)
-  using vector3 = storage::vector<3, value_type, std::array>;
-  /// Point in 3D space (does observe translations)
-  using point3 = vector3;
-  /// Point in 2D space
-  using point2 = storage::vector<2, value_type, std::array>;
-
-  /// 4x4 matrix type
-  using matrix44 = algebra::vc_soa::matrix44<scalar_type>;
-
-  /// Function (object) used for accessing a matrix element
-  using element_getter = algebra::vc_soa::element_getter<scalar_type>;
-
-  /// @}
-
-  /// @name Data objects
-  /// @{
-
-  matrix44 _data;
-  matrix44 _data_inv;
-
-  /// @}
-
-  /// Default constructor: identity
-  ALGEBRA_HOST_DEVICE
-  transform3() = default;
-
-  /// Contructor with arguments: t, x, y, z
-  ///
-  /// @param t the translation (or origin of the new frame)
-  /// @param x the x axis of the new frame
-  /// @param y the y axis of the new frame
-  /// @param z the z axis of the new frame, normal vector for planes
-  ALGEBRA_HOST_DEVICE
-  transform3(const vector3 &t, const vector3 &x, const vector3 &y,
-             const vector3 &z, [[maybe_unused]] bool get_inverse = true)
-      : _data{x, y, z, t}, _data_inv{invert(_data)} {}
-
-  /// Contructor with arguments: t, z, x
-  ///
-  /// @param t the translation (or origin of the new frame)
-  /// @param z the z axis of the new frame, normal vector for planes
-  /// @param x the x axis of the new frame
-  ///
-  /// @note y will be constructed by cross product
-  ALGEBRA_HOST_DEVICE
-  transform3(const vector3 &t, const vector3 &z, const vector3 &x,
-             bool get_inverse = true)
-      : transform3(t, x, cross(z, x), z, get_inverse) {}
-
-  /// Constructor with arguments: translation
-  ///
-  /// @param t is the transform
-  ALGEBRA_HOST_DEVICE
-  transform3(const vector3 &t) : _data{t}, _data_inv{invert(_data)} {}
-
-  /// Constructor with arguments: matrix
-  ///
-  /// @param m is the full 4x4 matrix with simd-vector elements
-  ALGEBRA_HOST_DEVICE
-  transform3(const matrix44 &m) : _data{m}, _data_inv{invert(_data)} {}
-
-  /// Defaults
-  transform3(const transform3 &rhs) = default;
-  ~transform3() = default;
-
-  /// Equality operator
-  ALGEBRA_HOST_DEVICE
-  inline constexpr bool operator==(const transform3 &rhs) const {
-    return (_data == rhs._data);
-  }
-
-  /// Matrix access operator
-  ALGEBRA_HOST_DEVICE
-  inline const value_type &operator()(std::size_t i, std::size_t j) const {
-    return element_getter{}(_data, i, j);
-  }
-  ALGEBRA_HOST_DEVICE
-  inline value_type &operator()(std::size_t i, std::size_t j) {
-    return element_getter{}(_data, i, j);
-  }
-
-  /// The determinant of a 4x4 matrix
-  ///
-  /// @param m is the matrix
-  ///
-  /// @return a sacalar determinant - no checking done
-  ALGEBRA_HOST_DEVICE
-  inline constexpr value_type determinant(const matrix44 &m) const {
-    return -m.z[0] * m.y[1] * m.x[2] + m.y[0] * m.z[1] * m.x[2] +
-           m.z[0] * m.x[1] * m.y[2] - m.x[0] * m.z[1] * m.y[2] -
-           m.y[0] * m.x[1] * m.z[2] + m.x[0] * m.y[1] * m.z[2];
-  }
-
-  /// The inverse of a 4x4 matrix
-  ///
-  /// @param m is the matrix
-  ///
-  /// @return an inverse matrix
-  ALGEBRA_HOST_DEVICE
-  inline constexpr matrix44 invert(const matrix44 &m) const {
-    matrix44 i;
-    i.x[0] = -m.z[1] * m.y[2] + m.y[1] * m.z[2];
-    i.x[1] = m.z[1] * m.x[2] - m.x[1] * m.z[2];
-    i.x[2] = -m.y[1] * m.x[2] + m.x[1] * m.y[2];
-    // i.x[3] = 0;
-    i.y[0] = m.z[0] * m.y[2] - m.y[0] * m.z[2];
-    i.y[1] = -m.z[0] * m.x[2] + m.x[0] * m.z[2];
-    i.y[2] = m.y[0] * m.x[2] - m.x[0] * m.y[2];
-    // i.y[3] = 0;
-    i.z[0] = -m.z[0] * m.y[1] + m.y[0] * m.z[1];
-    i.z[1] = m.z[0] * m.x[1] - m.x[0] * m.z[1];
-    i.z[2] = -m.y[0] * m.x[1] + m.x[0] * m.y[1];
-    // i.z[3] = 0;
-    i.t[0] = m.t[0] * m.z[1] * m.y[2] - m.z[0] * m.t[1] * m.y[2] -
-             m.t[0] * m.y[1] * m.z[2] + m.y[0] * m.t[1] * m.z[2] +
-             m.z[0] * m.y[1] * m.t[2] - m.y[0] * m.z[1] * m.t[2];
-    i.t[1] = m.z[0] * m.t[1] * m.x[2] - m.t[0] * m.z[1] * m.x[2] +
-             m.t[0] * m.x[1] * m.z[2] - m.x[0] * m.t[1] * m.z[2] -
-             m.z[0] * m.x[1] * m.t[2] + m.x[0] * m.z[1] * m.t[2];
-    i.t[2] = m.t[0] * m.y[1] * m.x[2] - m.y[0] * m.t[1] * m.x[2] -
-             m.t[0] * m.x[1] * m.y[2] + m.x[0] * m.t[1] * m.y[2] +
-             m.y[0] * m.x[1] * m.t[2] - m.x[0] * m.y[1] * m.t[2];
-    // i.t[3] = 1;
-    const value_type idet{value_type::One() / determinant(i)};
-
-    i.x = i.x * idet;
-    i.y = i.y * idet;
-    i.z = i.z * idet;
-    i.t = i.t * idet;
-
-    return i;
-  }
-
-  /// Rotate a vector into / from a frame
-  ///
-  /// @param m is the rotation matrix
-  /// @param v is the vector to be rotated
-  ALGEBRA_HOST_DEVICE
-  inline constexpr auto rotate(const matrix44 &m, const vector3 &v) const {
-
-    return m.x * v[0] + m.y * v[1] + m.z * v[2];
-  }
-
-  /// This method retrieves x axis
-  ALGEBRA_HOST_DEVICE
-  inline const auto &x() const { return _data.x; }
-
-  /// This method retrieves y axis
-  ALGEBRA_HOST_DEVICE
-  inline const auto &y() const { return _data.y; }
-
-  /// This method retrieves z axis
-  ALGEBRA_HOST_DEVICE
-  inline const auto &z() const { return _data.z; }
-
-  /// @returns the translation of the transform
-  ALGEBRA_HOST_DEVICE
-  inline point3 translation() const { return _data.t; }
-
-  /// @returns the 4x4 matrix of the transform
-  ALGEBRA_HOST_DEVICE
-  inline const matrix44 &matrix() const { return _data; }
-
-  /// @returns the 4x4 matrix of the inverse transform
-  ALGEBRA_HOST_DEVICE
-  inline const matrix44 &matrix_inverse() const { return _data_inv; }
-
-  /// This method transforms a point from the local 3D cartesian frame
-  /// to the global 3D cartesian frame
-  ///
-  /// @tparam point3_t 3D point
-  ///
-  /// @param v is the point to be transformed
-  ///
-  /// @return a global point
-  template <
-      typename point3_t,
-      std::enable_if_t<std::is_convertible_v<point3_t, point3>, bool> = true>
-  ALGEBRA_HOST_DEVICE inline point3 point_to_global(const point3_t &v) const {
-
-    return rotate(_data, v) + _data.t;
-  }
-
-  /// This method transforms a point from the global 3D cartesian frame
-  /// into the local 3D cartesian frame
-  ///
-  /// @tparam point3_t 3D point
-  ///
-  /// @param v is the point to be transformed
-  ///
-  /// @return a local point
-  template <
-      typename point3_t,
-      std::enable_if_t<std::is_convertible_v<point3_t, point3>, bool> = true>
-  ALGEBRA_HOST_DEVICE inline point3 point_to_local(const point3_t &v) const {
-
-    return rotate(_data_inv, v) + _data_inv.t;
-  }
-
-  /// This method transforms a vector from the local 3D cartesian frame
-  /// to the global 3D cartesian frame
-  ///
-  /// @tparam vector3_t 3D vector
-  ///
-  /// @param v is the vector to be transformed
-  ///
-  /// @return a vector in global coordinates
-  template <
-      typename vector3_t,
-      std::enable_if_t<std::is_convertible_v<vector3_t, vector3>, bool> = true>
-  ALGEBRA_HOST_DEVICE inline vector3 vector_to_global(
-      const vector3_t &v) const {
-
-    return rotate(_data, v);
-  }
-
-  /// This method transforms a vector from the global 3D cartesian frame
-  /// into the local 3D cartesian frame
-  ///
-  /// @tparam vector3_t 3D vector
-  ///
-  /// @param v is the vector to be transformed
-  ///
-  /// @return a vector in local coordinates
-  template <
-      typename vector3_t,
-      std::enable_if_t<std::is_convertible_v<vector3_t, vector3>, bool> = true>
-  ALGEBRA_HOST_DEVICE inline vector3 vector_to_local(const vector3_t &v) const {
-
-    return rotate(_data_inv, v);
-  }
-};  // struct transform3
-
-}  // namespace algebra::vc_soa::math
diff --git a/math/vc_soa/include/algebra/math/vc_soa.hpp b/math/vc_soa/include/algebra/math/vc_soa.hpp
index fc87b25b..4c4c15c3 100644
--- a/math/vc_soa/include/algebra/math/vc_soa.hpp
+++ b/math/vc_soa/include/algebra/math/vc_soa.hpp
@@ -9,5 +9,4 @@
 
 // Project include(s).
 #include "algebra/math/impl/vc_soa_getter.hpp"
-#include "algebra/math/impl/vc_soa_transform3.hpp"
 #include "algebra/math/impl/vc_soa_vector.hpp"
diff --git a/storage/common/CMakeLists.txt b/storage/common/CMakeLists.txt
index 025caf8c..1c084617 100644
--- a/storage/common/CMakeLists.txt
+++ b/storage/common/CMakeLists.txt
@@ -1,12 +1,13 @@
 # Algebra plugins library, part of the ACTS project (R&D line)
 #
-# (c) 2023 CERN for the benefit of the ACTS project
+# (c) 2023-2024 CERN for the benefit of the ACTS project
 #
 # Mozilla Public License Version 2.0
 
 # Set up the library.
 algebra_add_library(algebra_common_storage common_storage
    "include/algebra/storage/array_operators.hpp"
+   "include/algebra/storage/matrix44.hpp"
    "include/algebra/storage/vector.hpp")
 target_link_libraries(algebra_common_storage
    INTERFACE algebra::common)
diff --git a/storage/vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp b/storage/common/include/algebra/storage/matrix44.hpp
similarity index 54%
rename from storage/vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp
rename to storage/common/include/algebra/storage/matrix44.hpp
index 6c336cff..4cff0ba2 100644
--- a/storage/vc_soa/include/algebra/storage/impl/vc_soa_matrix44.hpp
+++ b/storage/common/include/algebra/storage/matrix44.hpp
@@ -6,6 +6,9 @@
  */
 #pragma once
 
+// Project include(s).
+#include "algebra/storage/vector.hpp"
+
 // Vc include(s).
 #ifdef _MSC_VER
 #pragma warning(push, 0)
@@ -20,16 +23,16 @@
 #include <cstddef>
 #include <limits>
 
-namespace algebra::vc_soa {
+namespace algebra::storage {
 
 /// 4x4 matrix type used by @c algebra::vc_soa::math::transform3 that has simd
 /// vectors as matrix elements
-template <typename scalar_t>
+template <template <typename, std::size_t> class array_t, typename value_t,
+          std::size_t N>
 struct matrix44 {
 
-  using vector_type = storage::vector<3, Vc::Vector<scalar_t>, std::array>;
-  using value_type = Vc::Vector<scalar_t>;
-  using scalar_type = scalar_t;
+  using vector_type = storage::vector<N, value_t, array_t>;
+  using value_type = value_t;
 
   /// Default constructor: Identity with no translation
   matrix44()
@@ -39,17 +42,17 @@ struct matrix44 {
         t{0.f, 0.f, 0.f} {}
 
   /// Identity rotation with translation @param translation
-  matrix44(const vector_type &v)
+  matrix44(const vector_type v)
       : x{1.f, 0.f, 0.f}, y{0.f, 1.f, 0.f}, z{0.f, 0.f, 1.f}, t{v} {}
 
   /// Construct from given column vectors @param x, @param y, @param z, @param t
-  matrix44(const vector_type &v_0, const vector_type &v_1,
-           const vector_type &v_2, const vector_type &v_3)
+  matrix44(const vector_type v_0, const vector_type v_1, const vector_type v_2,
+           const vector_type v_3)
       : x{v_0}, y{v_1}, z{v_2}, t{v_3} {}
 
   /// Identity rotation with translation from single elemenst @param t_0,
   /// @param t_1, @param t_2
-  matrix44(const value_type &t_0, const value_type &t_1, const value_type &t_2)
+  matrix44(const value_type t_0, const value_type t_1, const value_type t_2)
       : matrix44{{t_0, t_1, t_2}} {}
 
   /// Construct from elements (simd vector) of matrix column vectors:
@@ -57,33 +60,49 @@ struct matrix44 {
   /// - column 1: @param y_0, @param y_1, @param y_2
   /// - column 2: @param z_0, @param z_1, @param z_2
   /// - column 3: @param t_0, @param t_1, @param t_2
-  matrix44(const value_type &x_0, const value_type &x_1, const value_type &x_2,
-           const value_type &y_0, const value_type &y_1, const value_type &y_2,
-           const value_type &z_0, const value_type &z_1, const value_type &z_2,
-           const value_type &t_0, const value_type &t_1, const value_type &t_2)
+  matrix44(const value_type x_0, const value_type x_1, const value_type x_2,
+           const value_type y_0, const value_type y_1, const value_type y_2,
+           const value_type z_0, const value_type z_1, const value_type z_2,
+           const value_type t_0, const value_type t_1, const value_type t_2)
       : x{{x_0, x_1, x_2}},
         y{{y_0, y_1, y_2}},
         z{{z_0, z_1, z_2}},
         t{{t_0, t_1, t_2}} {}
 
   /// Equality operator between two matrices
+  /// @{
+  /// SoA
+  template <typename V = value_t,
+            typename std::enable_if_t<!std::is_scalar_v<V>, bool> = true>
   bool operator==(const matrix44 &rhs) const {
     return ((x == rhs.x) && (y == rhs.y) && (z == rhs.z) && (t == rhs.t));
   }
 
+  /// AoS
+  template <typename V = value_t,
+            typename std::enable_if_t<std::is_scalar_v<V>, bool> = true>
+  bool operator==(const matrix44 &rhs) const {
+    return ((x.get() == rhs.x.get()).isFull() &&
+            (y.get() == rhs.y.get()).isFull() &&
+            (z.get() == rhs.z.get()).isFull() &&
+            (t.get() == rhs.t.get()).isFull());
+  }
+  /// @}
+
   /// Data variables
   vector_type x, y, z, t;
 
 };  // struct matrix44
 
 /// Functor used to access elements of matrix44
-template <typename scalar_t>
 struct element_getter {
 
   /// Get const access to a matrix element
-  ALGEBRA_HOST inline const auto &operator()(const matrix44<scalar_t> &m,
-                                             std::size_t row,
-                                             std::size_t col) const {
+  template <template <typename, std::size_t> class array_t, typename value_t,
+            std::size_t N>
+  ALGEBRA_HOST inline decltype(auto) operator()(
+      const matrix44<array_t, value_t, N> &m, std::size_t row,
+      std::size_t col) const {
 
     // Make sure that the indices are valid.
     assert(row < 4u);
@@ -108,9 +127,12 @@ struct element_getter {
     }
   }
 
-  /// Get const access to a matrix element
-  ALGEBRA_HOST inline auto &operator()(matrix44<scalar_t> &m, std::size_t row,
-                                       std::size_t col) const {
+  /// Get non-const access to a matrix element
+  template <template <typename, std::size_t> class array_t, typename value_t,
+            std::size_t N>
+  ALGEBRA_HOST inline decltype(auto) operator()(
+      matrix44<array_t, value_t, N> &m, std::size_t row,
+      std::size_t col) const {
 
     // Make sure that the indices are valid.
     assert(row < 4u);
@@ -135,6 +157,18 @@ struct element_getter {
     }
   }
 
+  /// Get const access to a matrix element
+  template <template <typename, std::size_t> class array_t, typename value_t>
+  ALGEBRA_HOST inline decltype(auto) operator()(const array_t<value_t, 16> &m,
+                                                unsigned int row,
+                                                unsigned int col) const {
+    // Make sure that the indices are valid.
+    assert(row < 4);
+    assert(col < 4);
+    // Return the selected element.
+    return m[row * 4 + col];
+  }
+
 };  // struct element_getter
 
-}  // namespace algebra::vc_soa
+}  // namespace algebra::storage
diff --git a/storage/common/include/algebra/storage/vector.hpp b/storage/common/include/algebra/storage/vector.hpp
index f2eec8fd..84577ea8 100644
--- a/storage/common/include/algebra/storage/vector.hpp
+++ b/storage/common/include/algebra/storage/vector.hpp
@@ -84,6 +84,8 @@ class vector {
   constexpr operator const array_type &() const { return m_data; }
   /// @}
 
+  constexpr const auto &get() const { return m_data; }
+
   /// Subscript operator[]
   /// @{
   constexpr decltype(auto) operator[](std::size_t i) { return m_data[i]; }
diff --git a/storage/vc_aos/include/algebra/storage/vc_aos.hpp b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
index 3c9a2061..dc811e43 100644
--- a/storage/vc_aos/include/algebra/storage/vc_aos.hpp
+++ b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/storage/matrix44.hpp"
 #include "algebra/storage/vector.hpp"
 
 // System include(s).
diff --git a/storage/vc_soa/CMakeLists.txt b/storage/vc_soa/CMakeLists.txt
index cb1e761e..28e56a78 100644
--- a/storage/vc_soa/CMakeLists.txt
+++ b/storage/vc_soa/CMakeLists.txt
@@ -6,8 +6,7 @@
 
 # Set up the library.
 algebra_add_library( algebra_vc_soa_storage vc_soa_storage
-   "include/algebra/storage/vc_soa.hpp"
-   "include/algebra/storage/impl/vc_soa_matrix44.hpp"  )
+   "include/algebra/storage/vc_soa.hpp" )
 target_link_libraries( algebra_vc_soa_storage
    INTERFACE algebra::common algebra::common_storage Vc::Vc )
 algebra_test_public_headers( algebra_vc_soa_storage