diff --git a/benchmarks/array/array_getter.cpp b/benchmarks/array/array_getter.cpp
index b4d7890f..89f7364c 100644
--- a/benchmarks/array/array_getter.cpp
+++ b/benchmarks/array/array_getter.cpp
@@ -9,7 +9,6 @@
 #include "algebra/array_cmath.hpp"
 #include "benchmark/array/data_generator.hpp"
 #include "benchmark/common/benchmark_getter.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 
 // Benchmark include
 #include <benchmark/benchmark.h>
@@ -48,16 +47,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<phi_f_t>(cfg, "_single");
-  algebra::register_benchmark<phi_d_t>(cfg, "_double");
-  algebra::register_benchmark<theta_f_t>(cfg, "_single");
-  algebra::register_benchmark<theta_d_t>(cfg, "_double");
-  algebra::register_benchmark<perp_f_t>(cfg, "_single");
-  algebra::register_benchmark<perp_d_t>(cfg, "_double");
-  algebra::register_benchmark<norm_f_t>(cfg, "_single");
-  algebra::register_benchmark<norm_d_t>(cfg, "_double");
-  algebra::register_benchmark<eta_f_t>(cfg, "_single");
-  algebra::register_benchmark<eta_d_t>(cfg, "_double");
+  ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/array/array_matrix.cpp b/benchmarks/array/array_matrix.cpp
index 96835923..ac96ba1a 100644
--- a/benchmarks/array/array_matrix.cpp
+++ b/benchmarks/array/array_matrix.cpp
@@ -9,7 +9,6 @@
 #include "algebra/array_cmath.hpp"
 #include "benchmark/array/data_generator.hpp"
 #include "benchmark/common/benchmark_matrix.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 
 // Benchmark include
 #include <benchmark/benchmark.h>
@@ -75,26 +74,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<mat44_add_f_t>(cfg, "_4x4_add_single");
-  algebra::register_benchmark<mat44_add_d_t>(cfg, "_4x4_add_double");
-  algebra::register_benchmark<mat66_add_f_t>(cfg, "_6x6_add_single");
-  algebra::register_benchmark<mat66_add_d_t>(cfg, "_6x6_add_double");
-  algebra::register_benchmark<mat88_add_f_t>(cfg, "_8x8_add_single");
-  algebra::register_benchmark<mat88_add_d_t>(cfg, "_8x8_add_double");
-
-  algebra::register_benchmark<mat44_mul_f_t>(cfg, "_4x4_mul_single");
-  algebra::register_benchmark<mat44_mul_d_t>(cfg, "_4x4_mul_double");
-  algebra::register_benchmark<mat66_mul_f_t>(cfg, "_6x6_mul_single");
-  algebra::register_benchmark<mat66_mul_d_t>(cfg, "_6x6_mul_double");
-  algebra::register_benchmark<mat88_mul_f_t>(cfg, "_8x8_mul_single");
-  algebra::register_benchmark<mat88_mul_d_t>(cfg, "_8x8_mul_double");
-
-  algebra::register_benchmark<mat44_vec_f_t>(cfg, "_4x4_vec_single");
-  algebra::register_benchmark<mat44_vec_d_t>(cfg, "_4x4_vec_double");
-  algebra::register_benchmark<mat66_vec_f_t>(cfg, "_6x6_vec_single");
-  algebra::register_benchmark<mat66_vec_d_t>(cfg, "_6x6_vec_double");
-  algebra::register_benchmark<mat88_vec_f_t>(cfg, "_8x8_vec_single");
-  algebra::register_benchmark<mat88_vec_d_t>(cfg, "_8x8_vec_double");
+  ALGEBRA_PLUGINS_REGISTER_MATRIX_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/array/array_vector.cpp b/benchmarks/array/array_vector.cpp
index e680b481..5d21f9b0 100644
--- a/benchmarks/array/array_vector.cpp
+++ b/benchmarks/array/array_vector.cpp
@@ -9,7 +9,6 @@
 #include "algebra/array_cmath.hpp"
 #include "benchmark/array/data_generator.hpp"
 #include "benchmark/common/benchmark_vector.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 
 // Benchmark include
 #include <benchmark/benchmark.h>
@@ -50,16 +49,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<add_f_t>(cfg, "_single");
-  algebra::register_benchmark<add_d_t>(cfg, "_double");
-  algebra::register_benchmark<sub_f_t>(cfg, "_single");
-  algebra::register_benchmark<sub_d_t>(cfg, "_double");
-  algebra::register_benchmark<dot_f_t>(cfg, "_single");
-  algebra::register_benchmark<dot_d_t>(cfg, "_double");
-  algebra::register_benchmark<cross_f_t>(cfg, "_single");
-  algebra::register_benchmark<cross_d_t>(cfg, "_double");
-  algebra::register_benchmark<normlz_f_t>(cfg, "_single");
-  algebra::register_benchmark<normlz_d_t>(cfg, "_double");
+  ALGEBRA_PLUGINS_REGISTER_VECTOR_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/common/include/benchmark/common/benchmark_getter.hpp b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
index 916dd5b9..bb7bd205 100644
--- a/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
+++ b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
@@ -9,6 +9,7 @@
 
 // Project include(s)
 #include "benchmark_vector.hpp"
