diff --git a/benchmarks/array/include/benchmark/array/data_generator.hpp b/benchmarks/array/include/benchmark/array/data_generator.hpp
index b4bf5939..bfbfd847 100644
--- a/benchmarks/array/include/benchmark/array/data_generator.hpp
+++ b/benchmarks/array/include/benchmark/array/data_generator.hpp
@@ -9,6 +9,7 @@
 
 // Project include(s)
 #include "algebra/array_cmath.hpp"
+#include "algebra/concepts.hpp"
 
 // System include(s)
 #include <algorithm>
@@ -18,7 +19,7 @@
 namespace algebra {
 
 /// Fill an @c std::array based vector with random values
-template <typename vector_t>
+template <concepts::vector vector_t>
 inline void fill_random_vec(std::vector<vector_t> &collection) {
 
   // Generate a vector of the right type with random values
@@ -33,7 +34,7 @@ inline void fill_random_vec(std::vector<vector_t> &collection) {
 }
 
 /// Fill a @c std::array based transform3 with random values
-template <typename transform3_t>
+template <concepts::transform3D transform3_t>
 inline void fill_random_trf(std::vector<transform3_t> &collection) {
 
   using vector_t = typename transform3_t::vector3;
@@ -64,7 +65,7 @@ inline void fill_random_trf(std::vector<transform3_t> &collection) {
 }
 
 /// Fill a @c std::array based matrix with random values
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 inline void fill_random_matrix(std::vector<matrix_t> &collection) {
 
   using scalar_t = typename matrix_t::value_type::value_type;
diff --git a/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp b/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp
index e17ce107..a99ef315 100644
--- a/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp
+++ b/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp
@@ -17,14 +17,14 @@
 
 namespace algebra {
 
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 void fill_random_matrix(std::vector<matrix_t>&);
 
-template <typename vector_t>
+template <concepts::vector vector_t>
 void fill_random_vec(std::vector<vector_t>&);
 
 /// Benchmark for matrix operations
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 struct matrix_bm : public benchmark_base {
 
   /// Prefix for the benchmark name
@@ -59,8 +59,9 @@ struct matrix_bm : public benchmark_base {
 };
 
 /// Benchmark operations on a single matrix (transpose, inverse etc)
-template <typename matrix_t, typename unaryOP>
-struct matrix_unaryOP_bm : public matrix_bm<matrix_t> {
+template <concepts::matrix matrix_t, typename unaryOP>
+requires std::invocable<unaryOP, matrix_t> struct matrix_unaryOP_bm
+    : public matrix_bm<matrix_t> {
   using base_type = matrix_bm<matrix_t>;
 
   matrix_unaryOP_bm() = delete;
@@ -90,8 +91,9 @@ struct matrix_unaryOP_bm : public matrix_bm<matrix_t> {
 };
 
 /// Benchmark elementwise addition/subtraction/multiplication of matrices
-template <typename matrix_t, typename binaryOP>
-struct matrix_binaryOP_bm : public matrix_bm<matrix_t> {
+template <concepts::matrix matrix_t, typename binaryOP>
+requires std::invocable<binaryOP, matrix_t, matrix_t> struct matrix_binaryOP_bm
+    : public matrix_bm<matrix_t> {
   using base_type = matrix_bm<matrix_t>;
 
   matrix_binaryOP_bm() = delete;
@@ -121,7 +123,7 @@ struct matrix_binaryOP_bm : public matrix_bm<matrix_t> {
 };
 
 /// Benchmark matrix vector multiplication
-template <typename matrix_t, typename vector_t>
+template <concepts::matrix matrix_t, concepts::vector vector_t>
 struct matrix_vector_bm : public matrix_bm<matrix_t> {
   using base_type = matrix_bm<matrix_t>;
 
@@ -162,21 +164,21 @@ namespace bench_op {
 
 struct add {
   inline static const std::string name{"add"};
-  template <typename matrix_t>
+  template <concepts::matrix matrix_t>
   matrix_t operator()(const matrix_t& a, const matrix_t& b) const {
     return a + b;
   }
 };
 struct sub {
   inline static const std::string name{"sub"};
-  template <typename matrix_t>
+  template <concepts::matrix matrix_t>
   matrix_t operator()(const matrix_t& a, const matrix_t& b) const {
     return a - b;
   }
 };
 struct mul {
   inline static const std::string name{"mul"};
-  template <typename matrix_t>
+  template <concepts::matrix matrix_t>
   matrix_t operator()(const matrix_t& a, const matrix_t& b) const {
     return a * b;
   }
@@ -184,7 +186,7 @@ struct mul {
 
 struct transpose {
   inline static const std::string name{"transpose"};
-  template <typename matrix_t>
+  template <concepts::matrix matrix_t>
   auto operator()(const matrix_t& a) const {
     return algebra::matrix::transpose(a);
   }
diff --git a/benchmarks/common/include/benchmark/common/benchmark_transform3.hpp b/benchmarks/common/include/benchmark/common/benchmark_transform3.hpp
index 1c0e84d3..0442b987 100644
--- a/benchmarks/common/include/benchmark/common/benchmark_transform3.hpp
+++ b/benchmarks/common/include/benchmark/common/benchmark_transform3.hpp
@@ -19,11 +19,11 @@
 
 namespace algebra {
 
-template <typename transform3_t>
+template <concepts::transform3D transform3_t>
 void fill_random_trf(std::vector<transform3_t>&);
 
 /// Benchmark for vector operations
-template <typename transform3_t>
+template <concepts::transform3D transform3_t>
 struct transform3_bm : public vector_bm<typename transform3_t::vector3> {
  private:
   using base_type = vector_bm<typename transform3_t::vector3>;
diff --git a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
index 9d672055..a49f3a1c 100644
--- a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
+++ b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
@@ -20,11 +20,11 @@
 
 namespace algebra {
 
-template <typename vector_t>
+template <concepts::vector vector_t>
 void fill_random_vec(std::vector<vector_t> &);
 
 /// Benchmark for vector operations
-template <typename vector_t>
+template <concepts::vector vector_t>
 struct vector_bm : public benchmark_base {
 
   /// Prefix for the benchmark name
@@ -57,9 +57,10 @@ struct vector_bm : public benchmark_base {
 };
 
 /// Benchmark elementwise addition of vectors
-template <template <typename> class vector_t, typename scalar_t,
+template <template <typename> class vector_t, concepts::scalar scalar_t,
           typename unaryOP>
-struct vector_unaryOP_bm : public vector_bm<vector_t<scalar_t>> {
+requires std::invocable<unaryOP, vector_t<scalar_t>> struct vector_unaryOP_bm
+    : public vector_bm<vector_t<scalar_t>> {
   using base_type = vector_bm<vector_t<scalar_t>>;
 
   vector_unaryOP_bm() = delete;
@@ -89,9 +90,11 @@ struct vector_unaryOP_bm : public vector_bm<vector_t<scalar_t>> {
 };
 
 /// Benchmark elementwise addition of vectors
-template <template <typename> class vector_t, typename scalar_t,
+template <template <typename> class vector_t, concepts::scalar scalar_t,
           typename binaryOP>
-struct vector_binaryOP_bm : public vector_bm<vector_t<scalar_t>> {
+requires std::invocable<binaryOP, vector_t<scalar_t>,
+                        vector_t<scalar_t>> struct vector_binaryOP_bm
+    : public vector_bm<vector_t<scalar_t>> {
   using base_type = vector_bm<vector_t<scalar_t>>;
 
   vector_binaryOP_bm() = delete;
@@ -126,28 +129,28 @@ namespace bench_op {
 
 struct add {
   inline static const std::string name{"add"};
-  template <typename vector_t>
+  template <concepts::vector vector_t>
   auto operator()(const vector_t &a, const vector_t &b) const {
     return a + b;
   }
 };
 struct sub {
   inline static const std::string name{"sub"};
-  template <typename vector_t>
+  template <concepts::vector vector_t>
   auto operator()(const vector_t &a, const vector_t &b) const {
     return a - b;
   }
 };
 struct dot {
   inline static const std::string name{"dot"};
-  template <typename vector_t>
+  template <concepts::vector vector_t>
   auto operator()(const vector_t &a, const vector_t &b) const {
     return algebra::vector::dot(a, b);
   }
 };
 struct cross {
   inline static const std::string name{"cross"};
-  template <typename vector_t>
+  template <concepts::vector vector_t>
   auto operator()(const vector_t &a, const vector_t &b) const {
     return algebra::vector::cross(a, b);
   }
@@ -157,7 +160,7 @@ struct cross {
 #define ALGEBRA_PLUGINS_BENCH_VECTOR(OP)       \
   struct OP {                                  \
     inline static const std::string name{#OP}; \
-    template <typename vector_t>               \
+    template <concepts::vector vector_t>       \
     auto operator()(const vector_t &a) const { \
       return algebra::vector::OP(a);           \
     }                                          \
diff --git a/benchmarks/common/include/benchmark/common/register_benchmark.hpp b/benchmarks/common/include/benchmark/common/register_benchmark.hpp
index 6eec7a7b..89b329bd 100644
--- a/benchmarks/common/include/benchmark/common/register_benchmark.hpp
+++ b/benchmarks/common/include/benchmark/common/register_benchmark.hpp
@@ -7,6 +7,9 @@
 
 #pragma once
 
+// Project include(s)
+#include "benchmark_base.hpp"
+
 // Benchmark include
 #include <benchmark/benchmark.h>
 
@@ -15,8 +18,9 @@
 
 namespace algebra {
 
-template <typename benchmark_t, typename config_t>
-void register_benchmark(const config_t& cfg, const std::string& suffix) {
+template <typename benchmark_t>
+requires std::derived_from<benchmark_t, benchmark_base> void register_benchmark(
+    const benchmark_base::configuration& cfg, const std::string& suffix) {
   benchmark_t bench{cfg};
 
   ::benchmark::RegisterBenchmark((bench.name() + suffix).c_str(), bench)
diff --git a/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp
index e6963b5c..914d4d09 100644
--- a/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp
+++ b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp
@@ -17,7 +17,7 @@
 namespace algebra {
 
 /// Fill an @c Eigen3 based vector with random values
-template <typename vector_t>
+template <concepts::vector vector_t>
 inline void fill_random_vec(std::vector<vector_t> &collection) {
 
   auto rand_obj = []() { return vector_t::Random(); };
@@ -27,7 +27,7 @@ inline void fill_random_vec(std::vector<vector_t> &collection) {
 }
 
 /// Fill a @c Eigen3 based transform3 with random values
-template <typename transform3_t>
+template <concepts::transform3D transform3_t>
 inline void fill_random_trf(std::vector<transform3_t> &collection) {
 
   using vector_t = typename transform3_t::vector3;
@@ -52,7 +52,7 @@ inline void fill_random_trf(std::vector<transform3_t> &collection) {
 }
 
 /// Fill a @c Eigen3 based matrix with random values
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 inline void fill_random_matrix(std::vector<matrix_t> &collection) {
 
   auto rand_obj = []() { return matrix_t::Random(); };
diff --git a/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp b/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp
index 51b3ade3..b4ebad37 100644
--- a/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp
+++ b/benchmarks/vc_aos/include/benchmark/vc_aos/data_generator.hpp
@@ -17,7 +17,7 @@
 namespace algebra {
 
 /// Fill a @c Vc::SimdArray based vector with random values
-template <typename vector_aos_t>
+template <concepts::vector vector_aos_t>
 inline void fill_random_vec(std::vector<vector_aos_t> &collection) {
 
   // Generate a vector of the right type with random values
@@ -30,7 +30,7 @@ inline void fill_random_vec(std::vector<vector_aos_t> &collection) {
 }
 
 /// Fill a @c Vc::SimdArray based transform3 with random values
-template <typename transform3_t>
+template <concepts::transform3D transform3_t>
 inline void fill_random_trf(std::vector<transform3_t> &collection) {
   // Generate a random, but valid affine transformation
   auto rand_obj = []() {
@@ -59,7 +59,7 @@ inline void fill_random_trf(std::vector<transform3_t> &collection) {
 }
 
 /// Fill a @c Vc::SimdArray based matrix with random values
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 inline void fill_random_matrix(std::vector<matrix_t> &collection) {
   // Generate a random, but valid affine transformation
   auto rand_obj = []() {
diff --git a/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp b/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp
index baaff335..91871df0 100644
--- a/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp
+++ b/benchmarks/vc_soa/include/benchmark/vc_soa/data_generator.hpp
@@ -17,11 +17,11 @@
 namespace algebra {
 
 /// Fill a @c Vc::Vector based vector with random values
-template <typename vector_soa_t>
+template <concepts::vector vector_soa_t>
 inline void fill_random_vec(std::vector<vector_soa_t> &collection) {
   // Generate a vector of the right type with random values
   auto rand_obj = []() {
-    using simd_vector_t = typename vector_soa_t::value_type;
+    using simd_vector_t = typename vector_soa_t::scalar_type;
     vector_soa_t tmp{};
     tmp[0] = simd_vector_t::Random();
     tmp[1] = simd_vector_t::Random();
@@ -34,11 +34,11 @@ inline void fill_random_vec(std::vector<vector_soa_t> &collection) {
 }
 
 /// Fill a @c Vc::Vector based transform3 with random values
-template <typename transform3_t>
+template <concepts::transform3D transform3_t>
 inline void fill_random_trf(std::vector<transform3_t> &collection) {
   // Generate a random, but valid affine transformation
   auto rand_obj = []() {
-    using simd_vector_t = typename transform3_t::value_type;
+    using simd_vector_t = typename transform3_t::scalar_type;
     typename transform3_t::vector3 x_axis, z_axis, t;
     x_axis[0] = simd_vector_t::Random();
     x_axis[1] = simd_vector_t::Random();
@@ -66,11 +66,11 @@ inline void fill_random_trf(std::vector<transform3_t> &collection) {
 }
 
 /// Fill a @c Vc::Vector based matrix with random values
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 inline void fill_random_matrix(std::vector<matrix_t> &collection) {
   // Generate a random, but valid affine transformation
   auto rand_obj = []() {
-    using simd_vector_t = typename matrix_t::value_type;
+    using simd_vector_t = typename matrix_t::scalar_type;
 
     matrix_t m;
 
@@ -80,8 +80,9 @@ inline void fill_random_matrix(std::vector<matrix_t> &collection) {
 
       for (std::size_t i = 0u; i < matrix_t::rows(); ++i) {
         v[i] = simd_vector_t::Random();
-        m[j] = v;
       }
+
+      m[j] = v;
     }
 
     return m;
diff --git a/common/include/algebra/concepts.hpp b/common/include/algebra/concepts.hpp
new file mode 100644
index 00000000..b66948ec
--- /dev/null
+++ b/common/include/algebra/concepts.hpp
@@ -0,0 +1,107 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2024 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "algebra/type_traits.hpp"
+
+// System include(s).
+#include <concepts>
+
+namespace algebra::concepts {
+
+// Value concept: Single entry
+template <typename T>
+concept value = std::is_arithmetic_v<std::decay_t<T>>;
+
+/// Scalar concept: Elements of vectors/matrices (can be simd vectors)
+template <typename T>
+concept scalar = !algebra::traits::is_matrix<T> &&
+                 !algebra::traits::is_vector<T> && requires(T a, T b) {
+  { a + b }
+  ->std::convertible_to<T>;
+  { a - b }
+  ->std::convertible_to<T>;
+  { a* b }
+  ->std::convertible_to<T>;
+  { a / b }
+  ->std::convertible_to<T>;
+};
+
+/// Check if a scalar is simd
+template <typename T>
+concept simd_scalar = scalar<T> && !std::is_scalar_v<T>;
+
+/// Index concept to access vector/matrix elements
+template <typename T>
+concept index = std::is_integral_v<T> && !std::same_as<T, bool>;
+
+/// Vector concepts
+/// @{
+template <typename V>
+concept vector = algebra::traits::is_vector<V>;
+
+template <typename V>
+concept vector2D = vector<V> && (algebra::traits::size<V> == 2);
+
+template <typename V>
+concept vector3D = vector<V> && (algebra::traits::size<V> == 3);
+/// @}
+
+/// Point concepts
+/// @{
+template <typename V>
+concept point = vector<V>;
+
+template <typename V>
+concept point2D = point<V> && (algebra::traits::size<V> == 2);
+
+template <typename V>
+concept point3D = point<V> && (algebra::traits::size<V> == 3);
+/// @}
+
+/// Matrix concepts
+/// @{
+template <typename M>
+concept matrix = algebra::traits::is_matrix<M>;
+
+template <typename M>
+concept square_matrix = matrix<M>&& algebra::traits::is_square<M>;
+
+template <typename M>
+concept row_matrix = matrix<M> && (algebra::traits::rows<M> == 1);
+
+template <typename M>
+concept row_matrix3D = row_matrix<M> && (algebra::traits::rows<M> == 3);
+
+template <typename M>
+concept column_matrix = matrix<M> && (algebra::traits::columns<M> == 1);
+
+template <typename M>
+concept column_matrix3D = column_matrix<M> && (algebra::traits::rows<M> == 3);
+/// @}
+
+/// Transform concept
+template <typename T>
+concept transform3D = requires(T trf) {
+  // Local type definitions
+  requires scalar<typename T::scalar_type>;
+  requires vector3D<typename T::vector3>;
+  requires point2D<typename T::point2>;
+  requires point3D<typename T::point3>;
+
+  // Methods
+  trf.rotation();
+  trf.translation();
+  trf.point_to_global(typename T::vector3());
+  trf.point_to_local(typename T::vector3());
+  trf.vector_to_global(typename T::vector3());
+  trf.vector_to_local(typename T::vector3());
+};
+
+}  // namespace algebra::concepts
diff --git a/common/include/algebra/type_traits.hpp b/common/include/algebra/type_traits.hpp
index 5eac0c3e..3b338aad 100644
--- a/common/include/algebra/type_traits.hpp
+++ b/common/include/algebra/type_traits.hpp
@@ -12,7 +12,7 @@
 #include <cmath>
 #include <type_traits>
 
-namespace algebra::trait {
+namespace algebra::traits {
 
 /// Matrix traits
 /// @{
@@ -67,7 +67,15 @@ using get_matrix_t = typename matrix<M>::template other_type<T, ROWS, COLS>;
 /// Matrix dimensions
 /// @{
 template <typename M>
-struct dimensions {};
+struct dimensions {
+
+  using size_type = int;
+
+  // Error case
+  static constexpr size_type dim{-1};
+  static constexpr size_type rows{-1};
+  static constexpr size_type columns{-1};
+};
 
 /// Specilization for scalar types
 template <typename M>
@@ -81,13 +89,14 @@ requires std::is_fundamental_v<M> struct dimensions<M> {
 };
 
 template <class M>
-inline constexpr index_t<M> dim{dimensions<M>::sim};
+inline constexpr index_t<M> dim{dimensions<std::remove_cvref_t<M>>::sim};
 
 template <class M>
-inline constexpr index_t<M> rows{dimensions<M>::rows};
+inline constexpr index_t<M> rows{dimensions<std::remove_cvref_t<M>>::rows};
 
 template <class M>
-inline constexpr index_t<M> columns{dimensions<M>::columns};
+inline constexpr index_t<M> columns{
+    dimensions<std::remove_cvref_t<M>>::columns};
 
 template <class M>
 inline constexpr index_t<M> rank{std::min(rows<M>, columns<M>)};
@@ -95,11 +104,11 @@ inline constexpr index_t<M> rank{std::min(rows<M>, columns<M>)};
 template <class M>
 inline constexpr index_t<M> size{rows<M> * columns<M>};
 
-template <class M>
-inline constexpr bool is_matrix{dimensions<M>::dim == 2};
+template <class V>
+inline constexpr bool is_vector{dimensions<std::remove_cvref_t<V>>::dim == 1};
 
 template <class M>
-inline constexpr bool is_vector{dimensions<M>::dim == 1};
+inline constexpr bool is_matrix{dimensions<std::remove_cvref_t<M>>::dim == 2};
 
 template <class M>
 inline constexpr bool is_square{(rows<M> == columns<M>)};
@@ -120,13 +129,13 @@ template <class M>
 using block_getter_t = typename block_getter<M>::type;
 /// @}
 
-}  // namespace algebra::trait
+}  // namespace algebra::traits
 
 /// Default type trait specializations
 /// @{
 #define ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(A)                           \
                                                                         \
-  namespace trait {                                                     \
+  namespace traits {                                                    \
                                                                         \
   template <typename T, auto N>                                         \
   struct index<A::vector_type<T, N>> {                                  \
@@ -217,4 +226,4 @@ using block_getter_t = typename block_getter<M>::type;
     using type = A::block_getter;                                       \
   };                                                                    \
                                                                         \
-  }  // namespace algebra::trait
+  }  // namespace algebra::traits
diff --git a/frontend/array_cmath/include/algebra/array_cmath.hpp b/frontend/array_cmath/include/algebra/array_cmath.hpp
index 177787c3..ff66006c 100644
--- a/frontend/array_cmath/include/algebra/array_cmath.hpp
+++ b/frontend/array_cmath/include/algebra/array_cmath.hpp
@@ -58,6 +58,26 @@ using generic::math::theta;
 
 }  // namespace vector
 
+// Use special algorithms for 4 dimensional matrices
+namespace generic {
+
+// Determinant algorithms
+template <concepts::scalar T, auto ROWS, auto COLS>
+struct determinant_selector<4, array::matrix_type<T, ROWS, COLS>> {
+  using type =
+      matrix::determinant::hard_coded<array::matrix_type<T, ROWS, COLS>,
+                                      array::element_getter>;
+};
+
+// Inversion algorithms
+template <concepts::scalar T, auto ROWS, auto COLS>
+struct inversion_selector<4, array::matrix_type<T, ROWS, COLS>> {
+  using type = matrix::inverse::hard_coded<array::matrix_type<T, ROWS, COLS>,
+                                           array::element_getter>;
+};
+
+}  // namespace generic
+
 namespace matrix {
 
 /// @name Matrix functions on @c algebra::array::storage_type
@@ -81,7 +101,7 @@ namespace array {
 /// @name cmath based transforms on @c algebra::array
 /// @{
 
-template <typename T>
+template <concepts::scalar T>
 using transform3 =
     generic::math::transform3<array::size_type, T, array::matrix_type,
                               array::storage_type>;
diff --git a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp
index 235a2900..630d3323 100644
--- a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp
+++ b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp
@@ -75,7 +75,7 @@ using eigen::math::zero;
 
 namespace eigen {
 
-template <typename T>
+template <concepts::scalar T>
 using transform3 = math::transform3<T>;
 
 }  // namespace eigen
diff --git a/frontend/eigen_generic/include/algebra/eigen_generic.hpp b/frontend/eigen_generic/include/algebra/eigen_generic.hpp
index 4ab351c2..b652b1de 100644
--- a/frontend/eigen_generic/include/algebra/eigen_generic.hpp
+++ b/frontend/eigen_generic/include/algebra/eigen_generic.hpp
@@ -80,7 +80,7 @@ namespace eigen {
 /// @name generic based transforms on @c algebra::eigen
 /// @{
 
-template <typename T>
+template <concepts::scalar T>
 using transform3 =
     generic::math::transform3<eigen::size_type, T, eigen::matrix_type,
                               eigen::storage_type>;
diff --git a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp
index 39b2f21f..94cbaa2d 100644
--- a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp
+++ b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp
@@ -78,7 +78,7 @@ namespace fastor {
 /// @name Transform on @c algebra::fastor::storage_type
 /// @{
 
-template <typename T>
+template <concepts::scalar T>
 using transform3 = math::transform3<T>;
 
 /// @}
diff --git a/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp b/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp
index 9e579407..4ce28361 100644
--- a/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp
+++ b/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp
@@ -71,7 +71,7 @@ namespace smatrix {
 /// @name generic based transforms on @c algebra::smatrix
 /// @{
 
-template <typename T>
+template <concepts::scalar T>
 using transform3 =
     generic::math::transform3<smatrix::size_type, T, smatrix::matrix_type,
                               smatrix::storage_type>;
diff --git a/frontend/vc_aos/include/algebra/vc_aos.hpp b/frontend/vc_aos/include/algebra/vc_aos.hpp
index ef52aa6b..e6a05851 100644
--- a/frontend/vc_aos/include/algebra/vc_aos.hpp
+++ b/frontend/vc_aos/include/algebra/vc_aos.hpp
@@ -77,7 +77,7 @@ namespace vc_aos {
 /// @name Vc based transforms on @c algebra::vc_aos::storage_type
 /// @{
 
-template <typename T>
+template <concepts::value T>
 using transform3 = math::transform3<vc_aos::storage_type, T>;
 
 /// @}
diff --git a/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp b/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp
index 9a69860f..ac5f3431 100644
--- a/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp
+++ b/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp
@@ -79,7 +79,7 @@ namespace vc_aos {
 /// @name generic based transforms on @c algebra::vc_aos
 /// @{
 
-template <typename T>
+template <concepts::value T>
 using transform3 =
     generic::math::transform3<vc_aos::size_type, T, vc_aos::matrix_type,
                               vc_aos::storage_type>;
diff --git a/frontend/vc_soa/include/algebra/vc_soa.hpp b/frontend/vc_soa/include/algebra/vc_soa.hpp
index a4277882..90a0ff73 100644
--- a/frontend/vc_soa/include/algebra/vc_soa.hpp
+++ b/frontend/vc_soa/include/algebra/vc_soa.hpp
@@ -78,7 +78,7 @@ namespace vc_soa {
 /// @name Vc based transforms on @c algebra::vc_soa types
 /// @{
 
-template <typename T>
+template <concepts::value T>
 using transform3 =
     algebra::vc_aos::math::transform3<algebra::vc_soa::storage_type, T>;
 
diff --git a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp
index edf65628..dfd7f413 100644
--- a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp
+++ b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp
@@ -58,6 +58,26 @@ using generic::math::theta;
 
 }  // namespace vector
 
+// Use special algorithms for 4 dimensional matrices
+namespace generic {
+
+// Determinant algorithms
+template <concepts::scalar T, auto ROWS, auto COLS>
+struct determinant_selector<4, vecmem::matrix_type<T, ROWS, COLS>> {
+  using type =
+      matrix::determinant::hard_coded<vecmem::matrix_type<T, ROWS, COLS>,
+                                      vecmem::element_getter>;
+};
+
+// Inversion algorithms
+template <concepts::scalar T, auto ROWS, auto COLS>
+struct inversion_selector<4, vecmem::matrix_type<T, ROWS, COLS>> {
+  using type = matrix::inverse::hard_coded<vecmem::matrix_type<T, ROWS, COLS>,
+                                           vecmem::element_getter>;
+};
+
+}  // namespace generic
+
 namespace matrix {
 
 /// @name Matrix functions on @c algebra::vecmem::storage_type
@@ -81,7 +101,7 @@ namespace vecmem {
 /// @name cmath based transforms on @c algebra::vecmem
 /// @{
 
-template <typename T>
+template <concepts::scalar T>
 using transform3 =
     generic::math::transform3<std::size_t, T, vecmem::matrix_type,
                               vecmem::storage_type>;
diff --git a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp
index 67ea5317..69ca0353 100644
--- a/math/cmath/include/algebra/math/impl/cmath_matrix.hpp
+++ b/math/cmath/include/algebra/math/impl/cmath_matrix.hpp
@@ -8,26 +8,26 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
-#include "algebra/type_traits.hpp"
 
 namespace algebra::cmath {
 
 /// @returns zero matrix of type @tparam matrix_t
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 requires(std::is_scalar_v<typename matrix_t::value_type::value_type>)
     ALGEBRA_HOST_DEVICE inline matrix_t zero() {
   return matrix_t{};
 }
 
 /// @returns identity matrix of type @tparam matrix_t
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 requires(std::is_scalar_v<typename matrix_t::value_type::value_type>)
     ALGEBRA_HOST_DEVICE inline auto identity() {
   auto ret{zero<matrix_t>()};
 
-  for (std::size_t i = 0; i < algebra::trait::rank<matrix_t>; ++i) {
+  for (std::size_t i = 0; i < algebra::traits::rank<matrix_t>; ++i) {
     ret[i][i] = 1;
   }
 
@@ -35,7 +35,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 <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr void set_zero(
     array_t<array_t<scalar_t, ROWS>, COLS> &m) {
@@ -43,7 +43,7 @@ ALGEBRA_HOST_DEVICE constexpr void set_zero(
 }
 
 /// Set @param m as identity matrix
-template <std::size_t ROWS, std::size_t COLS, typename scalar_t,
+template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr void set_identity(
     array_t<array_t<scalar_t, ROWS>, COLS> &m) {
diff --git a/math/cmath/include/algebra/math/impl/cmath_operators.hpp b/math/cmath/include/algebra/math/impl/cmath_operators.hpp
index bb4a681f..83af29a5 100644
--- a/math/cmath/include/algebra/math/impl/cmath_operators.hpp
+++ b/math/cmath/include/algebra/math/impl/cmath_operators.hpp
@@ -19,48 +19,32 @@ namespace algebra::cmath {
 /// @name Operators on 2-element arrays
 /// @{
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 2> operator*(
-    const array_t<scalar_t, 2> &a, float s) {
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar1_t, concepts::scalar scalar2_t>
+ALGEBRA_HOST_DEVICE inline array_t<scalar1_t, 2> operator*(
+    const array_t<scalar1_t, 2> &a, scalar2_t s) {
 
-  return {a[0] * static_cast<scalar_t>(s), a[1] * static_cast<scalar_t>(s)};
+  return {a[0] * static_cast<scalar1_t>(s), a[1] * static_cast<scalar1_t>(s)};
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 2> operator*(
-    const array_t<scalar_t, 2> &a, double s) {
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar1_t, concepts::scalar scalar2_t>
+ALGEBRA_HOST_DEVICE inline array_t<scalar1_t, 2> operator*(
+    scalar2_t s, const array_t<scalar1_t, 2> &a) {
 
-  return {a[0] * static_cast<scalar_t>(s), a[1] * static_cast<scalar_t>(s)};
+  return {static_cast<scalar1_t>(s) * a[0], static_cast<scalar1_t>(s) * a[1]};
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 2> operator*(
-    float s, const array_t<scalar_t, 2> &a) {
-
-  return {static_cast<scalar_t>(s) * a[0], static_cast<scalar_t>(s) * a[1]};
-}
-
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 2> operator*(
-    double s, const array_t<scalar_t, 2> &a) {
-
-  return {static_cast<scalar_t>(s) * a[0], static_cast<scalar_t>(s) * a[1]};
-}
-
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 2> operator-(
     const array_t<scalar_t, 2> &a, const array_t<scalar_t, 2> &b) {
 
   return {a[0] - b[0], a[1] - b[1]};
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 2> operator+(
     const array_t<scalar_t, 2> &a, const array_t<scalar_t, 2> &b) {
 
@@ -72,52 +56,34 @@ ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 2> operator+(
 /// @name Operators on 3-element arrays
 /// @{
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 3> operator*(
-    const array_t<scalar_t, 3> &a, float s) {
-
-  return {a[0] * static_cast<scalar_t>(s), a[1] * static_cast<scalar_t>(s),
-          a[2] * static_cast<scalar_t>(s)};
-}
-
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 3> operator*(
-    const array_t<scalar_t, 3> &a, double s) {
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar1_t, concepts::scalar scalar2_t>
+ALGEBRA_HOST_DEVICE inline array_t<scalar1_t, 3> operator*(
+    const array_t<scalar1_t, 3> &a, scalar2_t s) {
 
-  return {a[0] * static_cast<scalar_t>(s), a[1] * static_cast<scalar_t>(s),
-          a[2] * static_cast<scalar_t>(s)};
+  return {a[0] * static_cast<scalar1_t>(s), a[1] * static_cast<scalar1_t>(s),
+          a[2] * static_cast<scalar1_t>(s)};
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 3> operator*(
-    float s, const array_t<scalar_t, 3> &a) {
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar1_t, concepts::scalar scalar2_t>
+ALGEBRA_HOST_DEVICE inline array_t<scalar1_t, 3> operator*(
+    scalar2_t s, const array_t<scalar1_t, 3> &a) {
 
-  return {static_cast<scalar_t>(s) * a[0], static_cast<scalar_t>(s) * a[1],
-          static_cast<scalar_t>(s) * a[2]};
+  return {static_cast<scalar1_t>(s) * a[0], static_cast<scalar1_t>(s) * a[1],
+          static_cast<scalar1_t>(s) * a[2]};
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 3> operator*(
-    double s, const array_t<scalar_t, 3> &a) {
-
-  return {static_cast<scalar_t>(s) * a[0], static_cast<scalar_t>(s) * a[1],
-          static_cast<scalar_t>(s) * a[2]};
-}
-
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 3> operator-(
     const array_t<scalar_t, 3> &a, const array_t<scalar_t, 3> &b) {
 
   return {a[0] - b[0], a[1] - b[1], a[2] - b[2]};
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 3> operator+(
     const array_t<scalar_t, 3> &a, const array_t<scalar_t, 3> &b) {
 
@@ -130,50 +96,16 @@ ALGEBRA_HOST_DEVICE inline array_t<scalar_t, 3> operator+(
 /// @{
 
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type ROWS, size_type COLS,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator*(
-    const array_t<array_t<scalar_t, ROWS>, COLS> &a, float s) {
-
-  array_t<array_t<scalar_t, ROWS>, COLS> ret;
-
-  for (size_type j = 0; j < COLS; ++j) {
-    for (size_type i = 0; i < ROWS; ++i) {
-      ret[j][i] = a[j][i] * static_cast<scalar_t>(s);
-    }
-  }
-
-  return ret;
-}
-
-template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type ROWS, size_type COLS,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator*(
-    const array_t<array_t<scalar_t, ROWS>, COLS> &a, double s) {
-
-  array_t<array_t<scalar_t, ROWS>, COLS> ret;
-
-  for (size_type j = 0; j < COLS; ++j) {
-    for (size_type i = 0; i < ROWS; ++i) {
-      ret[j][i] = a[j][i] * static_cast<scalar_t>(s);
-    }
-  }
-
-  return ret;
-}
-
-template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type ROWS, size_type COLS,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator*(
-    float s, const array_t<array_t<scalar_t, ROWS>, COLS> &a) {
+          concepts::scalar scalar1_t, concepts::scalar scalar2_t,
+          size_type ROWS, size_type COLS>
+ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar1_t, ROWS>, COLS> operator*(
+    const array_t<array_t<scalar1_t, ROWS>, COLS> &a, scalar2_t s) {
 
-  array_t<array_t<scalar_t, ROWS>, COLS> ret;
+  array_t<array_t<scalar1_t, ROWS>, COLS> ret;
 
   for (size_type j = 0; j < COLS; ++j) {
     for (size_type i = 0; i < ROWS; ++i) {
-      ret[j][i] = a[j][i] * static_cast<scalar_t>(s);
+      ret[j][i] = a[j][i] * static_cast<scalar1_t>(s);
     }
   }
 
@@ -181,16 +113,16 @@ ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator*(
 }
 
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type ROWS, size_type COLS,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
-ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator*(
-    double s, const array_t<array_t<scalar_t, ROWS>, COLS> &a) {
+          concepts::scalar scalar1_t, concepts::scalar scalar2_t,
+          size_type ROWS, size_type COLS>
+ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar1_t, ROWS>, COLS> operator*(
+    scalar2_t s, const array_t<array_t<scalar1_t, ROWS>, COLS> &a) {
 
-  array_t<array_t<scalar_t, ROWS>, COLS> ret;
+  array_t<array_t<scalar1_t, ROWS>, COLS> ret;
 
   for (size_type j = 0; j < COLS; ++j) {
     for (size_type i = 0; i < ROWS; ++i) {
-      ret[j][i] = a[j][i] * static_cast<scalar_t>(s);
+      ret[j][i] = a[j][i] * static_cast<scalar1_t>(s);
     }
   }
 
@@ -198,8 +130,7 @@ ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator*(
 }
 
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type M, size_type N, size_type O,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+          concepts::scalar scalar_t, size_type M, size_type N, size_type O>
 ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, M>, O> operator*(
     const array_t<array_t<scalar_t, M>, N> &A,
     const array_t<array_t<scalar_t, N>, O> &B) {
@@ -224,8 +155,7 @@ ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, M>, O> operator*(
 }
 
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type ROWS, size_type COLS,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+          concepts::scalar scalar_t, size_type ROWS, size_type COLS>
 ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator+(
     const array_t<array_t<scalar_t, ROWS>, COLS> &A,
     const array_t<array_t<scalar_t, ROWS>, COLS> &B) {
@@ -242,8 +172,7 @@ ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator+(
 }
 
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type ROWS, size_type COLS,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+          concepts::scalar scalar_t, size_type ROWS, size_type COLS>
 ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator-(
     const array_t<array_t<scalar_t, ROWS>, COLS> &A,
     const array_t<array_t<scalar_t, ROWS>, COLS> &B) {
@@ -265,8 +194,7 @@ ALGEBRA_HOST_DEVICE inline array_t<array_t<scalar_t, ROWS>, COLS> operator-(
 /// @{
 
 template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type ROWS, size_type COLS,
-          std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>
+          concepts::scalar scalar_t, size_type ROWS, size_type COLS>
 ALGEBRA_HOST_DEVICE inline array_t<scalar_t, ROWS> operator*(
     const array_t<array_t<scalar_t, ROWS>, COLS> &a,
     const array_t<scalar_t, COLS> &b) {
diff --git a/math/cmath/include/algebra/math/impl/cmath_vector.hpp b/math/cmath/include/algebra/math/impl/cmath_vector.hpp
index a4aad3e6..67e48608 100644
--- a/math/cmath/include/algebra/math/impl/cmath_vector.hpp
+++ b/math/cmath/include/algebra/math/impl/cmath_vector.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 
@@ -19,8 +20,9 @@ namespace algebra::cmath {
 /// @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>
+template <concepts::index size_type,
+          template <typename, size_type> class array_t,
+          concepts::scalar 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;
@@ -40,8 +42,9 @@ requires std::is_scalar_v<scalar_t> 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>
+template <concepts::index size_type,
+          template <typename, size_type> class array_t,
+          concepts::scalar 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) {
@@ -62,8 +65,9 @@ requires(COLS == 1 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
 /// @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>
+template <concepts::index size_type,
+          template <typename, size_type> class array_t,
+          concepts::scalar 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) {
@@ -84,8 +88,9 @@ requires(COLS == 1 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
 /// @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>
+template <concepts::index size_type,
+          template <typename, size_type> class array_t,
+          concepts::scalar 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) {
@@ -103,8 +108,9 @@ requires(COLS == 1 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
 /// This method retrieves the norm of a vector, no dimension restriction
 ///
 /// @param v the input vector
-template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N>
+template <concepts::index size_type,
+          template <typename, size_type> class array_t,
+          concepts::scalar 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) {
 
@@ -115,8 +121,9 @@ requires(N >= 2 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
 /// rows >= 3
 ///
 /// @param v the input vector
-template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N>
+template <concepts::index size_type,
+          template <typename, size_type> class array_t,
+          concepts::scalar 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 {
 
@@ -126,8 +133,9 @@ requires(N >= 3 && std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE
 /// Get a normalized version of the input vector
 ///
 /// @param v the input vector
-template <typename size_type, template <typename, size_type> class array_t,
-          typename scalar_t, size_type N>
+template <concepts::index size_type,
+          template <typename, size_type> class array_t,
+          concepts::scalar 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) {
diff --git a/math/common/include/algebra/math/common.hpp b/math/common/include/algebra/math/common.hpp
index de4731ce..a23aa357 100644
--- a/math/common/include/algebra/math/common.hpp
+++ b/math/common/include/algebra/math/common.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 
 // SYCL include(s).
@@ -29,43 +30,43 @@ namespace math_ns = std;
 #endif  // SYCL
 
 /// Absolute value of arg
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline auto fabs(scalar_t arg) {
   return math_ns::fabs(arg);
 }
 
 /// Fused multiply add
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline auto fma(scalar_t x, scalar_t y, scalar_t z) {
   return math_ns::fma(x, y, z);
 }
 
 /// Arc tangent of y/x
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline scalar_t atan2(scalar_t y, scalar_t x) {
   return math_ns::atan2(y, x);
 }
 
 /// Square root of arg
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline scalar_t sqrt(scalar_t arg) {
   return math_ns::sqrt(arg);
 }
 
 /// Inverse hyperbolic tangent of arg
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline scalar_t atanh(scalar_t arg) {
   return math_ns::atanh(arg);
 }
 
 /// Minimum of two values
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline scalar_t min(scalar_t a, scalar_t b) {
   return math_ns::min(a, b);
 }
 
 /// Maximum of two values
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline scalar_t max(scalar_t a, scalar_t b) {
   return math_ns::max(a, b);
 }
diff --git a/math/eigen/include/algebra/math/impl/eigen_matrix.hpp b/math/eigen/include/algebra/math/impl/eigen_matrix.hpp
index 35455ffe..014dd014 100644
--- a/math/eigen/include/algebra/math/impl/eigen_matrix.hpp
+++ b/math/eigen/include/algebra/math/impl/eigen_matrix.hpp
@@ -30,13 +30,13 @@
 namespace algebra::eigen::math {
 
 /// Create zero matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE inline matrix_t zero() {
   return matrix_t::Zero();
 }
 
 /// Create identity matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE inline matrix_t identity() {
   return matrix_t::Identity();
 }
diff --git a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp
index 0c155944..13c54d85 100644
--- a/math/eigen/include/algebra/math/impl/eigen_transform3.hpp
+++ b/math/eigen/include/algebra/math/impl/eigen_transform3.hpp
@@ -9,6 +9,7 @@
 
 // Project include(s).
 #include "algebra/qualifiers.hpp"
+#include "algebra/storage/impl/eigen_array.hpp"
 
 // Eigen include(s).
 #ifdef _MSC_VER
@@ -35,7 +36,7 @@
 namespace algebra::eigen::math {
 
 /// Transform wrapper class to ensure standard API within differnt plugins
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 struct transform3 {
 
   /// @name Type definitions for the struct
@@ -46,7 +47,7 @@ struct transform3 {
 
   /// Array type used by the transform
   template <int N>
-  using array_type = Eigen::Matrix<scalar_type, N, 1, 0, N, 1>;
+  using array_type = eigen::array<scalar_type, N>;
 
   /// 3-element "vector" type
   using vector3 = array_type<3>;
@@ -159,15 +160,13 @@ struct transform3 {
   ///
   /// @param m is the rotation matrix
   /// @param v is the vector to be rotated
-  template <
-      typename derived_type,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3,
-                       bool> = true,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::ColsAtCompileTime == 1,
-                       bool> = true>
-  ALGEBRA_HOST_DEVICE static inline auto rotate(
-      const Eigen::Transform<scalar_type, 3, Eigen::Affine> &m,
-      const Eigen::MatrixBase<derived_type> &v) {
+  template <typename derived_type>
+  requires(Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3 &&
+           Eigen::MatrixBase<derived_type>::ColsAtCompileTime ==
+               1) ALGEBRA_HOST_DEVICE
+      static inline auto rotate(
+          const Eigen::Transform<scalar_type, 3, Eigen::Affine> &m,
+          const Eigen::MatrixBase<derived_type> &v) {
 
     return m.matrix().template block<3, 3>(0, 0) * v;
   }
@@ -193,7 +192,7 @@ struct transform3 {
 
   /// This method retrieves the translation of a transform
   ALGEBRA_HOST_DEVICE
-  inline auto translation() const {
+  inline point3 translation() const {
 
     return _data.matrix().template block<3, 1>(0, 3);
   }
@@ -208,56 +207,48 @@ struct transform3 {
 
   /// This method transform from a point from the local 3D cartesian frame to
   /// the global 3D cartesian frame
-  template <
-      typename derived_type,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3,
-                       bool> = true,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::ColsAtCompileTime == 1,
-                       bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto point_to_global(
-      const Eigen::MatrixBase<derived_type> &v) const {
+  template <typename derived_type>
+  requires(Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3 &&
+           Eigen::MatrixBase<derived_type>::ColsAtCompileTime ==
+               1) ALGEBRA_HOST_DEVICE
+      inline auto point_to_global(
+          const Eigen::MatrixBase<derived_type> &v) const {
 
     return (_data * v);
   }
 
   /// This method transform from a vector from the global 3D cartesian frame
   /// into the local 3D cartesian frame
-  template <
-      typename derived_type,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3,
-                       bool> = true,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::ColsAtCompileTime == 1,
-                       bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto point_to_local(
-      const Eigen::MatrixBase<derived_type> &v) const {
+  template <typename derived_type>
+  requires(Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3 &&
+           Eigen::MatrixBase<derived_type>::ColsAtCompileTime ==
+               1) ALGEBRA_HOST_DEVICE
+      inline auto point_to_local(
+          const Eigen::MatrixBase<derived_type> &v) const {
 
     return (_data_inv * v);
   }
 
   /// This method transform from a vector from the local 3D cartesian frame to
   /// the global 3D cartesian frame
-  template <
-      typename derived_type,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3,
-                       bool> = true,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::ColsAtCompileTime == 1,
-                       bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto vector_to_global(
-      const Eigen::MatrixBase<derived_type> &v) const {
+  template <typename derived_type>
+  requires(Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3 &&
+           Eigen::MatrixBase<derived_type>::ColsAtCompileTime ==
+               1) ALGEBRA_HOST_DEVICE
+      inline auto vector_to_global(
+          const Eigen::MatrixBase<derived_type> &v) const {
 
     return (_data.linear() * v);
   }
 
   /// This method transform from a vector from the global 3D cartesian frame
   /// into the local 3D cartesian frame
-  template <
-      typename derived_type,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3,
-                       bool> = true,
-      std::enable_if_t<Eigen::MatrixBase<derived_type>::ColsAtCompileTime == 1,
-                       bool> = true>
-  ALGEBRA_HOST_DEVICE inline auto vector_to_local(
-      const Eigen::MatrixBase<derived_type> &v) const {
+  template <typename derived_type>
+  requires(Eigen::MatrixBase<derived_type>::RowsAtCompileTime == 3 &&
+           Eigen::MatrixBase<derived_type>::ColsAtCompileTime ==
+               1) ALGEBRA_HOST_DEVICE
+      inline auto vector_to_local(
+          const Eigen::MatrixBase<derived_type> &v) const {
 
     return (_data_inv.linear() * v);
   }
diff --git a/math/fastor/include/algebra/math/impl/fastor_matrix.hpp b/math/fastor/include/algebra/math/impl/fastor_matrix.hpp
index bfeff7fb..38a9ca76 100644
--- a/math/fastor/include/algebra/math/impl/fastor_matrix.hpp
+++ b/math/fastor/include/algebra/math/impl/fastor_matrix.hpp
@@ -23,17 +23,17 @@
 namespace algebra::fastor::math {
 
 /// Create zero matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE inline matrix_t zero() {
   return matrix_t(0);
 }
 
 /// Create identity matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE inline matrix_t identity() {
-  using scalar_t = algebra::trait::value_t<matrix_t>;
-  constexpr auto rows{algebra::trait::rows<matrix_t>};
-  constexpr auto cols{algebra::trait::columns<matrix_t>};
+  using scalar_t = algebra::traits::value_t<matrix_t>;
+  constexpr auto rows{algebra::traits::rows<matrix_t>};
+  constexpr auto cols{algebra::traits::columns<matrix_t>};
 
   if constexpr (rows >= cols) {
     matrix_type<scalar_t, rows, rows> identity_matrix;
@@ -49,13 +49,13 @@ ALGEBRA_HOST_DEVICE inline matrix_t identity() {
 }
 
 /// Set input matrix as zero matrix
-template <std::size_t ROWS, std::size_t COLS, typename scalar_t>
+template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type<scalar_t, ROWS, COLS> &m) {
   m.zeros();
 }
 
 /// Set input matrix as identity matrix
-template <std::size_t ROWS, std::size_t COLS, typename scalar_t>
+template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline void set_identity(
     matrix_type<scalar_t, ROWS, COLS> &m) {
 
@@ -63,7 +63,7 @@ ALGEBRA_HOST_DEVICE inline void set_identity(
 }
 
 /// Create transpose matrix
-template <std::size_t ROWS, std::size_t COLS, typename scalar_t>
+template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline matrix_type<scalar_t, COLS, ROWS> transpose(
     const matrix_type<scalar_t, ROWS, COLS> &m) {
 
@@ -71,7 +71,7 @@ ALGEBRA_HOST_DEVICE inline matrix_type<scalar_t, COLS, ROWS> transpose(
 }
 
 /// @returns the determinant of @param m
-template <std::size_t N, typename scalar_t>
+template <std::size_t N, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline scalar_t determinant(
     const matrix_type<scalar_t, N, N> &m) {
 
@@ -79,7 +79,7 @@ ALGEBRA_HOST_DEVICE inline scalar_t determinant(
 }
 
 /// @returns the inverse of @param m
-template <std::size_t ROWS, std::size_t COLS, typename scalar_t>
+template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline matrix_type<scalar_t, COLS, ROWS> inverse(
     const matrix_type<scalar_t, ROWS, COLS> &m) {
 
diff --git a/math/fastor/include/algebra/math/impl/fastor_transform3.hpp b/math/fastor/include/algebra/math/impl/fastor_transform3.hpp
index b5fedecf..b36f5d2d 100644
--- a/math/fastor/include/algebra/math/impl/fastor_transform3.hpp
+++ b/math/fastor/include/algebra/math/impl/fastor_transform3.hpp
@@ -25,7 +25,7 @@
 namespace algebra::fastor::math {
 
 /// Transform wrapper class to ensure standard API towards the Fastor library
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 struct transform3 {
   /// @name Type definitions for the struct
   /// @{
diff --git a/math/fastor/include/algebra/math/impl/fastor_vector.hpp b/math/fastor/include/algebra/math/impl/fastor_vector.hpp
index 08c2d603..3bf740e5 100644
--- a/math/fastor/include/algebra/math/impl/fastor_vector.hpp
+++ b/math/fastor/include/algebra/math/impl/fastor_vector.hpp
@@ -24,7 +24,7 @@ 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>
+template <concepts::scalar scalar_t, auto N>
 requires(N >= 3) ALGEBRA_HOST inline scalar_t
     theta(const Fastor::Tensor<scalar_t, N> &v) noexcept {
 
@@ -34,7 +34,7 @@ requires(N >= 3) ALGEBRA_HOST inline scalar_t
 /// This method retrieves the perpenticular magnitude of a vector with rows >= 2
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 requires(N >= 2) ALGEBRA_HOST inline scalar_t
     perp(const Fastor::Tensor<scalar_t, N> &v) noexcept {
 
@@ -45,7 +45,7 @@ requires(N >= 2) ALGEBRA_HOST inline scalar_t
 /// This method retrieves the norm of a vector, no dimension restriction
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t norm(const Fastor::Tensor<scalar_t, N> &v) {
 
   return Fastor::norm(v);
@@ -55,7 +55,7 @@ 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>
+template <concepts::scalar scalar_t, auto N>
 requires(N >= 3) ALGEBRA_HOST inline scalar_t
     eta(const Fastor::Tensor<scalar_t, N> &v) noexcept {
 
@@ -65,7 +65,7 @@ requires(N >= 3) ALGEBRA_HOST inline scalar_t
 /// Get a normalized version of the input vector
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline Fastor::Tensor<scalar_t, N> normalize(
     const Fastor::Tensor<scalar_t, N> &v) {
 
@@ -78,7 +78,7 @@ ALGEBRA_HOST inline Fastor::Tensor<scalar_t, N> normalize(
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST_DEVICE inline scalar_t dot(const Fastor::Tensor<scalar_t, N> &a,
                                         const Fastor::Tensor<scalar_t, N> &b) {
   return Fastor::inner(a, b);
@@ -90,7 +90,7 @@ ALGEBRA_HOST_DEVICE inline scalar_t dot(const Fastor::Tensor<scalar_t, N> &a,
 /// @param b the second input Tensor<scalar_t, N, 1>
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor<scalar_t, N> &a,
                                  const Fastor::Tensor<scalar_t, N, 1> &b) {
 
@@ -107,7 +107,7 @@ ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor<scalar_t, N> &a,
 /// @param b the first input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
                                  const Fastor::Tensor<scalar_t, N> &b) {
 
@@ -121,7 +121,7 @@ ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
 /// @param b the first input Tensor<scalar_t, 3, 1>
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
                                  const Fastor::Tensor<scalar_t, N, 1> &b) {
 
@@ -135,7 +135,7 @@ ALGEBRA_HOST inline scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
 /// @param b the second input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline Fastor::Tensor<scalar_t, 3> cross(
     const Fastor::Tensor<scalar_t, 3> &a,
     const Fastor::Tensor<scalar_t, 3> &b) {
@@ -148,7 +148,7 @@ ALGEBRA_HOST_DEVICE inline Fastor::Tensor<scalar_t, 3> cross(
 /// @param b the second input Tensor<scalar_t, 3, 1>
 ///
 /// @return a vector representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST inline Fastor::Tensor<scalar_t, 3> cross(
     const Fastor::Tensor<scalar_t, 3> &a,
     const Fastor::Tensor<scalar_t, 3, 1> &b) {
@@ -166,7 +166,7 @@ ALGEBRA_HOST inline Fastor::Tensor<scalar_t, 3> cross(
 /// @param b the first input vector
 ///
 /// @return a vector representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST inline Fastor::Tensor<scalar_t, 3> cross(
     const Fastor::Tensor<scalar_t, 3, 1> &a,
     const Fastor::Tensor<scalar_t, 3> &b) {
@@ -181,7 +181,7 @@ ALGEBRA_HOST inline Fastor::Tensor<scalar_t, 3> cross(
 /// @param b the first input Tensor<scalar_t, 3, 1>
 ///
 /// @return a vector representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST inline Fastor::Tensor<scalar_t, 3> cross(
     const Fastor::Tensor<scalar_t, 3, 1> &a,
     const Fastor::Tensor<scalar_t, 3, 1> &b) {
diff --git a/math/generic/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp b/math/generic/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp
index 9d929c42..3670c7bf 100644
--- a/math/generic/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp
+++ b/math/generic/include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 
@@ -17,12 +18,12 @@
 namespace algebra::generic::matrix::decomposition {
 
 /// "Partial Pivot LU Decomposition", assuming a N X N matrix
-template <class matrix_t, class element_getter_t>
+template <concepts::matrix matrix_t, class element_getter_t>
 struct partial_pivot_lud {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
-  using vector_type = algebra::trait::vector_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
+  using vector_type = algebra::traits::vector_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
@@ -40,10 +41,10 @@ struct partial_pivot_lud {
     int n_pivot = 0;
   };
 
-  ALGEBRA_HOST_DEVICE inline lud<algebra::trait::rank<matrix_t>> operator()(
+  ALGEBRA_HOST_DEVICE inline lud<algebra::traits::rank<matrix_t>> operator()(
       const matrix_t& m) const {
 
-    constexpr size_type N{algebra::trait::rank<matrix_t>};
+    constexpr size_type N{algebra::traits::rank<matrix_t>};
 
     // LU decomposition matrix
     matrix_t lu = m;
diff --git a/math/generic/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp b/math/generic/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp
index 461c5338..73224834 100644
--- a/math/generic/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp
+++ b/math/generic/include/algebra/math/algorithms/matrix/determinant/cofactor.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/type_traits.hpp"
 
@@ -17,17 +18,17 @@
 namespace algebra::generic::matrix::determinant {
 
 /// "Determinant getter", assuming a N X N matrix
-template <class matrix_t, class element_getter_t>
+template <concepts::square_matrix matrix_t, class element_getter_t>
 struct cofactor {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
 
   ALGEBRA_HOST_DEVICE inline scalar_type operator()(const matrix_t &m) const {
-    return determinant_getter_helper<algebra::trait::rank<matrix_t>>()(m);
+    return determinant_getter_helper<algebra::traits::rank<matrix_t>>()(m);
   }
 
   template <size_type N, typename Enable = void>
diff --git a/math/generic/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp b/math/generic/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp
index 14262203..c6b9acc8 100644
--- a/math/generic/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp
+++ b/math/generic/include/algebra/math/algorithms/matrix/determinant/hard_coded.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/type_traits.hpp"
 
@@ -17,18 +18,18 @@
 namespace algebra::generic::matrix::determinant {
 
 /// "Determinant getter", assuming a N X N matrix
-template <class matrix_t, class element_getter_t>
+template <concepts::square_matrix matrix_t, class element_getter_t>
 struct hard_coded {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
 
   // 2 X 2 matrix determinant
   template <typename M = matrix_t>
-  requires(algebra::trait::rank<M> == 2) ALGEBRA_HOST_DEVICE inline scalar_type
+  requires(algebra::traits::rank<M> == 2) ALGEBRA_HOST_DEVICE inline scalar_type
   operator()(const matrix_t &m) const {
 
     return element_getter()(m, 0, 0) * element_getter()(m, 1, 1) -
@@ -37,7 +38,7 @@ struct hard_coded {
 
   // 4 X 4 matrix determinant
   template <typename M = matrix_t>
-  requires(algebra::trait::rank<M> == 4) ALGEBRA_HOST_DEVICE inline scalar_type
+  requires(algebra::traits::rank<M> == 4) ALGEBRA_HOST_DEVICE inline scalar_type
   operator()(const matrix_t &m) const {
 
     return element_getter()(m, 0, 3) * element_getter()(m, 1, 2) *
@@ -89,6 +90,6 @@ struct hard_coded {
            element_getter()(m, 0, 0) * element_getter()(m, 1, 1) *
                element_getter()(m, 2, 2) * element_getter()(m, 3, 3);
   }
-};
+};  // namespace algebra::generic::matrix::determinant
 
 }  // namespace algebra::generic::matrix::determinant
diff --git a/math/generic/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp b/math/generic/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp
index a9a17f33..663db232 100644
--- a/math/generic/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp
+++ b/math/generic/include/algebra/math/algorithms/matrix/determinant/partial_pivot_lud.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/type_traits.hpp"
@@ -15,11 +16,11 @@
 namespace algebra::generic::matrix::determinant {
 
 /// "Partial Pivot LU Decomposition", assuming a N X N matrix
-template <class matrix_t, class element_getter_t>
+template <concepts::square_matrix matrix_t, class element_getter_t>
 struct partial_pivot_lud {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
@@ -30,9 +31,10 @@ struct partial_pivot_lud {
 
   ALGEBRA_HOST_DEVICE inline scalar_type operator()(const matrix_t& m) const {
 
-    constexpr size_type N{algebra::trait::rank<matrix_t>};
+    constexpr size_type N{algebra::traits::rank<matrix_t>};
 
-    const typename decomposition_t::template lud<algebra::trait::rank<matrix_t>>
+    const typename decomposition_t::template lud<
+        algebra::traits::rank<matrix_t>>
         decomp_res = decomposition_t()(m);
 
     // Get the LU decomposition matrix equal to (L - I) + U
diff --git a/math/generic/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp b/math/generic/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp
index 40eeed50..e9589167 100644
--- a/math/generic/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp
+++ b/math/generic/include/algebra/math/algorithms/matrix/inverse/cofactor.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/algorithms/matrix/determinant/cofactor.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/type_traits.hpp"
@@ -20,17 +21,17 @@ namespace algebra::generic::matrix {
 namespace adjoint {
 
 /// "Adjoint getter", assuming a N X N matrix
-template <class matrix_t, class element_getter_t>
+template <concepts::square_matrix matrix_t, class element_getter_t>
 struct cofactor {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
 
   ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t &m) const {
-    return adjoint_getter_helper<algebra::trait::rank<matrix_t>>()(m);
+    return adjoint_getter_helper<algebra::traits::rank<matrix_t>>()(m);
   }
 
   template <size_type N, typename Enable = void>
@@ -94,8 +95,8 @@ namespace inverse {
 template <class matrix_t, class element_getter_t>
 struct cofactor {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
@@ -106,7 +107,7 @@ struct cofactor {
 
   ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t &m) const {
 
-    constexpr size_type N{algebra::trait::rank<matrix_t>};
+    constexpr size_type N{algebra::traits::rank<matrix_t>};
 
     matrix_t ret;
 
diff --git a/math/generic/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp b/math/generic/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp
index ba817a1c..9f92f9a8 100644
--- a/math/generic/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp
+++ b/math/generic/include/algebra/math/algorithms/matrix/inverse/hard_coded.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/algorithms/matrix/determinant/hard_coded.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/type_traits.hpp"
@@ -18,11 +19,11 @@
 namespace algebra::generic::matrix::inverse {
 
 /// "inverse getter", assuming a N X N matrix
-template <class matrix_t, class element_getter_t>
+template <concepts::square_matrix matrix_t, class element_getter_t>
 struct hard_coded {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
@@ -32,7 +33,7 @@ struct hard_coded {
 
   // 2 X 2 matrix inverse
   template <typename M = matrix_t>
-  requires(algebra::trait::rank<M> == 2) ALGEBRA_HOST_DEVICE inline matrix_t
+  requires(algebra::traits::rank<M> == 2) ALGEBRA_HOST_DEVICE inline matrix_t
   operator()(const matrix_t &m) const {
 
     matrix_t ret;
@@ -49,7 +50,7 @@ struct hard_coded {
 
   // 4 X 4 matrix inverse
   template <typename M = matrix_t>
-  requires(algebra::trait::rank<M> == 4) ALGEBRA_HOST_DEVICE inline matrix_t
+  requires(algebra::traits::rank<M> == 4) ALGEBRA_HOST_DEVICE inline matrix_t
   operator()(const matrix_t &m) const {
 
     matrix_t ret;
@@ -270,6 +271,6 @@ struct hard_coded {
     }
     return ret;
   }
-};
+};  // namespace algebra::generic::matrix::inverse
 
 }  // namespace algebra::generic::matrix::inverse
diff --git a/math/generic/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp b/math/generic/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp
index aa29c4f1..6c901969 100644
--- a/math/generic/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp
+++ b/math/generic/include/algebra/math/algorithms/matrix/inverse/partial_pivot_lud.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/type_traits.hpp"
@@ -15,11 +16,11 @@
 namespace algebra::generic::matrix::inverse {
 
 /// "Partial Pivot LU Decomposition", assuming a N X N matrix
-template <class matrix_t, class element_getter_t>
+template <concepts::square_matrix matrix_t, class element_getter_t>
 struct partial_pivot_lud {
 
-  using scalar_type = algebra::trait::value_t<matrix_t>;
-  using size_type = algebra::trait::index_t<matrix_t>;
+  using scalar_type = algebra::traits::value_t<matrix_t>;
+  using size_type = algebra::traits::index_t<matrix_t>;
 
   /// Function (object) used for accessing a matrix element
   using element_getter = element_getter_t;
@@ -30,7 +31,7 @@ struct partial_pivot_lud {
 
   ALGEBRA_HOST_DEVICE inline matrix_t operator()(const matrix_t& m) const {
 
-    constexpr size_type N{algebra::trait::rank<matrix_t>};
+    constexpr size_type N{algebra::traits::rank<matrix_t>};
 
     const typename decomposition_t::template lud<N> decomp_res =
         decomposition_t()(m);
diff --git a/math/generic/include/algebra/math/algorithms/utils/algorithm_finder.hpp b/math/generic/include/algebra/math/algorithms/utils/algorithm_finder.hpp
index 918b7d6c..43bb450e 100644
--- a/math/generic/include/algebra/math/algorithms/utils/algorithm_finder.hpp
+++ b/math/generic/include/algebra/math/algorithms/utils/algorithm_finder.hpp
@@ -17,26 +17,24 @@ namespace algebra::generic {
 
 /// Get the type of determinant algorithm acording to matrix dimension
 /// @{
+
+// Default algorithm
 template <std::size_t N, typename... Args>
 struct determinant_selector {
   using type = matrix::determinant::partial_pivot_lud<Args...>;
 };
 
+/// Always use hard coded implementation for very small matrices
 template <typename... Args>
 struct determinant_selector<2, Args...> {
   using type = matrix::determinant::hard_coded<Args...>;
 };
 
-template <typename... Args>
-struct determinant_selector<4, Args...> {
-  using type = matrix::determinant::hard_coded<Args...>;
-};
-
 /// @tparam M matrix type
-template <class M>
+template <concepts::square_matrix M>
 using determinant_t =
-    typename determinant_selector<algebra::trait::rank<M>, M,
-                                  algebra::trait::element_getter_t<M>>::type;
+    typename determinant_selector<algebra::traits::rank<M>, M,
+                                  algebra::traits::element_getter_t<M>>::type;
 /// @}
 
 /// Get the type of inversion algorithm acording to matrix dimension
@@ -46,21 +44,17 @@ struct inversion_selector {
   using type = matrix::inverse::partial_pivot_lud<Args...>;
 };
 
+/// Always use hard coded implementation for very small matrices
 template <typename... Args>
 struct inversion_selector<2, Args...> {
   using type = matrix::inverse::hard_coded<Args...>;
 };
 
-template <typename... Args>
-struct inversion_selector<4, Args...> {
-  using type = matrix::inverse::hard_coded<Args...>;
-};
-
 /// @tparam M matrix type
-template <class M>
+template <concepts::square_matrix M>
 using inversion_t =
-    typename inversion_selector<algebra::trait::rank<M>, M,
-                                algebra::trait::element_getter_t<M>>::type;
+    typename inversion_selector<algebra::traits::rank<M>, M,
+                                algebra::traits::element_getter_t<M>>::type;
 /// @}generic
 
 }  // namespace algebra::generic
diff --git a/math/generic/include/algebra/math/impl/generic_matrix.hpp b/math/generic/include/algebra/math/impl/generic_matrix.hpp
index 41c49864..1bb4ce8d 100644
--- a/math/generic/include/algebra/math/impl/generic_matrix.hpp
+++ b/math/generic/include/algebra/math/impl/generic_matrix.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/algorithms/utils/algorithm_finder.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
@@ -15,16 +16,16 @@
 namespace algebra::generic::math {
 
 /// Create zero matrix - generic transform3
-template <class matrix_t>
-ALGEBRA_HOST_DEVICE inline matrix_t zero() {
+template <concepts::matrix M>
+ALGEBRA_HOST_DEVICE inline M zero() {
 
-  using index_t = algebra::trait::index_t<matrix_t>;
-  using element_getter_t = algebra::trait::element_getter_t<matrix_t>;
+  using index_t = algebra::traits::index_t<M>;
+  using element_getter_t = algebra::traits::element_getter_t<M>;
 
-  matrix_t ret;
+  M ret;
 
-  for (index_t j = 0; j < algebra::trait::columns<matrix_t>; ++j) {
-    for (index_t i = 0; i < algebra::trait::rows<matrix_t>; ++i) {
+  for (index_t j = 0; j < algebra::traits::columns<M>; ++j) {
+    for (index_t i = 0; i < algebra::traits::rows<M>; ++i) {
       element_getter_t{}(ret, i, j) = 0;
     }
   }
@@ -33,15 +34,15 @@ ALGEBRA_HOST_DEVICE inline matrix_t zero() {
 }
 
 /// Create identity matrix - generic transform3
-template <class matrix_t>
-ALGEBRA_HOST_DEVICE inline matrix_t identity() {
+template <concepts::matrix M>
+ALGEBRA_HOST_DEVICE inline M identity() {
 
-  using index_t = algebra::trait::index_t<matrix_t>;
-  using element_getter_t = algebra::trait::element_getter_t<matrix_t>;
+  using index_t = algebra::traits::index_t<M>;
+  using element_getter_t = algebra::traits::element_getter_t<M>;
 
-  auto ret{zero<matrix_t>()};
+  auto ret{zero<M>()};
 
-  for (index_t i = 0; i < algebra::trait::rank<matrix_t>; ++i) {
+  for (index_t i = 0; i < algebra::traits::rank<M>; ++i) {
     element_getter_t{}(ret, i, i) = 1;
   }
 
@@ -49,29 +50,29 @@ ALGEBRA_HOST_DEVICE inline matrix_t identity() {
 }
 
 /// Set @param m as zero matrix
-template <class matrix_t>
-ALGEBRA_HOST_DEVICE inline void set_zero(matrix_t &m) {
-  m = zero<matrix_t>();
+template <concepts::matrix M>
+ALGEBRA_HOST_DEVICE inline void set_zero(M &m) {
+  m = zero<M>();
 }
 
 /// Set @param m as identity matrix
-template <class matrix_t>
-ALGEBRA_HOST_DEVICE inline void set_identity(matrix_t &m) {
-  m = identity<matrix_t>();
+template <concepts::matrix M>
+ALGEBRA_HOST_DEVICE inline void set_identity(M &m) {
+  m = identity<M>();
 }
 
 /// @returns the transpose matrix of @param m
-template <class matrix_t>
-ALGEBRA_HOST_DEVICE inline auto transpose(const matrix_t &m) {
+template <concepts::matrix M>
+ALGEBRA_HOST_DEVICE inline auto transpose(const M &m) {
 
-  using index_t = algebra::trait::index_t<matrix_t>;
-  using value_t = algebra::trait::value_t<matrix_t>;
-  using element_getter_t = algebra::trait::element_getter_t<matrix_t>;
+  using index_t = algebra::traits::index_t<M>;
+  using value_t = algebra::traits::value_t<M>;
+  using element_getter_t = algebra::traits::element_getter_t<M>;
 
-  constexpr index_t rows{algebra::trait::rows<matrix_t>};
-  constexpr index_t columns{algebra::trait::columns<matrix_t>};
+  constexpr index_t rows{algebra::traits::rows<M>};
+  constexpr index_t columns{algebra::traits::columns<M>};
 
-  algebra::trait::get_matrix_t<matrix_t, columns, rows, value_t> ret;
+  algebra::traits::get_matrix_t<M, columns, rows, value_t> ret;
 
   for (index_t i = 0; i < rows; ++i) {
     for (index_t j = 0; j < columns; ++j) {
@@ -83,18 +84,18 @@ ALGEBRA_HOST_DEVICE inline auto transpose(const matrix_t &m) {
 }
 
 /// @returns the determinant of @param m
-template <class matrix_t>
-ALGEBRA_HOST_DEVICE inline algebra::trait::scalar_t<matrix_t> determinant(
-    const matrix_t &m) {
+template <concepts::square_matrix M>
+ALGEBRA_HOST_DEVICE inline algebra::traits::scalar_t<M> determinant(
+    const M &m) {
 
-  return determinant_t<matrix_t>{}(m);
+  return determinant_t<M>{}(m);
 }
 
 /// @returns the determinant of @param m
-template <class matrix_t>
-ALGEBRA_HOST_DEVICE inline matrix_t inverse(const matrix_t &m) {
+template <concepts::square_matrix M>
+ALGEBRA_HOST_DEVICE inline M inverse(const M &m) {
 
-  return inversion_t<matrix_t>{}(m);
+  return inversion_t<M>{}(m);
 }
 
 }  // namespace algebra::generic::math
diff --git a/math/generic/include/algebra/math/impl/generic_transform3.hpp b/math/generic/include/algebra/math/impl/generic_transform3.hpp
index ee740e3b..ce6c164c 100644
--- a/math/generic/include/algebra/math/impl/generic_transform3.hpp
+++ b/math/generic/include/algebra/math/impl/generic_transform3.hpp
@@ -17,7 +17,7 @@
 namespace algebra::generic::math {
 
 /// Transform wrapper class to ensure standard API within differnt plugins
-template <typename index_t, typename scalar_t,
+template <concepts::index index_t, concepts::scalar scalar_t,
           template <typename, index_t, index_t> class matrix_t,
           template <typename, index_t> class array_t>
 struct transform3 {
@@ -34,19 +34,20 @@ struct transform3 {
 
   // 4 x 4 Matrix (keep four rows for better alignment)
   using matrix44 = matrix_t<scalar_t, 4, 4>;
+  static_assert(concepts::square_matrix<matrix44>);
 
   /// 3-element "vector" type
-  using vector3 = algebra::trait::get_vector_t<matrix44, 3, scalar_t>;
+  using vector3 = algebra::traits::get_vector_t<matrix44, 3, scalar_t>;
   /// Point in 3D space
   using point3 = vector3;
   /// Point in 2D space
-  using point2 = algebra::trait::get_vector_t<matrix44, 2, scalar_t>;
+  using point2 = algebra::traits::get_vector_t<matrix44, 2, scalar_t>;
 
   /// Function (object) used for accessing a matrix element
-  using element_getter = algebra::trait::element_getter_t<matrix44>;
+  using element_getter = algebra::traits::element_getter_t<matrix44>;
 
   /// Function (object) used for accessing a matrix element
-  using block_getter = algebra::trait::block_getter_t<matrix44>;
+  using block_getter = algebra::traits::block_getter_t<matrix44>;
 
   /// Matrix inversion algorithm
   using matrix_inversion =
diff --git a/math/generic/include/algebra/math/impl/generic_vector.hpp b/math/generic/include/algebra/math/impl/generic_vector.hpp
index 51256f95..ba526a3e 100644
--- a/math/generic/include/algebra/math/impl/generic_vector.hpp
+++ b/math/generic/include/algebra/math/impl/generic_vector.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/type_traits.hpp"
@@ -17,12 +18,11 @@ namespace algebra::generic::math {
 /// This method retrieves phi from a vector with rows >= 2
 ///
 /// @param v the input vector
-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 {
+template <concepts::vector vector_t>
+ALGEBRA_HOST_DEVICE inline algebra::traits::scalar_t<vector_t> phi(
+    const vector_t &v) noexcept {
 
-  using element_getter_t = algebra::trait::element_getter_t<vector_t>;
+  using element_getter_t = algebra::traits::element_getter_t<vector_t>;
 
   return algebra::math::atan2(element_getter_t{}(v, 1),
                               element_getter_t{}(v, 0));
@@ -31,12 +31,11 @@ requires algebra::trait::is_vector<vector_t>
 /// This method retrieves the perpendicular magnitude of a vector with rows >= 2
 ///
 /// @param v the input vector
-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 {
+template <concepts::vector vector_t>
+ALGEBRA_HOST_DEVICE inline algebra::traits::scalar_t<vector_t> perp(
+    const vector_t &v) noexcept {
 
-  using element_getter_t = algebra::trait::element_getter_t<vector_t>;
+  using element_getter_t = algebra::traits::element_getter_t<vector_t>;
 
   return algebra::math::sqrt(
       algebra::math::fma(element_getter_t{}(v, 0), element_getter_t{}(v, 0),
@@ -46,12 +45,11 @@ requires algebra::trait::is_vector<vector_t>
 /// This method retrieves theta from a vector with rows >= 3
 ///
 /// @param v the input vector
-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 {
+template <concepts::vector vector_t>
+ALGEBRA_HOST_DEVICE inline algebra::traits::scalar_t<vector_t> theta(
+    const vector_t &v) noexcept {
 
-  using element_getter_t = algebra::trait::element_getter_t<vector_t>;
+  using element_getter_t = algebra::traits::element_getter_t<vector_t>;
 
   return algebra::math::atan2(perp(v), element_getter_t{}(v, 2));
 }
@@ -66,12 +64,14 @@ requires algebra::trait::is_vector<vector_t>
 ///
 /// @return a vector representing the cross product
 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) {
+requires(
+    (concepts::vector3D<vector1_t> || concepts::column_matrix3D<vector1_t>)&&(
+        concepts::vector3D<vector2_t> ||
+        concepts::column_matrix3D<vector2_t>)) ALGEBRA_HOST_DEVICE
+    inline algebra::traits::vector_t<vector1_t> cross(const vector1_t &a,
+                                                      const vector2_t &b) {
 
-  using element_getter_t = algebra::trait::element_getter_t<vector1_t>;
+  using element_getter_t = algebra::traits::element_getter_t<vector1_t>;
 
   return {
       algebra::math::fma(element_getter_t{}(a, 1), element_getter_t{}(b, 2),
@@ -92,18 +92,19 @@ requires(algebra::trait::columns<vector1_t> == 1 &&
 ///
 /// @return the scalar dot product value
 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) {
+requires((concepts::vector<vector1_t> || concepts::column_matrix<vector1_t>)&&(
+    concepts::vector<vector2_t> ||
+    concepts::column_matrix<vector2_t>)) ALGEBRA_HOST_DEVICE
+    inline algebra::traits::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>;
-  using element_getter_t = algebra::trait::element_getter_t<vector1_t>;
+  using scalar_t = algebra::traits::scalar_t<vector1_t>;
+  using index_t = algebra::traits::index_t<vector1_t>;
+  using element_getter_t = algebra::traits::element_getter_t<vector1_t>;
 
   scalar_t ret{element_getter_t{}(a, 0) * element_getter_t{}(b, 0)};
 
-  for (index_t i = 1; i < algebra::trait::size<vector1_t>; i++) {
+  for (index_t i = 1; i < algebra::traits::size<vector1_t>; i++) {
     ret = std::fma(element_getter_t{}(a, i), element_getter_t{}(b, i), ret);
   }
 
@@ -113,10 +114,9 @@ requires(algebra::trait::columns<vector1_t> == 1 &&
 /// This method retrieves the norm of a vector with rows >= 3
 ///
 /// @param v the input vector
-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) {
+template <concepts::vector vector_t>
+ALGEBRA_HOST_DEVICE inline algebra::traits::scalar_t<vector_t> norm(
+    const vector_t &v) {
 
   return algebra::math::sqrt(dot(v, v));
 }
@@ -125,12 +125,11 @@ requires algebra::trait::is_vector<vector_t>
 /// rows >= 3
 ///
 /// @param v the input vector
-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 {
+template <concepts::vector vector_t>
+ALGEBRA_HOST_DEVICE inline algebra::traits::scalar_t<vector_t> eta(
+    const vector_t &v) noexcept {
 
-  using element_getter_t = algebra::trait::element_getter_t<vector_t>;
+  using element_getter_t = algebra::traits::element_getter_t<vector_t>;
 
   return algebra::math::atanh(element_getter_t{}(v, 2) / norm(v));
 }
@@ -142,10 +141,10 @@ requires algebra::trait::is_vector<vector_t>
 /// @param v the input vector
 ///
 /// @returns the normalized vector
-template <typename vector_t>
+template <concepts::vector vector_t>
 ALGEBRA_HOST_DEVICE inline vector_t normalize(const vector_t &v) {
 
-  using scalar_t = algebra::trait::scalar_t<vector_t>;
+  using scalar_t = algebra::traits::scalar_t<vector_t>;
 
   return v * static_cast<scalar_t>(1.) / norm(v);
 }
diff --git a/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp b/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp
index 20fc20bf..b937a5ef 100644
--- a/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp
+++ b/math/smatrix/include/algebra/math/impl/smatrix_matrix.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/smatrix.hpp"
 
@@ -17,53 +18,40 @@
 namespace algebra::smatrix::math {
 
 /// Create zero matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE inline matrix_t zero() {
   return matrix_t();
 }
 
 /// Create identity matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE inline matrix_t identity() {
   return matrix_t(ROOT::Math::SMatrixIdentity());
 }
 
-// Set input matrix as zero matrix
-template <unsigned int ROWS, unsigned int COLS, typename scalar_t>
-ALGEBRA_HOST_DEVICE inline void set_zero(matrix_type<scalar_t, ROWS, COLS> &m) {
-
-  for (unsigned int i = 0; i < ROWS; ++i) {
-    for (unsigned int j = 0; j < COLS; ++j) {
-      m(i, j) = 0;
-    }
-  }
+/// Set input matrix as zero matrix
+template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t>
+ALGEBRA_HOST_DEVICE inline void set_zero(
+    ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m) {
+  m = zero<ROOT::Math::SMatrix<scalar_t, ROWS, COLS>>();
 }
 
-// Set input matrix as identity matrix
-template <unsigned int ROWS, unsigned int COLS, typename scalar_t>
+/// Set input matrix as identity matrix
+template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline void set_identity(
-    matrix_type<scalar_t, ROWS, COLS> &m) {
-
-  for (unsigned int i = 0; i < ROWS; ++i) {
-    for (unsigned int j = 0; j < COLS; ++j) {
-      if (i == j) {
-        m(i, j) = 1;
-      } else {
-        m(i, j) = 0;
-      }
-    }
-  }
+    ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m) {
+  m = identity<ROOT::Math::SMatrix<scalar_t, ROWS, COLS>>();
 }
 
-// Create transpose matrix
-template <unsigned int ROWS, unsigned int COLS, typename scalar_t>
+/// Create transpose matrix
+template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline matrix_type<scalar_t, COLS, ROWS> transpose(
-    const matrix_type<scalar_t, ROWS, COLS> &m) {
+    const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m) {
   return ROOT::Math::Transpose(m);
 }
 
 /// @returns the determinant of @param m
-template <unsigned int ROWS, unsigned int COLS, typename scalar_t>
+template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline scalar_t determinant(
     const matrix_type<scalar_t, ROWS, COLS> &m) {
 
@@ -74,7 +62,7 @@ ALGEBRA_HOST_DEVICE inline scalar_t determinant(
 }
 
 /// @returns the inverse of @param m
-template <unsigned int ROWS, unsigned int COLS, typename scalar_t>
+template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline matrix_type<scalar_t, COLS, ROWS> inverse(
     const matrix_type<scalar_t, ROWS, COLS> &m) {
 
diff --git a/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp b/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp
index 3bf1a7aa..1e6eb772 100644
--- a/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp
+++ b/math/smatrix/include/algebra/math/impl/smatrix_transform3.hpp
@@ -18,7 +18,7 @@
 namespace algebra::smatrix::math {
 
 /// Transform wrapper class to ensure standard API towards the ROOT::SMatrix lib
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 struct transform3 {
 
   /// @name Type definitions for the struct
diff --git a/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp b/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp
index bb5fca1d..194fbc26 100644
--- a/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp
+++ b/math/smatrix/include/algebra/math/impl/smatrix_vector.hpp
@@ -21,14 +21,14 @@ 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>
+template <concepts::scalar 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>
+template <concepts::scalar 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 {
 
@@ -38,7 +38,7 @@ requires(N >= 2) ALGEBRA_HOST inline scalar_t
 /// This method retrieves theta from a vector, vector base with rows >= 3
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 requires(N >= 3) ALGEBRA_HOST inline scalar_t
     theta(const ROOT::Math::SVector<scalar_t, N> &v) noexcept {
 
@@ -46,7 +46,7 @@ requires(N >= 3) ALGEBRA_HOST inline scalar_t
       TMath::ATan2(TMath::Sqrt(v[0] * v[0] + v[1] * v[1]), v[2]));
 }
 
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar 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 {
 
@@ -58,13 +58,13 @@ requires(N >= 3) ALGEBRA_HOST inline scalar_t
 /// This method retrieves the norm of a vector, no dimension restriction
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t norm(const ROOT::Math::SVector<scalar_t, N> &v) {
 
   return static_cast<scalar_t>(TMath::Sqrt(ROOT::Math::Dot(v, v)));
 }
 
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar scalar_t, class A, auto N>
 ALGEBRA_HOST inline scalar_t norm(
     const ROOT::Math::VecExpr<A, scalar_t, N> &v) {
 
@@ -75,14 +75,14 @@ ALGEBRA_HOST inline scalar_t norm(
 /// rows >= 3
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar 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>
+template <concepts::scalar 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 {
 
@@ -92,14 +92,14 @@ requires(N >= 3) ALGEBRA_HOST inline scalar_t
 /// This method retrieves the perpendicular magnitude of a vector with rows >= 2
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar 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>
+template <concepts::scalar 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 {
 
@@ -110,7 +110,7 @@ requires(N >= 2) ALGEBRA_HOST inline scalar_t
 /// Get a normalized version of the input vector
 ///
 /// @param v the input vector
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, N> normalize(
     const ROOT::Math::SVector<scalar_t, N> &v) {
 
@@ -120,7 +120,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, N> normalize(
 /// Get a normalized version of the input vector
 ///
 /// @param v the input vector
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar scalar_t, class A, auto N>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, N> normalize(
     const ROOT::Math::VecExpr<A, scalar_t, N> &v) {
 
@@ -133,7 +133,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, N> normalize(
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector<scalar_t, N> &a,
                                  const ROOT::Math::SVector<scalar_t, N> &b) {
 
@@ -146,7 +146,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector<scalar_t, N> &a,
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar scalar_t, class A, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector<scalar_t, N> &a,
                                  const ROOT::Math::VecExpr<A, scalar_t, N> &b) {
 
@@ -159,7 +159,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector<scalar_t, N> &a,
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar scalar_t, class A, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr<A, scalar_t, N> &a,
                                  const ROOT::Math::SVector<scalar_t, N> &b) {
 
@@ -172,7 +172,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr<A, scalar_t, N> &a,
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar scalar_t, class A, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr<A, scalar_t, N> &a,
                                  const ROOT::Math::VecExpr<A, scalar_t, N> &b) {
 
@@ -185,7 +185,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr<A, scalar_t, N> &a,
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar scalar_t, class A, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix<scalar_t, N, 1> &a,
                                  const ROOT::Math::VecExpr<A, scalar_t, N> &b) {
 
@@ -198,7 +198,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix<scalar_t, N, 1> &a,
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, class A, auto N>
+template <concepts::scalar scalar_t, class A, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr<A, scalar_t, N> &a,
                                  const ROOT::Math::SMatrix<scalar_t, N, 1> &b) {
   return dot(b, a);
@@ -210,7 +210,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::VecExpr<A, scalar_t, N> &a,
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix<scalar_t, N, 1> &a,
                                  const ROOT::Math::SVector<scalar_t, N> &b) {
 
@@ -223,7 +223,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SMatrix<scalar_t, N, 1> &a,
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename scalar_t, auto N>
+template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector<scalar_t, N> &a,
                                  const ROOT::Math::SMatrix<scalar_t, N, 1> &b) {
   return dot(b, a);
@@ -235,7 +235,7 @@ ALGEBRA_HOST inline scalar_t dot(const ROOT::Math::SVector<scalar_t, N> &a,
 /// @param b the second input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
     const ROOT::Math::SVector<scalar_t, 3> &a,
     const ROOT::Math::SVector<scalar_t, 3> &b) {
@@ -249,7 +249,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
 /// @param b the second input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t, class A>
+template <concepts::scalar scalar_t, class A>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
     const ROOT::Math::SVector<scalar_t, 3> &a,
     const ROOT::Math::VecExpr<A, scalar_t, 3> &b) {
@@ -263,7 +263,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
 /// @param b the second input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t, class A>
+template <concepts::scalar scalar_t, class A>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
     const ROOT::Math::VecExpr<A, scalar_t, 3> &a,
     const ROOT::Math::SVector<scalar_t, 3> &b) {
@@ -277,7 +277,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
 /// @param b the second input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t, class A>
+template <concepts::scalar scalar_t, class A>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
     const ROOT::Math::VecExpr<A, scalar_t, 3> &a,
     const ROOT::Math::VecExpr<A, scalar_t, 3> &b) {
@@ -291,7 +291,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
 /// @param b the second input matrix<3,1>
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
     const ROOT::Math::SVector<scalar_t, 3> &a,
     const ROOT::Math::SMatrix<scalar_t, 3, 1> &b) {
@@ -305,7 +305,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
 /// @param b the first input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
     const ROOT::Math::SMatrix<scalar_t, 3, 1> &a,
     const ROOT::Math::SVector<scalar_t, 3> &b) {
@@ -319,7 +319,7 @@ ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
 /// @param b the first input matrix<3,1>
 ///
 /// @return a vector (expression) representing the cross product
-template <typename scalar_t>
+template <concepts::scalar scalar_t>
 ALGEBRA_HOST inline ROOT::Math::SVector<scalar_t, 3> cross(
     const ROOT::Math::SMatrix<scalar_t, 3, 1> &a,
     const ROOT::Math::SMatrix<scalar_t, 3, 1> &b) {
diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp
index 1df40a6e..049b51f2 100644
--- a/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp
+++ b/math/vc_aos/include/algebra/math/impl/vc_aos_matrix.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/matrix.hpp"
 
@@ -22,7 +23,7 @@ using storage::zero;
 /// General case: Compute the determinant of a square matrix
 ///
 /// @returns the determinant
-template <std::size_t N, typename value_t,
+template <std::size_t N, concepts::value value_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr value_t determinant(
     const storage::matrix<array_t, value_t, N, N> &) noexcept {
@@ -30,7 +31,7 @@ ALGEBRA_HOST_DEVICE constexpr value_t determinant(
   return value_t(0);
 }
 
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::value value_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr storage::matrix<array_t, value_t, ROW, COL>
 inverse(const storage::matrix<array_t, value_t, ROW, COL> &m) noexcept {
diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp
index f31126e5..170fd41b 100644
--- a/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp
+++ b/math/vc_aos/include/algebra/math/impl/vc_aos_transform3.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/matrix.hpp"
@@ -34,7 +35,8 @@ using algebra::storage::operator-;
 using algebra::storage::operator+;
 
 /// Transform wrapper class to ensure standard API within differnt plugins
-template <template <typename, std::size_t> class array_t, typename scalar_t>
+template <template <typename, std::size_t> class array_t,
+          concepts::value value_t>
 struct transform3 {
 
  private:
@@ -51,24 +53,24 @@ struct transform3 {
   /// @{
 
   /// Scalar type used by the transform
-  using scalar_type = scalar_t;
+  using value_type = value_t;
   /// The type of the matrix elements (scalar for AoS, Vc::Vector for SoA)
-  using value_type =
-      std::conditional_t<Vc::is_simd_vector<array_t<scalar_t, 4>>::value,
-                         scalar_t, Vc::Vector<scalar_t>>;
+  using scalar_type =
+      std::conditional_t<Vc::is_simd_vector<array_t<value_t, 4>>::value,
+                         value_t, Vc::Vector<value_t>>;
 
   template <std::size_t N>
-  using array_type = array_t<value_type, N>;
+  using array_type = array_t<scalar_type, N>;
 
   /// 3-element "vector" type (does not observe translations)
-  using vector3 = storage::vector<3u, value_type, array_t>;
+  using vector3 = storage::vector<3u, scalar_type, array_t>;
   /// Point in 3D space (does observe translations)
   using point3 = vector3;
   /// Point in 2D space
-  using point2 = storage::vector<2u, value_type, array_t>;
+  using point2 = storage::vector<2u, scalar_type, array_t>;
 
   /// 4x4 matrix type (Last row is {0, 0, 0, 1} and can be omitted)
-  using matrix44 = storage::matrix<array_t, value_type, 3u, 4u>;
+  using matrix44 = storage::matrix<array_t, scalar_type, 3u, 4u>;
   using column_t = typename matrix44::vector_type;
 
   /// Function (object) used for accessing a matrix element
@@ -158,11 +160,11 @@ struct transform3 {
 
   /// Matrix access operator
   ALGEBRA_HOST_DEVICE
-  inline const value_type &operator()(std::size_t row, std::size_t col) const {
+  inline const scalar_type &operator()(std::size_t row, std::size_t col) const {
     return _data[col][row];
   }
   ALGEBRA_HOST_DEVICE
-  inline value_type &operator()(std::size_t row, std::size_t col) {
+  inline scalar_type &operator()(std::size_t row, std::size_t col) {
     return _data[col][row];
   }
 
@@ -172,7 +174,7 @@ struct transform3 {
   ///
   /// @return a sacalar determinant - no checking done
   ALGEBRA_HOST_DEVICE
-  constexpr value_type determinant(const matrix44 &m) const {
+  constexpr scalar_type determinant(const matrix44 &m) const {
     return -m[e_z][0] * m[e_y][1] * m[e_x][2] +
            m[e_y][0] * m[e_z][1] * m[e_x][2] +
            m[e_z][0] * m[e_x][1] * m[e_y][2] -
@@ -214,7 +216,7 @@ struct transform3 {
         m[e_t][0] * m[e_x][1] * m[e_y][2] + m[e_x][0] * m[e_t][1] * m[e_y][2] +
         m[e_y][0] * m[e_x][1] * m[e_t][2] - m[e_x][0] * m[e_y][1] * m[e_t][2];
     // i[e_t][3] = 1;
-    const value_type idet{value_type(1.f) / determinant(i)};
+    const scalar_type idet{scalar_type(1.f) / determinant(i)};
 
     i[e_x] = i[e_x] * idet;
     i[e_y] = i[e_y] * idet;
@@ -228,7 +230,7 @@ struct transform3 {
   ///
   /// @param m is the rotation matrix
   /// @param v is the vector to be rotated
-  template <typename vector3_type>
+  template <concepts::vector3D vector3_type>
   ALGEBRA_HOST_DEVICE constexpr auto rotate(const matrix44 &m,
                                             const vector3_type &v) const {
 
@@ -239,7 +241,7 @@ struct transform3 {
   ALGEBRA_HOST_DEVICE
   inline auto rotation() const {
 
-    using matrix_t = storage::matrix<array_t, value_type, 3u, 3u>;
+    using matrix_t = storage::matrix<array_t, scalar_type, 3u, 3u>;
 
     matrix_t submatrix;
     for (unsigned int icol = 0; icol < 3; ++icol) {
@@ -280,7 +282,7 @@ struct transform3 {
   /// @param v is the point to be transformed
   ///
   /// @return a global point
-  template <typename point3_type>
+  template <concepts::point3D point3_type>
   ALGEBRA_HOST_DEVICE inline auto point_to_global(const point3_type &p) const {
 
     return rotate(_data, p) + _data[e_t];
@@ -294,7 +296,7 @@ struct transform3 {
   /// @param v is the point to be transformed
   ///
   /// @return a local point
-  template <typename point3_type>
+  template <concepts::point3D point3_type>
   ALGEBRA_HOST_DEVICE inline auto point_to_local(const point3_type &p) const {
 
     return rotate(_data_inv, p) + _data_inv[e_t];
@@ -308,7 +310,7 @@ struct transform3 {
   /// @param v is the vector to be transformed
   ///
   /// @return a vector in global coordinates
-  template <typename vector3_type>
+  template <concepts::vector3D vector3_type>
   ALGEBRA_HOST_DEVICE inline auto vector_to_global(
       const vector3_type &v) const {
 
@@ -323,7 +325,7 @@ struct transform3 {
   /// @param v is the vector to be transformed
   ///
   /// @return a vector in global coordinates
-  template <typename vector3_type>
+  template <concepts::vector3D vector3_type>
   ALGEBRA_HOST_DEVICE inline auto vector_to_local(const vector3_type &v) const {
 
     return rotate(_data_inv, v);
diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp
index 112e9e3b..f2804c51 100644
--- a/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp
+++ b/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/vector.hpp"
@@ -50,11 +51,10 @@ requires(
 /// This method retrieves the norm of a vector, no dimension restriction
 ///
 /// @param v the input vector
-template <typename vector_t,
-          std::enable_if_t<(Vc::is_simd_vector<vector_t>::value ||
-                            algebra::detail::is_storage_vector_v<vector_t>),
-                           bool> = true>
-ALGEBRA_HOST_DEVICE inline auto norm(const vector_t &v) {
+template <typename vector_t>
+requires(Vc::is_simd_vector<vector_t>::value ||
+         algebra::detail::is_storage_vector_v<vector_t>) ALGEBRA_HOST_DEVICE
+    inline auto norm(const vector_t &v) {
 
   return algebra::math::sqrt(dot(v, v));
 }
@@ -76,11 +76,10 @@ requires(Vc::is_simd_vector<vector_type>::value ||
 /// rows >= 3
 ///
 /// @param v the input vector
-template <typename vector_t,
-          std::enable_if_t<(Vc::is_simd_vector<vector_t>::value ||
-                            algebra::detail::is_storage_vector_v<vector_t>),
-                           bool> = true>
-ALGEBRA_HOST_DEVICE inline auto eta(const vector_t &v) noexcept {
+template <typename vector_t>
+requires(Vc::is_simd_vector<vector_t>::value ||
+         algebra::detail::is_storage_vector_v<vector_t>) ALGEBRA_HOST_DEVICE
+    inline auto eta(const vector_t &v) noexcept {
 
   return algebra::math::atanh(v[2] / norm(v));
 }
diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp
index 5e30e421..a0afc6d8 100644
--- a/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp
+++ b/math/vc_soa/include/algebra/math/impl/vc_soa_matrix.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/impl/vc_soa_vector.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/matrix.hpp"
@@ -20,7 +21,7 @@ using storage::set_zero;
 using storage::transpose;
 using storage::zero;
 
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::value value_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr value_t determinant(
     const storage::matrix<array_t, value_t, ROW, COL> &) noexcept {
@@ -28,7 +29,7 @@ ALGEBRA_HOST_DEVICE constexpr value_t determinant(
   return value_t(0);
 }
 
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::value value_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr storage::matrix<array_t, value_t, ROW, COL>
 inverse(const storage::matrix<array_t, value_t, ROW, COL> &m) noexcept {
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 07fdf65c..364d0781 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
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/vector.hpp"
 
@@ -25,14 +26,14 @@ namespace algebra::vc_soa::math {
 /// This method retrieves phi from a vector, vector base with rows >= 2
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param v the input vector
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           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) {
+requires(N >= 2) ALGEBRA_HOST_DEVICE
+    inline auto phi(const storage::vector<N, Vc::Vector<value_t>, array_t> &v) {
 
   return Vc::atan2(v[1], v[0]);
 }
@@ -40,14 +41,14 @@ requires(N >= 2) ALGEBRA_HOST_DEVICE inline auto phi(
 /// This method retrieves the perpenticular magnitude of a vector with rows >= 2
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param v the input vector
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           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) {
+    const storage::vector<N, Vc::Vector<value_t>, array_t> &v) {
 
   return Vc::sqrt(Vc::fma(v[0], v[0], v[1] * v[1]));
 }
@@ -55,14 +56,14 @@ requires(N >= 2) ALGEBRA_HOST_DEVICE inline auto perp(
 /// This method retrieves theta from a vector, vector base with rows >= 3
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param v the input vector
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           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) {
+    const storage::vector<N, Vc::Vector<value_t>, array_t> &v) {
 
   return Vc::atan2(perp(v), v[2]);
 }
@@ -70,19 +71,19 @@ requires(N >= 3) ALGEBRA_HOST_DEVICE inline auto theta(
 /// Cross product between two input vectors - 3 Dim
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param a the first input vector
 /// @param b the second input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           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) {
+    inline storage::vector<N, Vc::Vector<value_t>, array_t> cross(
+        const storage::vector<N, Vc::Vector<value_t>, array_t> &a,
+        const storage::vector<N, Vc::Vector<value_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])};
@@ -91,18 +92,18 @@ requires(N == 3) ALGEBRA_HOST_DEVICE
 /// Dot product between two input vectors
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param a the first input vector
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline Vc::Vector<scalar_t> dot(
-    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &a,
-    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &b) {
+ALGEBRA_HOST_DEVICE inline Vc::Vector<value_t> dot(
+    const storage::vector<N, Vc::Vector<value_t>, array_t> &a,
+    const storage::vector<N, Vc::Vector<value_t>, array_t> &b) {
 
   auto ret = a[0] * b[0];
 
@@ -116,14 +117,14 @@ ALGEBRA_HOST_DEVICE inline Vc::Vector<scalar_t> dot(
 /// This method retrieves the norm of a vector, no dimension restriction
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param v the input vector
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline auto norm(
-    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
+    const storage::vector<N, Vc::Vector<value_t>, array_t> &v) {
 
   return Vc::sqrt(dot(v, v));
 }
@@ -131,16 +132,16 @@ ALGEBRA_HOST_DEVICE inline auto norm(
 /// Get a normalized version of the input vector
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param v the input vector
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline storage::vector<N, Vc::Vector<scalar_t>, array_t>
-normalize(const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
+ALGEBRA_HOST_DEVICE inline storage::vector<N, Vc::Vector<value_t>, array_t>
+normalize(const storage::vector<N, Vc::Vector<value_t>, array_t> &v) {
 
-  return (Vc::Vector<scalar_t>::One() / norm(v)) * v;
+  return (Vc::Vector<value_t>::One() / norm(v)) * v;
 
   // Less accurate, but faster
   // return Vc::reciprocal(norm(v)) * v;
@@ -153,17 +154,17 @@ normalize(const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
 /// rows >= 3
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param v the input vector
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           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) {
+requires(N >= 3) ALGEBRA_HOST_DEVICE
+    inline auto eta(const storage::vector<N, Vc::Vector<value_t>, array_t> &v) {
 
   // atanh does not exist in Vc
-  auto atanh_func = [](scalar_t e) { return std::atanh(e); };
+  auto atanh_func = [](value_t e) { return std::atanh(e); };
 
   return (v[2] / norm(v)).apply(atanh_func);
 
@@ -177,20 +178,20 @@ requires(N >= 3) ALGEBRA_HOST_DEVICE inline auto eta(
 /// Elementwise sum
 ///
 /// @tparam N dimension of the vector
-/// @tparam scalar_t scalar type
+/// @tparam value_t value type in the simd vectors
 /// @tparam array_t array type that holds the vector elements
 ///
 /// @param v the vector whose elements should be summed
 ///
 /// @return the sum of the elements
-template <std::size_t N, typename scalar_t,
+template <std::size_t N, concepts::value value_t,
           template <typename, std::size_t> class array_t>
-ALGEBRA_HOST_DEVICE inline Vc::Vector<scalar_t> sum(
-    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) {
+ALGEBRA_HOST_DEVICE inline Vc::Vector<value_t> sum(
+    const storage::vector<N, Vc::Vector<value_t>, array_t> &v) {
 
-  Vc::Vector<scalar_t> res{v[0]};
+  Vc::Vector<value_t> res{v[0]};
 
-  for (std::size_t i = 1u; i < Vc::Vector<scalar_t>::size(); ++i) {
+  for (std::size_t i = 1u; i < Vc::Vector<value_t>::size(); ++i) {
     res = res + v[i];
   }
 
diff --git a/storage/array/include/algebra/storage/array.hpp b/storage/array/include/algebra/storage/array.hpp
index ed79c674..ff064bd6 100644
--- a/storage/array/include/algebra/storage/array.hpp
+++ b/storage/array/include/algebra/storage/array.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s)
+#include "algebra/concepts.hpp"
 #include "algebra/storage/impl/cmath_getter.hpp"
 #include "algebra/type_traits.hpp"
 
@@ -25,23 +26,23 @@ using size_type = std::size_t;
 template <typename T, size_type N>
 using storage_type = std::array<T, N>;
 /// Vector type used in the Array storage model
-template <typename T, std::size_t N>
+template <concepts::scalar T, std::size_t N>
 using vector_type = storage_type<T, N>;
 /// Matrix type used in the Array storage model
-template <typename T, size_type ROWS, size_type COLS>
+template <concepts::scalar T, size_type ROWS, size_type COLS>
 using matrix_type = storage_type<storage_type<T, ROWS>, COLS>;
 
 /// 3-element "vector" type, using @c std::array
-template <typename T>
+template <concepts::scalar T>
 using vector3 = storage_type<T, 3>;
 /// Point in 3D space, using @c std::array
-template <typename T>
+template <concepts::scalar T>
 using point3 = vector3<T>;
 /// 2-element "vector" type, using @c std::array
-template <typename T>
+template <concepts::scalar T>
 using vector2 = storage_type<T, 2>;
 /// Point in 2D space, using @c std::array
-template <typename T>
+template <concepts::scalar T>
 using point2 = vector2<T>;
 
 /// Element Getter
diff --git a/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp b/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp
index 4ab385a5..1ee875a1 100644
--- a/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp
+++ b/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 
@@ -22,7 +23,7 @@ namespace algebra::cmath::storage {
 struct element_getter {
 
   /// Operator getting a reference to one element of a non-const matrix
-  template <std::size_t ROWS, std::size_t COLS, typename scalar_t,
+  template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
       array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
@@ -34,7 +35,7 @@ struct element_getter {
   }
 
   /// Operator getting one value of a const matrix
-  template <std::size_t ROWS, std::size_t COLS, typename scalar_t,
+  template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
@@ -46,7 +47,7 @@ struct element_getter {
   }
 
   /// Operator getting a reference to one element of a non-const matrix
-  template <std::size_t N, typename scalar_t,
+  template <std::size_t N, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
       array_t<array_t<scalar_t, N>, 1> &m, std::size_t row) const {
@@ -56,7 +57,7 @@ struct element_getter {
   }
 
   /// Operator getting a reference to one element of a const matrix
-  template <std::size_t N, typename scalar_t,
+  template <std::size_t N, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const array_t<array_t<scalar_t, N>, 1> &m, std::size_t row) const {
@@ -66,7 +67,7 @@ struct element_getter {
   }
 
   /// Operator getting a reference to one element of a non-const matrix
-  template <std::size_t N, typename scalar_t,
+  template <std::size_t N, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE inline scalar_t &operator()(array_t<scalar_t, N> &m,
                                                   std::size_t row) const {
@@ -76,7 +77,7 @@ struct element_getter {
   }
 
   /// Operator getting a reference to one element of a const matrix
-  template <std::size_t N, typename scalar_t,
+  template <std::size_t N, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(const array_t<scalar_t, N> &m,
                                                  std::size_t row) const {
@@ -88,8 +89,8 @@ struct element_getter {
 };  // struct element_getter
 
 /// Function extracting an element from a matrix (const)
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t ROWS, std::size_t COLS>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t ROWS, std::size_t COLS>
 ALGEBRA_HOST_DEVICE inline scalar_t element(
     const array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
     std::size_t col) {
@@ -98,8 +99,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t element(
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t ROWS, std::size_t COLS>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t ROWS, std::size_t COLS>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(
     array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
     std::size_t col) {
@@ -108,8 +109,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t &element(
 }
 
 /// Function extracting an element from a matrix (const)
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t element(
     const array_t<array_t<scalar_t, N>, 1> &m, std::size_t row) {
 
@@ -117,8 +118,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t element(
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(
     array_t<array_t<scalar_t, N>, 1> &m, std::size_t row) {
 
@@ -126,8 +127,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t &element(
 }
 
 /// Function extracting an element from a vector (const)
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t element(const array_t<scalar_t, N> &v,
                                             std::size_t row) {
 
@@ -135,8 +136,8 @@ ALGEBRA_HOST_DEVICE inline scalar_t element(const array_t<scalar_t, N> &v,
 }
 
 /// Function extracting an element from a vector (non-const)
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(array_t<scalar_t, N> &v,
                                              std::size_t row) {
 
@@ -148,7 +149,7 @@ struct block_getter {
 
   /// Operator producing a sub-matrix from a const matrix
   template <std::size_t ROWS, std::size_t COLS, std::size_t oROWS,
-            std::size_t oCOLS, typename scalar_t,
+            std::size_t oCOLS, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE auto operator()(
       const array_t<array_t<scalar_t, oROWS>, oCOLS> &m, std::size_t row,
@@ -166,7 +167,8 @@ struct block_getter {
 
   /// Operator producing a vector out of a const matrix
   template <std::size_t SIZE, std::size_t ROWS, std::size_t COLS,
-            typename scalar_t, template <typename, std::size_t> class array_t>
+            concepts::scalar scalar_t,
+            template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE inline array_t<scalar_t, SIZE> operator()(
       const array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
       std::size_t col) {
@@ -196,7 +198,8 @@ ALGEBRA_HOST_DEVICE decltype(auto) block(const input_matrix_type &m,
 
 /// Function extracting a vector from a matrix
 template <std::size_t SIZE, std::size_t ROWS, std::size_t COLS,
-          typename scalar_t, template <typename, std::size_t> class array_t>
+          concepts::scalar scalar_t,
+          template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline array_t<scalar_t, SIZE> vector(
     const array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
     std::size_t col) {
@@ -207,7 +210,8 @@ ALGEBRA_HOST_DEVICE inline array_t<scalar_t, SIZE> vector(
 /// Sets a matrix of dimension @tparam ROW and @tparam COL as submatrix of
 /// @param m beginning at row @param row and column @param col
 template <std::size_t ROWS, std::size_t COLS, class input_matrix_type,
-          typename scalar_t, template <typename, std::size_t> class array_t>
+          concepts::scalar scalar_t,
+          template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE void set_block(
     input_matrix_type &m, const array_t<array_t<scalar_t, ROWS>, COLS> &b,
     std::size_t row, std::size_t col) {
@@ -220,7 +224,7 @@ ALGEBRA_HOST_DEVICE void set_block(
 
 /// Sets a vector of length @tparam ROW as submatrix of
 /// @param m beginning at row @param row and column @param col
-template <std::size_t ROWS, typename scalar_t,
+template <std::size_t ROWS, concepts::scalar scalar_t,
           template <typename, std::size_t> class vector_t,
           class input_matrix_type>
 ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m,
diff --git a/storage/common/include/algebra/storage/array_operators.hpp b/storage/common/include/algebra/storage/array_operators.hpp
index cffc0cad..f9f8a59d 100644
--- a/storage/common/include/algebra/storage/array_operators.hpp
+++ b/storage/common/include/algebra/storage/array_operators.hpp
@@ -7,6 +7,7 @@
 #pragma once
 
 // Project include(s)
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 
 // System include(s).
@@ -17,10 +18,8 @@ namespace algebra::storage {
 
 /// @name Operators on N-element arrays
 /// @{
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N,
-          std::enable_if_t<!std::is_same_v<array_t<scalar_t, N>, scalar_t>,
-                           bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator*(
     const array_t<scalar_t, N> &a, scalar_t s) noexcept {
 
@@ -32,10 +31,8 @@ ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator*(
   return result;
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N,
-          std::enable_if_t<!std::is_same_v<array_t<scalar_t, N>, scalar_t>,
-                           bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator*(
     scalar_t s, const array_t<scalar_t, N> &a) noexcept {
 
@@ -47,9 +44,8 @@ ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator*(
   return result;
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N,
-          std::enable_if_t<std::is_object_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator-(
     const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
 
@@ -61,9 +57,8 @@ ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator-(
   return result;
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N,
-          std::enable_if_t<std::is_object_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator+(
     const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
 
@@ -75,9 +70,8 @@ ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator+(
   return result;
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N,
-          std::enable_if_t<std::is_object_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator*(
     const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
 
@@ -89,9 +83,8 @@ ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator*(
   return result;
 }
 
-template <template <typename, std::size_t> class array_t, typename scalar_t,
-          std::size_t N,
-          std::enable_if_t<std::is_object_v<scalar_t>, bool> = true>
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE constexpr array_t<scalar_t, N> operator/(
     const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
 
diff --git a/storage/common/include/algebra/storage/matrix.hpp b/storage/common/include/algebra/storage/matrix.hpp
index 74e614c8..a3d0e928 100644
--- a/storage/common/include/algebra/storage/matrix.hpp
+++ b/storage/common/include/algebra/storage/matrix.hpp
@@ -7,6 +7,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/storage/vector.hpp"
 #include "algebra/type_traits.hpp"
 
@@ -17,29 +18,29 @@
 namespace algebra::storage {
 
 /// Generic matrix type that can take vectors as columns
-template <template <typename, std::size_t> class array_t, typename value_t,
-          std::size_t ROW, std::size_t COL>
-struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
+template <template <typename, std::size_t> class array_t,
+          concepts::scalar scalar_t, std::size_t ROW, std::size_t COL>
+struct alignas(alignof(storage::vector<ROW, scalar_t, array_t>)) matrix {
 
   // The matrix consists of column vectors
-  using vector_type = storage::vector<ROW, value_t, array_t>;
+  using vector_type = storage::vector<ROW, scalar_t, array_t>;
   // Value type: Can be simd types
-  using value_type = value_t;
+  using scalar_type = scalar_t;
 
   /// Default constructor
   constexpr matrix() = default;
 
   /// Construct from given column vectors @param v
-  template <typename... vector_t>
+  template <concepts::vector... vector_t>
   ALGEBRA_HOST_DEVICE 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,
+  template <std::size_t R, std::size_t C, concepts::scalar S,
             template <typename, std::size_t> class A>
   ALGEBRA_HOST_DEVICE friend constexpr bool operator==(
-      const matrix<A, V, R, C> &lhs, const matrix<A, V, R, C> &rhs);
+      const matrix<A, S, R, C> &lhs, const matrix<A, S, R, C> &rhs);
 
   /// Subscript operator
   /// @{
@@ -75,14 +76,16 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
   /// @{
   // AoS
   template <std::size_t... I>
-  ALGEBRA_HOST_DEVICE requires(!std::is_scalar_v<value_t>) constexpr bool equal(
-      const matrix &rhs, std::index_sequence<I...>) const {
+  ALGEBRA_HOST_DEVICE requires(
+      !std::is_scalar_v<scalar_t>) constexpr bool equal(const matrix &rhs,
+                                                        std::index_sequence<
+                                                            I...>) const {
     return (... && (m_storage[I] == rhs[I]));
   }
 
   // SoA
   template <std::size_t... I>
-  ALGEBRA_HOST requires(std::is_scalar_v<value_t>) constexpr bool equal(
+  ALGEBRA_HOST requires(std::is_scalar_v<scalar_t>) constexpr bool equal(
       const matrix &rhs, std::index_sequence<I...>) const {
     return (... && ((m_storage[I].get() == rhs[I].get()).isFull()));
   }
@@ -90,41 +93,39 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
 
   /// Arithmetic operators
   /// @{
-  template <std::size_t R, std::size_t C, typename V,
+  template <std::size_t R, std::size_t C, typename S,
             template <typename, std::size_t> class A>
   ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator+(
-      const matrix<A, V, R, C> &lhs, const matrix<A, V, R, C> &rhs) noexcept;
+      const matrix<A, S, R, C> &lhs, const matrix<A, S, R, C> &rhs) noexcept;
 
-  template <std::size_t R, std::size_t C, typename V,
+  template <std::size_t R, std::size_t C, typename S,
             template <typename, std::size_t> class A>
   ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator-(
-      const matrix<A, V, R, C> &lhs, const matrix<A, V, R, C> &rhs) noexcept;
+      const matrix<A, S, R, C> &lhs, const matrix<A, S, R, C> &rhs) noexcept;
 
-  template <std::size_t R, std::size_t C, typename V, typename scalar_t,
+  template <std::size_t R, std::size_t C, typename S1, typename S2,
             template <typename, std::size_t> class A>
-  requires(std::is_scalar_v<scalar_t> ||
-           std::is_same_v<V, scalar_t>) ALGEBRA_HOST_DEVICE
+  requires(std::is_scalar_v<S2> || std::is_same_v<S1, S2>) ALGEBRA_HOST_DEVICE
       friend constexpr decltype(auto)
-      operator*(scalar_t a, const matrix<A, V, R, C> &rhs) noexcept;
+      operator*(S2 a, const matrix<A, S1, R, C> &rhs) noexcept;
 
-  template <std::size_t R, std::size_t C, typename V, typename scalar_t,
+  template <std::size_t R, std::size_t C, typename S1, typename S2,
             template <typename, std::size_t> class A>
-  requires(std::is_scalar_v<scalar_t> ||
-           std::is_same_v<V, scalar_t>) ALGEBRA_HOST_DEVICE
+  requires(std::is_scalar_v<S2> || std::is_same_v<S1, S2>) ALGEBRA_HOST_DEVICE
       friend constexpr decltype(auto)
-      operator*(const matrix<A, V, R, C> &lhs, scalar_t a) noexcept;
+      operator*(const matrix<A, S1, R, C> &lhs, S2 a) noexcept;
 
   /// Matrix-vector multiplication
-  template <std::size_t R, std::size_t C, typename V,
+  template <std::size_t R, std::size_t C, typename S,
             template <typename, std::size_t> class A>
   ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator*(
-      const matrix<A, V, R, C> &lhs, const vector<C, V, A> &v) noexcept;
+      const matrix<A, S, R, C> &lhs, const vector<C, S, A> &v) noexcept;
 
   /// Matrix-matrix multiplication
-  template <std::size_t LR, std::size_t C, std::size_t RC, typename V,
+  template <std::size_t LR, std::size_t C, std::size_t RC, typename S,
             template <typename, std::size_t> class A>
   ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator*(
-      const matrix<A, V, LR, C> &lhs, const matrix<A, V, C, RC> &rhs) noexcept;
+      const matrix<A, S, LR, C> &lhs, const matrix<A, S, C, RC> &rhs) noexcept;
   /// @}
 
   /// Matrix storage
@@ -133,7 +134,7 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
 };  // struct matrix
 
 /// Get a zero-initialized matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE constexpr matrix_t zero() noexcept {
 
   matrix_t m;
@@ -147,39 +148,39 @@ ALGEBRA_HOST_DEVICE constexpr matrix_t zero() noexcept {
 }
 
 /// Set a matrix to zero
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE constexpr void set_zero(matrix_t &m) noexcept {
   m = zero<matrix_t>();
 }
 
 /// Build an identity matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE constexpr matrix_t identity() noexcept {
 
   // Zero initialized
   matrix_t m{zero<matrix_t>()};
 
-  for (std::size_t i = 0u; i < algebra::trait::rank<matrix_t>; ++i) {
-    m[i][i] = typename matrix_t::value_type(1);
+  for (std::size_t i = 0u; i < algebra::traits::rank<matrix_t>; ++i) {
+    m[i][i] = typename matrix_t::scalar_type(1);
   }
 
   return m;
 }
 
 /// Set a matrix to zero
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE constexpr void set_identity(matrix_t &m) noexcept {
   m = identity<matrix_t>();
 }
 
 /// Transpose the matrix @param m
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t, std::size_t... I>
 ALGEBRA_HOST_DEVICE constexpr auto transpose(
-    const matrix<array_t, value_t, ROW, COL> &m,
+    const matrix<array_t, scalar_t, ROW, COL> &m,
     std::index_sequence<I...>) noexcept {
 
-  using matrix_T_t = matrix<array_t, value_t, COL, ROW>;
+  using matrix_T_t = matrix<array_t, scalar_t, COL, ROW>;
   using column_t = typename matrix_T_t::vector_type;
 
   matrix_T_t res_m;
@@ -192,33 +193,31 @@ ALGEBRA_HOST_DEVICE constexpr auto transpose(
 }
 
 /// Build an identity matrix
-template <typename matrix_t>
+template <concepts::matrix matrix_t>
 ALGEBRA_HOST_DEVICE constexpr auto transpose(const matrix_t &m) noexcept {
   return transpose(m, std::make_index_sequence<matrix_t::columns()>());
 }
 
 /// Equality operator between two matrices
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr bool operator==(
-    const matrix<array_t, value_t, ROW, COL> &lhs,
-    const matrix<array_t, value_t, ROW, COL> &rhs) {
+    const matrix<array_t, scalar_t, ROW, COL> &lhs,
+    const matrix<array_t, scalar_t, ROW, COL> &rhs) {
   return lhs.equal(rhs, std::make_index_sequence<COL>());
 }
 
 /// 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 {
+template <concepts::matrix matrix_t, concepts::scalar scalar_t,
+          std::size_t... J>
+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])...};
 }
 
 /// Matrix addition
-template <typename matrix_t, std::size_t... J>
+template <concepts::matrix matrix_t, std::size_t... J>
 ALGEBRA_HOST_DEVICE constexpr matrix_t matrix_add(
     const matrix_t &lhs, const matrix_t &rhs,
     std::index_sequence<J...>) noexcept {
@@ -226,7 +225,7 @@ ALGEBRA_HOST_DEVICE constexpr matrix_t matrix_add(
   return matrix_t{(lhs[J] + rhs[J])...};
 }
 
-template <typename matrix_t, std::size_t... J>
+template <concepts::matrix matrix_t, std::size_t... J>
 ALGEBRA_HOST_DEVICE constexpr decltype(auto) matrix_sub(
     const matrix_t &lhs, const matrix_t &rhs,
     std::index_sequence<J...>) noexcept {
@@ -236,61 +235,63 @@ ALGEBRA_HOST_DEVICE constexpr decltype(auto) matrix_sub(
 
 /// Arithmetic operators
 /// @{
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator+(
-    const matrix<array_t, value_t, ROW, COL> &lhs,
-    const matrix<array_t, value_t, ROW, COL> &rhs) noexcept {
+    const matrix<array_t, scalar_t, ROW, COL> &lhs,
+    const matrix<array_t, scalar_t, ROW, COL> &rhs) noexcept {
 
-  using matrix_t = matrix<array_t, value_t, ROW, COL>;
+  using matrix_t = matrix<array_t, scalar_t, ROW, COL>;
 
   return matrix_add(lhs, rhs, std::make_index_sequence<matrix_t::columns()>());
 }
 
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator-(
-    const matrix<array_t, value_t, ROW, COL> &lhs,
-    const matrix<array_t, value_t, ROW, COL> &rhs) noexcept {
+    const matrix<array_t, scalar_t, ROW, COL> &lhs,
+    const matrix<array_t, scalar_t, ROW, COL> &rhs) noexcept {
 
-  using matrix_t = matrix<array_t, value_t, ROW, COL>;
+  using matrix_t = matrix<array_t, scalar_t, ROW, COL>;
 
   return matrix_sub(lhs, rhs, std::make_index_sequence<matrix_t::columns()>());
 }
 
-template <std::size_t ROW, std::size_t COL, typename value_t, typename scalar_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar1_t,
+          concepts::scalar scalar2_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
+requires(std::is_scalar_v<scalar1_t> ||
+         std::is_same_v<scalar1_t, scalar2_t>) ALGEBRA_HOST_DEVICE
     constexpr decltype(auto)
-    operator*(scalar_t a,
-              const matrix<array_t, value_t, ROW, COL> &rhs) noexcept {
+    operator*(scalar1_t a,
+              const matrix<array_t, scalar2_t, ROW, COL> &rhs) noexcept {
 
-  using matrix_t = matrix<array_t, value_t, ROW, COL>;
+  using matrix_t = matrix<array_t, scalar2_t, ROW, COL>;
 
   return matrix_scalar_mul(a, rhs,
                            std::make_index_sequence<matrix_t::columns()>());
 }
 
-template <std::size_t ROW, std::size_t COL, typename value_t, typename scalar_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar1_t,
+          concepts::scalar scalar2_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
+requires(std::is_scalar_v<scalar1_t> ||
+         std::is_same_v<scalar1_t, scalar2_t>) ALGEBRA_HOST_DEVICE
     constexpr decltype(auto)
-    operator*(const matrix<array_t, value_t, ROW, COL> &lhs,
-              scalar_t a) noexcept {
+    operator*(const matrix<array_t, scalar1_t, ROW, COL> &lhs,
+              scalar2_t a) noexcept {
   return a * lhs;
 }
 
 /// Matrix-vector multiplication
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
-    const matrix<array_t, value_t, ROW, COL> &lhs,
-    const vector<COL, value_t, array_t> &v) noexcept {
+    const matrix<array_t, scalar_t, ROW, COL> &lhs,
+    const vector<COL, scalar_t, array_t> &v) noexcept {
 
   // Init vector
-  vector<ROW, value_t, array_t> res_v{v[0] * lhs[0]};
+  vector<ROW, scalar_t, array_t> res_v{v[0] * lhs[0]};
 
   // Add the rest per column
   for (std::size_t j = 1u; j < COL; ++j) {
@@ -302,13 +303,14 @@ ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
 }
 
 /// Matrix-matrix multiplication
-template <std::size_t LROW, std::size_t COL, std::size_t RCOL, typename value_t,
+template <std::size_t LROW, std::size_t COL, std::size_t RCOL,
+          concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
-    const matrix<array_t, value_t, LROW, COL> &lhs,
-    const matrix<array_t, value_t, COL, RCOL> &rhs) noexcept {
+    const matrix<array_t, scalar_t, LROW, COL> &lhs,
+    const matrix<array_t, scalar_t, COL, RCOL> &rhs) noexcept {
 
-  matrix<array_t, value_t, LROW, RCOL> res_m;
+  matrix<array_t, scalar_t, LROW, RCOL> res_m;
 
   for (std::size_t j = 0u; j < RCOL; ++j) {
     // Init column j
diff --git a/storage/common/include/algebra/storage/matrix_getter.hpp b/storage/common/include/algebra/storage/matrix_getter.hpp
index 522eac97..bdb586b6 100644
--- a/storage/common/include/algebra/storage/matrix_getter.hpp
+++ b/storage/common/include/algebra/storage/matrix_getter.hpp
@@ -7,9 +7,9 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/storage/matrix.hpp"
 #include "algebra/storage/vector.hpp"
-#include "algebra/type_traits.hpp"
 
 // System include(s).
 #include <cassert>
@@ -20,10 +20,10 @@ namespace algebra::storage {
 struct element_getter {
 
   /// Get const access to a matrix element
-  template <template <typename, std::size_t> class array_t, typename value_t,
-            std::size_t ROW, std::size_t COL>
+  template <template <typename, std::size_t> class array_t,
+            concepts::scalar scalar_t, std::size_t ROW, std::size_t COL>
   ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
-      const matrix<array_t, value_t, ROW, COL> &m, std::size_t row,
+      const matrix<array_t, scalar_t, ROW, COL> &m, std::size_t row,
       std::size_t col) const {
 
     // Make sure that the indices are valid.
@@ -35,10 +35,10 @@ struct element_getter {
   }
 
   /// Get non-const access to a matrix element
-  template <template <typename, std::size_t> class array_t, typename value_t,
-            std::size_t ROW, std::size_t COL>
+  template <template <typename, std::size_t> class array_t,
+            concepts::scalar scalar_t, std::size_t ROW, std::size_t COL>
   ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
-      matrix<array_t, value_t, ROW, COL> &m, std::size_t row,
+      matrix<array_t, scalar_t, ROW, COL> &m, std::size_t row,
       std::size_t col) const {
     assert(row < ROW);
     assert(col < COL);
@@ -47,40 +47,40 @@ struct element_getter {
   }
 
   /// Get const access to a matrix element
-  template <template <typename, std::size_t> class array_t, typename value_t,
-            std::size_t ROW>
+  template <template <typename, std::size_t> class array_t,
+            concepts::scalar scalar_t, std::size_t ROW>
   ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
-      const matrix<array_t, value_t, ROW, 1> &m, std::size_t row) const {
+      const matrix<array_t, scalar_t, ROW, 1> &m, std::size_t row) const {
 
     assert(row < ROW);
     return m[0][row];
   }
 
   /// Get non-const access to a matrix element
-  template <template <typename, std::size_t> class array_t, typename value_t,
-            std::size_t ROW>
+  template <template <typename, std::size_t> class array_t,
+            concepts::scalar scalar_t, std::size_t ROW>
   ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
-      matrix<array_t, value_t, ROW, 1> &m, std::size_t row) const {
+      matrix<array_t, scalar_t, ROW, 1> &m, std::size_t row) const {
 
     assert(row < ROW);
     return m[0][row];
   }
 
   /// Get const access to a vector element
-  template <template <typename, std::size_t> class array_t, typename value_t,
-            std::size_t N>
+  template <template <typename, std::size_t> class array_t,
+            concepts::scalar scalar_t, std::size_t N>
   ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
-      const vector<N, value_t, array_t> &v, std::size_t row) const {
+      const vector<N, scalar_t, array_t> &v, std::size_t row) const {
 
     assert(row < N);
     return v[row];
   }
 
   /// Get non-const access to a vector element
-  template <template <typename, std::size_t> class array_t, typename value_t,
-            std::size_t N>
+  template <template <typename, std::size_t> class array_t,
+            concepts::scalar scalar_t, std::size_t N>
   ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
-      vector<N, value_t, array_t> &v, std::size_t row) const {
+      vector<N, scalar_t, array_t> &v, std::size_t row) const {
 
     assert(row < N);
     return v[row];
@@ -89,51 +89,51 @@ struct element_getter {
 };  // struct element_getter
 
 /// Function extracting an element from a matrix (const)
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline decltype(auto) element(
-    const matrix<array_t, value_t, ROW, COL> &m, std::size_t row,
+    const matrix<array_t, scalar_t, ROW, COL> &m, std::size_t row,
     std::size_t col) {
   return element_getter{}(m, row, col);
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <std::size_t ROW, std::size_t COL, typename value_t,
+template <std::size_t ROW, std::size_t COL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline decltype(auto) element(
-    matrix<array_t, value_t, ROW, COL> &m, std::size_t row, std::size_t col) {
+    matrix<array_t, scalar_t, ROW, COL> &m, std::size_t row, std::size_t col) {
   return element_getter{}(m, row, col);
 }
 
 /// Function extracting an element from a 1D matrix (const)
-template <std::size_t ROW, typename value_t,
+template <std::size_t ROW, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline decltype(auto) element(
-    const matrix<array_t, value_t, ROW, 1> &m, std::size_t row) {
+    const matrix<array_t, scalar_t, ROW, 1> &m, std::size_t row) {
   return element_getter{}(m, row);
 }
 
 /// Function extracting an element from a 1D matrix (non-const)
-template <std::size_t ROW, typename value_t,
+template <std::size_t ROW, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline decltype(auto) element(
-    matrix<array_t, value_t, ROW, 1> &m, std::size_t row) {
+    matrix<array_t, scalar_t, ROW, 1> &m, std::size_t row) {
   return element_getter{}(m, row);
 }
 
 /// Function extracting an element from a vector (const)
-template <std::size_t N, typename value_t,
+template <std::size_t N, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline decltype(auto) element(
-    const vector<N, value_t, array_t> &v, std::size_t row) {
+    const vector<N, scalar_t, array_t> &v, std::size_t row) {
   return element_getter{}(v, row);
 }
 
 /// Function extracting an element from a vector (non-const)
-template <std::size_t N, typename value_t,
+template <std::size_t N, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE inline decltype(auto) element(
-    vector<N, value_t, array_t> &v, std::size_t row) {
+    vector<N, scalar_t, array_t> &v, std::size_t row) {
   return element_getter{}(v, row);
 }
 
@@ -142,10 +142,10 @@ struct block_getter {
 
   /// Get a block of a const matrix
   template <std::size_t ROWS, std::size_t COLS, std::size_t mROW,
-            std::size_t mCOL, typename value_t,
+            std::size_t mCOL, concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE constexpr auto operator()(
-      const matrix<array_t, value_t, mROW, mCOL> &m, const std::size_t row,
+      const matrix<array_t, scalar_t, mROW, mCOL> &m, const std::size_t row,
       const std::size_t col) noexcept {
 
     static_assert(ROWS <= mROW);
@@ -153,8 +153,8 @@ struct block_getter {
     assert(row + ROWS <= mROW);
     assert(col + COLS <= mCOL);
 
-    using input_matrix_t = matrix<array_t, value_t, mROW, mCOL>;
-    using matrix_t = matrix<array_t, value_t, ROWS, COLS>;
+    using input_matrix_t = matrix<array_t, scalar_t, mROW, mCOL>;
+    using matrix_t = matrix<array_t, scalar_t, ROWS, COLS>;
 
     matrix_t res_m;
 
@@ -180,9 +180,10 @@ 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>
+            concepts::scalar scalar_t,
+            template <typename, std::size_t> class array_t>
   ALGEBRA_HOST_DEVICE constexpr auto operator()(
-      const matrix<array_t, value_t, ROWS, COLS> &m, const std::size_t row,
+      const matrix<array_t, scalar_t, ROWS, COLS> &m, const std::size_t row,
       const std::size_t col) noexcept {
 
     static_assert(SIZE <= ROWS);
@@ -190,8 +191,8 @@ struct block_getter {
     assert(row + SIZE <= ROWS);
     assert(col <= COLS);
 
-    using input_matrix_t = matrix<array_t, value_t, ROWS, COLS>;
-    using vector_t = vector<SIZE, value_t, array_t>;
+    using input_matrix_t = matrix<array_t, scalar_t, ROWS, COLS>;
+    using vector_t = vector<SIZE, scalar_t, array_t>;
 
     vector_t res_v{};
 
@@ -212,29 +213,29 @@ struct block_getter {
 
 /// Get a block of a const matrix
 template <std::size_t ROWS, std::size_t COLS, std::size_t mROW,
-          std::size_t mCOL, typename value_t,
+          std::size_t mCOL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr auto block(
-    const matrix<array_t, value_t, mROW, mCOL> &m, const std::size_t row,
+    const matrix<array_t, scalar_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);
 }
 
 /// Get a block of a const matrix
 template <std::size_t ROWS, std::size_t COLS, std::size_t mROW,
-          std::size_t mCOL, typename value_t,
+          std::size_t mCOL, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 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,
+    matrix<array_t, scalar_t, mROW, mCOL> &m,
+    const matrix<array_t, scalar_t, ROWS, COLS> &b, const std::size_t row,
     const std::size_t col) noexcept {
   static_assert(ROWS <= mROW);
   static_assert(COLS <= mCOL);
   assert(row + ROWS <= mROW);
   assert(col + COLS <= mCOL);
 
-  using input_matrix_t = matrix<array_t, value_t, mROW, mCOL>;
-  using matrix_t = matrix<array_t, value_t, ROWS, COLS>;
+  using input_matrix_t = matrix<array_t, scalar_t, mROW, mCOL>;
+  using matrix_t = matrix<array_t, scalar_t, ROWS, COLS>;
 
   // Don't access single elements in underlying vectors unless necessary
   if constexpr (matrix_t::storage_rows() == input_matrix_t::storage_rows()) {
@@ -253,10 +254,10 @@ ALGEBRA_HOST_DEVICE constexpr void set_block(
 }
 
 /// Operator setting a block with a vector
-template <class matrix_t, std::size_t ROW, typename value_t,
+template <class matrix_t, std::size_t ROW, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr void set_block(
-    matrix_t &m, const storage::vector<ROW, value_t, array_t> &b,
+    matrix_t &m, const storage::vector<ROW, scalar_t, array_t> &b,
     std::size_t row, std::size_t col) noexcept {
   assert(row < ROW);
   assert(row < matrix_t::rows());
diff --git a/storage/common/include/algebra/storage/vector.hpp b/storage/common/include/algebra/storage/vector.hpp
index 81997721..8fdc8a09 100644
--- a/storage/common/include/algebra/storage/vector.hpp
+++ b/storage/common/include/algebra/storage/vector.hpp
@@ -12,10 +12,12 @@
 #endif
 
 // Project include(s)
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/array_operators.hpp"
 
 // System include(s).
+#include <array>
 #include <cstddef>
 #include <initializer_list>
 #include <type_traits>
@@ -41,65 +43,65 @@ consteval std::size_t nearest_power_of_two(std::size_t min_value,
 
 }  // namespace detail
 
-/// Vector wrapper for AoS vs interleaved SoA data. @c value_t can e.g. be a
+/// Vector wrapper for AoS vs interleaved SoA data. @c scalar_t can e.g. be a
 /// SIMD vector.
-template <std::size_t N, typename value_t,
+template <std::size_t N, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 class alignas(
-    alignof(array_t<value_t, detail::nearest_power_of_two(N, 2u)>)) vector {
+    alignof(array_t<scalar_t, detail::nearest_power_of_two(N, 2u)>)) vector {
 
  public:
   /// @returns the required size of the underlying array storage
   ALGEBRA_HOST_DEVICE
   static consteval std::size_t simd_size() {
-    return std::is_scalar_v<value_t> ? detail::nearest_power_of_two(N, 2u) : N;
+    return concepts::value<scalar_t> ? detail::nearest_power_of_two(N, 2u) : N;
   }
 
   // Value type is a simd vector in SoA and a scalar in AoS
-  using value_type = value_t;
+  using scalar_type = scalar_t;
   /// Underlying data array type
-  using array_type = array_t<value_t, simd_size()>;
+  using array_type = array_t<scalar_t, simd_size()>;
 
   /// Default contructor sets all entries to zero.
   ALGEBRA_HOST_DEVICE
   constexpr vector() {
-    if constexpr (std::is_scalar_v<value_type>) {
+    if constexpr (!concepts::simd_scalar<scalar_type>) {
       zero_fill(std::make_index_sequence<simd_size()>{});
     }
   }
 
   /// Construct from element values @param vals .
   ///
+  /// @{
+
+  /// Construct vector in SoA layout from simd scalars
+  template <concepts::simd_scalar... Scalars>
+  ALGEBRA_HOST_DEVICE requires(sizeof...(Scalars) ==
+                               N) constexpr vector(Scalars &&... scals)
+      : m_data{std::forward<Scalars>(scals)...} {}
+
   /// In order to avoid uninitialized values, which deteriorate the performance
   /// in explicitely vectorized code, the underlying data array is filled with
   /// zeroes if too few arguments are given.
-  template <typename... Values>
-  requires(std::conjunction_v<std::is_convertible<Values, value_type>...> &&
-           sizeof...(Values) <= N && !(N == 3 && simd_size() == 4) &&
-           !(N == 6 && simd_size() == 8)) ALGEBRA_HOST_DEVICE
-      constexpr vector(Values &&... vals)
-      : m_data{std::forward<Values>(vals)...} {
-    // Fill the uninitialized part of the vector register with zero
-    if constexpr ((sizeof...(Values) < simd_size()) &&
-                  std::is_scalar_v<value_type> &&
-                  (!std::conjunction_v<std::is_same<array_type, Values>...>)) {
+  template <concepts::value... Values>
+  ALGEBRA_HOST_DEVICE constexpr vector(Values &&... vals) {
+
+    static_assert(sizeof...(Values) <= N);
+
+    // Fill up last entries, if necessary (explicitly for now)
+    if constexpr ((simd_size() - N) == 1) {
+      m_data = {std::forward<Values>(vals)..., 0.f};
+    } else if constexpr ((simd_size() - N) == 2) {
+      m_data = {std::forward<Values>(vals)..., 0.f, 0.f};
+    } else if constexpr (sizeof...(Values) < simd_size()) {
+      // @TODO: Does this actually work, yet?
       zero_fill(std::make_index_sequence<simd_size() - sizeof...(Values)>{});
+    } else {
+      m_data = {std::forward<Values>(vals)...};
     }
   }
 
-  template <typename... Values>
-  requires(std::conjunction_v<std::is_convertible<Values, value_type>...> &&N ==
-               3 &&
-           simd_size() == 4) ALGEBRA_HOST_DEVICE
-      constexpr vector(Values &&... vals)
-      : m_data{std::forward<Values>(vals)..., 0.f} {}
-
-  template <typename... Values>
-  requires(std::conjunction_v<std::is_convertible<Values, value_type>...> &&N ==
-               6 &&
-           simd_size() == 8) ALGEBRA_HOST_DEVICE
-      constexpr vector(Values &&... vals)
-      : m_data{std::forward<Values>(vals)..., 0.f, 0.f} {}
+  /// @}
 
   /// Construct from existing array storage @param vals
   ALGEBRA_HOST_DEVICE
@@ -118,7 +120,7 @@ class alignas(
   ///
   /// @param list wrap an array of this data.
   ALGEBRA_HOST_DEVICE
-  constexpr vector &operator=(std::initializer_list<value_type> &list) {
+  constexpr vector &operator=(std::initializer_list<scalar_type> &list) {
     m_data = array_type(list);
     return *this;
   }
@@ -146,41 +148,41 @@ class alignas(
   ///
   /// @return Vector expression/return type according to the operation.
   ALGEBRA_HOST_DEVICE
-  constexpr decltype(auto) operator*=(value_type factor) noexcept {
+  constexpr decltype(auto) operator*=(scalar_type factor) noexcept {
     return m_data *= factor;
   }
 
   /// Equality operators
   /// @{
-  template <std::size_t M, typename o_value_t,
+  template <std::size_t M, concepts::scalar o_scalar_t,
             template <typename, std::size_t> class o_array_t,
             template <typename, std::size_t> class p_array_t>
   ALGEBRA_HOST_DEVICE friend constexpr bool operator==(
-      const vector<M, o_value_t, o_array_t> &,
-      const vector<M, o_value_t, p_array_t> &) noexcept;
+      const vector<M, o_scalar_t, o_array_t> &,
+      const vector<M, o_scalar_t, p_array_t> &) noexcept;
 
-  template <std::size_t M, typename o_value_t,
+  template <std::size_t M, concepts::scalar o_scalar_t,
             template <typename, std::size_t> class o_array_t,
             template <typename, std::size_t> class p_array_t, bool>
   ALGEBRA_HOST_DEVICE friend constexpr bool operator==(
-      const vector<M, o_value_t, o_array_t> &,
-      const p_array_t<o_value_t, M> &) noexcept;
+      const vector<M, o_scalar_t, o_array_t> &,
+      const p_array_t<o_scalar_t, M> &) noexcept;
   /// @}
 
   /// Inequality operator
-  template <typename other_type>
+  template <typename other_t>
   ALGEBRA_HOST_DEVICE constexpr bool operator!=(
-      const other_type &rhs) const noexcept {
+      const other_t &rhs) const noexcept {
     return ((*this == rhs) == false);
   }
 
   /// Elementwise comparison. Can result in a vector-of-masks for SoA vectors
-  template <typename other_type>
+  template <typename other_t>
   ALGEBRA_HOST_DEVICE constexpr auto compare(
-      const other_type &rhs) const noexcept {
+      const other_t &rhs) const noexcept {
     using result_t = decltype(m_data[0] == rhs[0]);
 
-    array_t<result_t, N> comp;
+    std::array<result_t, N> comp;
 
     for (unsigned int i{0u}; i < N; ++i) {
       comp[i] = (m_data[i] == rhs[i]);
@@ -190,26 +192,26 @@ class alignas(
   }
 
   /// Holds the data value for every vector element
-  array_t<value_t, simd_size()> m_data;
+  array_t<scalar_t, simd_size()> m_data;
 
  private:
   /// Sets the trailing uninitialized values to zero.
   template <std::size_t... Is>
   ALGEBRA_HOST_DEVICE constexpr void zero_fill(
       std::index_sequence<Is...>) noexcept {
-    ((m_data[simd_size() - sizeof...(Is) + Is] = value_t(0)), ...);
+    ((m_data[simd_size() - sizeof...(Is) + Is] = scalar_t(0)), ...);
   }
 };
 
 /// Friend operators
 /// @{
 
-template <std::size_t N, typename value_t,
+template <std::size_t N, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t,
           template <typename, std::size_t> class o_array_t>
-requires(std::is_scalar_v<value_t>) ALGEBRA_HOST_DEVICE constexpr bool
-operator==(const vector<N, value_t, array_t> &lhs,
-           const o_array_t<value_t, N> &rhs) noexcept {
+requires(std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE constexpr bool
+operator==(const vector<N, scalar_t, array_t> &lhs,
+           const o_array_t<scalar_t, N> &rhs) noexcept {
 
   const auto comp = lhs.compare(rhs);
   bool is_full = false;
@@ -221,12 +223,12 @@ operator==(const vector<N, value_t, array_t> &lhs,
   return is_full;
 }
 
-template <std::size_t N, typename value_t,
+template <std::size_t N, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t,
           template <typename, std::size_t> class o_array_t>
-requires(!std::is_scalar_v<value_t>) ALGEBRA_HOST_DEVICE constexpr bool
-operator==(const vector<N, value_t, array_t> &lhs,
-           const o_array_t<value_t, N> &rhs) noexcept {
+requires(!std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE constexpr bool
+operator==(const vector<N, scalar_t, array_t> &lhs,
+           const o_array_t<scalar_t, N> &rhs) noexcept {
 
   const auto comp = lhs.compare(rhs);
   bool is_full = false;
@@ -238,77 +240,66 @@ operator==(const vector<N, value_t, array_t> &lhs,
 
   return is_full;
 }
-template <std::size_t N, typename value_t,
+
+template <std::size_t N, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t,
           template <typename, std::size_t> class o_array_t>
 ALGEBRA_HOST_DEVICE constexpr bool operator==(
-    const vector<N, value_t, array_t> &lhs,
-    const vector<N, value_t, o_array_t> &rhs) noexcept {
+    const vector<N, scalar_t, array_t> &lhs,
+    const vector<N, scalar_t, o_array_t> &rhs) noexcept {
   return (lhs == rhs.m_data);
 }
 
 /// @}
 
 /// Macro declaring all instances of a specific arithmetic operator
-#define DECLARE_vector_OPERATORS(OP)                                           \
-  template <std::size_t N, typename value_t, typename scalar_t,                \
+#define DECLARE_VECTOR_OPERATORS(OP)                                           \
+  template <std::size_t N, concepts::scalar scalar_t, concepts::value value_t, \
             template <typename, std::size_t> class array_t>                    \
-  requires(std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE constexpr decltype( \
-      auto)                                                                    \
-  operator OP(const vector<N, value_t, array_t> &lhs, scalar_t rhs) noexcept { \
-    return lhs.m_data OP static_cast<value_t>(rhs);                            \
+  ALGEBRA_HOST_DEVICE inline constexpr decltype(auto) operator OP(             \
+      const vector<N, scalar_t, array_t> &lhs, value_t rhs) noexcept {         \
+    return lhs.m_data OP static_cast<scalar_t>(rhs);                           \
   }                                                                            \
-  template <std::size_t N, typename value_t, typename scalar_t,                \
+  template <std::size_t N, concepts::scalar scalar_t, concepts::value value_t, \
             template <typename, std::size_t> class array_t>                    \
-  requires(std::is_scalar_v<scalar_t>) ALGEBRA_HOST_DEVICE inline decltype(    \
-      auto)                                                                    \
-  operator OP(scalar_t lhs, const vector<N, value_t, array_t> &rhs) noexcept { \
-    return static_cast<value_t>(lhs) OP rhs.m_data;                            \
+  ALGEBRA_HOST_DEVICE inline decltype(auto) operator OP(                       \
+      value_t lhs, const vector<N, scalar_t, array_t> &rhs) noexcept {         \
+    return static_cast<scalar_t>(lhs) OP rhs.m_data;                           \
   }                                                                            \
-  template <std::size_t N, typename value_t,                                   \
+  template <std::size_t N, concepts::scalar scalar_t,                          \
             template <typename, std::size_t> class array_t>                    \
-  ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator OP(                    \
-      const vector<N, value_t, array_t> &lhs,                                  \
-      const vector<N, value_t, array_t> &rhs) noexcept {                       \
+  ALGEBRA_HOST_DEVICE inline constexpr decltype(auto) operator OP(             \
+      const vector<N, scalar_t, array_t> &lhs,                                 \
+      const vector<N, scalar_t, array_t> &rhs) noexcept {                      \
     return lhs.m_data OP rhs.m_data;                                           \
   }                                                                            \
-  template <std::size_t N, typename value_t,                                   \
-            template <typename, std::size_t> class array_t,                    \
-            typename other_type>                                               \
-  requires(                                                                    \
-      std::is_object<decltype(                                                 \
-          std::declval<typename vector<N, value_t, array_t>::array_type>()     \
-              OP std::declval<other_type>())>::value &&                        \
-      !std::is_scalar_v<other_type>)                                           \
-      ALGEBRA_HOST_DEVICE constexpr decltype(auto)                             \
-      operator OP(const vector<N, value_t, array_t> &lhs,                      \
-                  const other_type &rhs) noexcept {                            \
+  template <std::size_t N, concepts::scalar scalar_t,                          \
+            template <typename, std::size_t> class array_t, typename other_t>  \
+  requires(concepts::vector<other_t> || concepts::simd_scalar<other_t>)        \
+      ALGEBRA_HOST_DEVICE inline constexpr decltype(auto)                      \
+      operator OP(const vector<N, scalar_t, array_t> &lhs,                     \
+                  const other_t &rhs) noexcept {                               \
     return lhs.m_data OP rhs;                                                  \
   }                                                                            \
-  template <std::size_t N, typename value_t,                                   \
-            template <typename, std::size_t> class array_t,                    \
-            typename other_type>                                               \
-  requires(                                                                    \
-      std::is_object<decltype(                                                 \
-          std::declval<typename vector<N, value_t, array_t>::array_type>()     \
-              OP std::declval<other_type>())>::value &&                        \
-      !std::is_scalar_v<other_type>)                                           \
-      ALGEBRA_HOST_DEVICE constexpr decltype(auto)                             \
-      operator OP(const other_type &lhs,                                       \
-                  const vector<N, value_t, array_t> &rhs) noexcept {           \
+  template <std::size_t N, concepts::scalar scalar_t,                          \
+            template <typename, std::size_t> class array_t, typename other_t>  \
+  requires(concepts::vector<other_t> || concepts::simd_scalar<other_t>)        \
+      ALGEBRA_HOST_DEVICE inline constexpr decltype(auto)                      \
+      operator OP(const other_t &lhs,                                          \
+                  const vector<N, scalar_t, array_t> &rhs) noexcept {          \
     return lhs OP rhs.m_data;                                                  \
   }
 
 // Implement all arithmetic operations on top of @c vector.
 // clang-format off
-DECLARE_vector_OPERATORS(+)
-DECLARE_vector_OPERATORS(-)
-DECLARE_vector_OPERATORS(*)
-DECLARE_vector_OPERATORS(/)
+DECLARE_VECTOR_OPERATORS(+)
+DECLARE_VECTOR_OPERATORS(-)
+DECLARE_VECTOR_OPERATORS(*)
+DECLARE_VECTOR_OPERATORS(/)
 // clang-format on
 
 // Clean up.
-#undef DECLARE_vector_OPERATORS
+#undef DECLARE_VECTOR_OPERATORS
 
 }  // namespace storage
 
@@ -317,9 +308,9 @@ namespace detail {
 template <typename T>
 struct is_storage_vector : public std::false_type {};
 
-template <std::size_t N, typename value_t,
+template <std::size_t N, concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
-struct is_storage_vector<storage::vector<N, value_t, array_t>>
+struct is_storage_vector<storage::vector<N, scalar_t, array_t>>
     : public std::true_type {};
 
 template <typename T>
diff --git a/storage/eigen/include/algebra/storage/eigen.hpp b/storage/eigen/include/algebra/storage/eigen.hpp
index f6b8905c..efc698f1 100644
--- a/storage/eigen/include/algebra/storage/eigen.hpp
+++ b/storage/eigen/include/algebra/storage/eigen.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/storage/impl/eigen_array.hpp"
 #include "algebra/storage/impl/eigen_getter.hpp"
 #include "algebra/type_traits.hpp"
@@ -22,27 +23,27 @@ namespace eigen {
 /// size type for Eigen storage model
 using size_type = int;
 /// Array type used in the Eigen storage model
-template <typename T, size_type N>
+template <concepts::scalar T, size_type N>
 using storage_type = array<T, N>;
 /// Vector type used in the Eigen storage model
-template <typename T, size_type N>
+template <concepts::scalar T, size_type N>
 using vector_type = storage_type<T, N>;
 /// Matrix type used in the Eigen storage model
 /// If the number of rows is 1, make it RowMajor
-template <typename T, size_type ROWS, size_type COLS>
+template <concepts::scalar 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 eigen::vector_type
-template <typename T>
+template <concepts::scalar T>
 using vector3 = storage_type<T, 3>;
 /// Point in 3D space, using @c eigen::vector_type
-template <typename T>
+template <concepts::scalar T>
 using point3 = vector3<T>;
 /// 2-element "vector" type, using @c eigen::vector_type
-template <typename T>
+template <concepts::scalar T>
 using vector2 = storage_type<T, 2>;
 /// Point in 2D space, using @c eigen::vector_type
-template <typename T>
+template <concepts::scalar T>
 using point2 = vector2<T>;
 
 /// Element Getter
@@ -55,7 +56,7 @@ using block_getter = eigen::storage::block_getter;
 ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(eigen)
 
 // Extra specializations needed for some additional eigen types
-namespace trait {
+namespace traits {
 
 /// Index
 /// @{
@@ -172,6 +173,6 @@ struct block_getter<Eigen::MatrixBase<derived_t>> {
 };
 /// @}
 
-}  // namespace trait
+}  // namespace traits
 
 }  // namespace algebra
diff --git a/storage/eigen/include/algebra/storage/impl/eigen_array.hpp b/storage/eigen/include/algebra/storage/impl/eigen_array.hpp
index f36b243e..d86f100d 100644
--- a/storage/eigen/include/algebra/storage/impl/eigen_array.hpp
+++ b/storage/eigen/include/algebra/storage/impl/eigen_array.hpp
@@ -1,12 +1,15 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2020-2023 CERN for the benefit of the ACTS project
+ * (c) 2020-2024 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
 
 #pragma once
 
+// Project include(s).
+#include "algebra/concepts.hpp"
+
 // Eigen include(s).
 #ifdef _MSC_VER
 #pragma warning(push, 0)
@@ -29,7 +32,7 @@
 namespace algebra::eigen {
 
 /// Eigen array type
-template <typename T, int N>
+template <concepts::scalar T, int N>
 class array : public Eigen::Matrix<T, N, 1, 0, N, 1> {
 
  public:
diff --git a/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp b/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp
index 41073078..4c2dec98 100644
--- a/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp
+++ b/storage/eigen/include/algebra/storage/impl/eigen_getter.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 
@@ -35,44 +36,40 @@ 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>
-  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 {
+  template <typename derived_type, concepts::index size_type_1,
+            concepts::index size_type_2>
+  requires std::is_base_of_v<
+      Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
+      Eigen::MatrixBase<derived_type> >
+      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>
-  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 {
+  template <typename derived_type, concepts::index size_type_1,
+            concepts::index size_type_2>
+  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>
-  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 {
+  template <typename derived_type, concepts::index size_type>
+  requires std::is_base_of_v<
+      Eigen::DenseCoeffsBase<derived_type, Eigen::WriteAccessors>,
+      Eigen::MatrixBase<derived_type> >
+      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>
-  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 {
+  template <typename derived_type, concepts::index size_type>
+  ALGEBRA_HOST_DEVICE inline auto operator()(
+      const Eigen::MatrixBase<derived_type> &m, size_type row) const {
 
     return m(static_cast<Eigen::Index>(row));
   }
@@ -121,35 +118,43 @@ requires std::is_base_of_v<
 
 /// 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>
-  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 {
+
+  template <int kROWS, int kCOLS, typename derived_type,
+            concepts::index size_type_1, concepts::index size_type_2>
+  ALGEBRA_HOST_DEVICE decltype(auto) operator()(
+      const Eigen::MatrixBase<derived_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>
-  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 {
+  template <int kROWS, int kCOLS, typename derived_type,
+            concepts::index size_type_1, concepts::index size_type_2>
+  ALGEBRA_HOST_DEVICE decltype(auto) operator()(
+      Eigen::MatrixBase<derived_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>
-  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 {
+  template <int SIZE, typename derived_type, concepts::index size_type_1,
+            concepts::index size_type_2>
+  ALGEBRA_HOST_DEVICE decltype(auto) operator()(
+      Eigen::MatrixBase<derived_type> &m, size_type_1 row,
+      size_type_2 col) const {
+
+    return m.template block<SIZE, 1>(row, col);
+  }
+
+  template <int SIZE, typename derived_type, concepts::index size_type_1,
+            concepts::index size_type_2>
+  ALGEBRA_HOST_DEVICE decltype(auto) operator()(
+      const Eigen::MatrixBase<derived_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
diff --git a/storage/fastor/include/algebra/storage/fastor.hpp b/storage/fastor/include/algebra/storage/fastor.hpp
index 8b2ccb98..d63e47ce 100644
--- a/storage/fastor/include/algebra/storage/fastor.hpp
+++ b/storage/fastor/include/algebra/storage/fastor.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/storage/impl/fastor_getter.hpp"
 #include "algebra/storage/impl/fastor_matrix.hpp"
 #include "algebra/type_traits.hpp"
@@ -22,26 +23,26 @@ namespace fastor {
 /// size type for Fastor storage model
 using size_type = std::size_t;
 /// Array type used in the Fastor storage model
-template <typename T, size_type N>
+template <concepts::scalar 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>
+template <concepts::scalar 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>
+template <concepts::scalar T, size_type ROWS, size_type COLS>
 using matrix_type = algebra::fastor::Matrix<T, ROWS, COLS>;
 
 /// 3-element "vector" type, using @c Fastor::Tensor
-template <typename T>
+template <concepts::scalar T>
 using vector3 = storage_type<T, 3>;
 /// Point in 3D space, using @c Fastor::Tensor
-template <typename T>
+template <concepts::scalar T>
 using point3 = vector3<T>;
 /// 2-element "vector" type, using @c Fastor::Tensor
-template <typename T>
+template <concepts::scalar T>
 using vector2 = storage_type<T, 2>;
 /// Point in 2D space, using @c Fastor::Tensor
-template <typename T>
+template <concepts::scalar T>
 using point2 = vector2<T>;
 
 /// Element Getter
diff --git a/storage/fastor/include/algebra/storage/impl/fastor_getter.hpp b/storage/fastor/include/algebra/storage/impl/fastor_getter.hpp
index df9cdb7f..5c4013c0 100644
--- a/storage/fastor/include/algebra/storage/impl/fastor_getter.hpp
+++ b/storage/fastor/include/algebra/storage/impl/fastor_getter.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/impl/fastor_matrix.hpp"
@@ -30,7 +31,7 @@ namespace algebra::fastor::storage {
 /// Functor used to access elements of Fastor matrices
 struct element_getter {
 
-  template <std::size_t ROWS, std::size_t COLS, typename scalar_t>
+  template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
       Fastor::Tensor<scalar_t, ROWS, COLS> &m, std::size_t row,
       std::size_t col) const {
@@ -40,7 +41,7 @@ struct element_getter {
     return m(row, col);
   }
 
-  template <std::size_t ROWS, std::size_t COLS, typename scalar_t>
+  template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const Fastor::Tensor<scalar_t, ROWS, COLS> &m, std::size_t row,
       std::size_t col) const {
@@ -50,7 +51,7 @@ struct element_getter {
     return m(row, col);
   }
 
-  template <std::size_t N, typename scalar_t>
+  template <std::size_t N, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
       Fastor::Tensor<scalar_t, N, 1> &m, std::size_t row) const {
 
@@ -58,7 +59,7 @@ struct element_getter {
     return m(row);
   }
 
-  template <std::size_t N, typename scalar_t>
+  template <std::size_t N, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const Fastor::Tensor<scalar_t, N, 1> &m, std::size_t row) const {
 
@@ -66,7 +67,7 @@ struct element_getter {
     return m(row);
   }
 
-  template <std::size_t N, typename scalar_t>
+  template <std::size_t N, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
       Fastor::Tensor<scalar_t, N> &m, std::size_t row) const {
 
@@ -74,7 +75,7 @@ struct element_getter {
     return m(row);
   }
 
-  template <std::size_t N, typename scalar_t>
+  template <std::size_t N, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const Fastor::Tensor<scalar_t, N> &m, std::size_t row) const {
 
@@ -85,7 +86,7 @@ struct element_getter {
 };  // element_getter
 
 /// Function extracting an element from a matrix (const)
-template <typename scalar_t, std::size_t ROWS, std::size_t COLS>
+template <concepts::scalar scalar_t, std::size_t ROWS, std::size_t COLS>
 ALGEBRA_HOST_DEVICE inline scalar_t element(
     const Fastor::Tensor<scalar_t, ROWS, COLS> &m, std::size_t row,
     std::size_t col) {
@@ -93,35 +94,35 @@ ALGEBRA_HOST_DEVICE inline scalar_t element(
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <typename scalar_t, std::size_t ROWS, std::size_t COLS>
+template <concepts::scalar scalar_t, std::size_t ROWS, std::size_t COLS>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(
     Fastor::Tensor<scalar_t, ROWS, COLS> &m, std::size_t row, std::size_t col) {
   return element_getter()(m, row, col);
 }
 
 /// Function extracting an element from a matrix (const)
-template <typename scalar_t, std::size_t N>
+template <concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t element(
     const Fastor::Tensor<scalar_t, N, 1> &m, std::size_t row) {
   return element_getter()(m, row);
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <typename scalar_t, std::size_t N>
+template <concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(Fastor::Tensor<scalar_t, N, 1> &m,
                                              std::size_t row) {
   return element_getter()(m, row);
 }
 
 /// Function extracting an element from a vector (const)
-template <typename scalar_t, std::size_t N>
+template <concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t element(
     const Fastor::Tensor<scalar_t, N> &m, std::size_t row) {
   return element_getter()(m, row);
 }
 
 /// Function extracting an element from a vector (non-const)
-template <typename scalar_t, std::size_t N>
+template <concepts::scalar scalar_t, std::size_t N>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(Fastor::Tensor<scalar_t, N> &m,
                                              std::size_t row) {
   return element_getter()(m, row);
@@ -131,7 +132,7 @@ ALGEBRA_HOST_DEVICE inline scalar_t &element(Fastor::Tensor<scalar_t, N> &m,
 struct block_getter {
 
   template <std::size_t ROWS, std::size_t COLS, std::size_t oROWS,
-            std::size_t oCOLS, typename scalar_t>
+            std::size_t oCOLS, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE algebra::fastor::Matrix<scalar_t, ROWS, COLS> operator()(
       const Fastor::Tensor<scalar_t, oROWS, oCOLS> &m, std::size_t row,
       std::size_t col) const {
@@ -141,7 +142,7 @@ struct block_getter {
   }
 
   template <std::size_t SIZE, std::size_t oROWS, std::size_t oCOLS,
-            typename scalar_t>
+            concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE Fastor::Tensor<scalar_t, SIZE> operator()(
       const Fastor::Tensor<scalar_t, oROWS, oCOLS> &m, std::size_t row,
       std::size_t col) const {
@@ -155,7 +156,7 @@ struct block_getter {
 
 /// Operator getting a block of a const matrix
 template <std::size_t ROWS, std::size_t COLS, std::size_t oROWS,
-          std::size_t oCOLS, typename scalar_t>
+          std::size_t oCOLS, concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE decltype(auto) block(
     const Fastor::Tensor<scalar_t, oROWS, oCOLS> &m, std::size_t row,
     std::size_t col) {
@@ -164,7 +165,7 @@ ALGEBRA_HOST_DEVICE decltype(auto) block(
 
 /// Function extracting a slice from the matrix
 template <std::size_t SIZE, std::size_t ROWS, std::size_t COLS,
-          typename scalar_t>
+          concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline decltype(auto) vector(
     const Fastor::Tensor<scalar_t, ROWS, COLS> &m, std::size_t row,
     std::size_t col) {
@@ -173,8 +174,8 @@ ALGEBRA_HOST_DEVICE inline decltype(auto) vector(
 }
 
 /// Operator setting a block with a matrix
-template <std::size_t ROWS, std::size_t COLS, typename scalar_t,
-          class input_matrix_type>
+template <std::size_t ROWS, std::size_t COLS, concepts::scalar scalar_t,
+          concepts::matrix input_matrix_type>
 ALGEBRA_HOST_DEVICE void set_block(
     input_matrix_type &m, const Fastor::Tensor<scalar_t, ROWS, COLS> &b,
     std::size_t row, std::size_t col) {
@@ -183,7 +184,8 @@ ALGEBRA_HOST_DEVICE void set_block(
 }
 
 /// Operator setting a block with a vector
-template <std::size_t ROWS, typename scalar_t, class input_matrix_type>
+template <std::size_t ROWS, concepts::scalar scalar_t,
+          concepts::matrix input_matrix_type>
 ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m,
                                    const Fastor::Tensor<scalar_t, ROWS> &b,
                                    std::size_t row, std::size_t col) {
diff --git a/storage/fastor/include/algebra/storage/impl/fastor_matrix.hpp b/storage/fastor/include/algebra/storage/impl/fastor_matrix.hpp
index 15fa0d5f..a3d6a3e6 100644
--- a/storage/fastor/include/algebra/storage/impl/fastor_matrix.hpp
+++ b/storage/fastor/include/algebra/storage/impl/fastor_matrix.hpp
@@ -1,12 +1,15 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2023 CERN for the benefit of the ACTS project
+ * (c) 2023-2024 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
 
 #pragma once
 
+// Project include(s).
+#include "algebra/concepts.hpp"
+
 // Fastor include(s).
 #ifdef _MSC_VER
 #pragma warning(disable : 4244 4701 4702)
@@ -36,7 +39,7 @@ namespace algebra::fastor {
 /// class to behave the way a `Fastor::Tensor` would. Inheriting from
 /// `Fastor::Tensor` allows this class to reuse all the functions defined in the
 /// parent class (i.e. `Fastor::Tensor`).
-template <typename T, std::size_t M1, std::size_t N>
+template <concepts::scalar T, std::size_t M1, std::size_t N>
 class Matrix : public Fastor::Tensor<T, M1, N> {
 
  public:
@@ -53,8 +56,7 @@ class Matrix : public Fastor::Tensor<T, M1, N> {
   ///
   /// The `static_cast` is there to signal both to the compiler and the reader
   /// that we wish to interpret the `Matrix` object as a `Fastor::Tensor` here.
-  template <typename U, std::size_t M2,
-            std::enable_if_t<std::is_convertible_v<U, T>, bool> = true>
+  template <concepts::scalar U, std::size_t M2>
   inline Matrix<T, M1, M2> operator*(const Matrix<U, N, M2>& other) const {
     return Fastor::matmul(static_cast<Fastor::Tensor<T, M1, N>>(*this),
                           static_cast<Fastor::Tensor<T, N, M2>>(other));
@@ -71,8 +73,7 @@ class Matrix : public Fastor::Tensor<T, M1, N> {
   ///
   /// The `static_cast` is there to signal both to the compiler and the reader
   /// that we wish to interpret the `Matrix` object as a `Fastor::Tensor` here.
-  template <typename U,
-            std::enable_if_t<std::is_convertible_v<U, T>, bool> = true>
+  template <concepts::scalar U>
   inline Fastor::Tensor<T, M1> operator*(
       const Fastor::Tensor<U, N>& other) const {
     return Fastor::matmul(static_cast<Fastor::Tensor<T, M1, N>>(*this),
diff --git a/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp b/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp
index 1d32232e..cd4ded21 100644
--- a/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp
+++ b/storage/smatrix/include/algebra/storage/impl/smatrix_getter.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 
 // ROOT/Smatrix include(s).
@@ -25,7 +26,7 @@ namespace algebra::smatrix::storage {
 /// Functor used to access elements of SMatrix matrices
 struct element_getter {
 
-  template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
+  template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t &operator()(
       ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, unsigned int row,
       unsigned int col) const {
@@ -34,7 +35,7 @@ struct element_getter {
     return m(row, col);
   }
 
-  template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
+  template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, unsigned int row,
       unsigned int col) const {
@@ -43,28 +44,28 @@ struct element_getter {
     return m(row, col);
   }
 
-  template <typename scalar_t, unsigned int N>
+  template <unsigned int N, concepts::scalar scalar_t>
   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);
+    return m(row);
   }
 
-  template <typename scalar_t, unsigned int N>
+  template <unsigned int N, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const ROOT::Math::SMatrix<scalar_t, N, 1> &m, unsigned int row) const {
     assert(row < N);
     return m(row, 0);
   }
 
-  template <typename scalar_t, unsigned int N>
+  template <concepts::scalar 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 <typename scalar_t, unsigned int N>
+  template <unsigned int N, concepts::scalar scalar_t>
   ALGEBRA_HOST_DEVICE inline scalar_t operator()(
       const ROOT::Math::SVector<scalar_t, N> &m, unsigned int row) const {
     assert(row < N);
@@ -73,7 +74,7 @@ struct element_getter {
 };  // element_getter
 
 /// Function extracting an element from a matrix (const)
-template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
+template <concepts::scalar 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) {
@@ -82,7 +83,7 @@ ALGEBRA_HOST_DEVICE inline scalar_t element(
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <typename scalar_t, unsigned int ROWS, unsigned int COLS>
+template <concepts::scalar 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) {
@@ -91,28 +92,28 @@ ALGEBRA_HOST_DEVICE inline scalar_t &element(
 }
 
 /// Function extracting an element from a matrix (const)
-template <typename scalar_t, unsigned int N>
+template <concepts::scalar 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()(m, static_cast<unsigned int>(row));
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <typename scalar_t, unsigned int N>
+template <concepts::scalar 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()(m, static_cast<unsigned int>(row));
 }
 
 /// Function extracting an element from a matrix (const)
-template <typename scalar_t, unsigned int N>
+template <concepts::scalar 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()(m, static_cast<unsigned int>(row));
 }
 
 /// Function extracting an element from a matrix (non-const)
-template <typename scalar_t, unsigned int N>
+template <concepts::scalar scalar_t, unsigned int N>
 ALGEBRA_HOST_DEVICE inline scalar_t &element(
     ROOT::Math::SVector<scalar_t, N> &m, std::size_t row) {
   return element_getter()(m, static_cast<unsigned int>(row));
@@ -122,7 +123,7 @@ ALGEBRA_HOST_DEVICE inline scalar_t &element(
 struct block_getter {
 
   template <unsigned int ROWS, unsigned int COLS, unsigned int oROWS,
-            unsigned int oCOLS, typename scalar_t>
+            unsigned int oCOLS, concepts::scalar 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 {
@@ -131,7 +132,7 @@ struct block_getter {
   }
 
   template <unsigned int SIZE, unsigned int ROWS, unsigned int COLS,
-            typename scalar_t>
+            concepts::scalar 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 {
@@ -142,7 +143,7 @@ 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>
+          unsigned int oCOLS, concepts::scalar 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) {
@@ -153,7 +154,7 @@ ALGEBRA_HOST_DEVICE ROOT::Math::SMatrix<scalar_t, ROWS, COLS> block(
 /// Function extracting a slice from the matrix used by
 /// @c algebra::smatrix::transform3
 template <unsigned int SIZE, unsigned int ROWS, unsigned int COLS,
-          typename scalar_t>
+          concepts::scalar scalar_t>
 ALGEBRA_HOST_DEVICE inline auto vector(
     const ROOT::Math::SMatrix<scalar_t, ROWS, COLS> &m, std::size_t row,
     std::size_t col) {
@@ -163,8 +164,8 @@ ALGEBRA_HOST_DEVICE inline auto vector(
 }
 
 /// Operator setting a block with a matrix
-template <unsigned int ROWS, unsigned int COLS, typename scalar_t,
-          class input_matrix_type>
+template <unsigned int ROWS, unsigned int COLS, concepts::scalar scalar_t,
+          concepts::matrix 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) {
@@ -177,7 +178,8 @@ ALGEBRA_HOST_DEVICE void set_block(
 }
 
 /// Operator setting a block with a vector
-template <unsigned int ROWS, typename scalar_t, class input_matrix_type>
+template <unsigned int ROWS, concepts::scalar scalar_t,
+          concepts::matrix input_matrix_type>
 ALGEBRA_HOST_DEVICE void set_block(input_matrix_type &m,
                                    const ROOT::Math::SVector<scalar_t, ROWS> &b,
                                    unsigned int row, unsigned int col) {
diff --git a/storage/smatrix/include/algebra/storage/smatrix.hpp b/storage/smatrix/include/algebra/storage/smatrix.hpp
index 97169288..492fab83 100644
--- a/storage/smatrix/include/algebra/storage/smatrix.hpp
+++ b/storage/smatrix/include/algebra/storage/smatrix.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s)
+#include "algebra/concepts.hpp"
 #include "algebra/storage/impl/smatrix_getter.hpp"
 #include "algebra/type_traits.hpp"
 
@@ -25,26 +26,26 @@ namespace smatrix {
 /// size type for SMatrix storage model
 using size_type = unsigned int;
 /// Array type used in the SMatrix storage model
-template <typename T, size_type N>
+template <concepts::scalar 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>
+template <concepts::scalar 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>
+template <concepts::scalar T, size_type ROWS, size_type COLS>
 using matrix_type = ROOT::Math::SMatrix<T, ROWS, COLS>;
 
 /// 3-element "vector" type, using @c ROOT::Math::SVector
-template <typename T>
+template <concepts::scalar T>
 using vector3 = storage_type<T, 3>;
 /// Point in 3D space, using @c ROOT::Math::SVector
-template <typename T>
+template <concepts::scalar T>
 using point3 = vector3<T>;
 /// 2-element "vector" type, using @c ROOT::Math::SVector
-template <typename T>
+template <concepts::scalar T>
 using vector2 = storage_type<T, 2>;
 /// Point in 2D space, using @c ROOT::Math::SVector
-template <typename T>
+template <concepts::scalar T>
 using point2 = vector2<T>;
 
 /// Element Getter
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 aa54f67a..3faefe8b 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
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/matrix_getter.hpp"
 
@@ -18,10 +19,11 @@ using algebra::storage::element;
 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 <std::size_t SIZE, std::size_t ROW, std::size_t COL,
+          concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr auto vector(
-    const algebra::storage::matrix<array_t, value_t, ROW, COL> &m,
+    const algebra::storage::matrix<array_t, scalar_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,
                                                                     col);
diff --git a/storage/vc_aos/include/algebra/storage/vc_aos.hpp b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
index f3960d00..c93a0c2e 100644
--- a/storage/vc_aos/include/algebra/storage/vc_aos.hpp
+++ b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/storage/impl/vc_aos_getter.hpp"
 #include "algebra/storage/matrix.hpp"
 #include "algebra/storage/vector.hpp"
@@ -30,38 +31,41 @@ namespace algebra {
 
 namespace vc_aos {
 
-/// size type for Vc storage model
+/// Size type for Vc storage model
 using size_type = std::size_t;
 /// Array type used to store Vc::Vectors
-template <typename T, size_type N>
+template <concepts::value T, size_type N>
 using storage_type = Vc::SimdArray<T, N>;
-/// value type in a linear algebra vector: AoS layout
-template <typename T>
+/// Value type in a linear algebra vector: AoS layout
+template <concepts::value T>
 using value_type = T;
+/// Scalar type in a linear algebra vector: AoS layout
+template <concepts::value T>
+using scalar_type = T;
 /// Vector type used in the Vc AoS storage model
-template <typename T, std::size_t N>
+template <concepts::value T, std::size_t N>
 using vector_type = algebra::storage::vector<N, T, storage_type>;
 /// Matrix type used in the Vc AoS storage model
-template <typename T, size_type ROWS, size_type COLS>
+template <concepts::value T, size_type ROWS, size_type COLS>
 using matrix_type = algebra::storage::matrix<storage_type, T, ROWS, COLS>;
 
 /// 2-element "vector" type, using @c Vc::SimdArray
-template <typename T>
+template <concepts::value T>
 using vector2 = vector_type<T, 2>;
 /// Point in 2D space, using @c Vc::SimdArray
-template <typename T>
+template <concepts::value T>
 using point2 = vector2<T>;
 /// 3-element "vector" type, using @c Vc::SimdArray
-template <typename T>
+template <concepts::value T>
 using vector3 = vector_type<T, 3>;
 /// Point in 3D space, using @c Vc::SimdArray
-template <typename T>
+template <concepts::value T>
 using point3 = vector3<T>;
 /// 6-element "vector" type, using @c Vc::SimdArray
-template <typename T>
+template <concepts::value T>
 using vector6 = vector_type<T, 6>;
 /// 8-element "vector" type, using @c Vc::SimdArray
-template <typename T>
+template <concepts::value T>
 using vector8 = vector_type<T, 8>;
 
 /// Element Getter
@@ -73,4 +77,30 @@ using block_getter = algebra::storage::block_getter;
 
 ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(vc_aos)
 
+namespace traits {
+
+// Vector and storage types are different
+template <typename T, auto N>
+struct dimensions<vc_aos::storage_type<T, N>> {
+
+  using size_type = vc_aos::size_type;
+
+  static constexpr size_type dim{1};
+  static constexpr size_type rows{N};
+  static constexpr size_type columns{1};
+};
+
+// Vector and storage types are different
+template <typename T, auto N>
+struct dimensions<Vc_1::Vector<T, Vc_1::simd_abi::fixed_size<N>>> {
+
+  using size_type = vc_aos::size_type;
+
+  static constexpr size_type dim{1};
+  static constexpr size_type rows{N};
+  static constexpr size_type columns{1};
+};
+
+}  // namespace traits
+
 }  // 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 80418087..d16408a5 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
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/qualifiers.hpp"
 #include "algebra/storage/matrix_getter.hpp"
 
@@ -18,10 +19,11 @@ using algebra::storage::element;
 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 <std::size_t SIZE, std::size_t ROW, std::size_t COL,
+          concepts::scalar scalar_t,
           template <typename, std::size_t> class array_t>
 ALGEBRA_HOST_DEVICE constexpr auto vector(
-    const algebra::storage::matrix<array_t, value_t, ROW, COL> &m,
+    const algebra::storage::matrix<array_t, scalar_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,
                                                                     col);
diff --git a/storage/vc_soa/include/algebra/storage/vc_soa.hpp b/storage/vc_soa/include/algebra/storage/vc_soa.hpp
index c1821d45..4612226d 100644
--- a/storage/vc_soa/include/algebra/storage/vc_soa.hpp
+++ b/storage/vc_soa/include/algebra/storage/vc_soa.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/concepts.hpp"
 #include "algebra/storage/impl/vc_soa_getter.hpp"
 #include "algebra/storage/matrix.hpp"
 #include "algebra/storage/vector.hpp"
@@ -30,39 +31,42 @@ namespace algebra {
 
 namespace vc_soa {
 
-/// size type for Vc storage model
+/// Size type for Vc storage model
 using size_type = std::size_t;
-/// Array type used to store Vc::Vectors
-template <typename T, size_type N>
+/// Array type used to store Vc::Vectors or matrix columns
+template <concepts::scalar T, size_type N>
 using storage_type = std::array<T, N>;
-/// value type in a linear algebra vector: SoA layout
-template <typename T>
-using value_type = Vc::Vector<T>;
+/// Value type in a linear algebra vector: SoA layout
+template <concepts::value T>
+using value_type = T;
+/// Scalar type in a linear algebra vector: SoA layout
+template <concepts::value T>
+using scalar_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, Vc::Vector<T>, std::array>;
+template <concepts::value T, std::size_t N>
+using vector_type = algebra::storage::vector<N, Vc::Vector<T>, storage_type>;
 /// Matrix type used in the Vc SoA storage model
-template <typename T, size_type ROWS, size_type COLS>
+template <concepts::value T, size_type ROWS, size_type COLS>
 using matrix_type =
-    algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>;
+    algebra::storage::matrix<storage_type, Vc::Vector<T>, ROWS, COLS>;
 
 /// 2-element "vector" type, using @c Vc::Vector in every element
-template <typename T>
+template <concepts::value T>
 using vector2 = vector_type<T, 2>;
 /// Point in 2D space, using @c Vc::Vector in every element
-template <typename T>
+template <concepts::value T>
 using point2 = vector2<T>;
 /// 3-element "vector" type, using @c Vc::Vector in every element
-template <typename T>
+template <concepts::value T>
 using vector3 = vector_type<T, 3>;
 /// Point in 3D space, using @c Vc::Vector in every element
-template <typename T>
+template <concepts::value T>
 using point3 = vector3<T>;
 /// 6-element "vector" type, using @c Vc::Vector in every element
-template <typename T>
+template <concepts::value T>
 using vector6 = vector_type<T, 6>;
 /// 8-element "vector" type, using @c Vc::Vector in every element
-template <typename T>
+template <concepts::value T>
 using vector8 = vector_type<T, 8>;
 
 /// Element Getter
@@ -74,16 +78,28 @@ using block_getter = algebra::storage::block_getter;
 
 ALGEBRA_PLUGINS_DEFINE_TYPE_TRAITS(vc_soa)
 
-namespace trait {
+namespace traits {
 
 /// Make sure the simd scalar type is recognized correctly
 /// @{
-template <typename T, std::size_t ROWS, std::size_t COLS>
-struct scalar<algebra::storage::matrix<std::array, Vc::Vector<T>, ROWS, COLS>> {
+template <concepts::value T, std::size_t ROWS, std::size_t COLS>
+struct scalar<
+    algebra::storage::matrix<vc_soa::storage_type, Vc::Vector<T>, ROWS, COLS>> {
   using type = Vc::Vector<T>;
 };
 /// @}
 
-}  // namespace trait
+// Vector and storage types are different
+template <typename T, auto N>
+struct dimensions<vc_soa::storage_type<T, N>> {
+
+  using size_type = vc_soa::size_type;
+
+  static constexpr size_type dim{1};
+  static constexpr size_type rows{N};
+  static constexpr size_type columns{1};
+};
+
+}  // namespace traits
 
 }  // namespace algebra
diff --git a/storage/vecmem/include/algebra/storage/vecmem.hpp b/storage/vecmem/include/algebra/storage/vecmem.hpp
index ec76106a..07801edd 100644
--- a/storage/vecmem/include/algebra/storage/vecmem.hpp
+++ b/storage/vecmem/include/algebra/storage/vecmem.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s)
+#include "algebra/concepts.hpp"
 #include "algebra/storage/impl/cmath_getter.hpp"
 #include "algebra/type_traits.hpp"
 
@@ -27,23 +28,23 @@ using size_type = std::size_t;
 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>
+template <concepts::scalar 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>
+template <concepts::scalar T, std::size_t ROWS, std::size_t COLS>
 using matrix_type = storage_type<storage_type<T, ROWS>, COLS>;
 
 /// 3-element "vector" type, using @c vecmem::static_array
-template <typename T>
+template <concepts::scalar T>
 using vector3 = storage_type<T, 3>;
 /// Point in 3D space, using @c vecmem::static_array
-template <typename T>
+template <concepts::scalar T>
 using point3 = vector3<T>;
 /// 2-element "vector" type, using @c vecmem::static_array
-template <typename T>
+template <concepts::scalar T>
 using vector2 = storage_type<T, 2>;
 /// Point in 2D space, using @c vecmem::static_array
-template <typename T>
+template <concepts::scalar T>
 using point2 = vector2<T>;
 
 /// Element Getter
diff --git a/tests/common/test_base.hpp b/tests/common/test_base.hpp
index 6f433e42..dc1c942f 100644
--- a/tests/common/test_base.hpp
+++ b/tests/common/test_base.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Local include(s).
+#include "algebra/concepts.hpp"
 #include "test_types.hpp"
 
 // System include(s).
@@ -18,10 +19,13 @@ template <class T>
 class test_base {};
 
 /// Test base class, using a @c test_types type argument
-template <typename scalar_t, typename point2_t, typename point3_t,
-          typename vector2_t, typename vector3_t, typename transform3_t,
-          typename size_ty,
-          template <typename, size_ty, size_ty> class matrix_t>
+template <
+    algebra::concepts::scalar scalar_t, algebra::concepts::point2D point2_t,
+    algebra::concepts::point3D point3_t, algebra::concepts::vector2D vector2_t,
+    algebra::concepts::vector3D vector3_t,
+    algebra::concepts::transform3D transform3_t,
+    algebra::concepts::index size_ty,
+    template <typename, size_ty, size_ty> class matrix_t>
 class test_base<test_types<scalar_t, point2_t, point3_t, vector2_t, vector3_t,
                            transform3_t, size_ty, matrix_t> > {
 
diff --git a/tests/common/test_host_basics.hpp b/tests/common/test_host_basics.hpp
index c7a478c5..e69f54dd 100644
--- a/tests/common/test_host_basics.hpp
+++ b/tests/common/test_host_basics.hpp
@@ -75,6 +75,20 @@ TYPED_TEST_P(test_host_basics_vector, local_vectors) {
 // This defines the vector3 test suite
 TYPED_TEST_P(test_host_basics_vector, vector3) {
 
+  // Test concepts
+  static_assert(algebra::concepts::scalar<typename TypeParam::scalar>);
+  static_assert(!algebra::concepts::vector<typename TypeParam::scalar>);
+
+  static_assert(!algebra::concepts::scalar<typename TypeParam::vector3>);
+  static_assert(algebra::concepts::vector<typename TypeParam::vector3>);
+  static_assert(algebra::concepts::vector3D<typename TypeParam::vector3>);
+  static_assert(!algebra::concepts::vector2D<typename TypeParam::vector3>);
+
+  static_assert(!algebra::concepts::scalar<typename TypeParam::vector2>);
+  static_assert(algebra::concepts::vector<typename TypeParam::vector2>);
+  static_assert(algebra::concepts::vector2D<typename TypeParam::vector2>);
+  static_assert(!algebra::concepts::vector3D<typename TypeParam::vector2>);
+
   // Construction
   typename TypeParam::vector3 vA{0.f, 1.f, 2.f};
   ASSERT_EQ(vA[0], 0.f);
@@ -147,24 +161,39 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_2x3) {
   using matrix_2x3_t = typename TypeParam::template matrix<2, 3>;
 
   // Test type traits
-  static_assert(std::is_same_v<algebra::trait::index_t<matrix_2x3_t>,
+  static_assert(std::is_same_v<algebra::traits::index_t<matrix_2x3_t>,
                                typename TypeParam::size_type>);
-  static_assert(std::is_same_v<algebra::trait::value_t<matrix_2x3_t>,
+  static_assert(std::is_same_v<algebra::traits::value_t<matrix_2x3_t>,
                                typename TypeParam::scalar>);
-  static_assert(std::is_same_v<algebra::trait::scalar_t<matrix_2x3_t>,
+  static_assert(std::is_same_v<algebra::traits::scalar_t<matrix_2x3_t>,
                                typename TypeParam::scalar>);
-  static_assert(std::is_same_v<algebra::trait::vector_t<matrix_2x3_t>,
+  static_assert(std::is_same_v<algebra::traits::vector_t<matrix_2x3_t>,
                                typename TypeParam::vector2>);
 
-  static_assert(algebra::trait::rows<matrix_2x3_t> == 2);
-  static_assert(algebra::trait::columns<matrix_2x3_t> == 3);
-  static_assert(algebra::trait::rank<matrix_2x3_t> == 2);
-  static_assert(algebra::trait::size<matrix_2x3_t> == 6);
-  static_assert(!algebra::trait::is_square<matrix_2x3_t>);
+  static_assert(algebra::traits::rows<matrix_2x3_t> == 2);
+  static_assert(algebra::traits::columns<matrix_2x3_t> == 3);
+  static_assert(algebra::traits::rank<matrix_2x3_t> == 2);
+  static_assert(algebra::traits::size<matrix_2x3_t> == 6);
+  static_assert(!algebra::traits::is_square<matrix_2x3_t>);
   static_assert(
-      algebra::trait::is_square<typename TypeParam::template matrix<2, 2>>);
+      algebra::traits::is_square<typename TypeParam::template matrix<2, 2>>);
   static_assert(
-      algebra::trait::is_square<typename TypeParam::template matrix<3, 3>>);
+      algebra::traits::is_square<typename TypeParam::template matrix<3, 3>>);
+
+  // Test concepts
+  static_assert(algebra::concepts::matrix<matrix_2x3_t>);
+  static_assert(!algebra::concepts::scalar<matrix_2x3_t>);
+  static_assert(!algebra::concepts::vector<matrix_2x3_t>);
+  static_assert(!algebra::concepts::square_matrix<matrix_2x3_t>);
+
+  static_assert(
+      algebra::concepts::index<algebra::traits::index_t<matrix_2x3_t>>);
+  static_assert(
+      algebra::concepts::value<algebra::traits::value_t<matrix_2x3_t>>);
+  static_assert(
+      algebra::concepts::scalar<algebra::traits::scalar_t<matrix_2x3_t>>);
+  static_assert(
+      algebra::concepts::vector<algebra::traits::vector_t<matrix_2x3_t>>);
 
   // Test on matrix - vector operations
   typename TypeParam::vector3 vE{1.f, 2.f, 3.f};
@@ -211,22 +240,28 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_6x4) {
   matrix_6x4_t m;
 
   // Test type traits
-  static_assert(std::is_same_v<algebra::trait::index_t<matrix_6x4_t>,
+  static_assert(std::is_same_v<algebra::traits::index_t<matrix_6x4_t>,
                                typename TypeParam::size_type>);
-  static_assert(std::is_same_v<algebra::trait::value_t<matrix_6x4_t>,
+  static_assert(std::is_same_v<algebra::traits::value_t<matrix_6x4_t>,
                                typename TypeParam::scalar>);
-  static_assert(std::is_same_v<algebra::trait::scalar_t<matrix_6x4_t>,
+  static_assert(std::is_same_v<algebra::traits::scalar_t<matrix_6x4_t>,
                                typename TypeParam::scalar>);
 
-  static_assert(algebra::trait::rows<matrix_6x4_t> == 6);
-  static_assert(algebra::trait::columns<matrix_6x4_t> == 4);
-  static_assert(algebra::trait::rank<matrix_6x4_t> == 4);
-  static_assert(algebra::trait::size<matrix_6x4_t> == 24);
-  static_assert(!algebra::trait::is_square<matrix_6x4_t>);
+  static_assert(algebra::traits::rows<matrix_6x4_t> == 6);
+  static_assert(algebra::traits::columns<matrix_6x4_t> == 4);
+  static_assert(algebra::traits::rank<matrix_6x4_t> == 4);
+  static_assert(algebra::traits::size<matrix_6x4_t> == 24);
+  static_assert(!algebra::traits::is_square<matrix_6x4_t>);
   static_assert(
-      algebra::trait::is_square<typename TypeParam::template matrix<4, 4>>);
+      algebra::traits::is_square<typename TypeParam::template matrix<4, 4>>);
   static_assert(
-      algebra::trait::is_square<typename TypeParam::template matrix<6, 6>>);
+      algebra::traits::is_square<typename TypeParam::template matrix<6, 6>>);
+
+  // Test concepts
+  static_assert(algebra::concepts::matrix<matrix_6x4_t>);
+  static_assert(!algebra::concepts::scalar<matrix_6x4_t>);
+  static_assert(!algebra::concepts::vector<matrix_6x4_t>);
+  static_assert(!algebra::concepts::square_matrix<matrix_6x4_t>);
 
   // Fill it.
   for (typename TypeParam::size_type i = 0; i < ROWS; ++i) {
@@ -619,6 +654,8 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_small_mixed) {
 // This defines the transform3 test suite
 TYPED_TEST_P(test_host_basics_transform, transform3) {
 
+  static_assert(algebra::concepts::transform3D<typename TypeParam::transform3>);
+
   // Preparation work
   typename TypeParam::vector3 z =
       algebra::vector::normalize(typename TypeParam::vector3{3.f, 2.f, 1.f});
diff --git a/tests/common/test_types.hpp b/tests/common/test_types.hpp
index 305d25ef..1342f590 100644
--- a/tests/common/test_types.hpp
+++ b/tests/common/test_types.hpp
@@ -7,11 +7,17 @@
 
 #pragma once
 
+// Local include(s).
+#include "algebra/concepts.hpp"
+
 /// Simple struct holding the types that describe a given plugin
-template <typename scalar_t, typename point2_t, typename point3_t,
-          typename vector2_t, typename vector3_t, typename transform3_t,
-          typename size_ty,
-          template <typename, size_ty, size_ty> class matrix_t>
+template <
+    algebra::concepts::scalar scalar_t, algebra::concepts::point2D point2_t,
+    algebra::concepts::point3D point3_t, algebra::concepts::vector2D vector2_t,
+    algebra::concepts::vector3D vector3_t,
+    algebra::concepts::transform3D transform3_t,
+    algebra::concepts::index size_ty,
+    template <typename, size_ty, size_ty> class matrix_t>
 struct test_types {
 
   using scalar = scalar_t;
diff --git a/tests/vc_soa/vc_soa.cpp b/tests/vc_soa/vc_soa.cpp
index b81cddbe..15bff51b 100644
--- a/tests/vc_soa/vc_soa.cpp
+++ b/tests/vc_soa/vc_soa.cpp
@@ -15,23 +15,23 @@
 
 using namespace algebra;
 
-using scalar_t = float;
+using value_t = float;
 
 constexpr float tol{1e-5f};
 
 /// This test the vector functions on an SoA (Vc::Vector) based vector
 TEST(test_vc_host, vc_soa_vector) {
 
-  using vector3_v = vc_soa::vector3<scalar_t>;
+  using vector3_v = vc_soa::vector3<value_t>;
   // Value type is Vc::Vector<float>
-  using value_t = typename vector3_v::value_type;
+  using scalar_t = typename vector3_v::scalar_type;
 
   vector3_v a{1.f, 2.f, 3.f};
   vector3_v b{4.f, 5.f, 6.f};
 
-  EXPECT_TRUE((a[0] == value_t(1.f)).isFull());
-  EXPECT_TRUE((a[1] == value_t(2.f)).isFull());
-  EXPECT_TRUE((a[2] == value_t(3.f)).isFull());
+  EXPECT_TRUE((a[0] == scalar_t(1.f)).isFull());
+  EXPECT_TRUE((a[1] == scalar_t(2.f)).isFull());
+  EXPECT_TRUE((a[2] == scalar_t(3.f)).isFull());
 
   // Masked comparison
   auto m = a.compare(a);
@@ -50,54 +50,54 @@ TEST(test_vc_host, vc_soa_vector) {
 
   // Addition
   auto v_add = a + b;
-  EXPECT_TRUE((v_add[0] == value_t(5.f)).isFull());
-  EXPECT_TRUE((v_add[1] == value_t(7.f)).isFull());
-  EXPECT_TRUE((v_add[2] == value_t(9.f)).isFull());
+  EXPECT_TRUE((v_add[0] == scalar_t(5.f)).isFull());
+  EXPECT_TRUE((v_add[1] == scalar_t(7.f)).isFull());
+  EXPECT_TRUE((v_add[2] == scalar_t(9.f)).isFull());
 
   // Subration
   auto v_sub = a - b;
-  EXPECT_TRUE((v_sub[0] == value_t(-3.f)).isFull());
-  EXPECT_TRUE((v_sub[1] == value_t(-3.f)).isFull());
-  EXPECT_TRUE((v_sub[2] == value_t(-3.f)).isFull());
+  EXPECT_TRUE((v_sub[0] == scalar_t(-3.f)).isFull());
+  EXPECT_TRUE((v_sub[1] == scalar_t(-3.f)).isFull());
+  EXPECT_TRUE((v_sub[2] == scalar_t(-3.f)).isFull());
 
   // Multiplication
   auto v_mul = a * b;
-  EXPECT_TRUE((v_mul[0] == value_t(4.f)).isFull());
-  EXPECT_TRUE((v_mul[1] == value_t(10.f)).isFull());
-  EXPECT_TRUE((v_mul[2] == value_t(18.f)).isFull());
+  EXPECT_TRUE((v_mul[0] == scalar_t(4.f)).isFull());
+  EXPECT_TRUE((v_mul[1] == scalar_t(10.f)).isFull());
+  EXPECT_TRUE((v_mul[2] == scalar_t(18.f)).isFull());
 
   // Division
   auto v_div = a / b;
-  EXPECT_TRUE((v_div[0] == value_t(0.25f)).isFull());
-  EXPECT_TRUE((v_div[1] == value_t(0.4f)).isFull());
-  EXPECT_TRUE((v_div[2] == value_t(0.5f)).isFull());
+  EXPECT_TRUE((v_div[0] == scalar_t(0.25f)).isFull());
+  EXPECT_TRUE((v_div[1] == scalar_t(0.4f)).isFull());
+  EXPECT_TRUE((v_div[2] == scalar_t(0.5f)).isFull());
 
   // Scalar multiplication
   auto v_smul = 2.f * b;
-  EXPECT_TRUE((v_smul[0] == value_t(8.f)).isFull());
-  EXPECT_TRUE((v_smul[1] == value_t(10.f)).isFull());
-  EXPECT_TRUE((v_smul[2] == value_t(12.f)).isFull());
+  EXPECT_TRUE((v_smul[0] == scalar_t(8.f)).isFull());
+  EXPECT_TRUE((v_smul[1] == scalar_t(10.f)).isFull());
+  EXPECT_TRUE((v_smul[2] == scalar_t(12.f)).isFull());
 
   // Expression
   auto v_expr = (b / a) - (2.5f * b) + vector3_v{};
-  EXPECT_TRUE((v_expr[0] == value_t(-6.f)).isFull());
-  EXPECT_TRUE((v_expr[1] == value_t(-10.f)).isFull());
-  EXPECT_TRUE((v_expr[2] == value_t(-13.f)).isFull());
+  EXPECT_TRUE((v_expr[0] == scalar_t(-6.f)).isFull());
+  EXPECT_TRUE((v_expr[1] == scalar_t(-10.f)).isFull());
+  EXPECT_TRUE((v_expr[2] == scalar_t(-13.f)).isFull());
 
   auto d{vector::dot(a, b)};
-  EXPECT_TRUE((d == value_t(32.f)).isFull());
+  EXPECT_TRUE((d == scalar_t(32.f)).isFull());
 
-  value_t norms_a{vector::norm(vector::normalize(a))};
-  value_t norms_b{vector::norm(vector::normalize(b))};
+  scalar_t norms_a{vector::norm(vector::normalize(a))};
+  scalar_t norms_b{vector::norm(vector::normalize(b))};
   for (unsigned int i{0u}; i < norms_a.size(); ++i) {
     EXPECT_NEAR(norms_a[i], 1.f, tol);
     EXPECT_NEAR(norms_b[i], 1.f, tol);
   }
 
   auto cr{vector::cross(a, b)};
-  EXPECT_TRUE((cr[0] == value_t(-3.f)).isFull());
-  EXPECT_TRUE((cr[1] == value_t(6.f)).isFull());
-  EXPECT_TRUE((cr[2] == value_t(-3.f)).isFull());
+  EXPECT_TRUE((cr[0] == scalar_t(-3.f)).isFull());
+  EXPECT_TRUE((cr[1] == scalar_t(6.f)).isFull());
+  EXPECT_TRUE((cr[2] == scalar_t(-3.f)).isFull());
 
   static_assert(std::is_convertible_v<decltype(v_expr), vector3_v>,
                 "expression type not convertible");
@@ -106,7 +106,7 @@ TEST(test_vc_host, vc_soa_vector) {
 /// This test the getter functions on an SoA (Vc::Vector) based vector
 TEST(test_vc_host, vc_soa_getter) {
 
-  using vector3_v = vc_soa::vector3<scalar_t>;
+  using vector3_v = vc_soa::vector3<value_t>;
 
   vector3_v a{1.f, 2.f, 3.f};
 
@@ -114,7 +114,7 @@ TEST(test_vc_host, vc_soa_getter) {
 
   // Phi angle
   auto v_phi = vector::phi(a);
-  EXPECT_NEAR(v_phi[0], static_cast<scalar_t>(std::atan2(2., 1.)), tol);
+  EXPECT_NEAR(v_phi[0], static_cast<value_t>(std::atan2(2., 1.)), tol);
 
   // Perpendicular projection
   auto v_perp = vector::perp(a);
@@ -122,7 +122,7 @@ TEST(test_vc_host, vc_soa_getter) {
 
   // Theta angle
   auto v_theta = vector::theta(a);
-  EXPECT_NEAR(v_theta[0], static_cast<scalar_t>(std::atan2(std::sqrt(5.), 3.)),
+  EXPECT_NEAR(v_theta[0], static_cast<value_t>(std::atan2(std::sqrt(5.), 3.)),
               tol);
 
   // Norm of the vector
@@ -132,32 +132,32 @@ TEST(test_vc_host, vc_soa_getter) {
   // Eta of the vector
   auto v_eta = vector::eta(a);
   EXPECT_NEAR(v_eta[0],
-              static_cast<scalar_t>(std::atanh(1. / std::sqrt(14.) * 3.)), tol);
+              static_cast<value_t>(std::atanh(1. / std::sqrt(14.) * 3.)), tol);
 }
 
 /// This test an SoA (Vc::Vector) based affine transform3
 TEST(test_vc_host, vc_soa_transform3) {
 
-  using vector3 = vc_soa::vector3<scalar_t>;
-  using point3 = vc_soa::point3<scalar_t>;
-  // Value type is Vc::Vector<float>
-  using value_t = typename vector3::value_type;
-  using transform3 = vc_soa::transform3<scalar_t>;
+  using vector3 = vc_soa::vector3<value_t>;
+  using point3 = vc_soa::point3<value_t>;
+  // Scalar type is Vc::Vector<float>
+  using scalar_t = typename vector3::scalar_type;
+  using transform3 = vc_soa::transform3<value_t>;
 
   transform3 idty{};
 
-  EXPECT_TRUE((idty(0, 0) == value_t::One()).isFull());
-  EXPECT_TRUE((idty(1, 0) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(2, 0) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(0, 1) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(1, 1) == value_t::One()).isFull());
-  EXPECT_TRUE((idty(2, 1) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(0, 2) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(1, 2) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(2, 2) == value_t::One()).isFull());
-  EXPECT_TRUE((idty(0, 3) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(1, 3) == value_t::Zero()).isFull());
-  EXPECT_TRUE((idty(2, 3) == value_t::Zero()).isFull());
+  EXPECT_TRUE((idty(0, 0) == scalar_t::One()).isFull());
+  EXPECT_TRUE((idty(1, 0) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(2, 0) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(0, 1) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(1, 1) == scalar_t::One()).isFull());
+  EXPECT_TRUE((idty(2, 1) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(0, 2) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(1, 2) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(2, 2) == scalar_t::One()).isFull());
+  EXPECT_TRUE((idty(0, 3) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(1, 3) == scalar_t::Zero()).isFull());
+  EXPECT_TRUE((idty(2, 3) == scalar_t::Zero()).isFull());
 
   // Preparatioon work
   vector3 z = vector::normalize(vector3{3.f, 2.f, 1.f});
@@ -180,9 +180,9 @@ TEST(test_vc_host, vc_soa_transform3) {
   EXPECT_TRUE((trf2(0, 2) == z[0]).isFull());
   EXPECT_TRUE((trf2(1, 2) == z[1]).isFull());
   EXPECT_TRUE((trf2(2, 2) == z[2]).isFull());
-  EXPECT_TRUE((trf2(0, 3) == 2.f * value_t::One()).isFull());
-  EXPECT_TRUE((trf2(1, 3) == 3.f * value_t::One()).isFull());
-  EXPECT_TRUE((trf2(2, 3) == 4.f * value_t::One()).isFull());
+  EXPECT_TRUE((trf2(0, 3) == 2.f * scalar_t::One()).isFull());
+  EXPECT_TRUE((trf2(1, 3) == 3.f * scalar_t::One()).isFull());
+  EXPECT_TRUE((trf2(2, 3) == 4.f * scalar_t::One()).isFull());
 
   // Check that local origin translates into global translation
   point3 lzero = {0.f, 0.f, 0.f};
@@ -221,45 +221,45 @@ TEST(test_vc_host, vc_soa_transform3) {
 /// This test an SoA (Vc::Vector) based 2x3 matrix
 TEST(test_vc_host, vc_soa_matrix3) {
 
-  using matrix_2x3_t = vc_soa::matrix_type<scalar_t, 2, 3>;
+  using matrix_2x3_t = vc_soa::matrix_type<value_t, 2, 3>;
 
   // Test type traits
   static_assert(
-      std::is_same_v<algebra::trait::index_t<matrix_2x3_t>, std::size_t>);
-  static_assert(std::is_same_v<algebra::trait::value_t<matrix_2x3_t>, float>);
-  static_assert(std::is_same_v<algebra::trait::scalar_t<matrix_2x3_t>,
+      std::is_same_v<algebra::traits::index_t<matrix_2x3_t>, std::size_t>);
+  static_assert(std::is_same_v<algebra::traits::value_t<matrix_2x3_t>, float>);
+  static_assert(std::is_same_v<algebra::traits::scalar_t<matrix_2x3_t>,
                                Vc::Vector<float>>);
-  static_assert(std::is_same_v<algebra::trait::vector_t<matrix_2x3_t>,
-                               vc_soa::vector2<scalar_t>>);
-
-  static_assert(algebra::trait::rows<matrix_2x3_t> == 2);
-  static_assert(algebra::trait::columns<matrix_2x3_t> == 3);
-  static_assert(algebra::trait::rank<matrix_2x3_t> == 2);
-  static_assert(algebra::trait::size<matrix_2x3_t> == 6);
-  static_assert(!algebra::trait::is_square<matrix_2x3_t>);
-  static_assert(algebra::trait::is_square<vc_soa::matrix_type<scalar_t, 2, 2>>);
-  static_assert(algebra::trait::is_square<vc_soa::matrix_type<scalar_t, 3, 3>>);
+  static_assert(std::is_same_v<algebra::traits::vector_t<matrix_2x3_t>,
+                               vc_soa::vector2<value_t>>);
+
+  static_assert(algebra::traits::rows<matrix_2x3_t> == 2);
+  static_assert(algebra::traits::columns<matrix_2x3_t> == 3);
+  static_assert(algebra::traits::rank<matrix_2x3_t> == 2);
+  static_assert(algebra::traits::size<matrix_2x3_t> == 6);
+  static_assert(!algebra::traits::is_square<matrix_2x3_t>);
+  static_assert(algebra::traits::is_square<vc_soa::matrix_type<value_t, 2, 2>>);
+  static_assert(algebra::traits::is_square<vc_soa::matrix_type<value_t, 3, 3>>);
 }
 
 /// This test an SoA (Vc::Vector) based 6x4 matrix
 TEST(test_vc_host, vc_soa_matrix64) {
 
   // Create the matrix.
-  using matrix_6x4_t = vc_soa::matrix_type<scalar_t, 6, 4>;
+  using matrix_6x4_t = vc_soa::matrix_type<value_t, 6, 4>;
   matrix_6x4_t m;
 
   // Test type traits
   static_assert(
-      std::is_same_v<algebra::trait::index_t<matrix_6x4_t>, std::size_t>);
-  static_assert(std::is_same_v<algebra::trait::value_t<matrix_6x4_t>, float>);
-  static_assert(std::is_same_v<algebra::trait::scalar_t<matrix_6x4_t>,
+      std::is_same_v<algebra::traits::index_t<matrix_6x4_t>, std::size_t>);
+  static_assert(std::is_same_v<algebra::traits::value_t<matrix_6x4_t>, float>);
+  static_assert(std::is_same_v<algebra::traits::scalar_t<matrix_6x4_t>,
                                Vc::Vector<float>>);
 
-  static_assert(algebra::trait::rows<matrix_6x4_t> == 6);
-  static_assert(algebra::trait::columns<matrix_6x4_t> == 4);
-  static_assert(algebra::trait::rank<matrix_6x4_t> == 4);
-  static_assert(algebra::trait::size<matrix_6x4_t> == 24);
-  static_assert(!algebra::trait::is_square<matrix_6x4_t>);
-  static_assert(algebra::trait::is_square<vc_soa::matrix_type<scalar_t, 4, 4>>);
-  static_assert(algebra::trait::is_square<vc_soa::matrix_type<scalar_t, 6, 6>>);
+  static_assert(algebra::traits::rows<matrix_6x4_t> == 6);
+  static_assert(algebra::traits::columns<matrix_6x4_t> == 4);
+  static_assert(algebra::traits::rank<matrix_6x4_t> == 4);
+  static_assert(algebra::traits::size<matrix_6x4_t> == 24);
+  static_assert(!algebra::traits::is_square<matrix_6x4_t>);
+  static_assert(algebra::traits::is_square<vc_soa::matrix_type<value_t, 4, 4>>);
+  static_assert(algebra::traits::is_square<vc_soa::matrix_type<value_t, 6, 6>>);
 }