diff --git a/benchmarks/common/include/benchmark/common/benchmark_getter.hpp b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
index 30db7232..c6c9c9ae 100644
--- a/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
+++ b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
@@ -11,4 +11,5 @@
 #include "benchmark_vector.hpp"
 #include "register_benchmark.hpp"
 
-namespace algebra::bench_op {}  // namespace algebra::bench_op
+namespace algebra::bench_op { /* TODO: Implement */
+}  // namespace algebra::bench_op
diff --git a/common/include/algebra/type_traits.hpp b/common/include/algebra/type_traits.hpp
index 53adf683..5eac0c3e 100644
--- a/common/include/algebra/type_traits.hpp
+++ b/common/include/algebra/type_traits.hpp
@@ -61,7 +61,7 @@ template <class M>
 struct matrix {};
 
 template <class M, index_t<M> ROWS, index_t<M> COLS, typename T>
-using other_matrix_t = typename matrix<M>::template other_type<T, ROWS, COLS>;
+using get_matrix_t = typename matrix<M>::template other_type<T, ROWS, COLS>;
 /// @}
 
 /// Matrix dimensions
@@ -71,8 +71,7 @@ struct dimensions {};
 
 /// Specilization for scalar types
 template <typename M>
-  requires std::is_fundamental_v<M>
-struct dimensions<M> {
+requires std::is_fundamental_v<M> struct dimensions<M> {
 
   using size_type = std::size_t;
 
@@ -121,6 +120,101 @@ template <class M>
 using block_getter_t = typename block_getter<M>::type;
 /// @}
 
-/// @}
-
 }  // namespace algebra::trait
+
+/// Default type trait specializations
+/// @{
+#define ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(A)                           \
+                                                                        \
+  namespace trait {                                                     \
+                                                                        \
+  template <typename T, auto N>                                         \
+  struct index<A::vector_type<T, N>> {                                  \
+    using type = algebra::A::size_type;                                 \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto ROWS, auto COLS>                           \
+  struct index<A::matrix_type<T, ROWS, COLS>> {                         \
+    using type = algebra::A::size_type;                                 \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto N>                                         \
+  struct dimensions<A::vector_type<T, N>> {                             \
+                                                                        \
+    using size_type = index_t<A::vector_type<T, N>>;                    \
+                                                                        \
+    static constexpr size_type dim{1};                                  \
+    static constexpr size_type rows{N};                                 \
+    static constexpr size_type columns{1};                              \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto ROWS, auto COLS>                           \
+  struct dimensions<A::matrix_type<T, ROWS, COLS>> {                    \
+                                                                        \
+    using size_type = index_t<A::matrix_type<T, ROWS, COLS>>;           \
+                                                                        \
+    static constexpr size_type dim{2};                                  \
+    static constexpr size_type rows{ROWS};                              \
+    static constexpr size_type columns{COLS};                           \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto N>                                         \
+  struct value<A::vector_type<T, N>> {                                  \
+    using type = T;                                                     \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto ROWS, auto COLS>                           \
+  struct value<A::matrix_type<T, ROWS, COLS>> {                         \
+    using type = T;                                                     \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto N>                                         \
+  struct vector<A::vector_type<T, N>> {                                 \
+                                                                        \
+    template <typename other_T, auto other_N>                           \
+    using other_type = A::vector_type<other_T, other_N>;                \
+                                                                        \
+    using type = other_type<T, N>;                                      \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto ROWS, auto COLS>                           \
+  struct vector<A::matrix_type<T, ROWS, COLS>> {                        \
+                                                                        \
+    template <typename other_T, auto other_N>                           \
+    using other_type = A::vector_type<other_T, other_N>;                \
+                                                                        \
+    using type = other_type<T, ROWS>;                                   \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto ROWS, auto COLS>                           \
+  struct matrix<A::matrix_type<T, ROWS, COLS>> {                        \
+    template <typename other_T, auto other_ROWS, auto other_COLS>       \
+    using other_type = A::matrix_type<other_T, other_ROWS, other_COLS>; \
+                                                                        \
+    using type = A::matrix_type<T, ROWS, COLS>;                         \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto N>                                         \
+  struct matrix<A::vector_type<T, N>> {                                 \
+    template <typename other_T, auto other_ROWS, auto other_COLS>       \
+    using other_type = A::matrix_type<other_T, other_ROWS, other_COLS>; \
+                                                                        \
+    using type = other_type<T, N, 1>;                                   \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto N>                                         \
+  struct element_getter<A::vector_type<T, N>> {                         \
+    using type = A::element_getter;                                     \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto ROWS, auto COLS>                           \
+  struct element_getter<A::matrix_type<T, ROWS, COLS>> {                \
+    using type = A::element_getter;                                     \
+  };                                                                    \
+                                                                        \
+  template <typename T, auto ROWS, auto COLS>                           \
+  struct block_getter<A::matrix_type<T, ROWS, COLS>> {                  \
+    using type = A::block_getter;                                       \
+  };                                                                    \
+                                                                        \
+  }  // namespace algebra::trait
diff --git a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp
index a5df35f4..67ea5317 100644
--- a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp
+++ b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp
@@ -37,7 +37,7 @@ requires(std::is_scalar_v<typename matrix_t::value_type::value_type>)
 /// Set @param m as zero matrix
 template <std::size_t ROWS, std::size_t COLS, typename scalar_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline constexpr void set_zero(
+ALGEBRA_HOST_DEVICE constexpr void set_zero(
     array_t<array_t<scalar_t, ROWS>, COLS> &m) {
   m = zero<array_t<array_t<scalar_t, ROWS>, COLS>>();
 }
@@ -45,7 +45,7 @@ ALGEBRA_HOST_DEVICE inline constexpr void set_zero(
 /// Set @param m as identity matrix
 template <std::size_t ROWS, std::size_t COLS, typename scalar_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline constexpr void set_identity(
+ALGEBRA_HOST_DEVICE constexpr void set_identity(
     array_t<array_t<scalar_t, ROWS>, COLS> &m) {
   m = identity<array_t<array_t<scalar_t, ROWS>, COLS>>();
 }
diff --git a/math/cmath/include/algebra/math/impl/cmath_vector.hpp b/math/cmath/include/algebra/math/impl/cmath_vector.hpp
index 12a6866e..a4aad3e6 100644
--- a/math/cmath/include/algebra/math/impl/cmath_vector.hpp
+++ b/math/cmath/include/algebra/math/impl/cmath_vector.hpp
@@ -20,10 +20,9 @@ namespace algebra::cmath {
 ///
 /// @return the scalar dot product value
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline scalar_t dot(const array_t<scalar_t, N> &a,
-                                        const array_t<scalar_t, N> &b) {
+          typename scalar_t, size_type N>
+requires std::is_scalar_v<scalar_t> ALGEBRA_HOST_DEVICE inline scalar_t dot(
+    const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) {
   array_t<scalar_t, N> tmp;
   for (size_type i = 0; i < N; i++) {
     tmp[i] = a[i] * b[i];
@@ -41,13 +40,11 @@ ALGEBRA_HOST_DEVICE inline scalar_t dot(const array_t<scalar_t, N> &a,
 /// @param b the second input matrix with single column
 ///
 /// @return the scalar dot product value
-template <
-    typename size_type, template <typename, size_type> class array_t,
-    typename scalar_t, size_type N, size_type COLS,
-    std::enable_if_t<COLS == 1 && std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline scalar_t dot(
-    const array_t<scalar_t, N> &a,
-    const array_t<array_t<scalar_t, N>, COLS> &b) {
+template <typename size_type, template <typename, size_type> class array_t,
+          typename scalar_t, size_type N, size_type COLS>
+requires(COLS == 1 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
+    inline scalar_t dot(const array_t<scalar_t, N> &a,
+                        const array_t<array_t<scalar_t, N>, COLS> &b) {
   array_t<scalar_t, N> tmp;
   for (size_type i = 0; i < N; i++) {
     tmp[i] = a[i] * b[0][i];
@@ -65,13 +62,11 @@ ALGEBRA_HOST_DEVICE inline scalar_t dot(
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <
-    typename size_type, template <typename, size_type> class array_t,
-    typename scalar_t, size_type N, size_type COLS,
-    std::enable_if_t<COLS == 1 && std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline scalar_t dot(
-    const array_t<array_t<scalar_t, N>, COLS> &a,
-    const array_t<scalar_t, N> &b) {
+template <typename size_type, template <typename, size_type> class array_t,
+          typename scalar_t, size_type N, size_type COLS>
+requires(COLS == 1 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
+    inline scalar_t dot(const array_t<array_t<scalar_t, N>, COLS> &a,
+                        const array_t<scalar_t, N> &b) {
   array_t<scalar_t, N> tmp;
   for (size_type i = 0; i < N; i++) {
     tmp[i] = a[0][i] * b[i];
@@ -89,13 +84,11 @@ ALGEBRA_HOST_DEVICE inline scalar_t dot(
 /// @param b the second input matrix with single column
 ///
 /// @return the scalar dot product value
-template <
-    typename size_type, template <typename, size_type> class array_t,
-    typename scalar_t, size_type N, size_type COLS,
-    std::enable_if_t<COLS == 1 && std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline scalar_t dot(
-    const array_t<array_t<scalar_t, N>, COLS> &a,
-    const array_t<array_t<scalar_t, N>, COLS> &b) {
+template <typename size_type, template <typename, size_type> class array_t,
+          typename scalar_t, size_type N, size_type COLS>
+requires(COLS == 1 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
+    inline scalar_t dot(const array_t<array_t<scalar_t, N>, COLS> &a,
+                        const array_t<array_t<scalar_t, N>, COLS> &b) {
   array_t<scalar_t, N> tmp;
   for (size_type i = 0; i < N; i++) {
     tmp[i] = a[0][i] * b[0][i];
@@ -111,9 +104,9 @@ ALGEBRA_HOST_DEVICE inline scalar_t dot(
 ///
 /// @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) && std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline scalar_t norm(const array_t<scalar_t, N> &v) {
+          typename scalar_t, size_type N>
+requires(N >= 2 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
+    inline scalar_t norm(const array_t<scalar_t, N> &v) {
 
   return algebra::math::sqrt(dot(v, v));
 }
@@ -123,10 +116,9 @@ 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) && std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline scalar_t eta(
-    const array_t<scalar_t, N> &v) noexcept {
+          typename scalar_t, size_type N>
+requires(N >= 3 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
+    inline scalar_t eta(const array_t<scalar_t, N> &v) noexcept {
 
   return algebra::math::atanh(v[2] / norm(v));
 }
@@ -135,10 +127,10 @@ ALGEBRA_HOST_DEVICE inline scalar_t eta(
 ///
 /// @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<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, N> normalize(
-    const array_t<scalar_t, N> &v) {
+          typename scalar_t, size_type N>
+requires std::is_scalar_v<scalar_t>
+    ALGEBRA_HOST_DEVICE inline array_t<scalar_t, N> normalize(
+        const array_t<scalar_t, N> &v) {
 
   return (static_cast<scalar_t>(1.) / norm(v)) * v;
 }
diff --git a/math/eigen/include/algebra/math/impl/eigen_vector.hpp b/math/eigen/include/algebra/math/impl/eigen_vector.hpp
index bf8cbd41..3486cf81 100644
--- a/math/eigen/include/algebra/math/impl/eigen_vector.hpp
+++ b/math/eigen/include/algebra/math/impl/eigen_vector.hpp
@@ -41,12 +41,10 @@ ALGEBRA_HOST_DEVICE inline auto norm(const Eigen::MatrixBase<derived_type> &v) {
 /// rows >= 3
 ///
 /// @param v the input vector
-template <
-    typename derived_type,
-    std::enable_if_t<Eigen::MatrixBase<derived_type>::RowsAtCompileTime >= 3,
-                     bool> = true>
-ALGEBRA_HOST_DEVICE inline auto eta(
-    const Eigen::MatrixBase<derived_type> &v) noexcept {
+template <typename derived_type>
+requires(Eigen::MatrixBase<derived_type>::RowsAtCompileTime >=
+         3) ALGEBRA_HOST_DEVICE
+    inline auto eta(const Eigen::MatrixBase<derived_type> &v) noexcept {
 
   return algebra::math::atanh(v[2] / v.norm());
 }
diff --git a/math/fastor/include/algebra/math/impl/fastor_vector.hpp b/math/fastor/include/algebra/math/impl/fastor_vector.hpp
index b6f4df2c..08c2d603 100644
--- a/math/fastor/include/algebra/math/impl/fastor_vector.hpp
+++ b/math/fastor/include/algebra/math/impl/fastor_vector.hpp
@@ -24,9 +24,9 @@ namespace algebra::fastor::math {
 /// This method retrieves theta from a vector, vector base with rows >= 3
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N, std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST inline scalar_t theta(
-    const Fastor::Tensor<scalar_t, N> &v) noexcept {
+template <typename scalar_t, auto N>
+requires(N >= 3) ALGEBRA_HOST inline scalar_t
+    theta(const Fastor::Tensor<scalar_t, N> &v) noexcept {
 
   return algebra::math::atan2(Fastor::norm(v(Fastor::fseq<0, 2>())), v[2]);
 }
@@ -34,9 +34,9 @@ ALGEBRA_HOST inline scalar_t theta(
 /// This method retrieves the perpenticular magnitude of a vector with rows >= 2
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N, std::enable_if_t<N >= 2, bool> = true>
-ALGEBRA_HOST inline scalar_t perp(
-    const Fastor::Tensor<scalar_t, N> &v) noexcept {
+template <typename scalar_t, auto N>
+requires(N >= 2) ALGEBRA_HOST inline scalar_t
+    perp(const Fastor::Tensor<scalar_t, N> &v) noexcept {
 
   return algebra::math::sqrt(
       Fastor::inner(v(Fastor::fseq<0, 2>()), v(Fastor::fseq<0, 2>())));
@@ -55,9 +55,9 @@ ALGEBRA_HOST inline scalar_t norm(const Fastor::Tensor<scalar_t, N> &v) {
 /// rows >= 3
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N, std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST inline scalar_t eta(
-    const Fastor::Tensor<scalar_t, N> &v) noexcept {
+template <typename scalar_t, auto N>
+requires(N >= 3) ALGEBRA_HOST inline scalar_t
+    eta(const Fastor::Tensor<scalar_t, N> &v) noexcept {
 
   return algebra::math::atanh(v[2] / Fastor::norm(v));
 }
diff --git a/math/generic/include/algebra/math/impl/generic_matrix.hpp b/math/generic/include/algebra/math/impl/generic_matrix.hpp
index c460e820..41c49864 100644
--- a/math/generic/include/algebra/math/impl/generic_matrix.hpp
+++ b/math/generic/include/algebra/math/impl/generic_matrix.hpp
@@ -71,7 +71,7 @@ ALGEBRA_HOST_DEVICE inline auto transpose(const matrix_t &m) {
   constexpr index_t rows{algebra::trait::rows<matrix_t>};
   constexpr index_t columns{algebra::trait::columns<matrix_t>};
 
-  algebra::trait::other_matrix_t<matrix_t, columns, rows, value_t> ret;
+  algebra::trait::get_matrix_t<matrix_t, columns, rows, value_t> ret;
 
   for (index_t i = 0; i < rows; ++i) {
     for (index_t j = 0; j < columns; ++j) {
diff --git a/math/generic/include/algebra/math/impl/generic_transform3.hpp b/math/generic/include/algebra/math/impl/generic_transform3.hpp
index 3b14fd82..ee740e3b 100644
--- a/math/generic/include/algebra/math/impl/generic_transform3.hpp
+++ b/math/generic/include/algebra/math/impl/generic_transform3.hpp
@@ -138,9 +138,7 @@ struct transform3 {
   ///
   /// @param m is the full 4x4 matrix
   ALGEBRA_HOST_DEVICE
-  explicit transform3(const matrix44 &m) {
-
-    _data = m;
+  explicit transform3(const matrix44 &m) : _data{m} {
     _data_inv = matrix_inversion{}(_data);
   }
 
diff --git a/math/generic/include/algebra/math/impl/generic_vector.hpp b/math/generic/include/algebra/math/impl/generic_vector.hpp
index 7261ec54..51256f95 100644
--- a/math/generic/include/algebra/math/impl/generic_vector.hpp
+++ b/math/generic/include/algebra/math/impl/generic_vector.hpp
@@ -17,10 +17,10 @@ namespace algebra::generic::math {
 /// This method retrieves phi from a vector with rows >= 2
 ///
 /// @param v the input vector
-template <typename vector_t,
-          std::enable_if_t<algebra::trait::is_vector<vector_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> phi(
-    const vector_t &v) noexcept {
+template <typename vector_t>
+requires algebra::trait::is_vector<vector_t>
+    ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> phi(
+        const vector_t &v) noexcept {
 
   using element_getter_t = algebra::trait::element_getter_t<vector_t>;
 
@@ -31,10 +31,10 @@ ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> phi(
 /// This method retrieves the perpendicular magnitude of a vector with rows >= 2
 ///
 /// @param v the input vector
-template <typename vector_t,
-          std::enable_if_t<algebra::trait::is_vector<vector_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> perp(
-    const vector_t &v) noexcept {
+template <typename vector_t>
+requires algebra::trait::is_vector<vector_t>
+    ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> perp(
+        const vector_t &v) noexcept {
 
   using element_getter_t = algebra::trait::element_getter_t<vector_t>;
 
@@ -46,10 +46,10 @@ ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> perp(
 /// This method retrieves theta from a vector with rows >= 3
 ///
 /// @param v the input vector
-template <typename vector_t,
-          std::enable_if_t<algebra::trait::is_vector<vector_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> theta(
-    const vector_t &v) noexcept {
+template <typename vector_t>
+requires algebra::trait::is_vector<vector_t>
+    ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> theta(
+        const vector_t &v) noexcept {
 
   using element_getter_t = algebra::trait::element_getter_t<vector_t>;
 
@@ -65,12 +65,11 @@ ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> theta(
 /// @param b the second input vector
 ///
 /// @return a vector representing the cross product
-template <typename vector1_t, typename vector2_t,
-          std::enable_if_t<(algebra::trait::columns<vector1_t> == 1 &&
-                            algebra::trait::columns<vector2_t> == 1),
-                           bool> = true>
-ALGEBRA_HOST_DEVICE inline algebra::trait::vector_t<vector1_t> cross(
-    const vector1_t &a, const vector2_t &b) {
+template <typename vector1_t, typename vector2_t>
+requires(algebra::trait::columns<vector1_t> == 1 &&
+         algebra::trait::columns<vector2_t> == 1) ALGEBRA_HOST_DEVICE
+    inline algebra::trait::vector_t<vector1_t> cross(const vector1_t &a,
+                                                     const vector2_t &b) {
 
   using element_getter_t = algebra::trait::element_getter_t<vector1_t>;
 
@@ -92,12 +91,11 @@ ALGEBRA_HOST_DEVICE inline algebra::trait::vector_t<vector1_t> cross(
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename vector1_t, typename vector2_t,
-          std::enable_if_t<(algebra::trait::columns<vector1_t> == 1 &&
-                            algebra::trait::columns<vector2_t> == 1),
-                           bool> = true>
-ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector1_t> dot(
-    const vector1_t &a, const vector2_t &b) {
+template <typename vector1_t, typename vector2_t>
+requires(algebra::trait::columns<vector1_t> == 1 &&
+         algebra::trait::columns<vector2_t> == 1) ALGEBRA_HOST_DEVICE
+    inline algebra::trait::scalar_t<vector1_t> dot(const vector1_t &a,
+                                                   const vector2_t &b) {
 
   using scalar_t = algebra::trait::scalar_t<vector1_t>;
   using index_t = algebra::trait::index_t<vector1_t>;
@@ -115,10 +113,10 @@ ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector1_t> dot(
 /// This method retrieves the norm of a vector with rows >= 3
 ///
 /// @param v the input vector
-template <typename vector_t,
-          std::enable_if_t<algebra::trait::is_vector<vector_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> norm(
-    const vector_t &v) {
+template <typename vector_t>
+requires algebra::trait::is_vector<vector_t>
+    ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> norm(
+        const vector_t &v) {
 
   return algebra::math::sqrt(dot(v, v));
 }
@@ -127,10 +125,10 @@ ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> norm(
 /// rows >= 3
 ///
 /// @param v the input vector
-template <typename vector_t,
-          std::enable_if_t<algebra::trait::is_vector<vector_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> eta(
-    const vector_t &v) noexcept {
+template <typename vector_t>
+requires algebra::trait::is_vector<vector_t>
+    ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<vector_t> eta(
+        const vector_t &v) noexcept {
 
   using element_getter_t = algebra::trait::element_getter_t<vector_t>;
 
diff --git a/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp b/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp
index 7989fc52..bb5fca1d 100644
--- a/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp
+++ b/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp
@@ -21,17 +21,16 @@ namespace algebra::smatrix::math {
 /// This method retrieves phi from a vector, vector base with rows >= 2
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N, std::enable_if_t<N >= 2, bool> = true>
-ALGEBRA_HOST inline scalar_t phi(
-    const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
+template <typename scalar_t, auto N>
+requires(N >= 2) ALGEBRA_HOST inline scalar_t
+    phi(const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(TMath::ATan2(v[1], v[0]));
 }
 
-template <typename scalar_t, class A, auto N,
-          std::enable_if_t<N >= 2, bool> = true>
-ALGEBRA_HOST inline scalar_t phi(
-    const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
+template <typename scalar_t, class A, auto N>
+requires(N >= 2) ALGEBRA_HOST inline scalar_t
+    phi(const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(TMath::ATan2(v.apply(1), v.apply(0)));
 }
@@ -39,18 +38,17 @@ ALGEBRA_HOST inline scalar_t phi(
 /// This method retrieves theta from a vector, vector base with rows >= 3
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N, std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST inline scalar_t theta(
-    const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
+template <typename scalar_t, auto N>
+requires(N >= 3) ALGEBRA_HOST inline scalar_t
+    theta(const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(
       TMath::ATan2(TMath::Sqrt(v[0] * v[0] + v[1] * v[1]), v[2]));
 }
 
-template <typename scalar_t, class A, auto N,
-          std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST inline scalar_t theta(
-    const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
+template <typename scalar_t, class A, auto N>
+requires(N >= 3) ALGEBRA_HOST inline scalar_t
+    theta(const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(TMath::ATan2(
       TMath::Sqrt(v.apply(0) * v.apply(0) + v.apply(1) * v.apply(1)),
@@ -77,17 +75,16 @@ ALGEBRA_HOST inline scalar_t norm(
 /// rows >= 3
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N, std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST inline scalar_t eta(
-    const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
+template <typename scalar_t, auto N>
+requires(N >= 3) ALGEBRA_HOST inline scalar_t
+    eta(const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(TMath::ATanH(v[2] / norm(v)));
 }
 
-template <typename scalar_t, class A, auto N,
-          std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST inline scalar_t eta(
-    const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
+template <typename scalar_t, class A, auto N>
+requires(N >= 3) ALGEBRA_HOST inline scalar_t
+    eta(const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(TMath::ATanH(v.apply(2) / norm(v)));
 }
@@ -95,17 +92,16 @@ ALGEBRA_HOST inline scalar_t eta(
 /// This method retrieves the perpendicular magnitude of a vector with rows >= 2
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N, std::enable_if_t<N >= 2, bool> = true>
-ALGEBRA_HOST inline scalar_t perp(
-    const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
+template <typename scalar_t, auto N>
+requires(N >= 2) ALGEBRA_HOST inline scalar_t
+    perp(const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(TMath::Sqrt(v[0] * v[0] + v[1] * v[1]));
 }
 
-template <typename scalar_t, class A, auto N,
-          std::enable_if_t<N >= 2, bool> = true>
-ALGEBRA_HOST inline scalar_t perp(
-    const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
+template <typename scalar_t, class A, auto N>
+requires(N >= 2) ALGEBRA_HOST inline scalar_t
+    perp(const ROOT::Math::VecExpr<A, scalar_t, N> &v) noexcept {
 
   return static_cast<scalar_t>(
       TMath::Sqrt(v.apply(0) * v.apply(0) + v.apply(1) * v.apply(1)));
diff --git a/math/smatrix/include/algebra/math/smatrix.hpp b/math/smatrix/include/algebra/math/smatrix.hpp
index 418bb931..e55507f4 100644
--- a/math/smatrix/include/algebra/math/smatrix.hpp
+++ b/math/smatrix/include/algebra/math/smatrix.hpp
@@ -8,7 +8,6 @@
 #pragma once
 
 // Project include(s).
-#include "algebra/math/impl/smatrix_getter.hpp"
 #include "algebra/math/impl/smatrix_matrix.hpp"
 #include "algebra/math/impl/smatrix_transform3.hpp"
 #include "algebra/math/impl/smatrix_vector.hpp"
diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp
index e8473f75..07fdf65c 100644
--- a/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp
+++ b/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp
@@ -30,9 +30,8 @@ namespace algebra::vc_soa::math {
 ///
 /// @param v the input vector
 template <std::size_t N, typename scalar_t,
-          template <typename, std::size_t> class array_t,
-          std::enable_if_t<N >= 2, bool> = true>
-ALGEBRA_HOST_DEVICE inline auto phi(
+          template <typename, std::size_t> class array_t>
+requires(N >= 2) ALGEBRA_HOST_DEVICE inline auto phi(
     const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
 
   return Vc::atan2(v[1], v[0]);
@@ -46,9 +45,8 @@ ALGEBRA_HOST_DEVICE inline auto phi(
 ///
 /// @param v the input vector
 template <std::size_t N, typename scalar_t,
-          template <typename, std::size_t> class array_t,
-          std::enable_if_t<N >= 2, bool> = true>
-ALGEBRA_HOST_DEVICE inline auto perp(
+          template <typename, std::size_t> class array_t>
+requires(N >= 2) ALGEBRA_HOST_DEVICE inline auto perp(
     const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
 
   return Vc::sqrt(Vc::fma(v[0], v[0], v[1] * v[1]));
@@ -62,9 +60,8 @@ ALGEBRA_HOST_DEVICE inline auto perp(
 ///
 /// @param v the input vector
 template <std::size_t N, typename scalar_t,
-          template <typename, std::size_t> class array_t,
-          std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST_DEVICE inline auto theta(
+          template <typename, std::size_t> class array_t>
+requires(N >= 3) ALGEBRA_HOST_DEVICE inline auto theta(
     const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
 
   return Vc::atan2(perp(v), v[2]);
@@ -81,11 +78,11 @@ ALGEBRA_HOST_DEVICE inline auto theta(
 ///
 /// @return a vector (expression) representing the cross product
 template <std::size_t N, typename scalar_t,
-          template <typename, std::size_t> class array_t,
-          std::enable_if_t<N == 3, bool> = true>
-ALGEBRA_HOST_DEVICE inline storage::vector<N, Vc::Vector<scalar_t>, array_t>
-cross(const storage::vector<N, Vc::Vector<scalar_t>, array_t> &a,
-      const storage::vector<N, Vc::Vector<scalar_t>, array_t> &b) {
+          template <typename, std::size_t> class array_t>
+requires(N == 3) ALGEBRA_HOST_DEVICE
+    inline storage::vector<N, Vc::Vector<scalar_t>, array_t> cross(
+        const storage::vector<N, Vc::Vector<scalar_t>, array_t> &a,
+        const storage::vector<N, Vc::Vector<scalar_t>, array_t> &b) {
 
   return {Vc::fma(a[1], b[2], -b[1] * a[2]), Vc::fma(a[2], b[0], -b[2] * a[0]),
           Vc::fma(a[0], b[1], -b[0] * a[1])};
@@ -161,9 +158,8 @@ normalize(const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
 ///
 /// @param v the input vector
 template <std::size_t N, typename scalar_t,
-          template <typename, std::size_t> class array_t,
-          std::enable_if_t<N >= 3, bool> = true>
-ALGEBRA_HOST_DEVICE inline auto eta(
+          template <typename, std::size_t> class array_t>
+requires(N >= 3) ALGEBRA_HOST_DEVICE inline auto eta(
     const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
 
   // atanh does not exist in Vc
diff --git a/storage/array/include/algebra/storage/array.hpp b/storage/array/include/algebra/storage/array.hpp
index 49226ecf..ed79c674 100644
--- a/storage/array/include/algebra/storage/array.hpp
+++ b/storage/array/include/algebra/storage/array.hpp
@@ -44,120 +44,13 @@ using vector2 = storage_type<T, 2>;
 template <typename T>
 using point2 = vector2<T>;
 
-}  // namespace array
-
-namespace trait {
-
-/// Type trait specializations
-/// @{
-
-/// Index
-/// @{
-template <typename T, std::size_t N>
-struct index<std::array<T, N>> {
-  using type = algebra::array::size_type;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct index<std::array<std::array<T, ROWS>, COLS>> {
-  using type = algebra::array::size_type;
-};
-/// @}
-
-/// Dimension
-/// @{
-template <typename T, std::size_t N>
-struct dimensions<std::array<T, N>> {
-
-  using size_type = index_t<std::array<T, N>>;
-
-  static constexpr size_type dim{1};
-  static constexpr size_type rows{N};
-  static constexpr size_type columns{1};
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct dimensions<std::array<std::array<T, ROWS>, COLS>> {
-
-  using size_type = index_t<std::array<std::array<T, ROWS>, COLS>>;
-
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-/// @}
-
-/// Value
-/// @{
-template <typename T, std::size_t N>
-struct value<std::array<T, N>> {
-  using type = T;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct value<std::array<std::array<T, ROWS>, COLS>> {
-  using type = T;
-};
-/// @}
+/// Element Getter
+using element_getter = cmath::storage::element_getter;
+/// Block Getter
+using block_getter = cmath::storage::block_getter;
 
-/// Vector
-/// @{
-template <typename T, std::size_t N>
-struct vector<std::array<T, N>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type = std::array<other_T, other_N>;
-
-  using type = other_type<T, N>;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct vector<std::array<std::array<T, ROWS>, COLS>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type = std::array<other_T, other_N>;
-
-  using type = other_type<T, ROWS>;
-};
-/// @}
-
-/// Matrix
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct matrix<std::array<std::array<T, ROWS>, COLS>> {
-  template <typename other_T, std::size_t other_ROWS, std::size_t other_COLS>
-  using other_type = std::array<std::array<other_T, other_ROWS>, other_COLS>;
-
-  using type = std::array<std::array<T, ROWS>, COLS>;
-};
-
-template <typename T, int N>
-struct matrix<std::array<T, N>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type = std::array<std::array<other_T, other_ROWS>, other_COLS>;
-
-  using type = other_type<T, N, 1>;
-};
-/// @}
-
-/// Elemet/Block Getter
-/// @{
-template <typename T, std::size_t N>
-struct element_getter<std::array<T, N>> {
-  using type = cmath::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct element_getter<std::array<std::array<T, ROWS>, COLS>> {
-  using type = cmath::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct block_getter<std::array<std::array<T, ROWS>, COLS>> {
-  using type = cmath::storage::block_getter;
-};
-/// @}
+}  // namespace array
 
-}  // namespace trait
+ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(array)
 
 }  // namespace algebra
diff --git a/storage/common/include/algebra/storage/matrix.hpp b/storage/common/include/algebra/storage/matrix.hpp
index 58db1629..721d5a1d 100644
--- a/storage/common/include/algebra/storage/matrix.hpp
+++ b/storage/common/include/algebra/storage/matrix.hpp
@@ -31,8 +31,8 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
 
   /// Construct from given column vectors @param v
   template <typename... vector_t>
-    requires(sizeof...(vector_t) == COL)
-  explicit matrix(vector_t &&...v) : m_storage{std::forward<vector_t>(v)...} {}
+  requires(sizeof...(vector_t) == COL) explicit matrix(vector_t &&... v)
+      : m_storage{std::forward<vector_t>(v)...} {}
 
   /// Equality operator between two matrices
   template <std::size_t R, std::size_t C, typename V,
@@ -69,15 +69,15 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
   /// @{
   // AoS
   template <std::size_t... I>
-    requires(!std::is_scalar_v<value_t>)
-  constexpr bool equal(const matrix &rhs, std::index_sequence<I...>) const {
+  requires(!std::is_scalar_v<value_t>) constexpr bool equal(
+      const matrix &rhs, std::index_sequence<I...>) const {
     return (... && (m_storage[I] == rhs[I]));
   }
 
   // SoA
   template <std::size_t... I>
-    requires(std::is_scalar_v<value_t>)
-  constexpr bool equal(const matrix &rhs, std::index_sequence<I...>) const {
+  requires(std::is_scalar_v<value_t>) constexpr bool equal(
+      const matrix &rhs, std::index_sequence<I...>) const {
     return (... && ((m_storage[I].get() == rhs[I].get()).isFull()));
   }
   /// @}
@@ -96,15 +96,16 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
 
   template <std::size_t R, std::size_t C, typename V, typename scalar_t,
             template <typename, std::size_t> class A>
-    requires(std::is_scalar_v<scalar_t> || std::is_same_v<V, scalar_t>)
-  ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator*(
-      scalar_t a, const matrix<A, V, R, C> &rhs) noexcept;
+  requires(std::is_scalar_v<scalar_t> ||
+           std::is_same_v<V, scalar_t>) ALGEBRA_HOST_DEVICE
+      friend constexpr decltype(auto)
+      operator*(scalar_t a, const matrix<A, V, R, C> &rhs) noexcept;
 
   template <std::size_t R, std::size_t C, typename V, typename scalar_t,
             template <typename, std::size_t> class A>
-    requires(std::is_scalar_v<scalar_t> || std::is_same_v<V, scalar_t>)
-  friend constexpr decltype(auto) operator*(const matrix<A, V, R, C> &lhs,
-                                            scalar_t a) noexcept;
+  requires(std::is_scalar_v<scalar_t> ||
+           std::is_same_v<V, scalar_t>) friend constexpr decltype(auto)
+  operator*(const matrix<A, V, R, C> &lhs, scalar_t a) noexcept;
 
   /// Matrix-vector multiplication
   template <std::size_t R, std::size_t C, typename V,
@@ -199,10 +200,11 @@ constexpr bool operator==(const matrix<array_t, value_t, ROW, COL> &lhs,
 
 /// Scalar multiplication
 template <typename matrix_t, typename scalar_t, std::size_t... J>
-  requires(std::is_scalar_v<scalar_t> ||
-           std::is_same_v<typename matrix_t::value_type, scalar_t>)
-ALGEBRA_HOST_DEVICE constexpr matrix_t matrix_scalar_mul(
-    scalar_t a, const matrix_t &rhs, std::index_sequence<J...>) noexcept {
+requires(std::is_scalar_v<scalar_t> ||
+         std::is_same_v<typename matrix_t::value_type, scalar_t>)
+    ALGEBRA_HOST_DEVICE constexpr matrix_t
+    matrix_scalar_mul(scalar_t a, const matrix_t &rhs,
+                      std::index_sequence<J...>) noexcept {
 
   return matrix_t{(a * rhs[J])...};
 }
@@ -250,9 +252,11 @@ ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator-(
 
 template <std::size_t ROW, std::size_t COL, typename value_t, typename scalar_t,
           template <typename, std::size_t> class array_t>
-  requires(std::is_scalar_v<scalar_t> || std::is_same_v<value_t, scalar_t>)
-ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
-    scalar_t a, const matrix<array_t, value_t, ROW, COL> &rhs) noexcept {
+requires(std::is_scalar_v<scalar_t> ||
+         std::is_same_v<value_t, scalar_t>) ALGEBRA_HOST_DEVICE
+    constexpr decltype(auto)
+    operator*(scalar_t a,
+              const matrix<array_t, value_t, ROW, COL> &rhs) noexcept {
 
   using matrix_t = matrix<array_t, value_t, ROW, COL>;
 
@@ -262,9 +266,9 @@ ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
 
 template <std::size_t ROW, std::size_t COL, typename value_t, typename scalar_t,
           template <typename, std::size_t> class array_t>
-  requires(std::is_scalar_v<scalar_t> || std::is_same_v<value_t, scalar_t>)
-constexpr decltype(auto) operator*(
-    const matrix<array_t, value_t, ROW, COL> &lhs, scalar_t a) noexcept {
+requires(std::is_scalar_v<scalar_t> ||
+         std::is_same_v<value_t, scalar_t>) constexpr decltype(auto)
+operator*(const matrix<array_t, value_t, ROW, COL> &lhs, scalar_t a) noexcept {
   return a * lhs;
 }
 
diff --git a/storage/common/include/algebra/storage/matrix_getter.hpp b/storage/common/include/algebra/storage/matrix_getter.hpp
index cd3a0717..bd6db011 100644
--- a/storage/common/include/algebra/storage/matrix_getter.hpp
+++ b/storage/common/include/algebra/storage/matrix_getter.hpp
@@ -144,7 +144,7 @@ struct block_getter {
   template <std::size_t ROWS, std::size_t COLS, std::size_t mROW,
             std::size_t mCOL, typename value_t,
             template <typename, std::size_t> class array_t>
-  ALGEBRA_HOST_DEVICE inline constexpr auto operator()(
+  ALGEBRA_HOST_DEVICE constexpr auto operator()(
       const matrix<array_t, value_t, mROW, mCOL> &m, const std::size_t row,
       const std::size_t col) noexcept {
 
@@ -181,7 +181,7 @@ struct block_getter {
   /// Get a vector of a const matrix
   template <std::size_t SIZE, std::size_t ROWS, std::size_t COLS,
             typename value_t, template <typename, std::size_t> class array_t>
-  ALGEBRA_HOST_DEVICE inline constexpr auto operator()(
+  ALGEBRA_HOST_DEVICE constexpr auto operator()(
       const matrix<array_t, value_t, ROWS, COLS> &m, const std::size_t row,
       const std::size_t col) noexcept {
 
@@ -214,7 +214,7 @@ struct block_getter {
 template <std::size_t ROWS, std::size_t COLS, std::size_t mROW,
           std::size_t mCOL, typename value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline constexpr auto block(
+ALGEBRA_HOST_DEVICE constexpr auto block(
     const matrix<array_t, value_t, mROW, mCOL> &m, const std::size_t row,
     const std::size_t col) noexcept {
   return block_getter{}.template operator()<ROWS, COLS>(m, row, col);
@@ -224,7 +224,7 @@ ALGEBRA_HOST_DEVICE inline constexpr auto block(
 template <std::size_t ROWS, std::size_t COLS, std::size_t mROW,
           std::size_t mCOL, typename value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline constexpr void set_block(
+ALGEBRA_HOST_DEVICE constexpr void set_block(
     matrix<array_t, value_t, mROW, mCOL> &m,
     const matrix<array_t, value_t, ROWS, COLS> &b, const std::size_t row,
     const std::size_t col) noexcept {
@@ -255,7 +255,7 @@ ALGEBRA_HOST_DEVICE inline constexpr void set_block(
 /// Operator setting a block with a vector
 template <class matrix_t, std::size_t ROW, typename value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline constexpr void set_block(
+ALGEBRA_HOST_DEVICE constexpr void set_block(
     matrix_t &m, const storage::vector<ROW, value_t, array_t> &b,
     std::size_t row, std::size_t col) noexcept {
   assert(row < ROW);
diff --git a/storage/eigen/include/algebra/storage/eigen.hpp b/storage/eigen/include/algebra/storage/eigen.hpp
index cb193743..f6b8905c 100644
--- a/storage/eigen/include/algebra/storage/eigen.hpp
+++ b/storage/eigen/include/algebra/storage/eigen.hpp
@@ -32,43 +32,38 @@ using vector_type = storage_type<T, N>;
 template <typename T, size_type ROWS, size_type COLS>
 using matrix_type = Eigen::Matrix<T, ROWS, COLS, (ROWS == 1), ROWS, COLS>;
 
-/// 3-element "vector" type, using @c algebra::eigen::array
+/// 3-element "vector" type, using @c eigen::vector_type
 template <typename T>
 using vector3 = storage_type<T, 3>;
-/// Point in 3D space, using @c algebra::eigen::array
+/// Point in 3D space, using @c eigen::vector_type
 template <typename T>
 using point3 = vector3<T>;
-/// 2-element "vector" type, using @c algebra::eigen::array
+/// 2-element "vector" type, using @c eigen::vector_type
 template <typename T>
 using vector2 = storage_type<T, 2>;
-/// Point in 2D space, using @c algebra::eigen::array
+/// Point in 2D space, using @c eigen::vector_type
 template <typename T>
 using point2 = vector2<T>;
 
+/// Element Getter
+using element_getter = eigen::storage::element_getter;
+/// Block Getter
+using block_getter = eigen::storage::block_getter;
+
 }  // namespace eigen
 
-namespace trait {
+ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(eigen)
 
-/// Type trait specializations
-/// @{
+// Extra specializations needed for some additional eigen types
+namespace trait {
 
 /// Index
 /// @{
-template <typename T, int N>
-struct index<algebra::eigen::array<T, N>> {
-  using type = algebra::eigen::size_type;
-};
-
 template <typename derived_t>
 struct index<Eigen::MatrixBase<derived_t>> {
   using type = algebra::eigen::size_type;
 };
 
-template <typename T, int ROWS, int COLS, int OPT>
-struct index<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>> {
-  using type = algebra::eigen::size_type;
-};
-
 template <typename T, int ROWS, int COLS, int OPT, int bROWS, int bCOLS>
 struct index<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>, bROWS,
                           bCOLS, false>> {
@@ -78,16 +73,6 @@ struct index<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>, bROWS,
 
 /// Dimensions
 /// @{
-template <typename T, int N>
-struct dimensions<algebra::eigen::array<T, N>> {
-
-  using size_type = index_t<algebra::eigen::array<T, N>>;
-
-  static constexpr size_type dim{1};
-  static constexpr size_type rows{N};
-  static constexpr size_type columns{1};
-};
-
 template <typename derived_t>
 struct dimensions<Eigen::MatrixBase<derived_t>> {
 
@@ -100,16 +85,6 @@ struct dimensions<Eigen::MatrixBase<derived_t>> {
       Eigen::MatrixBase<derived_t>::ColsAtCompileTime};
 };
 
-template <typename T, int ROWS, int COLS, int OPT>
-struct dimensions<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>> {
-
-  using size_type = index_t<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>>;
-
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-
 template <typename T, int ROWS, int COLS, int OPT, int bROWS, int bCOLS>
 struct dimensions<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>,
                                bROWS, bCOLS, false>> {
@@ -126,21 +101,11 @@ struct dimensions<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>,
 
 /// Value
 /// @{
-template <typename T, int N>
-struct value<algebra::eigen::array<T, N>> {
-  using type = T;
-};
-
 template <typename derived_t>
 struct value<Eigen::MatrixBase<derived_t>> {
   using type = typename Eigen::MatrixBase<derived_t>::value_type;
 };
 
-template <typename T, int ROWS, int COLS, int OPT>
-struct value<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>> {
-  using type = T;
-};
-
 template <typename T, int ROWS, int COLS, int OPT, int bROWS, int bCOLS>
 struct value<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>, bROWS,
                           bCOLS, false>> {
@@ -150,33 +115,15 @@ struct value<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>, bROWS,
 
 /// Vector
 /// @{
-template <typename T, int N>
-struct vector<algebra::eigen::array<T, N>> {
-
-  template <typename other_T, int other_N>
-  using other_type = algebra::eigen::array<other_T, other_N>;
-
-  using type = other_type<T, N>;
-};
-
 template <typename derived_t>
 struct vector<Eigen::MatrixBase<derived_t>> {
 
   template <typename other_T, int other_N>
-  using other_type = algebra::eigen::array<other_T, other_N>;
+  using other_type = eigen::vector_type<other_T, other_N>;
 
   using type = other_type<value_t<Eigen::MatrixBase<derived_t>>,
                           rows<Eigen::MatrixBase<derived_t>>>;
 };
-
-template <typename T, int ROWS, int COLS, int OPT>
-struct vector<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>> {
-
-  template <typename other_T, int other_N>
-  using other_type = algebra::eigen::array<other_T, other_N>;
-
-  using type = other_type<T, ROWS>;
-};
 /// @}
 
 /// Matrix
@@ -192,22 +139,6 @@ struct matrix<Eigen::MatrixBase<derived_t>> {
   using type = Eigen::MatrixBase<derived_t>;
 };
 
-template <typename T, int ROWS, int COLS, int OPT>
-struct matrix<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type = eigen::matrix_type<other_T, other_ROWS, other_COLS>;
-
-  using type = Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>;
-};
-
-template <typename T, int N>
-struct matrix<algebra::eigen::array<T, N>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type = eigen::matrix_type<other_T, other_ROWS, other_COLS>;
-
-  using type = other_type<T, N, 1>;
-};
-
 template <typename T, int ROWS, int COLS, int OPT, int bROWS, int bCOLS>
 struct matrix<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>, bROWS,
                            bCOLS, false>> {
@@ -221,21 +152,11 @@ struct matrix<Eigen::Block<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>, bROWS,
 
 /// Element getter
 /// @{
-template <typename T, int N>
-struct element_getter<algebra::eigen::array<T, N>> {
-  using type = eigen::storage::element_getter;
-};
-
 template <typename derived_t>
 struct element_getter<Eigen::MatrixBase<derived_t>> {
   using type = eigen::storage::element_getter;
 };
 
-template <typename T, int ROWS, int COLS, int OPT>
-struct element_getter<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>> {
-  using type = eigen::storage::element_getter;
-};
-
 template <typename T, int ROWS, int COLS, int OPT, int bROWS, int bCOLS>
 struct element_getter<Eigen::Block<
     Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>, bROWS, bCOLS, false>> {
@@ -249,12 +170,6 @@ template <typename derived_t>
 struct block_getter<Eigen::MatrixBase<derived_t>> {
   using type = eigen::storage::block_getter;
 };
-
-template <typename T, int ROWS, int COLS, int OPT>
-struct block_getter<Eigen::Matrix<T, ROWS, COLS, OPT, ROWS, COLS>> {
-  using type = eigen::storage::block_getter;
-};
-/// @}
 /// @}
 
 }  // namespace trait
diff --git a/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp b/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp
index 970a6f85..41073078 100644
--- a/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp
+++ b/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp
@@ -35,83 +35,73 @@ namespace algebra::eigen::storage {
 /// Functor used to access elements of Eigen matrices
 struct element_getter {
   /// Get non-const access to a matrix element
-  template <typename derived_type, typename size_type_1, typename size_type_2,
-            std::enable_if_t<
-                std::is_base_of<
-                    Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
-                    Eigen::MatrixBase<derived_type> >::value &&
-                    std::is_convertible<size_type_1, Eigen::Index>::value &&
-                    std::is_convertible<size_type_2, Eigen::Index>::value,
-                bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto &operator()(
-      Eigen::MatrixBase<derived_type> &m, size_type_1 row,
-      size_type_2 col) const {
+  template <typename derived_type, typename size_type_1, typename size_type_2>
+  requires(std::is_base_of_v<
+           Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
+           Eigen::MatrixBase<derived_type> >
+               &&std::is_convertible_v<size_type_1, Eigen::Index>
+                   &&std::is_convertible_v<size_type_2, Eigen::Index>)
+      ALGEBRA_HOST_DEVICE inline auto &
+      operator()(Eigen::MatrixBase<derived_type> &m, size_type_1 row,
+                 size_type_2 col) const {
 
     return m(static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
   }
   /// Get const access to a matrix element
-  template <typename derived_type, typename size_type_1, typename size_type_2,
-            std::enable_if_t<
-                std::is_convertible<size_type_1, Eigen::Index>::value &&
-                    std::is_convertible<size_type_2, Eigen::Index>::value,
-                bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto operator()(
-      const Eigen::MatrixBase<derived_type> &m, size_type_1 row,
-      size_type_2 col) const {
+  template <typename derived_type, typename size_type_1, typename size_type_2>
+  requires(std::is_convertible_v<size_type_1, Eigen::Index>
+               &&std::is_convertible_v<size_type_2, Eigen::Index>)
+      ALGEBRA_HOST_DEVICE inline auto
+      operator()(const Eigen::MatrixBase<derived_type> &m, size_type_1 row,
+                 size_type_2 col) const {
 
     return m(static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
   }
   /// Get non-const access to a matrix element
-  template <typename derived_type, typename size_type,
-            std::enable_if_t<
-                std::is_base_of<
-                    Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
-                    Eigen::MatrixBase<derived_type> >::value &&
-                    std::is_convertible<size_type, Eigen::Index>::value,
-                bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto &operator()(
-      Eigen::MatrixBase<derived_type> &m, size_type row) const {
+  template <typename derived_type, typename size_type>
+  requires(std::is_base_of_v<
+           Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
+           Eigen::MatrixBase<derived_type> >
+               &&std::is_convertible_v<size_type, Eigen::Index>)
+      ALGEBRA_HOST_DEVICE inline auto &
+      operator()(Eigen::MatrixBase<derived_type> &m, size_type row) const {
 
     return m(static_cast<Eigen::Index>(row));
   }
   /// Get const access to a matrix element
-  template <typename derived_type, typename size_type,
-            std::enable_if_t<
-                std::is_convertible<size_type, Eigen::Index>::value,
-                bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto operator()(
-      const Eigen::MatrixBase<derived_type> &m, size_type row) const {
+  template <typename derived_type, typename size_type>
+  requires std::is_convertible_v<size_type, Eigen::Index>
+      ALGEBRA_HOST_DEVICE inline auto operator()(
+          const Eigen::MatrixBase<derived_type> &m, size_type row) const {
 
     return m(static_cast<Eigen::Index>(row));
   }
 };  // struct element_getter
 
 /// Function extracting an element from a matrix (const)
-template <
-    typename derived_type>
+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, static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
+  return element_getter()(m, static_cast<Eigen::Index>(row),
+                          static_cast<Eigen::Index>(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, static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
+template <typename derived_type>
+requires std::is_base_of_v<
+    Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
+    Eigen::MatrixBase<derived_type> >
+    ALGEBRA_HOST_DEVICE inline auto &element(Eigen::MatrixBase<derived_type> &m,
+                                             std::size_t row, std::size_t col) {
+
+  return element_getter()(m, static_cast<Eigen::Index>(row),
+                          static_cast<Eigen::Index>(col));
 }
 
 /// Function extracting an element from a matrix (const)
-template <
-    typename derived_type>
+template <typename derived_type>
 ALGEBRA_HOST_DEVICE inline auto element(
     const Eigen::MatrixBase<derived_type> &m, std::size_t row) {
 
@@ -119,14 +109,12 @@ ALGEBRA_HOST_DEVICE inline auto element(
 }
 
 /// 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) {
+template <typename derived_type>
+requires std::is_base_of_v<
+    Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
+    Eigen::MatrixBase<derived_type> >
+    ALGEBRA_HOST_DEVICE inline auto &element(Eigen::MatrixBase<derived_type> &m,
+                                             std::size_t row) {
 
   return element_getter()(m, static_cast<Eigen::Index>(row));
 }
@@ -134,70 +122,65 @@ ALGEBRA_HOST_DEVICE inline auto &element(Eigen::MatrixBase<derived_type> &m,
 /// Functor used to extract a block from Eigen matrices
 struct block_getter {
   template <int kROWS, int kCOLS, typename matrix_type, typename size_type_1,
-            typename size_type_2,
-            std::enable_if_t<
-                std::is_convertible<size_type_1, Eigen::Index>::value &&
-                    std::is_convertible<size_type_2, Eigen::Index>::value,
-                bool> = true>
-  ALGEBRA_HOST_DEVICE decltype(auto) operator()(const matrix_type &m, size_type_1 row,
-                                      size_type_2 col) const {
+            typename size_type_2>
+  requires(std::is_convertible_v<size_type_1, Eigen::Index>
+               &&std::is_convertible_v<size_type_2, Eigen::Index>)
+      ALGEBRA_HOST_DEVICE decltype(auto)
+      operator()(const matrix_type &m, size_type_1 row, size_type_2 col) const {
 
     return m.template block<kROWS, kCOLS>(row, col);
   }
 
   template <int kROWS, int kCOLS, typename matrix_type, typename size_type_1,
-            typename size_type_2,
-            std::enable_if_t<
-                std::is_convertible<size_type_1, Eigen::Index>::value &&
-                    std::is_convertible<size_type_2, Eigen::Index>::value,
-                bool> = true>
-  ALGEBRA_HOST_DEVICE decltype(auto) operator()(matrix_type &m, size_type_1 row,
-                                      size_type_2 col) const {
+            typename size_type_2>
+  requires(std::is_convertible_v<size_type_1, Eigen::Index>
+               &&std::is_convertible_v<size_type_2, Eigen::Index>)
+      ALGEBRA_HOST_DEVICE decltype(auto)
+      operator()(matrix_type &m, size_type_1 row, size_type_2 col) const {
 
     return m.template block<kROWS, kCOLS>(row, col);
   }
 
   template <int SIZE, typename matrix_type, typename size_type_1,
-            typename size_type_2,
-            std::enable_if_t<
-                std::is_convertible<size_type_1, Eigen::Index>::value &&
-                    std::is_convertible<size_type_2, Eigen::Index>::value,
-                bool> = true>
-  ALGEBRA_HOST_DEVICE decltype(auto) operator()(matrix_type &m, size_type_1 row,
-                                      size_type_2 col) const {
+            typename size_type_2>
+  requires(std::is_convertible_v<size_type_1, Eigen::Index>
+               &&std::is_convertible_v<size_type_2, Eigen::Index>)
+      ALGEBRA_HOST_DEVICE decltype(auto)
+      operator()(matrix_type &m, size_type_1 row, size_type_2 col) const {
 
     return m.template block<SIZE, 1>(row, col);
   }
 };  // struct block_getter
 
 /// Operator getting a block of a const matrix
-template <
-    int ROWS, int COLS, class derived_type>
-ALGEBRA_HOST_DEVICE decltype(auto) block(const Eigen::MatrixBase<derived_type> &m,
-                               std::size_t row, std::size_t col) {
-  return block_getter{}.template operator()<ROWS, COLS>(m, static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
+template <int ROWS, int COLS, class derived_type>
+ALGEBRA_HOST_DEVICE decltype(auto) block(
+    const Eigen::MatrixBase<derived_type> &m, std::size_t row,
+    std::size_t col) {
+  return block_getter{}.template operator()<ROWS, COLS>(
+      m, static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
 }
 
 /// Operator getting a block of a const matrix
-template <
-    int ROWS, int COLS, class derived_type>
+template <int ROWS, int COLS, class derived_type>
 ALGEBRA_HOST_DEVICE decltype(auto) block(Eigen::MatrixBase<derived_type> &m,
-                               std::size_t row, std::size_t col) {
-  return block_getter{}.template operator()<ROWS, COLS>(m, static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
+                                         std::size_t row, std::size_t col) {
+  return block_getter{}.template operator()<ROWS, COLS>(
+      m, static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
 }
 
 /// Function extracting a slice from the matrix
 template <int SIZE, typename derived_type>
-ALGEBRA_HOST_DEVICE inline decltype(auto) vector(const Eigen::MatrixBase<derived_type>& m,
-                                       std::size_t row, std::size_t col) {
+ALGEBRA_HOST_DEVICE inline decltype(auto) vector(
+    const Eigen::MatrixBase<derived_type> &m, std::size_t row,
+    std::size_t col) {
 
-  return block_getter{}.template operator()<SIZE>(m, static_cast<Eigen::Index>(row),
-                                   static_cast<Eigen::Index>(col));
+  return block_getter{}.template operator()<SIZE>(
+      m, static_cast<Eigen::Index>(row), static_cast<Eigen::Index>(col));
 }
 
 /// Operator setting a block
-template <
-    typename derived_type1, typename derived_type2>
+template <typename derived_type1, typename derived_type2>
 ALGEBRA_HOST_DEVICE void set_block(Eigen::MatrixBase<derived_type1> &m,
                                    const Eigen::MatrixBase<derived_type2> &b,
                                    std::size_t row, std::size_t col) {
diff --git a/storage/fastor/include/algebra/storage/fastor.hpp b/storage/fastor/include/algebra/storage/fastor.hpp
index c58271f9..8b2ccb98 100644
--- a/storage/fastor/include/algebra/storage/fastor.hpp
+++ b/storage/fastor/include/algebra/storage/fastor.hpp
@@ -24,6 +24,9 @@ using size_type = std::size_t;
 /// Array type used in the Fastor storage model
 template <typename T, size_type N>
 using storage_type = Fastor::Tensor<T, N>;
+/// Vector type used in the Fastor storage model
+template <typename T, size_type N>
+using vector_type = storage_type<T, N>;
 /// Matrix type used in the Fastor storage model
 template <typename T, size_type ROWS, size_type COLS>
 using matrix_type = algebra::fastor::Matrix<T, ROWS, COLS>;
@@ -41,172 +44,13 @@ using vector2 = storage_type<T, 2>;
 template <typename T>
 using point2 = vector2<T>;
 
-}  // namespace fastor
-
-namespace trait {
-
-/// Type trait specializations
-/// @{
-
-/// Index
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct index<Fastor::Tensor<T, ROWS, COLS>> {
-  using type = std::size_t;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct index<algebra::fastor::Matrix<T, ROWS, COLS>> {
-  using type = std::size_t;
-};
-
-template <typename T, std::size_t N>
-struct index<Fastor::Tensor<T, N>> {
-  using type = std::size_t;
-};
-/// @}
-
-/// Dimension
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct dimensions<Fastor::Tensor<T, ROWS, COLS>> {
-
-  using size_type = std::size_t;
-
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct dimensions<algebra::fastor::Matrix<T, ROWS, COLS>> {
-
-  using size_type = std::size_t;
-
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-
-template <typename T, std::size_t N>
-struct dimensions<Fastor::Tensor<T, N>> {
-
-  using size_type = std::size_t;
-
-  static constexpr size_type dim{1};
-  static constexpr size_type rows{N};
-  static constexpr size_type columns{1};
-};
-/// @}
-
-/// Value
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct value<Fastor::Tensor<T, ROWS, COLS>> {
-  using type = T;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct value<algebra::fastor::Matrix<T, ROWS, COLS>> {
-  using type = T;
-};
-
-template <typename T, std::size_t N>
-struct value<Fastor::Tensor<T, N>> {
-  using type = T;
-};
-/// @}
-
-/// Vector
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct vector<Fastor::Tensor<T, ROWS, COLS>> {
+/// Element Getter
+using element_getter = fastor::storage::element_getter;
+/// Block Getter
+using block_getter = fastor::storage::block_getter;
 
-  template <typename other_T, std::size_t other_N>
-  using other_type = Fastor::Tensor<other_T, other_N>;
-
-  using type = other_type<T, ROWS>;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct vector<algebra::fastor::Matrix<T, ROWS, COLS>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type = Fastor::Tensor<other_T, other_N>;
-
-  using type = Fastor::Tensor<T, ROWS>;
-};
-
-template <typename T, std::size_t N>
-struct vector<Fastor::Tensor<T, N>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type = Fastor::Tensor<other_T, other_N>;
-
-  using type = Fastor::Tensor<T, N>;
-};
-/// @}
-
-/// Matrix
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct matrix<Fastor::Tensor<T, ROWS, COLS>> {
-  template <typename other_T, std::size_t other_ROWS, std::size_t other_COLS>
-  using other_type = Fastor::Tensor<other_T, other_ROWS, other_COLS>;
-
-  using type = Fastor::Tensor<T, ROWS, COLS>;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct matrix<algebra::fastor::Matrix<T, ROWS, COLS>> {
-  template <typename other_T, std::size_t other_ROWS, std::size_t other_COLS>
-  using other_type = algebra::fastor::Matrix<other_T, other_ROWS, other_COLS>;
-
-  using type = algebra::fastor::Matrix<T, ROWS, COLS>;
-};
-
-template <typename T, int N>
-struct matrix<Fastor::Tensor<T, N>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type = algebra::fastor::Matrix<other_T, other_ROWS, other_COLS>;
-
-  using type = other_type<T, N, 1>;
-};
-/// @}
-
-/// Element/Block Getter
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct element_getter<Fastor::Tensor<T, ROWS, COLS>> {
-  using type = fastor::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct element_getter<algebra::fastor::Matrix<T, ROWS, COLS>> {
-  using type = fastor::storage::element_getter;
-};
-
-template <typename T, std::size_t N>
-struct element_getter<Fastor::Tensor<T, N>> {
-  using type = fastor::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct block_getter<Fastor::Tensor<T, ROWS, COLS>> {
-  using type = fastor::storage::block_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct block_getter<algebra::fastor::Matrix<T, ROWS, COLS>> {
-  using type = fastor::storage::block_getter;
-};
-
-template <typename T, std::size_t N>
-struct block_getter<Fastor::Tensor<T, N>> {
-  using type = fastor::storage::block_getter;
-};
-/// @}
+}  // namespace fastor
 
-}  // namespace trait
+ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(fastor)
 
 }  // namespace algebra
diff --git a/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp b/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp
index 5184e0ba..1d32232e 100644
--- a/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp
+++ b/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp
@@ -23,19 +23,18 @@
 namespace algebra::smatrix::storage {
 
 /// Functor used to access elements of SMatrix matrices
-template <typename scalar_t>
 struct element_getter {
 
-  template <unsigned int ROWS, unsigned int COLS, typename scalar_t>
-  ALGEBRA_HOST_DEVICE inline scalar_t &operator()(ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m,
-                                                  unsigned int row,
-                                                  unsigned int col) const {
+  template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
+  ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
+      ROOT::Math::SMatrix<scalar_t, 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, typename scalar_t>
+  template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, unsigned int row,
       unsigned int col) const {
@@ -44,31 +43,31 @@ struct element_getter {
     return m(row, col);
   }
 
-  template <unsigned int N, typename scalar_t>
-  ALGEBRA_HOST_DEVICE inline scalar_t &operator()(ROOT::Math::SMatrix<scalar_t, N, 1> &m,
-                                                  unsigned int row) const {
-    assert(row < ROWS);
-    return m(row);
+  template <typename scalar_t, unsigned int N>
+  ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
+      ROOT::Math::SMatrix<scalar_t, N, 1> &m, unsigned int row) const {
+    assert(row < N);
+    return m(row, 0);
   }
 
-  template <unsigned int N, typename scalar_t>
+  template <typename scalar_t, unsigned int N>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const ROOT::Math::SMatrix<scalar_t, N, 1> &m, unsigned int row) const {
-    assert(row < ROWS);
-    return m(row);
+    assert(row < N);
+    return m(row, 0);
   }
 
-  template <unsigned int N, typename scalar_t>
-  ALGEBRA_HOST_DEVICE inline scalar_t &operator()(ROOT::Math::SVector<scalar_t, N> &m,
-                                                  unsigned int row) const {
-    assert(row < ROWS);
+  template <typename scalar_t, unsigned int N>
+  ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
+      ROOT::Math::SVector<scalar_t, N> &m, unsigned int row) const {
+    assert(row < N);
     return m(row);
   }
 
-  template <unsigned int N, typename scalar_t>
+  template <typename scalar_t, unsigned int N>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const ROOT::Math::SVector<scalar_t, N> &m, unsigned int row) const {
-    assert(row < ROWS);
+    assert(row < N);
     return m(row);
   }
 };  // element_getter
@@ -78,7 +77,8 @@ 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, static_cast<unsigned int>(row), static_cast<unsigned int>(col));
+  return element_getter()(m, static_cast<unsigned int>(row),
+                          static_cast<unsigned int>(col));
 }
 
 /// Function extracting an element from a matrix (non-const)
@@ -86,62 +86,68 @@ 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, static_cast<unsigned int>(row), static_cast<unsigned int>(col));
+  return element_getter()(m, static_cast<unsigned int>(row),
+                          static_cast<unsigned int>(col));
 }
 
 /// Function extracting an element from a matrix (const)
 template <typename scalar_t, unsigned int N>
 ALGEBRA_HOST_DEVICE inline scalar_t element(
     const ROOT::Math::SMatrix<scalar_t, N, 1> &m, std::size_t row) {
-  return element_getter<scalar_t>()(m, static_cast<unsigned int>(row));
+  return element_getter()(m, static_cast<unsigned int>(row));
 }
 
 /// Function extracting an element from a matrix (non-const)
 template <typename scalar_t, unsigned int N>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(
     ROOT::Math::SMatrix<scalar_t, N, 1> &m, std::size_t row) {
-  return element_getter<scalar_t>()(m, static_cast<unsigned int>(row));
+  return element_getter()(m, static_cast<unsigned int>(row));
 }
 
 /// Function extracting an element from a matrix (const)
 template <typename scalar_t, unsigned int N>
 ALGEBRA_HOST_DEVICE inline scalar_t element(
     const ROOT::Math::SVector<scalar_t, N> &m, std::size_t row) {
-  return element_getter<scalar_t>()(m, static_cast<unsigned int>(row));
+  return element_getter()(m, static_cast<unsigned int>(row));
 }
 
 /// Function extracting an element from a matrix (non-const)
 template <typename scalar_t, unsigned int N>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(
-    ROOT::Math::SVector<scalar_t, N, 1> &m, std::size_t row) {
-  return element_getter<scalar_t>()(m, static_cast<unsigned int>(row));
+    ROOT::Math::SVector<scalar_t, N> &m, std::size_t row) {
+  return element_getter()(m, static_cast<unsigned int>(row));
 }
 
 /// Functor used to extract a block from SMatrix matrices
 struct block_getter {
 
-  template <unsigned int ROWS, unsigned int COLS, unsigned int oROWS, unsigned int oCOLS, typename scalar_t>
+  template <unsigned int ROWS, unsigned int COLS, unsigned int oROWS,
+            unsigned int oCOLS, typename scalar_t>
   ALGEBRA_HOST_DEVICE ROOT::Math::SMatrix<scalar_t, ROWS, COLS> operator()(
-      const ROOT::Math::SMatrix<scalar_t, oROWS, oCOLS> &m, unsigned int row, unsigned int col) const {
+      const ROOT::Math::SMatrix<scalar_t, oROWS, oCOLS> &m, unsigned int row,
+      unsigned int col) const {
 
     return m.template Sub<ROOT::Math::SMatrix<scalar_t, ROWS, COLS>>(row, col);
   }
 
-  template <unsigned int SIZE, unsigned int ROWS, unsigned int COLS, typename scalar_t>
+  template <unsigned int SIZE, unsigned int ROWS, unsigned int COLS,
+            typename scalar_t>
   ALGEBRA_HOST_DEVICE ROOT::Math::SVector<scalar_t, SIZE> operator()(
-      const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, unsigned int row, unsigned int col) const {
+      const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, unsigned int row,
+      unsigned int col) const {
 
-  return m.template SubCol<ROOT::Math::SVector<scalar_t, SIZE>>(col, row);
+    return m.template SubCol<ROOT::Math::SVector<scalar_t, SIZE>>(col, row);
   }
 };  // struct block_getter
 
 /// Operator getting a block of a const matrix
-template <unsigned int ROWS, unsigned int COLS, unsigned int oROWS, unsigned int oCOLS,
-          typename scalar_t>
-ALGEBRA_HOST_DEVICE matrix_type<ROWS, COLS> block(const ROOT::Math::SMatrix<scalar_t, oROWS, oCOLS>&m,
-                                                  std::size_t row,
-                                                  std::size_t col) {
-  return block_getter{}.template operator()<ROWS, COLS>(m, static_cast<unsigned int>(row), static_cast<unsigned int>(col));
+template <unsigned int ROWS, unsigned int COLS, unsigned int oROWS,
+          unsigned int oCOLS, typename scalar_t>
+ALGEBRA_HOST_DEVICE ROOT::Math::SMatrix<scalar_t, ROWS, COLS> block(
+    const ROOT::Math::SMatrix<scalar_t, oROWS, oCOLS> &m, std::size_t row,
+    std::size_t col) {
+  return block_getter{}.template operator()<ROWS, COLS>(
+      m, static_cast<unsigned int>(row), static_cast<unsigned int>(col));
 }
 
 /// Function extracting a slice from the matrix used by
@@ -149,32 +155,35 @@ ALGEBRA_HOST_DEVICE matrix_type<ROWS, COLS> block(const ROOT::Math::SMatrix<scal
 template <unsigned int SIZE, unsigned int ROWS, unsigned int COLS,
           typename scalar_t>
 ALGEBRA_HOST_DEVICE inline auto vector(
-    const ROOT::Math::SMatrix<scalar_t, ROWS, COLS>& m, std::size_t row,
+    const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, std::size_t row,
     std::size_t col) {
 
-  return block_getter{}.template operator()<SIZE>(m, static_cast<unsigned int>(row), static_cast<unsigned int>(col));
+  return block_getter{}.template operator()<SIZE>(
+      m, static_cast<unsigned int>(row), static_cast<unsigned int>(col));
 }
 
 /// Operator setting a block with a matrix
-template <unsigned int ROWS, unsigned int COLS, class input_matrix_type>
-ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m,
-                                   const matrix_type<ROWS, COLS> &b,
-                                   std::size_t row, std::size_t col) {
+template <unsigned int ROWS, unsigned int COLS, typename scalar_t,
+          class input_matrix_type>
+ALGEBRA_HOST_DEVICE void set_block(
+    input_matrix_type &m, const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &b,
+    std::size_t row, std::size_t col) {
   for (unsigned int i = 0; i < ROWS; ++i) {
     for (unsigned int j = 0; j < COLS; ++j) {
-      m(i + static_cast<unsigned int>(row), j + static_cast<unsigned int>(col)) = b(i, j);
+      m(i + static_cast<unsigned int>(row),
+        j + static_cast<unsigned int>(col)) = b(i, j);
     }
   }
 }
 
 /// Operator setting a block with a vector
-template <unsigned int ROWS, class input_matrix_type>
+template <unsigned int ROWS, typename scalar_t, class input_matrix_type>
 ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m,
-                                   const vector_type<ROWS> &b, unsigned int row,
-                                   unsigned int col) {
+                                   const ROOT::Math::SVector<scalar_t, ROWS> &b,
+                                   unsigned int row, unsigned int col) {
   for (unsigned int i = 0; i < ROWS; ++i) {
     m(i + row, col) = b[i];
   }
 }
 
-}  // namespace algebra::smatrix
+}  // namespace algebra::smatrix::storage
diff --git a/storage/smatrix/include/algebra/storage/smatrix.hpp b/storage/smatrix/include/algebra/storage/smatrix.hpp
index 8527c10f..97169288 100644
--- a/storage/smatrix/include/algebra/storage/smatrix.hpp
+++ b/storage/smatrix/include/algebra/storage/smatrix.hpp
@@ -27,6 +27,9 @@ using size_type = unsigned int;
 /// Array type used in the SMatrix storage model
 template <typename T, size_type N>
 using storage_type = ROOT::Math::SVector<T, N>;
+/// Vector type used in the SMatrix storage model
+template <typename T, size_type N>
+using vector_type = storage_type<T, N>;
 /// Matrix type used in the SMatrix storage model
 template <typename T, size_type ROWS, size_type COLS>
 using matrix_type = ROOT::Math::SMatrix<T, ROWS, COLS>;
@@ -44,101 +47,13 @@ using vector2 = storage_type<T, 2>;
 template <typename T>
 using point2 = vector2<T>;
 
-}  // namespace smatrix
-
-namespace trait {
-
-/// Type trait specializations
-/// @{
-
-/// Index
-/// @{
-template <typename T, unsigned int ROWS, unsigned int COLS, typename expr>
-struct index<ROOT::Math::SMatrix<T, ROWS, COLS, expr>> {
-  using type = unsigned int;
-};
-/// @}
-
-/// Dimension
-/// @{
-template <typename T, unsigned int ROWS, unsigned int COLS, typename expr>
-struct dimensions<ROOT::Math::SMatrix<T, ROWS, COLS, expr>> {
-
-  using size_type = unsigned int;
-
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-/// @}
-
-/// Value
-/// @{
-template <typename T, unsigned int ROWS, unsigned int COLS, typename expr>
-struct value<ROOT::Math::SMatrix<T, ROWS, COLS, expr>> {
-  using type = T;
-};
-/// @}
-
-/// Vector
-/// @{
-template <typename T, unsigned int ROWS, unsigned int COLS, typename expr>
-struct vector<ROOT::Math::SMatrix<T, ROWS, COLS, expr>> {
-
-  template <typename other_T, unsigned int other_N>
-  using other_type = ROOT::Math::SVector<other_T, other_N>;
+/// Element Getter
+using element_getter = smatrix::storage::element_getter;
+/// Block Getter
+using block_getter = smatrix::storage::block_getter;
 
-  using type = other_type<T, ROWS>;
-};
-
-template <typename T, unsigned int N>
-struct vector<ROOT::Math::SVector<T, N>> {
-  using type = ROOT::Math::SVector<T, N>;
-
-  template <typename other_T, unsigned int other_N>
-  using other_type = ROOT::Math::SVector<other_T, other_N>;
-
-  using type = other_type<T, N>;
-};
-/// @}
-
-/// Matrix
-/// @{
-template <typename T, unsigned int ROWS, unsigned int COLS, typename expr>
-struct matrix<ROOT::Math::SMatrix<T, ROWS, COLS, expr>> {
-  template <typename other_T, unsigned int other_ROWS, unsigned int other_COLS>
-  using other_type = ROOT::Math::SMatrix<other_T, other_ROWS, other_COLS, expr>;
-
-  using type = ROOT::Math::SMatrix<T, ROWS, COLS, expr>;
-};
-
-template <typename T, int N>
-struct matrix<ROOT::Math::SVector<T, N>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type = ROOT::Math::SMatrix<other_T, other_ROWS, other_COLS>;
-
-  using type = other_type<T, N, 1>;
-};
-/// @}
-
-/// Elemet/Block Getter
-/// @{
-template <typename T, unsigned int ROWS, unsigned int COLS, typename expr>
-struct element_getter<ROOT::Math::SMatrix<T, ROWS, COLS, expr>> {
-  using type = smatrix::storage::element_getter;
-};
-
-template <typename T, unsigned int N>
-struct element_getter<ROOT::Math::SVector<T, N>> {
-  using type = smatrix::storage::element_getter;
-};
-
-template <typename T, unsigned int ROWS, unsigned int COLS, typename expr>
-struct block_getter<ROOT::Math::SMatrix<T, ROWS, COLS, expr>> {
-  using type = smatrix::storage::block_getter;
-};
-/// @}
+}  // namespace smatrix
 
-}  // namespace trait
+ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(smatrix)
 
 }  // namespace algebra
diff --git a/storage/vc_aos/include/algebra/storage/impl/vc_aos_getter.hpp b/storage/vc_aos/include/algebra/storage/impl/vc_aos_getter.hpp
index 459655a3..aa54f67a 100644
--- a/storage/vc_aos/include/algebra/storage/impl/vc_aos_getter.hpp
+++ b/storage/vc_aos/include/algebra/storage/impl/vc_aos_getter.hpp
@@ -20,7 +20,7 @@ using algebra::storage::set_block;
 /// Get a vector of a const matrix
 template <std::size_t SIZE, std::size_t ROW, std::size_t COL, typename value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline constexpr auto vector(
+ALGEBRA_HOST_DEVICE constexpr auto vector(
     const algebra::storage::matrix<array_t, value_t, ROW, COL> &m,
     const std::size_t row, const std::size_t col) noexcept {
   return algebra::storage::block_getter{}.template operator()<SIZE>(m, row,
diff --git a/storage/vc_aos/include/algebra/storage/vc_aos.hpp b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
index 4a5d334f..f3960d00 100644
--- a/storage/vc_aos/include/algebra/storage/vc_aos.hpp
+++ b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
@@ -64,128 +64,13 @@ using vector6 = vector_type<T, 6>;
 template <typename T>
 using vector8 = vector_type<T, 8>;
 
-}  // namespace vc_aos
-
-namespace trait {
-
-/// Type trait specializations
-/// @{
-
-/// Index
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct index<algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>> {
-  using type = algebra::vc_aos::size_type;
-};
-
-template <typename T, std::size_t N>
-struct index<algebra::storage::vector<N, T, vc_aos::storage_type>> {
-  using type = algebra::vc_aos::size_type;
-};
-/// @}
-
-/// Dimension
-/// @{
-template <typename T, std::size_t N>
-struct dimensions<algebra::storage::vector<N, T, vc_aos::storage_type>> {
-
-  using size_type =
-      index_t<algebra::storage::vector<N, T, vc_aos::storage_type>>;
-
-  static constexpr size_type dim{1};
-  static constexpr size_type rows{N};
-  static constexpr size_type columns{1};
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct dimensions<
-    algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>> {
-
-  using size_type =
-      index_t<algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>>;
-
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-/// @}
-
-/// Value
-/// @{
-template <typename T, std::size_t N>
-struct value<algebra::storage::vector<N, T, vc_aos::storage_type>> {
-  using type = T;
-};
+/// Element Getter
+using element_getter = algebra::storage::element_getter;
+/// Block Getter
+using block_getter = algebra::storage::block_getter;
 
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct value<algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>> {
-  using type = T;
-};
-/// @}
-
-/// Matrix
-/// @{
-template <typename T, std::size_t N>
-struct vector<algebra::storage::vector<N, T, vc_aos::storage_type>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type =
-      algebra::storage::vector<other_N, other_T, vc_aos::storage_type>;
-
-  using type = other_type<T, N>;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct vector<algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type =
-      algebra::storage::vector<other_N, other_T, vc_aos::storage_type>;
-
-  using type = other_type<T, ROWS>;
-};
-/// @}
-
-/// Matrix
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct matrix<algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>> {
-  template <typename other_T, std::size_t other_ROWS, std::size_t other_COLS>
-  using other_type = algebra::storage::matrix<vc_aos::storage_type, other_T,
-                                              other_ROWS, other_COLS>;
-
-  using type = algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>;
-};
-
-template <typename T, int N>
-struct matrix<algebra::storage::vector<N, T, vc_aos::storage_type>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type = algebra::storage::matrix<vc_aos::storage_type, other_T,
-                                              other_ROWS, other_COLS>;
-
-  using type = other_type<T, N, 1>;
-};
-/// @}
-
-/// Elemet/Block Getter
-/// @{
-template <typename T, std::size_t N>
-struct element_getter<algebra::storage::vector<N, T, vc_aos::storage_type>> {
-  using type = algebra::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct element_getter<
-    algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>> {
-  using type = algebra::storage::element_getter;
-};
+}  // namespace vc_aos
 
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct block_getter<
-    algebra::storage::matrix<vc_aos::storage_type, T, ROWS, COLS>> {
-  using type = algebra::storage::block_getter;
-};
-/// @}
-}  // namespace trait
+ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(vc_aos)
 
 }  // namespace algebra
diff --git a/storage/vc_soa/include/algebra/storage/impl/vc_soa_getter.hpp b/storage/vc_soa/include/algebra/storage/impl/vc_soa_getter.hpp
index 41a3e6fb..80418087 100644
--- a/storage/vc_soa/include/algebra/storage/impl/vc_soa_getter.hpp
+++ b/storage/vc_soa/include/algebra/storage/impl/vc_soa_getter.hpp
@@ -20,7 +20,7 @@ using algebra::storage::set_block;
 /// Get a vector of a const matrix
 template <std::size_t SIZE, std::size_t ROW, std::size_t COL, typename value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline constexpr auto vector(
+ALGEBRA_HOST_DEVICE constexpr auto vector(
     const algebra::storage::matrix<array_t, value_t, ROW, COL> &m,
     const std::size_t row, const std::size_t col) noexcept {
   return algebra::storage::block_getter{}.template operator()<SIZE>(m, row,
diff --git a/storage/vc_soa/include/algebra/storage/vc_soa.hpp b/storage/vc_soa/include/algebra/storage/vc_soa.hpp
index 905062be..c1821d45 100644
--- a/storage/vc_soa/include/algebra/storage/vc_soa.hpp
+++ b/storage/vc_soa/include/algebra/storage/vc_soa.hpp
@@ -40,11 +40,11 @@ template <typename T>
 using value_type = Vc::Vector<T>;
 /// Vector type used in the Vc SoA storage model
 template <typename T, std::size_t N>
-using vector_type = algebra::storage::vector<N, value_type<T>, storage_type>;
+using vector_type = algebra::storage::vector<N, Vc::Vector<T>, std::array>;
 /// Matrix type used in the Vc SoA storage model
 template <typename T, size_type ROWS, size_type COLS>
 using matrix_type =
-    algebra::storage::matrix<storage_type, value_type<T>, ROWS, COLS>;
+    algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>;
 
 /// 2-element "vector" type, using @c Vc::Vector in every element
 template <typename T>
@@ -65,113 +65,25 @@ using vector6 = vector_type<T, 6>;
 template <typename T>
 using vector8 = vector_type<T, 8>;
 
-}  // namespace vc_soa
-
-namespace trait {
+/// Element Getter
+using element_getter = algebra::storage::element_getter;
+/// Block Getter
+using block_getter = algebra::storage::block_getter;
 
-/// Type trait specializations
-/// @{
-
-/// Index
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct index<algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
-  using type = algebra::vc_soa::size_type;
-};
-/// @}
-
-/// Dimension
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct dimensions<
-    algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
+}  // namespace vc_soa
 
-  using size_type =
-      index_t<algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>>;
+ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(vc_soa)
 
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-/// @}
+namespace trait {
 
-/// Value/Scalar
+/// Make sure the simd scalar type is recognized correctly
 /// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct value<algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
-  using type = T;
-};
-
 template <typename T, std::size_t ROWS, std::size_t COLS>
 struct scalar<algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
   using type = Vc::Vector<T>;
 };
 /// @}
 
-/// Vector
-/// @{
-template <typename T, std::size_t N>
-struct vector<algebra::storage::vector<N, Vc::Vector<T>, std::array>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type =
-      algebra::storage::vector<other_N, Vc::Vector<other_T>, std::array>;
-
-  using type = other_type<T, N>;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct vector<algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type =
-      algebra::storage::vector<other_N, Vc::Vector<other_T>, std::array>;
-
-  using type = other_type<T, ROWS>;
-};
-/// @}
-
-/// Matrix
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct matrix<algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
-  template <typename other_T, std::size_t other_ROWS, std::size_t other_COLS>
-  using other_type = algebra::storage::matrix<std::array, Vc::Vector<other_T>,
-                                              other_ROWS, other_COLS>;
-
-  using type = algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>;
-};
-
-template <typename T, int N>
-struct matrix<algebra::storage::vector<N, Vc::Vector<T>, std::array>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type = algebra::storage::matrix<std::array, Vc::Vector<other_T>,
-                                              other_ROWS, other_COLS>;
-
-  using type = other_type<T, N, 1>;
-};
-/// @}
-
-/// Elemet/Block Getter
-/// @{
-template <typename T, std::size_t N>
-struct element_getter<algebra::storage::vector<N, Vc::Vector<T>, std::array>> {
-  using type = algebra::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct element_getter<
-    algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
-  using type = algebra::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct block_getter<
-    algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
-  using type = algebra::storage::block_getter;
-};
-/// @}
-
 }  // namespace trait
 
 }  // namespace algebra
diff --git a/storage/vecmem/include/algebra/storage/vecmem.hpp b/storage/vecmem/include/algebra/storage/vecmem.hpp
index 77ef938e..ec76106a 100644
--- a/storage/vecmem/include/algebra/storage/vecmem.hpp
+++ b/storage/vecmem/include/algebra/storage/vecmem.hpp
@@ -21,9 +21,14 @@ namespace algebra {
 
 namespace vecmem {
 
+/// size type for VecMem storage model
+using size_type = std::size_t;
 /// Array type used in the VecMem storage model
 template <typename T, std::size_t N>
 using storage_type = ::vecmem::static_array<T, N>;
+/// Vector type used in the VecMem storage model
+template <typename T, std::size_t N>
+using vector_type = storage_type<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>;
@@ -41,128 +46,13 @@ using vector2 = storage_type<T, 2>;
 template <typename T>
 using point2 = vector2<T>;
 
-}  // namespace vecmem
-
-namespace trait {
-
-/// Type trait specializations
-/// @{
-
-/// Index
-/// @{
-template <typename T, std::size_t N>
-struct index<::vecmem::static_array<T, N>> {
-  using type = std::size_t;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct index<::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>> {
-  using type = std::size_t;
-};
-/// @}
-
-/// Dimension
-/// @{
-template <typename T, std::size_t N>
-struct dimensions<::vecmem::static_array<T, N>> {
-
-  using size_type = index_t<::vecmem::static_array<T, N>>;
-
-  static constexpr size_type dim{1};
-  static constexpr size_type rows{N};
-  static constexpr size_type columns{1};
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct dimensions<
-    ::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>> {
-
-  using size_type =
-      index_t<::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>>;
-
-  static constexpr size_type dim{2};
-  static constexpr size_type rows{ROWS};
-  static constexpr size_type columns{COLS};
-};
-/// @}
-
-/// Value
-/// @{
-template <typename T, std::size_t N>
-struct value<::vecmem::static_array<T, N>> {
-  using type = T;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct value<::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>> {
-  using type = T;
-};
-/// @}
+/// Element Getter
+using element_getter = cmath::storage::element_getter;
+/// Block Getter
+using block_getter = cmath::storage::block_getter;
 
-/// Vector
-/// @{
-template <typename T, std::size_t N>
-struct vector<::vecmem::static_array<T, N>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type = ::vecmem::static_array<other_T, other_N>;
-
-  using type = other_type<T, N>;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct vector<::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>> {
-
-  template <typename other_T, std::size_t other_N>
-  using other_type = ::vecmem::static_array<other_T, other_N>;
-
-  using type = other_type<T, ROWS>;
-};
-/// @}
-
-/// Matrix
-/// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct matrix<::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>> {
-  template <typename other_T, std::size_t other_ROWS, std::size_t other_COLS>
-  using other_type =
-      ::vecmem::static_array<::vecmem::static_array<other_T, other_ROWS>,
-                             other_COLS>;
-
-  using type = ::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>;
-};
-
-template <typename T, int N>
-struct matrix<::vecmem::static_array<T, N>> {
-  template <typename other_T, int other_ROWS, int other_COLS>
-  using other_type =
-      ::vecmem::static_array<::vecmem::static_array<other_T, other_ROWS>,
-                             other_COLS>;
-
-  using type = other_type<T, N, 1>;
-};
-/// @}
-
-/// Elemet/Block Getter
-/// @{
-template <typename T, std::size_t N>
-struct element_getter<::vecmem::static_array<T, N>> {
-  using type = cmath::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct element_getter<
-    ::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>> {
-  using type = cmath::storage::element_getter;
-};
-
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct block_getter<
-    ::vecmem::static_array<::vecmem::static_array<T, ROWS>, COLS>> {
-  using type = cmath::storage::block_getter;
-};
-/// @}
+}  // namespace vecmem
 
-}  // namespace trait
+ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(vecmem)
 
 }  // namespace algebra
diff --git a/tests/accelerator/sycl/CMakeLists.txt b/tests/accelerator/sycl/CMakeLists.txt
index b8faa892..e8df7445 100644
--- a/tests/accelerator/sycl/CMakeLists.txt
+++ b/tests/accelerator/sycl/CMakeLists.txt
@@ -24,8 +24,8 @@ target_link_libraries( algebra_tests_sycl_common
 add_library( algebra::tests_sycl_common ALIAS algebra_tests_sycl_common )
 
 # Set up all of the (available) SYCL tests.
-algebra_add_test( array_sycl "array_generic.sycl"
-   LINK_LIBRARIES algebra::array_generic algebra::tests_sycl_common
+algebra_add_test( array_sycl "array_cmath.sycl"
+   LINK_LIBRARIES algebra::array_cmath algebra::tests_sycl_common
                   GTest::gtest_main )
 
 if( ALGEBRA_PLUGINS_INCLUDE_EIGEN )
diff --git a/tests/accelerator/sycl/array_cmath.sycl b/tests/accelerator/sycl/array_cmath.sycl
index bca3bcbe..c353adf2 100644
--- a/tests/accelerator/sycl/array_cmath.sycl
+++ b/tests/accelerator/sycl/array_cmath.sycl
@@ -6,7 +6,7 @@
  */
 
 // Project include(s).
-#include "algebra/array_generic.hpp"
+#include "algebra/array_cmath.hpp"
 
 // Local include(s).
 #include "test_sycl_basics.hpp"
@@ -23,9 +23,9 @@ struct test_specialisation_name {
   static std::string GetName(int i) {
     switch (i) {
       case 0:
-        return "sycl_array_generic<float>";
+        return "sycl_array_cmath<float>";
       case 1:
-        return "sycl_array_generic<double>";
+        return "sycl_array_cmath<double>";
       default:
         return "unknown";
     }
@@ -44,6 +44,6 @@ typedef testing::Types<
                algebra::array::vector3<double>,
                algebra::array::transform3<double>, std::size_t,
                algebra::array::matrix_type>>
-    array_generic_types;
+    array_cmath_types;
 INSTANTIATE_TYPED_TEST_SUITE_P(algebra_plugins, test_sycl_basics,
-                               array_generic_types, test_specialisation_name);
+                               array_cmath_types, test_specialisation_name);
diff --git a/tests/accelerator/sycl/eigen_cmath.sycl b/tests/accelerator/sycl/eigen_generic.sycl
similarity index 100%
rename from tests/accelerator/sycl/eigen_cmath.sycl
rename to tests/accelerator/sycl/eigen_generic.sycl