+#include "register_benchmark.hpp"
 
 namespace algebra::bench_op {
 
@@ -29,4 +30,17 @@ ALGEBRA_PLUGINS_BENCH_GETTER(perp)
 ALGEBRA_PLUGINS_BENCH_GETTER(norm)
 ALGEBRA_PLUGINS_BENCH_GETTER(eta)
 
+// Macro for registering all getter benchmarks
+#define ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(CFG)        \
+  algebra::register_benchmark<phi_f_t>(CFG, "_single");   \
+  algebra::register_benchmark<phi_d_t>(CFG, "_double");   \
+  algebra::register_benchmark<theta_f_t>(CFG, "_single"); \
+  algebra::register_benchmark<theta_d_t>(CFG, "_double"); \
+  algebra::register_benchmark<perp_f_t>(CFG, "_single");  \
+  algebra::register_benchmark<perp_d_t>(CFG, "_double");  \
+  algebra::register_benchmark<norm_f_t>(CFG, "_single");  \
+  algebra::register_benchmark<norm_d_t>(CFG, "_double");  \
+  algebra::register_benchmark<eta_f_t>(CFG, "_single");   \
+  algebra::register_benchmark<eta_d_t>(CFG, "_double");
+
 }  // namespace algebra::bench_op
diff --git a/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp b/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp
index cb9d7f50..3da3333b 100644
--- a/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp
+++ b/benchmarks/common/include/benchmark/common/benchmark_matrix.hpp
@@ -9,6 +9,7 @@
 
 // Project include(s)
 #include "benchmark_base.hpp"
+#include "register_benchmark.hpp"
 
 // System include(s)
 #include <string>
@@ -183,4 +184,27 @@ struct mul {
 
 }  // namespace bench_op
 
+// Macro for registering all vector benchmarks
+#define ALGEBRA_PLUGINS_REGISTER_MATRIX_BENCH(CFG)                    \
+  algebra::register_benchmark<mat44_add_f_t>(CFG, "_4x4_add_single"); \
+  algebra::register_benchmark<mat44_add_d_t>(CFG, "_4x4_add_double"); \
+  algebra::register_benchmark<mat66_add_f_t>(CFG, "_6x6_add_single"); \
+  algebra::register_benchmark<mat66_add_d_t>(CFG, "_6x6_add_double"); \
+  algebra::register_benchmark<mat88_add_f_t>(CFG, "_8x8_add_single"); \
+  algebra::register_benchmark<mat88_add_d_t>(CFG, "_8x8_add_double"); \
+                                                                      \
+  algebra::register_benchmark<mat44_mul_f_t>(CFG, "_4x4_mul_single"); \
+  algebra::register_benchmark<mat44_mul_d_t>(CFG, "_4x4_mul_double"); \
+  algebra::register_benchmark<mat66_mul_f_t>(CFG, "_6x6_mul_single"); \
+  algebra::register_benchmark<mat66_mul_d_t>(CFG, "_6x6_mul_double"); \
+  algebra::register_benchmark<mat88_mul_f_t>(CFG, "_8x8_mul_single"); \
+  algebra::register_benchmark<mat88_mul_d_t>(CFG, "_8x8_mul_double"); \
+                                                                      \
+  algebra::register_benchmark<mat44_vec_f_t>(CFG, "_4x4_vec_single"); \
+  algebra::register_benchmark<mat44_vec_d_t>(CFG, "_4x4_vec_double"); \
+  algebra::register_benchmark<mat66_vec_f_t>(CFG, "_6x6_vec_single"); \
+  algebra::register_benchmark<mat66_vec_d_t>(CFG, "_6x6_vec_double"); \
+  algebra::register_benchmark<mat88_vec_f_t>(CFG, "_8x8_vec_single"); \
+  algebra::register_benchmark<mat88_vec_d_t>(CFG, "_8x8_vec_double");
+
 }  // namespace algebra
diff --git a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
index cb128012..93de4c5c 100644
--- a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
+++ b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
@@ -9,6 +9,7 @@
 
 // Project include(s)
 #include "benchmark_base.hpp"
+#include "register_benchmark.hpp"
 
 // System include(s)
 #include <chrono>
@@ -161,4 +162,17 @@ struct normalize {
 
 }  // namespace bench_op
 
+// Macro for registering all vector benchmarks
+#define ALGEBRA_PLUGINS_REGISTER_VECTOR_BENCH(CFG)         \
+  algebra::register_benchmark<add_f_t>(CFG, "_single");    \
+  algebra::register_benchmark<add_d_t>(CFG, "_double");    \
+  algebra::register_benchmark<sub_f_t>(CFG, "_single");    \
+  algebra::register_benchmark<sub_d_t>(CFG, "_double");    \
+  algebra::register_benchmark<dot_f_t>(CFG, "_single");    \
+  algebra::register_benchmark<dot_d_t>(CFG, "_double");    \
+  algebra::register_benchmark<cross_f_t>(CFG, "_single");  \
+  algebra::register_benchmark<cross_d_t>(CFG, "_double");  \
+  algebra::register_benchmark<normlz_f_t>(CFG, "_single"); \
+  algebra::register_benchmark<normlz_d_t>(CFG, "_double");
+
 }  // namespace algebra
