diff --git a/benchmarks/common/python/plot_benchmarks.py b/benchmarks/common/python/plot_benchmarks.py
new file mode 100644
index 00000000..05de102e
--- /dev/null
+++ b/benchmarks/common/python/plot_benchmarks.py
@@ -0,0 +1,167 @@
+# Algebra plugins library, part of the ACTS project
+#
+# (c) 2024 CERN for the benefit of the ACTS project
+#
+# Mozilla Public License Version 2.0 
+
+# python includes
+import argparse
+import logging
+import json
+import pandas as pd
+import os
+import sys
+from datetime import datetime
+import pathlib
+
+""" Parse commandline arguments """
+def parse_arguments():
+#----------------------------------------------------------------arg parsing
+
+    parser = argparse.ArgumentParser(description = "Compare benchmarks")
+    parser.add_argument("--debug", "-d",
+                        help=("Enables debug output"), 
+                        action="store_true")
+    parser.add_argument("--logfile",
+                        help=("Write log in file"), 
+                        default = "", type=str)
+    parser.add_argument("--input", "-i", nargs='+',
+                        help=("Input material scan data file."),
+                        default = "", type=str)
+    parser.add_argument("--outdir", "-o",
+                        help=("Output directory for plots."),
+                        default = "./material_plots/", type=str)
+    parser.add_argument("--output_format", "-of",
+                        help=("Format of the plot files (svg|png|pdf)."),
+                        default = "png", type=str)
+
+    args = parser.parse_args()
+
+#---------------------------------------------------------------------config
+
+    # Check output path
+    if not os.path.isdir(args.outdir):
+        os.mkdir(args.outdir, 0o755)
+    outdir = args.outdir
+
+    # Check input data files
+    if len(args.input) == 0:
+        logging.error(f"Please specify an input data file!")
+        sys.exit(1)
+
+    for input_file in args.input:
+        if not os.path.isfile(input_file):
+            logging.error(f"Data file does not exist! ({input_file})")
+            sys.exit(1)
+
+    if not args.output_format in ["svg", "png", "pdf"]:
+        logging.error(f"Unknown output file format: {args.output_format}")
+        sys.exit(1)
+
+    # Set log level
+    logLevel = logging.INFO
+    if args.debug:
+        logLevel = logging.DEBUG
+
+    # Check logfile path
+    if args.logfile != "":
+        logDirName  = os.path.dirname(args.logfile)
+
+        if logDirName != "" and not os.path.isdir(logDirName):
+            os.mkdir(logDirName, 0o755)
+
+        if not os.path.isfile(args.logfile):
+            with open(args.logfile, 'x'): pass
+
+        # Write log in logfile
+        logging.basicConfig(filename=args.logfile, 
+                            format=("%(levelname)s (%(module)s):"
+                                    " %(message)s"), level=logLevel)
+    else:
+        # Write log to terminal
+        logging.basicConfig(format=("%(levelname)s (%(module)s):"
+                                    " %(message)s"), level=logLevel)
+
+    logging.info("\n--------------------------------------------------------\n"
+                 "Plotting benchmarks - "+\
+                 str(datetime.now().strftime("%d/%m/%Y %H:%M"))+\
+                 "\n->files: " + str(len(args.input)) +\
+                 "\n--------------------------------------------------------\n")
+
+    return args
+
+#----------------------------------------------------------------prepare data
+
+"""Read benchmark data into pandas frame"""
+def read_data(args):
+
+    benchmarks = {}
+    for input_file in args.input:
+        with open(input_file, 'r') as file:
+            extension = pathlib.Path(input_file).suffix
+            name = input_file.removesuffix(extension)
+            try:
+                json_data = json.load(file)
+                benchmarks[name] = pd.DataFrame(json_data["benchmarks"])
+            except ValueError:
+                logging.error('Could not parse the benchmark data')
+                exit(1)
+
+        benchmarks[name]['label'] = benchmarks[name]['name'].apply(lambda x: x.split("/")[0].split("_")[1:])
+
+    return benchmarks
+
+
+"""Display the benchmark data"""
+def plot_benchmarks(data, args):
+    x = np.arange(len(data[0] / 3))  # Always three entries per benchmark case
+    width = 0.15  # the width of the bars
+
+    # Reorganize into benchmark families
+    families = {}
+    for plugin in data:
+        for index, row in df.sort_values(by=['family_index']).iterrows():
+            print(row['familiy_index'])
+
+
+        bars_eigen  = ax.bar(x - 1.5*width, eigen, width, label='eigen', yerr=eigen_err)
+        bars_AoS    = ax.bar(x - 0.5*width, AoS, width, label='AoS', yerr=AoS_err)
+        bars_hybrid = ax.bar(x + 0.5*width, hybrid, width, label='hybrid', yerr=hybrid_err)
+        bars_horiz  = ax.bar(x + 1.5*width, horiz, width, label='horiz', yerr=horiz_err)
+
+    
+    fig, ax = plt.subplots()
+
+    # Add some text for labels, title and custom x-axis tick labels, etc.
+    ax.set_ylabel('speedup')
+    ax.set_title('speed up relative to Eigen sse(double) - preliminary')
+    ax.set_xticks(x)
+    ax.set_xticklabels(labels)
+    ax.legend()
+
+    autolabel(bars_eigen, ax)
+    autolabel(bars_AoS, ax)
+    autolabel(bars_hybrid, ax)
+    autolabel(bars_horiz, ax)
+
+    fig.tight_layout()
+
+    plot_file_name = "./aggregate_plots/speedup_comparision.png"
+    fig.savefig(plot_file_name, dpi=100)
+    fig.clf()
+    plt.close(fig)
+
+
+def main():
+    args = parse_arguments()
+    data = read_data(args)
+    plot_benchmarks(data, args)
+
+    print(data['vc_aos_vector_bench'])
+
+#-------------------------------------------------------------------------------
+
+if __name__ == "__main__":
+    main()
+
+#------------------------------------------------------------------------------- 
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