diff --git a/benchmarks/eigen/eigen_getter.cpp b/benchmarks/eigen/eigen_getter.cpp
index 2761a5b0..aa0fd427 100644
--- a/benchmarks/eigen/eigen_getter.cpp
+++ b/benchmarks/eigen/eigen_getter.cpp
@@ -8,7 +8,6 @@
 // Project include(s)
 #include "algebra/eigen_eigen.hpp"
 #include "benchmark/common/benchmark_getter.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 #include "benchmark/eigen/data_generator.hpp"
 
 // Benchmark include
@@ -48,16 +47,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<phi_f_t>(cfg, "_single");
-  algebra::register_benchmark<phi_d_t>(cfg, "_double");
-  algebra::register_benchmark<theta_f_t>(cfg, "_single");
-  algebra::register_benchmark<theta_d_t>(cfg, "_double");
-  algebra::register_benchmark<perp_f_t>(cfg, "_single");
-  algebra::register_benchmark<perp_d_t>(cfg, "_double");
-  algebra::register_benchmark<norm_f_t>(cfg, "_single");
-  algebra::register_benchmark<norm_d_t>(cfg, "_double");
-  algebra::register_benchmark<eta_f_t>(cfg, "_single");
-  algebra::register_benchmark<eta_d_t>(cfg, "_double");
+  ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/eigen/eigen_matrix.cpp b/benchmarks/eigen/eigen_matrix.cpp
index f66eb43c..b7d30ada 100644
--- a/benchmarks/eigen/eigen_matrix.cpp
+++ b/benchmarks/eigen/eigen_matrix.cpp
@@ -8,7 +8,6 @@
 // Project include(s)
 #include "algebra/eigen_eigen.hpp"
 #include "benchmark/common/benchmark_matrix.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 #include "benchmark/eigen/data_generator.hpp"
 
 // Benchmark include
@@ -75,26 +74,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<mat44_add_f_t>(cfg, "_4x4_add_single");
-  algebra::register_benchmark<mat44_add_d_t>(cfg, "_4x4_add_double");
-  algebra::register_benchmark<mat66_add_f_t>(cfg, "_6x6_add_single");
-  algebra::register_benchmark<mat66_add_d_t>(cfg, "_6x6_add_double");
-  algebra::register_benchmark<mat88_add_f_t>(cfg, "_8x8_add_single");
-  algebra::register_benchmark<mat88_add_d_t>(cfg, "_8x8_add_double");
-
-  algebra::register_benchmark<mat44_mul_f_t>(cfg, "_4x4_mul_single");
-  algebra::register_benchmark<mat44_mul_d_t>(cfg, "_4x4_mul_double");
-  algebra::register_benchmark<mat66_mul_f_t>(cfg, "_6x6_mul_single");
-  algebra::register_benchmark<mat66_mul_d_t>(cfg, "_6x6_mul_double");
-  algebra::register_benchmark<mat88_mul_f_t>(cfg, "_8x8_mul_single");
-  algebra::register_benchmark<mat88_mul_d_t>(cfg, "_8x8_mul_double");
-
-  algebra::register_benchmark<mat44_vec_f_t>(cfg, "_4x4_vec_single");
-  algebra::register_benchmark<mat44_vec_d_t>(cfg, "_4x4_vec_double");
-  algebra::register_benchmark<mat66_vec_f_t>(cfg, "_6x6_vec_single");
-  algebra::register_benchmark<mat66_vec_d_t>(cfg, "_6x6_vec_double");
-  algebra::register_benchmark<mat88_vec_f_t>(cfg, "_8x8_vec_single");
-  algebra::register_benchmark<mat88_vec_d_t>(cfg, "_8x8_vec_double");
+  ALGEBRA_PLUGINS_REGISTER_MATRIX_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/eigen/eigen_vector.cpp b/benchmarks/eigen/eigen_vector.cpp
index e36fe0e7..de791f86 100644
--- a/benchmarks/eigen/eigen_vector.cpp
+++ b/benchmarks/eigen/eigen_vector.cpp
@@ -8,7 +8,6 @@
 // Project include(s)
 #include "algebra/eigen_eigen.hpp"
 #include "benchmark/common/benchmark_vector.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 #include "benchmark/eigen/data_generator.hpp"
 
 // Benchmark include
@@ -50,16 +49,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<add_f_t>(cfg, "_single");
-  algebra::register_benchmark<add_d_t>(cfg, "_double");
-  algebra::register_benchmark<sub_f_t>(cfg, "_single");
-  algebra::register_benchmark<sub_d_t>(cfg, "_double");
-  algebra::register_benchmark<dot_f_t>(cfg, "_single");
-  algebra::register_benchmark<dot_d_t>(cfg, "_double");
-  algebra::register_benchmark<cross_f_t>(cfg, "_single");
-  algebra::register_benchmark<cross_d_t>(cfg, "_double");
-  algebra::register_benchmark<normlz_f_t>(cfg, "_single");
-  algebra::register_benchmark<normlz_d_t>(cfg, "_double");
+  ALGEBRA_PLUGINS_REGISTER_VECTOR_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/vc_aos/vc_aos_getter.cpp b/benchmarks/vc_aos/vc_aos_getter.cpp
index f27bca7d..8f391020 100644
--- a/benchmarks/vc_aos/vc_aos_getter.cpp
+++ b/benchmarks/vc_aos/vc_aos_getter.cpp
@@ -8,7 +8,6 @@
 // Project include(s)
 #include "algebra/vc_aos.hpp"
 #include "benchmark/common/benchmark_getter.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 #include "benchmark/vc_aos/data_generator.hpp"
 
 // Benchmark include
@@ -48,16 +47,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<phi_f_t>(cfg, "_single");
-  algebra::register_benchmark<phi_d_t>(cfg, "_double");
-  algebra::register_benchmark<theta_f_t>(cfg, "_single");
-  algebra::register_benchmark<theta_d_t>(cfg, "_double");
-  algebra::register_benchmark<perp_f_t>(cfg, "_single");
-  algebra::register_benchmark<perp_d_t>(cfg, "_double");
-  algebra::register_benchmark<norm_f_t>(cfg, "_single");
-  algebra::register_benchmark<norm_d_t>(cfg, "_double");
-  algebra::register_benchmark<eta_f_t>(cfg, "_single");
-  algebra::register_benchmark<eta_d_t>(cfg, "_double");
+  ALGEBRA_PLUGINS_REGISTER_GETTER_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/vc_aos/vc_aos_matrix.cpp b/benchmarks/vc_aos/vc_aos_matrix.cpp
index d728bb1b..95be41d7 100644
--- a/benchmarks/vc_aos/vc_aos_matrix.cpp
+++ b/benchmarks/vc_aos/vc_aos_matrix.cpp
@@ -8,7 +8,6 @@
 // Project include(s)
 #include "algebra/vc_aos.hpp"
 #include "benchmark/common/benchmark_matrix.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 #include "benchmark/vc_aos/data_generator.hpp"
 
 // Benchmark include
@@ -75,26 +74,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<mat44_add_f_t>(cfg, "_4x4_add_single");
-  algebra::register_benchmark<mat44_add_d_t>(cfg, "_4x4_add_double");
-  algebra::register_benchmark<mat66_add_f_t>(cfg, "_6x6_add_single");
-  algebra::register_benchmark<mat66_add_d_t>(cfg, "_6x6_add_double");
-  algebra::register_benchmark<mat88_add_f_t>(cfg, "_8x8_add_single");
-  algebra::register_benchmark<mat88_add_d_t>(cfg, "_8x8_add_double");
-
-  algebra::register_benchmark<mat44_mul_f_t>(cfg, "_4x4_mul_single");
-  algebra::register_benchmark<mat44_mul_d_t>(cfg, "_4x4_mul_double");
-  algebra::register_benchmark<mat66_mul_f_t>(cfg, "_6x6_mul_single");
-  algebra::register_benchmark<mat66_mul_d_t>(cfg, "_6x6_mul_double");
-  algebra::register_benchmark<mat88_mul_f_t>(cfg, "_8x8_mul_single");
-  algebra::register_benchmark<mat88_mul_d_t>(cfg, "_8x8_mul_double");
-
-  algebra::register_benchmark<mat44_vec_f_t>(cfg, "_4x4_vec_single");
-  algebra::register_benchmark<mat44_vec_d_t>(cfg, "_4x4_vec_double");
-  algebra::register_benchmark<mat66_vec_f_t>(cfg, "_6x6_vec_single");
-  algebra::register_benchmark<mat66_vec_d_t>(cfg, "_6x6_vec_double");
-  algebra::register_benchmark<mat88_vec_f_t>(cfg, "_8x8_vec_single");
-  algebra::register_benchmark<mat88_vec_d_t>(cfg, "_8x8_vec_double");
+  ALGEBRA_PLUGINS_REGISTER_MATRIX_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/benchmarks/vc_aos/vc_aos_vector.cpp b/benchmarks/vc_aos/vc_aos_vector.cpp
index 39f0e028..fbe8261c 100644
--- a/benchmarks/vc_aos/vc_aos_vector.cpp
+++ b/benchmarks/vc_aos/vc_aos_vector.cpp
@@ -8,7 +8,6 @@
 // Project include(s)
 #include "algebra/vc_aos.hpp"
 #include "benchmark/common/benchmark_vector.hpp"
-#include "benchmark/common/register_benchmark.hpp"
 #include "benchmark/vc_aos/data_generator.hpp"
 
 // Benchmark include
@@ -48,16 +47,7 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<add_f_t>(cfg, "_single");
-  algebra::register_benchmark<add_d_t>(cfg, "_double");
-  algebra::register_benchmark<sub_f_t>(cfg, "_single");
-  algebra::register_benchmark<sub_d_t>(cfg, "_double");
-  algebra::register_benchmark<dot_f_t>(cfg, "_single");
-  algebra::register_benchmark<dot_d_t>(cfg, "_double");
-  algebra::register_benchmark<cross_f_t>(cfg, "_single");
-  algebra::register_benchmark<cross_d_t>(cfg, "_double");
-  algebra::register_benchmark<normlz_f_t>(cfg, "_single");
-  algebra::register_benchmark<normlz_d_t>(cfg, "_double");
+  ALGEBRA_PLUGINS_REGISTER_VECTOR_BENCH(cfg)
 
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
diff --git a/extern/vc/CMakeLists.txt b/extern/vc/CMakeLists.txt
index adb221bb..050ecd72 100644
--- a/extern/vc/CMakeLists.txt
+++ b/extern/vc/CMakeLists.txt
@@ -52,7 +52,7 @@ target_include_directories( Vc
 # Disable a warning for GCC which the Vc build turns on explicitly for that
 # compiler, just to then go and trigger it... :-/
 if( "${CMAKE_CXX_COMPILER_ID}" MATCHES "GNU" )
-   target_compile_options( Vc PRIVATE "-Wno-old-style-cast -Wdeprecated-enum-enum-conversion -Wdeprecated-declarations" )
+   target_compile_options( Vc PRIVATE "-Wno-old-style-cast )
 endif()
 
 # Set up an alias for the Vc target, with the same name that it will have
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 17e617d8..d61521d2 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
@@ -104,22 +104,6 @@ struct actor {
       const matrix_type<ROWS, COLS> &m) {
     return storage::transpose(m);
   }
-
-  /// Get determinant using a specific algorithm @tparam determinant_t
-  /*template <std::size_t N, typename determinant_t = determinant_t>
-  ALGEBRA_HOST_DEVICE static constexpr constexpr value_type determinant(
-      const matrix_type<N, N> &m) {
-
-    return determinant_t{}(m);
-  }
-
-  /// Get inverse using a specific algorithm @tparam inverse_t
-  template <std::size_t N, typename inverse_t = inverse_t>
-  ALGEBRA_HOST_DEVICE static constexpr constexpr matrix_type<N, N> inverse(
-      const matrix_type<N, N> &m) {
-
-    return inverse_t{}(m);
-  }*/
 };
 
 }  // namespace matrix
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 db6be7ad..277053f7 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
@@ -36,6 +36,7 @@ using algebra::storage::operator+;
 template <template <typename, std::size_t> class array_t, typename scalar_t>
 struct transform3 {
 
+ private:
   // How to address matrix columns
   enum elem : std::size_t {
     e_x = 0u,
@@ -44,6 +45,7 @@ struct transform3 {
     e_t = 3u,
   };
 
+ public:
   /// @name Type definitions for the struct
   /// @{
 
diff --git a/storage/common/include/algebra/storage/matrix.hpp b/storage/common/include/algebra/storage/matrix.hpp
index 9d4dd7a4..592e8676 100644
--- a/storage/common/include/algebra/storage/matrix.hpp
+++ b/storage/common/include/algebra/storage/matrix.hpp
@@ -35,21 +35,24 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
 
   /// Construct from given column vectors @param v
   template <typename... vector_t>
-  requires(sizeof...(vector_t) == COL) explicit matrix(vector_t &&... v)
-      : m_storage{std::forward<vector_t>(v)...} {}
+  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 <typename, std::size_t> class A>
-  friend constexpr bool operator==(const matrix<A, V, R, C> &lhs,
-                                   const matrix<A, V, R, C> &rhs);
+  ALGEBRA_HOST_DEVICE friend constexpr bool operator==(
+      const matrix<A, V, R, C> &lhs, const matrix<A, V, R, C> &rhs);
 
   /// Subscript operator
   /// @{
+  ALGEBRA_HOST_DEVICE
   constexpr const vector_type &operator[](const std::size_t i) const {
     assert(i < COL);
     return m_storage[i];
   }
+  ALGEBRA_HOST_DEVICE
   constexpr vector_type &operator[](const std::size_t i) {
     assert(i < COL);
     return m_storage[i];
@@ -57,15 +60,18 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
   /// @}
 
   /// @returns the number of rows
+  ALGEBRA_HOST_DEVICE
   static consteval std::size_t rows() { return ROW; }
 
   /// @returns the number of rows of the underlying vector storage
   /// @note can be different from the matrix rows due to padding
+  ALGEBRA_HOST_DEVICE
   static consteval std::size_t storage_rows() {
     return vector_type::simd_size();
   }
 
   /// @returns the number of columns
+  ALGEBRA_HOST_DEVICE
   static consteval std::size_t columns() { return COL; }
 
  private:
@@ -73,15 +79,17 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
   /// @{
   // AoS
   template <std::size_t... I>
-  requires(!std::is_scalar_v<value_t>) constexpr bool equal(
-      const matrix &rhs, std::index_sequence<I...>) const {
+  ALGEBRA_HOST_DEVICE
+    requires(!std::is_scalar_v<value_t>)
+  constexpr bool equal(const matrix &rhs, std::index_sequence<I...>) const {
     return (... && (m_storage[I] == rhs[I]));
   }
 
   // SoA
   template <std::size_t... I>
-  requires(std::is_scalar_v<value_t>) constexpr bool equal(
-      const matrix &rhs, std::index_sequence<I...>) const {
+  ALGEBRA_HOST
+    requires(std::is_scalar_v<value_t>)
+  constexpr bool equal(const matrix &rhs, std::index_sequence<I...>) const {
     return (... && ((m_storage[I].get() == rhs[I].get()).isFull()));
   }
   /// @}
@@ -100,16 +108,15 @@ struct alignas(alignof(storage::vector<ROW, value_t, array_t>)) matrix {
 
   template <std::size_t R, std::size_t C, typename V, typename scalar_t,
             template <typename, std::size_t> class A>
-  requires(std::is_scalar_v<scalar_t> ||
-           std::is_same_v<V, scalar_t>) ALGEBRA_HOST_DEVICE
-      friend constexpr decltype(auto)
-      operator*(scalar_t a, const matrix<A, V, R, C> &rhs) noexcept;
+    requires(std::is_scalar_v<scalar_t> || std::is_same_v<V, scalar_t>)
+  ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator*(
+      scalar_t a, const matrix<A, V, R, C> &rhs) noexcept;
 
   template <std::size_t R, std::size_t C, typename V, typename scalar_t,
             template <typename, std::size_t> class A>
-  requires(std::is_scalar_v<scalar_t> ||
-           std::is_same_v<V, scalar_t>) friend constexpr decltype(auto)
-  operator*(const matrix<A, V, R, C> &lhs, scalar_t a) noexcept;
+    requires(std::is_scalar_v<scalar_t> || std::is_same_v<V, scalar_t>)
+  ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator*(
+      const matrix<A, V, R, C> &lhs, scalar_t a) noexcept;
 
   /// Matrix-vector multiplication
   template <std::size_t R, std::size_t C, typename V,
@@ -135,7 +142,7 @@ 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>
-  ALGEBRA_HOST inline decltype(auto) operator()(
+  ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
       const matrix<array_t, value_t, ROW, COL> &m, std::size_t row,
       std::size_t col) const {
 
@@ -150,7 +157,7 @@ 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>
-  ALGEBRA_HOST inline decltype(auto) operator()(
+  ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
       matrix<array_t, value_t, ROW, COL> &m, std::size_t row,
       std::size_t col) const {
 
@@ -164,9 +171,8 @@ struct element_getter {
 
   /// Get const access to a matrix44 element from a flat array
   template <template <typename, std::size_t> class array_t, typename value_t>
-  ALGEBRA_HOST inline decltype(auto) operator()(const array_t<value_t, 16> &m,
-                                                unsigned int row,
-                                                unsigned int col) const {
+  ALGEBRA_HOST_DEVICE inline decltype(auto) operator()(
+      const array_t<value_t, 16> &m, unsigned int row, unsigned int col) const {
     // Make sure that the indices are valid.
     assert(row < 4);
     assert(col < 4);
@@ -263,11 +269,10 @@ ALGEBRA_HOST_DEVICE constexpr auto block(
       }
       return;
     }
-  } else {
-    for (std::size_t j = col; j < col + COLS; ++j) {
-      for (std::size_t i = row; i < row + ROWS; ++i) {
-        res_m[j - col][i - row] = m[j][i];
-      }
+  }
+  for (std::size_t j = col; j < col + COLS; ++j) {
+    for (std::size_t i = row; i < row + ROWS; ++i) {
+      res_m[j - col][i - row] = m[j][i];
     }
   }
 
@@ -330,18 +335,18 @@ ALGEBRA_HOST_DEVICE constexpr void set_block(
 /// Equality operator between two matrices
 template <std::size_t ROW, std::size_t COL, typename value_t,
           template <typename, std::size_t> class array_t>
-constexpr bool operator==(const matrix<array_t, value_t, ROW, COL> &lhs,
-                          const matrix<array_t, value_t, ROW, COL> &rhs) {
+ALGEBRA_HOST_DEVICE constexpr bool operator==(
+    const matrix<array_t, value_t, ROW, COL> &lhs,
+    const matrix<array_t, value_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 {
+  requires(std::is_scalar_v<scalar_t> ||
+           std::is_same_v<typename matrix_t::value_type, scalar_t>)
+ALGEBRA_HOST_DEVICE constexpr matrix_t matrix_scalar_mul(
+    scalar_t a, const matrix_t &rhs, std::index_sequence<J...>) noexcept {
 
   return matrix_t{(a * rhs[J])...};
 }
@@ -389,11 +394,9 @@ ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator-(
 
 template <std::size_t ROW, std::size_t COL, typename value_t, typename scalar_t,
           template <typename, std::size_t> class array_t>
-requires(std::is_scalar_v<scalar_t> ||
-         std::is_same_v<value_t, scalar_t>) ALGEBRA_HOST_DEVICE
-    constexpr decltype(auto)
-    operator*(scalar_t a,
-              const matrix<array_t, value_t, ROW, COL> &rhs) noexcept {
+  requires(std::is_scalar_v<scalar_t> || std::is_same_v<value_t, scalar_t>)
+ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
+    scalar_t a, const matrix<array_t, value_t, ROW, COL> &rhs) noexcept {
 
   using matrix_t = matrix<array_t, value_t, ROW, COL>;
 
@@ -403,9 +406,9 @@ requires(std::is_scalar_v<scalar_t> ||
 
 template <std::size_t ROW, std::size_t COL, typename value_t, typename scalar_t,
           template <typename, std::size_t> class array_t>
-requires(std::is_scalar_v<scalar_t> ||
-         std::is_same_v<value_t, scalar_t>) constexpr decltype(auto)
-operator*(const matrix<array_t, value_t, ROW, COL> &lhs, scalar_t a) noexcept {
+  requires(std::is_scalar_v<scalar_t> || std::is_same_v<value_t, scalar_t>)
+ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
+    const matrix<array_t, value_t, ROW, COL> &lhs, scalar_t a) noexcept {
   return a * lhs;
 }
 
diff --git a/storage/common/include/algebra/storage/vector.hpp b/storage/common/include/algebra/storage/vector.hpp
index 39118063..742bb84e 100644
--- a/storage/common/include/algebra/storage/vector.hpp
+++ b/storage/common/include/algebra/storage/vector.hpp
@@ -29,6 +29,7 @@ namespace detail {
 /// Make sure the vector/matrix dimension aligns with simd sizes
 /// @see
 /// https://gitlab.in2p3.fr/CodeursIntensifs/Fast5x5/-/blob/master/fast5x5.hpp?ref_type=heads
+ALGEBRA_HOST_DEVICE
 consteval std::size_t nearest_power_of_two(std::size_t min_value,
                                            std::size_t current_value) {
   // Computes the nearest power of two relative to `min_value` starting from the
@@ -49,6 +50,7 @@ class alignas(
 
  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;
   }
@@ -59,6 +61,7 @@ class alignas(
   using array_type = array_t<value_t, simd_size()>;
 
   /// Default contructor sets all entries to zero.
+  ALGEBRA_HOST_DEVICE
   constexpr vector() {
     if constexpr (std::is_scalar_v<value_type>) {
       zero_fill(std::make_index_sequence<simd_size()>{});
@@ -70,24 +73,10 @@ class alignas(
   /// 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,
-            std::enable_if_t<std::conjunction_v<
-                                 std::is_convertible<Values, value_type>...> &&
-                                 sizeof...(Values) <= N,
-                             bool> = true>
-  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>...>)) {
-      zero_fill(
-          std::make_index_sequence<simd_size() - sizeof...(Values)>{});
-    }
-  }*/
   template <typename... Values>
-  requires(std::conjunction_v<std::is_convertible<Values, value_type>...> &&
-           sizeof...(Values) <= N &&
-           simd_size() != 4) constexpr vector(Values &&... vals)
+    requires(std::conjunction_v<std::is_convertible<Values, value_type>...> &&
+             sizeof...(Values) <= N && simd_size() != 4)
+  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()) &&
@@ -98,23 +87,25 @@ class alignas(
   }
 
   template <typename... Values>
-  requires(std::conjunction_v<std::is_convertible<Values, value_type>...> &&N ==
-               3 &&
-           simd_size() == 4) constexpr vector(Values &&... vals)
+    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) constexpr vector(Values &&... vals)
+    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 .
+  /// Construct from existing array storage @param vals
+  ALGEBRA_HOST_DEVICE
   constexpr vector(const array_type &vals) : m_data{vals} {}
 
   /// Assignment operator from wrapped data.
   ///
   /// @param lhs wrap a copy of this data.
+  ALGEBRA_HOST_DEVICE
   constexpr const vector &operator=(const array_type &lhs) {
     m_data = lhs;
     return *this;
@@ -123,6 +114,7 @@ class alignas(
   /// Assignment operator from @c std::initializer_list .
   ///
   /// @param list wrap an array of this data.
+  ALGEBRA_HOST_DEVICE
   constexpr vector &operator=(std::initializer_list<value_type> &list) {
     m_data = array_type(list);
     return *this;
@@ -130,21 +122,27 @@ class alignas(
 
   /// Conversion operator from wrapper to underlying data array.
   /// @{
+  ALGEBRA_HOST_DEVICE
   constexpr operator array_type &() { return m_data; }
+  ALGEBRA_HOST_DEVICE
   constexpr operator const array_type &() const { return m_data; }
   /// @}
 
+  ALGEBRA_HOST_DEVICE
   constexpr const auto &get() const { return m_data; }
 
   /// Subscript operator[]
   /// @{
+  ALGEBRA_HOST_DEVICE
   constexpr decltype(auto) operator[](std::size_t i) { return m_data[i]; }
+  ALGEBRA_HOST_DEVICE
   constexpr decltype(auto) operator[](std::size_t i) const { return m_data[i]; }
   /// @}
 
   /// Operator*=
   ///
   /// @return Vector expression/return type according to the operation.
+  ALGEBRA_HOST_DEVICE
   constexpr decltype(auto) operator*=(value_type factor) noexcept {
     return m_data *= factor;
   }
@@ -154,26 +152,29 @@ class alignas(
   template <std::size_t M, typename o_value_t,
             template <typename, std::size_t> class o_array_t,
             template <typename, std::size_t> class p_array_t>
-  friend constexpr bool operator==(
+  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;
 
   template <std::size_t M, typename o_value_t,
             template <typename, std::size_t> class o_array_t,
             template <typename, std::size_t> class p_array_t, bool>
-  friend constexpr bool operator==(const vector<M, o_value_t, o_array_t> &,
-                                   const p_array_t<o_value_t, M> &) noexcept;
+  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;
   /// @}
 
   /// Inequality operator
   template <typename other_type>
-  constexpr bool operator!=(const other_type &rhs) const noexcept {
+  ALGEBRA_HOST_DEVICE constexpr bool operator!=(
+      const other_type &rhs) const noexcept {
     return ((*this == rhs) == false);
   }
 
   /// Elementwise comparison. Can result in a vector-of-masks for SoA vectors
   template <typename other_type>
-  constexpr auto compare(const other_type &rhs) const noexcept {
+  ALGEBRA_HOST_DEVICE constexpr auto compare(
+      const other_type &rhs) const noexcept {
     using result_t = decltype(m_data[0] == rhs[0]);
 
     array_t<result_t, N> comp;
@@ -191,7 +192,8 @@ class alignas(
  private:
   /// Sets the trailing uninitialized values to zero.
   template <std::size_t... Is>
-  constexpr void zero_fill(std::index_sequence<Is...>) noexcept {
+  ALGEBRA_HOST_DEVICE constexpr void zero_fill(
+      std::index_sequence<Is...>) noexcept {
     ((m_data[simd_size() - sizeof...(Is) + Is] = value_t(0)), ...);
   }
 };
@@ -202,7 +204,8 @@ class alignas(
 template <std::size_t N, typename value_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>) constexpr bool operator==(
+  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 {
 
@@ -219,7 +222,8 @@ requires(std::is_scalar_v<value_t>) constexpr bool operator==(
 template <std::size_t N, typename value_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>) constexpr bool operator==(
+  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 {
 
@@ -236,57 +240,62 @@ requires(!std::is_scalar_v<value_t>) constexpr bool operator==(
 template <std::size_t N, typename value_t,
           template <typename, std::size_t> class array_t,
           template <typename, std::size_t> class o_array_t>
-constexpr bool operator==(const vector<N, value_t, array_t> &lhs,
-                          const vector<N, value_t, o_array_t> &rhs) noexcept {
+ALGEBRA_HOST_DEVICE constexpr bool operator==(
+    const vector<N, value_t, array_t> &lhs,
+    const vector<N, value_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,              \
-            template <typename, std::size_t> class array_t>                  \
-  requires(std::is_scalar_v<scalar_t>) 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);                          \
-  }                                                                          \
-  template <std::size_t N, typename value_t, typename scalar_t,              \
-            template <typename, std::size_t> class array_t>                  \
-  requires(std::is_scalar_v<scalar_t>) 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;                          \
-  }                                                                          \
-  template <std::size_t N, typename value_t,                                 \
-            template <typename, std::size_t> class array_t>                  \
-  constexpr decltype(auto) operator OP(                                      \
-      const vector<N, value_t, array_t> &lhs,                                \
-      const vector<N, value_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>) constexpr decltype(auto)                \
-  operator OP(const vector<N, value_t, array_t> &lhs,                        \
-              const other_type &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>) constexpr decltype(auto)                \
-  operator OP(const other_type &lhs,                                         \
-              const vector<N, value_t, array_t> &rhs) noexcept {             \
-    return lhs OP rhs.m_data;                                                \
+#define DECLARE_vector_OPERATORS(OP)                                          \
+  template <std::size_t N, typename value_t, typename scalar_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);                           \
+  }                                                                           \
+  template <std::size_t N, typename value_t, typename scalar_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;                           \
+  }                                                                           \
+  template <std::size_t N, typename value_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 {                      \
+    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 {                                       \
+    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 {                      \
+    return lhs OP rhs.m_data;                                                 \
   }
 
 // Implement all arithmetic operations on top of @c vector.
diff --git a/storage/vc_aos/include/algebra/storage/vc_aos.hpp b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
index 6286e212..44220505 100644
--- a/storage/vc_aos/include/algebra/storage/vc_aos.hpp
+++ b/storage/vc_aos/include/algebra/storage/vc_aos.hpp
@@ -55,7 +55,7 @@ template <typename T>
 using point3 = vector3<T>;
 /// 6-element "vector" type, using @c Vc::SimdArray
 template <typename T>
-using vector6 = vector_type<T, 8>;
+using vector6 = vector_type<T, 6>;
 /// 8-element "vector" type, using @c Vc::SimdArray
 template <typename T>
 using vector8 = vector_type<T, 8>;
diff --git a/tests/accelerator/sycl/common/test_sycl_basics.hpp b/tests/accelerator/sycl/common/test_sycl_basics.hpp
index 32933cc1..630185f2 100644
--- a/tests/accelerator/sycl/common/test_sycl_basics.hpp
+++ b/tests/accelerator/sycl/common/test_sycl_basics.hpp
@@ -40,7 +40,7 @@ class test_sycl_basics : public test_basics_base<T> {
 TYPED_TEST_SUITE_P(test_sycl_basics);
 
 /// Test for some basic 2D "vector operations"
-/*TYPED_TEST_P(test_sycl_basics, vector_2d_ops) {
+TYPED_TEST_P(test_sycl_basics, vector_2d_ops) {
 
   // Don't run the test at double precision, if the SYCL device doesn't
   // support it.
@@ -137,7 +137,7 @@ TYPED_TEST_P(test_sycl_basics, matrix22_ops) {
 
   // Compare the outputs.
   this->compareOutputs();
-}*/
+}
 
 /// Test for some operations with @c transform3
 TYPED_TEST_P(test_sycl_basics, transform3) {
@@ -165,6 +165,5 @@ TYPED_TEST_P(test_sycl_basics, transform3) {
   this->compareOutputs();
 }
 
-REGISTER_TYPED_TEST_SUITE_P(test_sycl_basics, /*vector_2d_ops, vector_3d_ops,
-                            matrix64_ops, matrix22_ops,*/
-                            transform3);
+REGISTER_TYPED_TEST_SUITE_P(test_sycl_basics, vector_2d_ops, vector_3d_ops,
+                            matrix64_ops, matrix22_ops, transform3);