diff --git a/.gitignore b/.gitignore
index 259148fa..8baa945a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,3 +30,7 @@
 *.exe
 *.out
 *.app
+
+# Local folders
+/bin
+/build
\ No newline at end of file
diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt
new file mode 100644
index 00000000..49dd801e
--- /dev/null
+++ b/benchmarks/CMakeLists.txt
@@ -0,0 +1,72 @@
+# Algebra plugins library, part of the ACTS project (R&D line)
+#
+# (c) 2023 CERN for the benefit of the ACTS project
+#
+# Mozilla Public License Version 2.0
+
+# Set the default C++ compiler flags.
+include( algebra-plugins-compiler-options-cpp )
+
+# Set up an INTERFACE library for the common header(s).
+add_library( algebra_bench_common INTERFACE )
+target_include_directories( algebra_bench_common INTERFACE
+   "${CMAKE_CURRENT_SOURCE_DIR}/common/include" )
+if( "${CMAKE_CXX_COMPILER_ID}" MATCHES "MSVC" )
+   target_compile_definitions( algebra_bench_common INTERFACE
+      -D_USE_MATH_DEFINES )
+endif()
+add_library( algebra::bench_common ALIAS algebra_bench_common )
+
+# Set up all of the (available) "host" benchmarks.
+add_library( algebra_bench_array INTERFACE )
+target_include_directories( algebra_bench_array INTERFACE 
+   "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/array/include>" )
+target_link_libraries(algebra_bench_array INTERFACE algebra::array_cmath 
+                                                    algebra::common_storage)
+add_library( algebra::bench_array ALIAS algebra_bench_array )
+
+algebra_add_benchmark( array_getter
+   "array/array_getter.cpp"
+   LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                   algebra::bench_array algebra::array_cmath )
+algebra_add_benchmark( array_vector
+   "array/array_vector.cpp"
+   LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                   algebra::bench_array algebra::array_cmath )
+
+if( ALGEBRA_PLUGINS_INCLUDE_EIGEN )
+   add_library( algebra_bench_eigen INTERFACE )
+   target_include_directories( algebra_bench_eigen INTERFACE 
+      "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/eigen/include>" )
+   target_link_libraries(algebra_bench_eigen INTERFACE algebra::eigen_eigen
+                                                    algebra::common_storage)
+   add_library( algebra::bench_eigen ALIAS algebra_bench_eigen )
+
+   algebra_add_benchmark( eigen_getter
+      "eigen/eigen_getter.cpp"
+      LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                     algebra::bench_eigen algebra::eigen_eigen )
+   algebra_add_benchmark( eigen_vector
+      "eigen/eigen_vector.cpp"
+      LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                     algebra::bench_eigen algebra::eigen_eigen )
+endif()
+
+if( ALGEBRA_PLUGINS_INCLUDE_VC )
+   add_library( algebra_bench_vc INTERFACE )
+   target_include_directories( algebra_bench_vc INTERFACE 
+      "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/vc/include>" )
+   target_link_libraries(algebra_bench_vc INTERFACE algebra::vc_vc
+                                                    algebra::vc_soa
+                                                    algebra::common_storage)
+   add_library( algebra::bench_vc ALIAS algebra_bench_vc )
+
+   algebra_add_benchmark( vc_getter
+      "vc/vc_soa_getter.cpp"
+      LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                     algebra::bench_vc algebra::vc_vc algebra::vc_soa )
+   algebra_add_benchmark( vc_vector
+      "vc/vc_soa_vector.cpp"
+      LINK_LIBRARIES benchmark::benchmark algebra::bench_common
+                     algebra::bench_vc algebra::vc_vc algebra::vc_soa )
+endif()
diff --git a/benchmarks/array/array_getter.cpp b/benchmarks/array/array_getter.cpp
new file mode 100644
index 00000000..e53cb24a
--- /dev/null
+++ b/benchmarks/array/array_getter.cpp
@@ -0,0 +1,93 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/array_cmath.hpp"
+#include "benchmark/array/data_generator.hpp"
+#include "benchmark/common/benchmark_getter.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+// System include(s)
+#include <string>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+  constexpr std::size_t n_warmup{static_cast<std::size_t>(0.1 * n_samples)};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg{};
+  cfg.n_samples(n_samples).n_warmup(n_warmup);
+  cfg.do_sleep(false);
+
+  vector_unaryOP_bm<array::vector3, float, bench_op::phi> v_phi_s{cfg};
+  vector_unaryOP_bm<array::vector3, float, bench_op::theta> v_theta_s{cfg};
+  vector_unaryOP_bm<array::vector3, float, bench_op::perp> v_perp_s{cfg};
+  vector_unaryOP_bm<array::vector3, float, bench_op::norm> v_norm_s{cfg};
+  vector_unaryOP_bm<array::vector3, float, bench_op::eta> v_eta_s{cfg};
+
+  vector_unaryOP_bm<array::vector3, double, bench_op::phi> v_phi_d{cfg};
+  vector_unaryOP_bm<array::vector3, double, bench_op::theta> v_theta_d{cfg};
+  vector_unaryOP_bm<array::vector3, double, bench_op::perp> v_perp_d{cfg};
+  vector_unaryOP_bm<array::vector3, double, bench_op::norm> v_norm_d{cfg};
+  vector_unaryOP_bm<array::vector3, double, bench_op::eta> v_eta_d{cfg};
+
+  std::cout << "Algebra-Plugins 'getter' benchmark (std::array)\n"
+            << "-----------------------------------------------\n\n"
+            << cfg;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_phi_s.name() + "_single").c_str(), v_phi_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_phi_d.name() + "_double").c_str(), v_phi_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_theta_s.name() + "_single").c_str(),
+                                 v_theta_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_theta_d.name() + "_double").c_str(),
+                                 v_theta_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_perp_s.name() + "_single").c_str(),
+                                 v_perp_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_perp_d.name() + "_double").c_str(),
+                                 v_perp_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_norm_s.name() + "_single").c_str(),
+                                 v_norm_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_norm_d.name() + "_double").c_str(),
+                                 v_norm_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_eta_s.name() + "_single").c_str(), v_eta_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_eta_d.name() + "_double").c_str(), v_eta_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}
diff --git a/benchmarks/array/array_vector.cpp b/benchmarks/array/array_vector.cpp
new file mode 100644
index 00000000..20b7ae05
--- /dev/null
+++ b/benchmarks/array/array_vector.cpp
@@ -0,0 +1,93 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/array_cmath.hpp"
+#include "benchmark/array/data_generator.hpp"
+#include "benchmark/common/benchmark_vector.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+// System include(s)
+#include <string>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+  constexpr std::size_t n_warmup{static_cast<std::size_t>(0.1 * n_samples)};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg{};
+  cfg.n_samples(n_samples).n_warmup(n_warmup);
+  cfg.do_sleep(false);
+
+  vector_binaryOP_bm<array::vector3, float, bench_op::add> v_add_s{cfg};
+  vector_binaryOP_bm<array::vector3, float, bench_op::sub> v_sub_s{cfg};
+  vector_binaryOP_bm<array::vector3, float, bench_op::dot> v_dot_s{cfg};
+  vector_binaryOP_bm<array::vector3, float, bench_op::cross> v_cross_s{cfg};
+  vector_unaryOP_bm<array::vector3, float, bench_op::normalize> v_normalize_s{
+      cfg};
+
+  vector_binaryOP_bm<array::vector3, double, bench_op::add> v_add_d{cfg};
+  vector_binaryOP_bm<array::vector3, double, bench_op::sub> v_sub_d{cfg};
+  vector_binaryOP_bm<array::vector3, double, bench_op::dot> v_dot_d{cfg};
+  vector_binaryOP_bm<array::vector3, double, bench_op::cross> v_cross_d{cfg};
+  vector_unaryOP_bm<array::vector3, double, bench_op::normalize> v_normalize_d{
+      cfg};
+
+  std::cout << "Algebra-Plugins 'vector' benchmark (std::array)\n"
+            << "-----------------------------------------------\n\n"
+            << cfg;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_add_s.name() + "_single").c_str(), v_add_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_add_d.name() + "_double").c_str(), v_add_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_sub_s.name() + "_single").c_str(), v_sub_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_sub_d.name() + "_double").c_str(), v_sub_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_dot_s.name() + "_single").c_str(), v_dot_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_dot_d.name() + "_double").c_str(), v_dot_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_cross_s.name() + "_single").c_str(),
+                                 v_cross_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_cross_d.name() + "_double").c_str(),
+                                 v_cross_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_normalize_s.name() + "_single").c_str(),
+                                 v_normalize_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_normalize_d.name() + "_double").c_str(),
+                                 v_normalize_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}
diff --git a/benchmarks/array/include/benchmark/array/data_generator.hpp b/benchmarks/array/include/benchmark/array/data_generator.hpp
new file mode 100644
index 00000000..7f4ac7ca
--- /dev/null
+++ b/benchmarks/array/include/benchmark/array/data_generator.hpp
@@ -0,0 +1,36 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "algebra/storage/array.hpp"
+#include "algebra/storage/vector.hpp"
+
+// System include(s)
+#include <algorithm>
+#include <random>
+#include <vector>
+
+namespace algebra {
+
+/// Fill an @c std::array based vector with random values
+template <typename vector_t>
+inline void fill_random(std::vector<vector_t> &collection) {
+
+  // Generate a vector of the right type with random values
+  std::random_device rd;
+  std::mt19937 mt(rd());
+  std::uniform_real_distribution<typename vector_t::value_type> dist(0.f, 1.f);
+
+  auto rand_obj = [&]() { return vector_t{dist(mt), dist(mt), dist(mt)}; };
+
+  collection.resize(collection.capacity());
+  std::generate(collection.begin(), collection.end(), rand_obj);
+}
+
+}  // namespace algebra
\ No newline at end of file
diff --git a/benchmarks/common/include/benchmark/common/benchmark_base.hpp b/benchmarks/common/include/benchmark/common/benchmark_base.hpp
new file mode 100644
index 00000000..b3ca9fdd
--- /dev/null
+++ b/benchmarks/common/include/benchmark/common/benchmark_base.hpp
@@ -0,0 +1,108 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+// System include(s)
+#include <ostream>
+#include <string>
+
+namespace algebra {
+
+/// Base type for linear algebra benchmarks with google benchmark
+struct benchmark_base {
+  /// Local configuration type
+  struct configuration {
+    /// Size of data sample to be used in benchmark
+    std::size_t m_samples{100u};
+    /// Run a number of operations before the benchmark
+    bool m_warmup = true;
+    // Sleep after building data sample
+    bool m_sleep = false;
+    // Size of data in warm-up round
+    std::size_t m_n_warmup{static_cast<std::size_t>(0.1f * m_samples)};
+    // Size of data in warm-up round
+    std::size_t m_n_sleep{1u};
+
+    /// Setters
+    /// @{
+    configuration& n_samples(std::size_t n) {
+      m_samples = n;
+      return *this;
+    }
+    configuration& do_warmup(bool b) {
+      m_warmup = b;
+      return *this;
+    }
+    configuration& n_warmup(std::size_t n) {
+      m_n_warmup = n;
+      m_warmup = true;
+      return *this;
+    }
+    configuration& do_sleep(bool b) {
+      m_sleep = b;
+      return *this;
+    }
+    configuration& n_sleep(std::size_t n) {
+      m_n_sleep = n;
+      m_sleep = true;
+      return *this;
+    }
+    /// @}
+
+    /// Getters
+    /// @{
+    std::size_t n_samples() const { return m_samples; }
+    constexpr bool do_warmup() const { return m_warmup; }
+    constexpr std::size_t n_warmup() const { return m_n_warmup; }
+    constexpr bool do_sleep() const { return m_sleep; }
+    constexpr std::size_t n_sleep() const { return m_n_sleep; }
+    /// @}
+
+    /// Print configuration
+    friend std::ostream& operator<<(std::ostream& os, const configuration& c);
+  };
+
+  /// The benchmark configuration
+  configuration m_cfg{};
+
+  /// Default construction
+  benchmark_base() = default;
+
+  /// Construct from an externally provided configuration @param cfg
+  benchmark_base(configuration cfg) : m_cfg{cfg} {}
+
+  /// @returns the benchmark configuration
+  configuration& config() { return m_cfg; }
+
+  /// Default destructor
+  virtual ~benchmark_base() = default;
+
+  /// @returns the benchmark name
+  virtual std::string name() const = 0;
+
+  /// Benchmark case
+  virtual void operator()(::benchmark::State&) = 0;
+};
+
+std::ostream& operator<<(std::ostream& os,
+                         const benchmark_base::configuration& cfg) {
+  os << " -> running:\t " << cfg.n_samples() << " samples" << std::endl;
+  if (cfg.do_warmup()) {
+    os << " -> warmup: \t " << cfg.n_warmup() << " samples" << std::endl;
+  }
+  if (cfg.do_sleep()) {
+    os << " -> cool down:\t " << cfg.n_sleep() << "s" << std::endl;
+  }
+  os << std::endl;
+  return os;
+}
+
+}  // namespace algebra
\ No newline at end of file
diff --git a/benchmarks/common/include/benchmark/common/benchmark_getter.hpp b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
new file mode 100644
index 00000000..f697e2f7
--- /dev/null
+++ b/benchmarks/common/include/benchmark/common/benchmark_getter.hpp
@@ -0,0 +1,53 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "benchmark_vector.hpp"
+
+namespace algebra::bench_op {
+
+// Functions to be benchmarked
+
+struct phi {
+  inline static const std::string name{"phi"};
+  template <typename vector_t>
+  decltype(auto) operator()(const vector_t &a) const {
+    return algebra::getter::phi(a);
+  }
+};
+struct theta {
+  inline static const std::string name{"theta"};
+  template <typename vector_t>
+  decltype(auto) operator()(const vector_t &a) const {
+    return algebra::getter::theta(a);
+  }
+};
+struct perp {
+  inline static const std::string name{"perp"};
+  template <typename vector_t>
+  decltype(auto) operator()(const vector_t &a) const {
+    return algebra::getter::perp(a);
+  }
+};
+struct norm {
+  inline static const std::string name{"norm"};
+  template <typename vector_t>
+  decltype(auto) operator()(const vector_t &a) const {
+    return algebra::getter::norm(a);
+  }
+};
+struct eta {
+  inline static const std::string name{"eta"};
+  template <typename vector_t>
+  decltype(auto) operator()(const vector_t &a) const {
+    return algebra::getter::eta(a);
+  }
+};
+
+}  // namespace algebra::bench_op
\ No newline at end of file
diff --git a/benchmarks/common/include/benchmark/common/benchmark_vector.hpp b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
new file mode 100644
index 00000000..59572195
--- /dev/null
+++ b/benchmarks/common/include/benchmark/common/benchmark_vector.hpp
@@ -0,0 +1,181 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "benchmark_base.hpp"
+
+// System include(s)
+#include <chrono>
+#include <iostream>
+#include <string>
+#include <thread>
+#include <vector>
+
+namespace algebra {
+
+template <typename vector_t>
+void fill_random(std::vector<vector_t> &);
+
+/// Benchmark for vector operations
+template <typename vector_t>
+struct vector_bm : public benchmark_base {
+
+  /// Prefix for the benchmark name
+  inline static const std::string name{"vector"};
+
+  std::vector<vector_t> a, b, results;
+
+  /// No default construction: Cannot prepare data
+  vector_bm() = delete;
+
+  /// Construct from an externally provided configuration @param cfg
+  vector_bm(benchmark_base::configuration cfg) : benchmark_base{cfg} {
+
+    const std::size_t n_data{this->m_cfg.n_samples() + this->m_cfg.n_warmup()};
+
+    a.reserve(n_data);
+    b.reserve(n_data);
+
+    fill_random(a);
+    fill_random(b);
+  }
+
+  /// Clear state
+  virtual ~vector_bm() {
+    a.clear();
+    b.clear();
+  }
+};
+
+/// Benchmark elementwise addition of vectors
+template <template <typename> class vector_t, typename scalar_t,
+          typename unaryOP>
+struct vector_unaryOP_bm : public vector_bm<vector_t<scalar_t>> {
+  using base_type = vector_bm<vector_t<scalar_t>>;
+
+  vector_unaryOP_bm() = delete;
+  vector_unaryOP_bm(benchmark_base::configuration cfg) : base_type{cfg} {}
+  std::string name() const override {
+    return base_type::name + "_" + unaryOP::name;
+  }
+
+  void operator()(::benchmark::State &state) override {
+
+    const std::size_t n_samples{this->m_cfg.n_samples()};
+    const std::size_t n_warmup{this->m_cfg.n_warmup()};
+
+    // Spin down before benchmark (Thread zero is counting the clock)
+    if (state.thread_index() == 0 and this->m_cfg.do_sleep()) {
+      std::this_thread::sleep_for(std::chrono::seconds(this->m_cfg.n_sleep()));
+    }
+
+    // Run the benchmark
+    for (auto _ : state) {
+      // Warm-up
+      state.PauseTiming();
+      if (this->m_cfg.do_warmup()) {
+        for (std::size_t i{0u}; i < n_warmup; ++i) {
+          ::benchmark::DoNotOptimize(unaryOP{}(this->a[i]));
+          benchmark::ClobberMemory();
+        }
+      }
+      state.ResumeTiming();
+
+      for (std::size_t i{n_warmup}; i < n_samples + n_warmup; ++i) {
+        ::benchmark::DoNotOptimize(unaryOP{}(this->a[i]));
+        benchmark::ClobberMemory();
+      }
+    }
+  }
+};
+
+/// Benchmark elementwise addition of vectors
+template <template <typename> class vector_t, typename scalar_t,
+          typename binaryOP>
+struct vector_binaryOP_bm : public vector_bm<vector_t<scalar_t>> {
+  using base_type = vector_bm<vector_t<scalar_t>>;
+
+  vector_binaryOP_bm() = delete;
+  vector_binaryOP_bm(benchmark_base::configuration cfg) : base_type{cfg} {}
+  std::string name() const override {
+    return base_type::name + "_" + binaryOP::name;
+  }
+
+  void operator()(::benchmark::State &state) override {
+
+    const std::size_t n_samples{this->m_cfg.n_samples()};
+    const std::size_t n_warmup{this->m_cfg.n_warmup()};
+
+    // Spin down before benchmark (Thread zero is counting the clock)
+    if (state.thread_index() == 0 and this->m_cfg.do_sleep()) {
+      std::this_thread::sleep_for(std::chrono::seconds(this->m_cfg.n_sleep()));
+    }
+
+    // Run the benchmark
+    for (auto _ : state) {
+      // Warm-up
+      state.PauseTiming();
+      if (this->m_cfg.do_warmup()) {
+        for (std::size_t i{0u}; i < n_warmup; ++i) {
+          ::benchmark::DoNotOptimize(binaryOP{}(this->a[i], this->b[i]));
+          benchmark::ClobberMemory();
+        }
+      }
+      state.ResumeTiming();
+
+      for (std::size_t i{n_warmup}; i < n_samples + n_warmup; ++i) {
+        ::benchmark::DoNotOptimize(binaryOP{}(this->a[i], this->b[i]));
+        benchmark::ClobberMemory();
+      }
+    }
+  }
+};
+
+// Functions to be benchmarked
+namespace bench_op {
+
+struct add {
+  inline static const std::string name{"add"};
+  template <typename vector_t>
+  decltype(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>
+  decltype(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>
+  decltype(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>
+  decltype(auto) operator()(const vector_t &a, const vector_t &b) const {
+    return algebra::vector::cross(a, b);
+  }
+};
+struct normalize {
+  inline static const std::string name{"normalize"};
+  template <typename vector_t>
+  decltype(auto) operator()(const vector_t &a) const {
+    return algebra::vector::normalize(a);
+  }
+};
+
+}  // namespace bench_op
+
+}  // namespace algebra
\ No newline at end of file
diff --git a/benchmarks/eigen/eigen_getter.cpp b/benchmarks/eigen/eigen_getter.cpp
new file mode 100644
index 00000000..14ce10eb
--- /dev/null
+++ b/benchmarks/eigen/eigen_getter.cpp
@@ -0,0 +1,93 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/eigen_eigen.hpp"
+#include "benchmark/common/benchmark_getter.hpp"
+#include "benchmark/eigen/data_generator.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+// System include(s)
+#include <string>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+  constexpr std::size_t n_warmup{static_cast<std::size_t>(0.1 * n_samples)};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg{};
+  cfg.n_samples(n_samples).n_warmup(n_warmup);
+  cfg.do_sleep(false);
+
+  vector_unaryOP_bm<eigen::vector3, float, bench_op::phi> v_phi_s{cfg};
+  vector_unaryOP_bm<eigen::vector3, float, bench_op::theta> v_theta_s{cfg};
+  vector_unaryOP_bm<eigen::vector3, float, bench_op::perp> v_perp_s{cfg};
+  vector_unaryOP_bm<eigen::vector3, float, bench_op::norm> v_norm_s{cfg};
+  vector_unaryOP_bm<eigen::vector3, float, bench_op::eta> v_eta_s{cfg};
+
+  vector_unaryOP_bm<eigen::vector3, double, bench_op::phi> v_phi_d{cfg};
+  vector_unaryOP_bm<eigen::vector3, double, bench_op::theta> v_theta_d{cfg};
+  vector_unaryOP_bm<eigen::vector3, double, bench_op::perp> v_perp_d{cfg};
+  vector_unaryOP_bm<eigen::vector3, double, bench_op::norm> v_norm_d{cfg};
+  vector_unaryOP_bm<eigen::vector3, double, bench_op::eta> v_eta_d{cfg};
+
+  std::cout << "Algebra-Plugins 'getter' benchmark (Eigen)\n"
+            << "------------------------------------------\n\n"
+            << cfg;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_phi_s.name() + "_single").c_str(), v_phi_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_phi_d.name() + "_double").c_str(), v_phi_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_theta_s.name() + "_single").c_str(),
+                                 v_theta_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_theta_d.name() + "_double").c_str(),
+                                 v_theta_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_perp_s.name() + "_single").c_str(),
+                                 v_perp_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_perp_d.name() + "_double").c_str(),
+                                 v_perp_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_norm_s.name() + "_single").c_str(),
+                                 v_norm_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_norm_d.name() + "_double").c_str(),
+                                 v_norm_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_eta_s.name() + "_single").c_str(), v_eta_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_eta_d.name() + "_double").c_str(), v_eta_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}
diff --git a/benchmarks/eigen/eigen_vector.cpp b/benchmarks/eigen/eigen_vector.cpp
new file mode 100644
index 00000000..ba3c03a0
--- /dev/null
+++ b/benchmarks/eigen/eigen_vector.cpp
@@ -0,0 +1,93 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/eigen_eigen.hpp"
+#include "benchmark/common/benchmark_vector.hpp"
+#include "benchmark/eigen/data_generator.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+// System include(s)
+#include <string>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+  constexpr std::size_t n_warmup{static_cast<std::size_t>(0.1 * n_samples)};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg{};
+  cfg.n_samples(n_samples).n_warmup(n_warmup);
+  cfg.do_sleep(false);
+
+  vector_binaryOP_bm<eigen::vector3, float, bench_op::add> v_add_s{cfg};
+  vector_binaryOP_bm<eigen::vector3, float, bench_op::sub> v_sub_s{cfg};
+  vector_binaryOP_bm<eigen::vector3, float, bench_op::dot> v_dot_s{cfg};
+  vector_binaryOP_bm<eigen::vector3, float, bench_op::cross> v_cross_s{cfg};
+  vector_unaryOP_bm<eigen::vector3, float, bench_op::normalize> v_normalize_s{
+      cfg};
+
+  vector_binaryOP_bm<eigen::vector3, double, bench_op::add> v_add_d{cfg};
+  vector_binaryOP_bm<eigen::vector3, double, bench_op::sub> v_sub_d{cfg};
+  vector_binaryOP_bm<eigen::vector3, double, bench_op::dot> v_dot_d{cfg};
+  vector_binaryOP_bm<eigen::vector3, double, bench_op::cross> v_cross_d{cfg};
+  vector_unaryOP_bm<eigen::vector3, double, bench_op::normalize> v_normalize_d{
+      cfg};
+
+  std::cout << "Algebra-Plugins 'vector' benchmark (Eigen)\n"
+            << "------------------------------------------\n\n"
+            << cfg;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_add_s.name() + "_single").c_str(), v_add_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_add_d.name() + "_double").c_str(), v_add_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_sub_s.name() + "_single").c_str(), v_sub_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_sub_d.name() + "_double").c_str(), v_sub_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_dot_s.name() + "_single").c_str(), v_dot_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_dot_d.name() + "_double").c_str(), v_dot_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_cross_s.name() + "_single").c_str(),
+                                 v_cross_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_cross_d.name() + "_double").c_str(),
+                                 v_cross_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_normalize_s.name() + "_single").c_str(),
+                                 v_normalize_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_normalize_d.name() + "_double").c_str(),
+                                 v_normalize_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}
diff --git a/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp
new file mode 100644
index 00000000..17cbad7b
--- /dev/null
+++ b/benchmarks/eigen/include/benchmark/eigen/data_generator.hpp
@@ -0,0 +1,30 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "algebra/storage/eigen.hpp"
+#include "algebra/storage/vector.hpp"
+
+// System include(s)
+#include <algorithm>
+#include <random>
+
+namespace algebra {
+
+/// Fill a @c Eigen3 based vector with random values
+template <typename vector_t>
+inline void fill_random(std::vector<vector_t> &collection) {
+
+  auto rand_obj = []() { return vector_t::Random(); };
+
+  collection.resize(collection.capacity());
+  std::generate(collection.begin(), collection.end(), rand_obj);
+}
+
+}  // namespace algebra
\ No newline at end of file
diff --git a/benchmarks/vc/include/benchmark/vc/data_generator.hpp b/benchmarks/vc/include/benchmark/vc/data_generator.hpp
new file mode 100644
index 00000000..3820e20f
--- /dev/null
+++ b/benchmarks/vc/include/benchmark/vc/data_generator.hpp
@@ -0,0 +1,53 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "algebra/storage/vc.hpp"
+#include "algebra/storage/vc_soa.hpp"
+#include "algebra/storage/vector.hpp"
+
+// System include(s)
+#include <algorithm>
+#include <vector>
+
+namespace algebra {
+
+// @todo Leave for later
+
+/// Fill a @c Vc::SimdArray based vector with random values
+/*template <typename vector_t>
+inline void fill_random(
+    std::vector<vector3_s<value_t>, allocator_t<vector3_s<value_t>>>
+        &collection) {
+
+  using vector_t = vector3_s<value_t>;
+  // Generate a vector of the right type with random values
+  auto rand_obj = [&]() { return vector_t{vector_t::array_type::Random()}; };
+
+  std::generate(collection.begin(), collection.end(), rand_obj);
+}*/
+
+/// Fill a @c Vc::Vector based vector with random values
+template <typename vector_soa_t>
+inline void fill_random(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;
+    vector_soa_t tmp{};
+    tmp[0] = simd_vector_t::Random();
+    tmp[1] = simd_vector_t::Random();
+    tmp[2] = simd_vector_t::Random();
+    return tmp;
+  };
+
+  collection.resize(collection.capacity());
+  std::generate(collection.begin(), collection.end(), rand_obj);
+}
+
+}  // namespace algebra
\ No newline at end of file
diff --git a/benchmarks/vc/vc_soa_getter.cpp b/benchmarks/vc/vc_soa_getter.cpp
new file mode 100644
index 00000000..95da5c94
--- /dev/null
+++ b/benchmarks/vc/vc_soa_getter.cpp
@@ -0,0 +1,102 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/vc_soa.hpp"
+#include "benchmark/common/benchmark_getter.hpp"
+#include "benchmark/vc/data_generator.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+// System include(s)
+#include <string>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg_s{};
+  // Reduce the number of samples, since a single SoA struct contains multiple
+  // vectors
+  cfg_s.n_samples(n_samples / Vc::float_v::Size)
+      .n_warmup(static_cast<std::size_t>(0.1 * cfg_s.n_samples()));
+  cfg_s.do_sleep(false);
+
+  // For double precision we need more samples (less vectors per SoA)
+  algebra::benchmark_base::configuration cfg_d{cfg_s};
+  cfg_d.n_samples(n_samples / Vc::double_v::Size)
+      .n_warmup(static_cast<std::size_t>(0.1 * cfg_d.n_samples()));
+
+  vector_unaryOP_bm<vc_soa::vector3, float, bench_op::phi> v_phi_s{cfg_s};
+  vector_unaryOP_bm<vc_soa::vector3, float, bench_op::theta> v_theta_s{cfg_s};
+  vector_unaryOP_bm<vc_soa::vector3, float, bench_op::perp> v_perp_s{cfg_s};
+  vector_unaryOP_bm<vc_soa::vector3, float, bench_op::norm> v_norm_s{cfg_s};
+  vector_unaryOP_bm<vc_soa::vector3, float, bench_op::eta> v_eta_s{cfg_s};
+
+  vector_unaryOP_bm<vc_soa::vector3, double, bench_op::phi> v_phi_d{cfg_d};
+  vector_unaryOP_bm<vc_soa::vector3, double, bench_op::theta> v_theta_d{cfg_d};
+  vector_unaryOP_bm<vc_soa::vector3, double, bench_op::perp> v_perp_d{cfg_d};
+  vector_unaryOP_bm<vc_soa::vector3, double, bench_op::norm> v_norm_d{cfg_d};
+  vector_unaryOP_bm<vc_soa::vector3, double, bench_op::eta> v_eta_d{cfg_d};
+
+  std::cout << "Algebra-Plugins 'getter' benchmark (Vc SoA)\n"
+            << "-------------------------------------------\n\n"
+            << "(single)\n"
+            << cfg_s << "(double)\n"
+            << cfg_d;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_phi_s.name() + "_single").c_str(), v_phi_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_phi_d.name() + "_double").c_str(), v_phi_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_theta_s.name() + "_single").c_str(),
+                                 v_theta_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_theta_d.name() + "_double").c_str(),
+                                 v_theta_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_perp_s.name() + "_single").c_str(),
+                                 v_perp_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_perp_d.name() + "_double").c_str(),
+                                 v_perp_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_norm_s.name() + "_single").c_str(),
+                                 v_norm_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_norm_d.name() + "_double").c_str(),
+                                 v_norm_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_eta_s.name() + "_single").c_str(), v_eta_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_eta_d.name() + "_double").c_str(), v_eta_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}
diff --git a/benchmarks/vc/vc_soa_vector.cpp b/benchmarks/vc/vc_soa_vector.cpp
new file mode 100644
index 00000000..b283e09f
--- /dev/null
+++ b/benchmarks/vc/vc_soa_vector.cpp
@@ -0,0 +1,99 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s)
+#include "algebra/vc_soa.hpp"
+#include "benchmark/common/benchmark_vector.hpp"
+#include "benchmark/vc/data_generator.hpp"
+
+// Benchmark include
+#include <benchmark/benchmark.h>
+
+using namespace algebra;
+
+/// Run vector benchmarks
+int main(int argc, char** argv) {
+
+  constexpr std::size_t n_samples{160000};
+
+  //
+  // Prepare benchmarks
+  //
+  algebra::benchmark_base::configuration cfg_s{};
+  // Reduce the number of samples, since a single SoA struct contains multiple
+  // vectors
+  cfg_s.n_samples(n_samples / Vc::float_v::Size)
+      .n_warmup(static_cast<std::size_t>(0.1 * cfg_s.n_samples()));
+  cfg_s.do_sleep(false);
+
+  // For double precision we need more samples (less vectors per SoA)
+  algebra::benchmark_base::configuration cfg_d{cfg_s};
+  cfg_d.n_samples(n_samples / Vc::double_v::Size)
+      .n_warmup(static_cast<std::size_t>(0.1 * cfg_d.n_samples()));
+
+  vector_binaryOP_bm<vc_soa::vector3, float, bench_op::add> v_add_s{cfg_s};
+  vector_binaryOP_bm<vc_soa::vector3, float, bench_op::sub> v_sub_s{cfg_s};
+  vector_binaryOP_bm<vc_soa::vector3, float, bench_op::dot> v_dot_s{cfg_s};
+  vector_binaryOP_bm<vc_soa::vector3, float, bench_op::cross> v_cross_s{cfg_s};
+  vector_unaryOP_bm<vc_soa::vector3, float, bench_op::normalize> v_normalize_s{
+      cfg_s};
+
+  vector_binaryOP_bm<vc_soa::vector3, double, bench_op::add> v_add_d{cfg_d};
+  vector_binaryOP_bm<vc_soa::vector3, double, bench_op::sub> v_sub_d{cfg_d};
+  vector_binaryOP_bm<vc_soa::vector3, double, bench_op::dot> v_dot_d{cfg_d};
+  vector_binaryOP_bm<vc_soa::vector3, double, bench_op::cross> v_cross_d{cfg_d};
+  vector_unaryOP_bm<vc_soa::vector3, double, bench_op::normalize> v_normalize_d{
+      cfg_d};
+
+  std::cout << "Algebra-Plugins 'vector' benchmark (Vc SoA)\n"
+            << "-------------------------------------------\n\n"
+            << "(single)\n"
+            << cfg_s << "(double)\n"
+            << cfg_d;
+
+  //
+  // Register all benchmarks
+  //
+  ::benchmark::RegisterBenchmark((v_add_s.name() + "_single").c_str(), v_add_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_add_d.name() + "_double").c_str(), v_add_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_sub_s.name() + "_single").c_str(), v_sub_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_sub_d.name() + "_double").c_str(), v_sub_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_dot_s.name() + "_single").c_str(), v_dot_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_dot_d.name() + "_double").c_str(), v_dot_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_cross_s.name() + "_single").c_str(),
+                                 v_cross_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_cross_d.name() + "_double").c_str(),
+                                 v_cross_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_normalize_s.name() + "_single").c_str(),
+                                 v_normalize_s)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+  ::benchmark::RegisterBenchmark((v_normalize_d.name() + "_double").c_str(),
+                                 v_normalize_d)
+      ->MeasureProcessCPUTime()
+      ->ThreadPerCpu();
+
+  ::benchmark::Initialize(&argc, argv);
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+}
\ No newline at end of file
diff --git a/cmake/algebra-plugins-functions.cmake b/cmake/algebra-plugins-functions.cmake
index 0fac495c..a6f01f7d 100644
--- a/cmake/algebra-plugins-functions.cmake
+++ b/cmake/algebra-plugins-functions.cmake
@@ -63,6 +63,25 @@ function( algebra_add_test name )
 
 endfunction( algebra_add_test )
 
+# Helper function for setting up the algebra plugin benchmarks.
+#
+# Usage: algebra_add_benchmark( array source1.cpp source2.cpp
+#                               LINK_LIBRARIES algebra::array )
+#
+function( algebra_add_benchmark name )
+
+   # Parse the function's options.
+   cmake_parse_arguments( ARG "" "" "LINK_LIBRARIES" ${ARGN} )
+
+   # Create the test executable.
+   set( bench_exe_name "algebra_benchmark_${name}" )
+   add_executable( ${bench_exe_name} ${ARG_UNPARSED_ARGUMENTS} )
+   if( ARG_LINK_LIBRARIES )
+      target_link_libraries( ${bench_exe_name} PRIVATE ${ARG_LINK_LIBRARIES} )
+   endif()
+
+endfunction( algebra_add_benchmark )
+
 # Helper function for adding individual flags to "flag variables".
 #
 # Usage: algebra_add_flag( CMAKE_CXX_FLAGS "-Wall" )
diff --git a/frontend/CMakeLists.txt b/frontend/CMakeLists.txt
index 0f214b08..8da36578 100644
--- a/frontend/CMakeLists.txt
+++ b/frontend/CMakeLists.txt
@@ -20,6 +20,7 @@ endif()
 if( ALGEBRA_PLUGINS_INCLUDE_VC )
    add_subdirectory( vc_cmath )
    add_subdirectory( vc_vc )
+   add_subdirectory( vc_soa )
 endif()
 
 if( ALGEBRA_PLUGINS_INCLUDE_FASTOR )
diff --git a/frontend/vc_soa/CMakeLists.txt b/frontend/vc_soa/CMakeLists.txt
new file mode 100644
index 00000000..a41cc132
--- /dev/null
+++ b/frontend/vc_soa/CMakeLists.txt
@@ -0,0 +1,12 @@
+# Algebra plugins library, part of the ACTS project (R&D line)
+#
+# (c) 2023 CERN for the benefit of the ACTS project
+#
+# Mozilla Public License Version 2.0
+
+# Set up the library.
+algebra_add_library( algebra_vc_soa vc_soa
+   "include/algebra/vc_soa.hpp" )
+target_link_libraries( algebra_vc_soa
+   INTERFACE algebra::common algebra::vc_storage algebra::cmath_math
+             algebra::vc_soa_math )
diff --git a/frontend/vc_soa/include/algebra/vc_soa.hpp b/frontend/vc_soa/include/algebra/vc_soa.hpp
new file mode 100644
index 00000000..08b8052d
--- /dev/null
+++ b/frontend/vc_soa/include/algebra/vc_soa.hpp
@@ -0,0 +1,80 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+#include "algebra/math/vc_soa.hpp"
+#include "algebra/storage/vc_soa.hpp"
+
+// System include(s).
+#include <cassert>
+#include <type_traits>
+
+/// @name Operators on @c algebra::storage::vector types
+/// @{
+
+using algebra::storage::operator*;
+using algebra::storage::operator/;
+using algebra::storage::operator-;
+using algebra::storage::operator+;
+
+/// @}
+
+namespace algebra {
+namespace vc_soa {
+
+/// @name Vc based transforms on @c algebra::vc::storage_type
+/// @{
+
+//@todo
+
+/// @}
+
+}  // namespace vc_soa
+
+namespace getter {
+
+/// @name Getter functions on @c algebra::vc_soa types
+/// @{
+
+using vc_soa::math::eta;
+using vc_soa::math::norm;
+using vc_soa::math::perp;
+using vc_soa::math::phi;
+using vc_soa::math::theta;
+
+/// @}
+
+}  // namespace getter
+
+namespace vector {
+
+/// @name Vector functions on @c algebra::vc_soa types
+/// @{
+
+using vc_soa::math::cross;
+using vc_soa::math::dot;
+using vc_soa::math::normalize;
+
+/// @}
+
+}  // namespace vector
+
+namespace matrix {
+
+using size_type = vc_soa::size_type;
+
+template <typename T, size_type N>
+using array_type = vc_soa::storage_type<T, N>;
+
+template <typename T, size_type ROWS, size_type COLS>
+using matrix_type = vc_soa::matrix_type<T, ROWS, COLS>;
+
+}  // namespace matrix
+
+}  // namespace algebra
\ No newline at end of file
diff --git a/math/CMakeLists.txt b/math/CMakeLists.txt
index 97045cee..eef16b01 100644
--- a/math/CMakeLists.txt
+++ b/math/CMakeLists.txt
@@ -18,4 +18,5 @@ if( ALGEBRA_PLUGINS_INCLUDE_SMATRIX )
 endif()
 if( ALGEBRA_PLUGINS_INCLUDE_VC )
    add_subdirectory( vc )
+   add_subdirectory( vc_soa )
 endif()
diff --git a/math/cmath/CMakeLists.txt b/math/cmath/CMakeLists.txt
index b31d7c89..bae4ff20 100644
--- a/math/cmath/CMakeLists.txt
+++ b/math/cmath/CMakeLists.txt
@@ -14,6 +14,7 @@ algebra_add_library(algebra_cmath_math cmath_math
    "include/algebra/math/impl/cmath_operators.hpp"
    "include/algebra/math/impl/cmath_transform3.hpp"
    "include/algebra/math/impl/cmath_vector.hpp"
+   "include/algebra/math/impl/cmath_vector_soa.hpp"
    # algorithms include
    "include/algebra/math/algorithms/matrix/decomposition/partial_pivot_lud.hpp"
    "include/algebra/math/algorithms/matrix/determinant/cofactor.hpp"
diff --git a/math/cmath/include/algebra/math/impl/cmath_operators.hpp b/math/cmath/include/algebra/math/impl/cmath_operators.hpp
index ef624da7..0c56df91 100644
--- a/math/cmath/include/algebra/math/impl/cmath_operators.hpp
+++ b/math/cmath/include/algebra/math/impl/cmath_operators.hpp
@@ -1,6 +1,6 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2020-2022 CERN for the benefit of the ACTS project
+ * (c) 2020-2023 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -12,6 +12,7 @@
 
 // System include(s).
 #include <cstddef>
+#include <type_traits>
 
 namespace algebra::cmath {
 
diff --git a/math/vc_soa/CMakeLists.txt b/math/vc_soa/CMakeLists.txt
new file mode 100644
index 00000000..d3418d8c
--- /dev/null
+++ b/math/vc_soa/CMakeLists.txt
@@ -0,0 +1,13 @@
+# Algebra plugins library, part of the ACTS project (R&D line)
+#
+# (c) 2021-2022 CERN for the benefit of the ACTS project
+#
+# Mozilla Public License Version 2.0
+
+# Set up the library.
+algebra_add_library( algebra_vc_soa_math vc_soa_math
+   "include/algebra/math/vc_soa.hpp"
+   "include/algebra/math/impl/vc_soa_getter.hpp"
+   "include/algebra/math/impl/vc_soa_vector.hpp" )
+target_link_libraries( algebra_vc_soa_math
+   INTERFACE algebra::common algebra::common_math algebra::common_storage Vc::Vc )
diff --git a/math/vc_soa/include/algebra/math/impl/vc_soa_getter.hpp b/math/vc_soa/include/algebra/math/impl/vc_soa_getter.hpp
new file mode 100644
index 00000000..53a88835
--- /dev/null
+++ b/math/vc_soa/include/algebra/math/impl/vc_soa_getter.hpp
@@ -0,0 +1,122 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2022-2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+#include "algebra/qualifiers.hpp"
+#include "algebra/storage/vector.hpp"
+
+// Vc include(s).
+#ifdef _MSC_VER
+#pragma warning(push, 0)
+#endif  // MSVC
+#include <Vc/Vc>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif  // MSVC
+
+// System include(s)
+#include <cmath>
+#include <type_traits>
+
+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 array_t array type that holds the vector elements
+///
+/// @param v the input vector
+template <std::size_t N, typename scalar_t,
+          template <typename, std::size_t> class array_t,
+          std::enable_if_t<N >= 2, bool> = true>
+ALGEBRA_HOST_DEVICE inline auto phi(
+    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) noexcept {
+
+  return Vc::atan2(v[1], v[0]);
+}
+
+/// This method retrieves the perpenticular magnitude of a vector with rows >= 2
+///
+/// @tparam N dimension of the vector
+/// @tparam scalar_t scalar type
+/// @tparam array_t array type that holds the vector elements
+///
+/// @param v the input vector
+template <std::size_t N, typename scalar_t,
+          template <typename, std::size_t> class array_t,
+          std::enable_if_t<N >= 2, bool> = true>
+ALGEBRA_HOST_DEVICE inline auto perp(
+    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) noexcept {
+
+  auto tmp = Vc::Vector<scalar_t>::Zero();
+  tmp = Vc::fma(v[0], v[0], tmp);
+  tmp = Vc::fma(v[1], v[1], tmp);
+
+  return Vc::sqrt(tmp);
+}
+
+/// This method retrieves theta from a vector, vector base with rows >= 3
+///
+/// @tparam N dimension of the vector
+/// @tparam scalar_t scalar type
+/// @tparam array_t array type that holds the vector elements
+///
+/// @param v the input vector
+template <std::size_t N, typename scalar_t,
+          template <typename, std::size_t> class array_t,
+          std::enable_if_t<N >= 3, bool> = true>
+ALGEBRA_HOST_DEVICE inline auto theta(
+    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) noexcept {
+
+  return Vc::atan2(perp(v), v[2]);
+}
+
+/// This method retrieves the norm of a vector, no dimension restriction
+///
+/// @tparam N dimension of the vector
+/// @tparam scalar_t scalar type
+/// @tparam array_t array type that holds the vector elements
+///
+/// @param v the input vector
+template <std::size_t N, typename scalar_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) {
+
+  return Vc::sqrt(dot(v, v));
+}
+
+/// This method retrieves the pseudo-rapidity from a vector or vector base with
+/// rows >= 3
+///
+/// @tparam N dimension of the vector
+/// @tparam scalar_t scalar type
+/// @tparam array_t array type that holds the vector elements
+///
+/// @param v the input vector
+template <std::size_t N, typename scalar_t,
+          template <typename, std::size_t> class array_t,
+          std::enable_if_t<N >= 3, bool> = true>
+ALGEBRA_HOST_DEVICE inline auto eta(
+    const storage::vector<N, Vc::Vector<scalar_t>, array_t> &v) noexcept {
+
+  // atanh does not exist in Vc
+  auto atanh_func = [](scalar_t e) { return std::atanh(e); };
+
+  return (v[2] / norm(v)).apply(atanh_func);
+
+  // Faster, but less accurate
+  // return (Vc::reciprocal(norm(v)) * v[2]).apply(atanh_func);
+
+  // Even faster, but even less accurate
+  // return (Vc::rsqrt(dot(v, v)) * v[2]).apply(atanh_func);
+}
+
+}  // namespace algebra::vc_soa::math
\ No newline at end of file
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
new file mode 100644
index 00000000..c81b80d9
--- /dev/null
+++ b/math/vc_soa/include/algebra/math/impl/vc_soa_vector.hpp
@@ -0,0 +1,93 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+#include "algebra/qualifiers.hpp"
+#include "algebra/storage/vector.hpp"
+
+// Vc include(s).
+#ifdef _MSC_VER
+#pragma warning(push, 0)
+#endif  // MSVC
+#include <Vc/Vc>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif  // MSVC
+
+// System include(s)
+#include <cmath>
+
+namespace algebra::vc_soa::math {
+
+/// Cross product between two input vectors - 3 Dim
+///
+/// @tparam N dimension of the vector
+/// @tparam scalar_t scalar type
+/// @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 <typename, std::size_t> class array_t,
+          std::enable_if_t<N == 3, bool> = true>
+ALGEBRA_HOST_DEVICE inline storage::vector<N, Vc::Vector<scalar_t>, array_t>
+cross(const storage::vector<N, Vc::Vector<scalar_t>, array_t> &a,
+      const storage::vector<N, Vc::Vector<scalar_t>, array_t> &b) {
+
+  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])};
+}
+
+/// Dot product between two input vectors
+///
+/// @tparam N dimension of the vector
+/// @tparam scalar_t scalar type
+/// @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 <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) {
+
+  auto ret = Vc::Vector<scalar_t>::Zero();
+  for (unsigned int i{0u}; i < N; i++) {
+    ret = Vc::fma(a[i], b[i], ret);
+  }
+  return ret;
+}
+
+/// Get a normalized version of the input vector
+///
+/// @tparam N dimension of the vector
+/// @tparam scalar_t scalar type
+/// @tparam array_t array type that holds the vector elements
+///
+/// @param v the input vector
+template <std::size_t N, typename scalar_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) {
+
+  return (Vc::Vector<scalar_t>::One() / Vc::sqrt(dot(v, v))) * v;
+
+  // Less accurate, but faster
+  // return Vc::reciprocal(Vc::sqrt(dot(v, v))) * v;
+
+  // Even less accurate, but even faster
+  // return Vc::rsqrt(dot(v, v)) * v;
+}
+
+}  // namespace algebra::vc_soa::math
\ No newline at end of file
diff --git a/math/vc_soa/include/algebra/math/vc_soa.hpp b/math/vc_soa/include/algebra/math/vc_soa.hpp
new file mode 100644
index 00000000..df6c689a
--- /dev/null
+++ b/math/vc_soa/include/algebra/math/vc_soa.hpp
@@ -0,0 +1,14 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+// clang-format off
+#include "algebra/math/impl/vc_soa_vector.hpp"
+#include "algebra/math/impl/vc_soa_getter.hpp"
+// clang-format on
\ No newline at end of file
diff --git a/storage/CMakeLists.txt b/storage/CMakeLists.txt
index affd4c48..cc170aa2 100644
--- a/storage/CMakeLists.txt
+++ b/storage/CMakeLists.txt
@@ -6,6 +6,7 @@
 
 # Set up all enabled libraries.
 add_subdirectory( array )
+add_subdirectory( common )
 if( ALGEBRA_PLUGINS_INCLUDE_EIGEN )
    add_subdirectory( eigen )
 endif()
diff --git a/storage/common/CMakeLists.txt b/storage/common/CMakeLists.txt
new file mode 100644
index 00000000..025caf8c
--- /dev/null
+++ b/storage/common/CMakeLists.txt
@@ -0,0 +1,12 @@
+# Algebra plugins library, part of the ACTS project (R&D line)
+#
+# (c) 2023 CERN for the benefit of the ACTS project
+#
+# Mozilla Public License Version 2.0
+
+# Set up the library.
+algebra_add_library(algebra_common_storage common_storage
+   "include/algebra/storage/array_operators.hpp"
+   "include/algebra/storage/vector.hpp")
+target_link_libraries(algebra_common_storage
+   INTERFACE algebra::common)
diff --git a/storage/common/include/algebra/storage/array_operators.hpp b/storage/common/include/algebra/storage/array_operators.hpp
new file mode 100644
index 00000000..792c850b
--- /dev/null
+++ b/storage/common/include/algebra/storage/array_operators.hpp
@@ -0,0 +1,107 @@
+/** Algebra plugins, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+#pragma once
+
+// Project include(s)
+#include "algebra/qualifiers.hpp"
+
+// System include(s).
+#include <cstddef>
+#include <type_traits>
+
+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>
+ALGEBRA_HOST_DEVICE inline constexpr array_t<scalar_t, N> operator*(
+    const array_t<scalar_t, N> &a, scalar_t s) noexcept {
+
+  array_t<scalar_t, N> result;
+  for (std::size_t i{0u}; i < N; ++i) {
+    result[i] = a[i] * s;
+  }
+
+  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>
+ALGEBRA_HOST_DEVICE inline constexpr array_t<scalar_t, N> operator*(
+    scalar_t s, const array_t<scalar_t, N> &a) noexcept {
+
+  array_t<scalar_t, N> result;
+  for (std::size_t i{0u}; i < N; ++i) {
+    result[i] = s * a[i];
+  }
+
+  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>
+ALGEBRA_HOST_DEVICE inline constexpr array_t<scalar_t, N> operator-(
+    const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
+
+  array_t<scalar_t, N> result;
+  for (std::size_t i{0u}; i < N; ++i) {
+    result[i] = a[i] - b[i];
+  }
+
+  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>
+ALGEBRA_HOST_DEVICE inline constexpr array_t<scalar_t, N> operator+(
+    const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
+
+  array_t<scalar_t, N> result;
+  for (std::size_t i{0u}; i < N; ++i) {
+    result[i] = a[i] + b[i];
+  }
+
+  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>
+ALGEBRA_HOST_DEVICE inline constexpr array_t<scalar_t, N> operator*(
+    const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
+
+  array_t<scalar_t, N> result;
+  for (std::size_t i{0u}; i < N; ++i) {
+    result[i] = a[i] * b[i];
+  }
+
+  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>
+ALGEBRA_HOST_DEVICE inline constexpr array_t<scalar_t, N> operator/(
+    const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) noexcept {
+
+  array_t<scalar_t, N> result;
+  for (std::size_t i{0u}; i < N; ++i) {
+    result[i] = a[i] / b[i];
+  }
+
+  return result;
+}
+/// @}
+
+}  // namespace algebra::storage
diff --git a/storage/common/include/algebra/storage/vector.hpp b/storage/common/include/algebra/storage/vector.hpp
new file mode 100644
index 00000000..f3b31aca
--- /dev/null
+++ b/storage/common/include/algebra/storage/vector.hpp
@@ -0,0 +1,243 @@
+/** Algebra plugins, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+#pragma once
+
+// Project include(s)
+#include "algebra/qualifiers.hpp"
+#include "algebra/storage/array_operators.hpp"
+
+// System include(s).
+#include <cstddef>
+#include <initializer_list>
+#include <type_traits>
+#include <utility>
+
+namespace algebra::storage {
+
+/// Vector wrapper for AoS vs interleaved SoA data. @c value_t can e.g. be a
+/// SIMD vector.
+template <std::size_t N, typename value_t,
+          template <typename, std::size_t> class array_t>
+class vector {
+ public:
+  // Value type is a simd vector in SoA and a scalar in AoS
+  using value_type = value_t;
+  /// Underlying data array type
+  using array_type = array_t<value_t, N>;
+
+  /// Default contructor sets all entries to zero.
+  constexpr vector() { zero_fill(std::make_index_sequence<N>{}); }
+
+  /// Construct from element values @param vals .
+  ///
+  /// 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>
+  constexpr vector(Values &&...vals) : m_data{vals...} {
+    if constexpr ((sizeof...(Values) < N) &&
+                  (!std::conjunction_v<std::is_same<array_type, Values>...>)) {
+      zero_fill(std::make_index_sequence<N - sizeof...(Values)>{});
+    }
+  }
+
+  /// Construct from existing array storage @param vals .
+  constexpr vector(const array_type &vals) : m_data{vals} {}
+
+  /// Assignment operator from wrapped data.
+  ///
+  /// @param lhs wrap a copy of this data.
+  constexpr const vector &operator=(const array_type &lhs) {
+    m_data = lhs;
+    return *this;
+  }
+
+  /// Assignment operator from @c std::initializer_list .
+  ///
+  /// @param list wrap an array of this data.
+  constexpr vector &operator=(std::initializer_list<value_type> &list) {
+    m_data = array_type(list);
+    return *this;
+  }
+
+  /// Conversion operator from wrapper to underlying data array.
+  /// @{
+  constexpr operator array_type &() { return m_data; }
+  constexpr operator const array_type &() const { return m_data; }
+  /// @}
+
+  /// Subscript operator[]
+  /// @{
+  constexpr decltype(auto) operator[](std::size_t i) { return m_data[i]; }
+  constexpr decltype(auto) operator[](std::size_t i) const { return m_data[i]; }
+  /// @}
+
+  /// Operator*=
+  ///
+  /// @return Vector expression/return type according to the operation.
+  constexpr decltype(auto) operator*=(value_type factor) noexcept {
+    return m_data *= factor;
+  }
+
+  /// Equality operators
+  /// @{
+  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==(
+      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;
+  /// @}
+
+  /// Inequality operator
+  template <typename other_type>
+  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 {
+    using result_t = decltype(m_data[0] == rhs[0]);
+
+    array_t<result_t, N> comp;
+
+    for (unsigned int i{0u}; i < N; ++i) {
+      comp[i] = (m_data[i] == rhs[i]);
+    }
+
+    return comp;
+  }
+
+  /// Holds the data value for every vector element
+  array_t<value_t, N> m_data;
+
+ private:
+  /// Sets the trailing uninitialized values to zero.
+  template <std::size_t... Is>
+  constexpr void zero_fill(std::index_sequence<Is...>) noexcept {
+    if constexpr (sizeof...(Is) > 0) {
+      //((m_data[N - sizeof...(Is) + Is] = value_t(0)), ...);
+    }
+  }
+};
+
+/// Friend operators
+/// @{
+
+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,
+          std::enable_if_t<std::is_scalar_v<value_t>, bool> = true>
+constexpr bool operator==(const vector<N, value_t, array_t> &lhs,
+                          const o_array_t<value_t, N> &rhs) noexcept {
+
+  const auto comp = lhs.compare(rhs);
+  bool is_full = false;
+
+  for (unsigned int i{0u}; i < N; ++i) {
+    is_full |= comp[i];
+  }
+
+  return is_full;
+}
+
+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,
+          std::enable_if_t<!std::is_scalar_v<value_t>, bool> = true>
+constexpr bool operator==(const vector<N, value_t, array_t> &lhs,
+                          const o_array_t<value_t, N> &rhs) noexcept {
+
+  const auto comp = lhs.compare(rhs);
+  bool is_full = false;
+
+  for (unsigned int i{0u}; i < N; ++i) {
+    // Ducktyping the Vc::Vector::MaskType
+    is_full |= comp[i].isFull();
+  }
+
+  return is_full;
+}
+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 {
+  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,                   \
+            std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>        \
+  inline 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,                   \
+            std::enable_if_t<std::is_scalar_v<scalar_t>, bool> = true>        \
+  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>                   \
+  inline 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,    \
+      std::enable_if_t<                                                       \
+          std::is_object<decltype(std::declval<typename vector<               \
+                                      N, value_t, array_t>::array_type>()     \
+                                      OP std::declval<other_type>())>::value, \
+          bool> = true>                                                       \
+  inline 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,    \
+      std::enable_if_t<                                                       \
+          std::is_object<decltype(std::declval<typename vector<               \
+                                      N, value_t, array_t>::array_type>()     \
+                                      OP std::declval<other_type>())>::value, \
+          bool> = true>                                                       \
+  inline 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.
+// clang-format off
+DECLARE_vector_OPERATORS(+)
+DECLARE_vector_OPERATORS(-)
+DECLARE_vector_OPERATORS(*)
+DECLARE_vector_OPERATORS(/)
+// clang-format on
+
+// Clean up.
+#undef DECLARE_vector_OPERATORS
+
+}  // namespace algebra::storage
\ No newline at end of file
diff --git a/storage/vc/CMakeLists.txt b/storage/vc/CMakeLists.txt
index b9761199..070d1fc2 100644
--- a/storage/vc/CMakeLists.txt
+++ b/storage/vc/CMakeLists.txt
@@ -7,6 +7,7 @@
 # Set up the library.
 algebra_add_library( algebra_vc_storage vc_storage
    "include/algebra/storage/vc.hpp"
+   "include/algebra/storage/vc_soa.hpp"
    "include/algebra/storage/impl/vc_array4_wrapper.hpp" )
 target_link_libraries( algebra_vc_storage
-   INTERFACE algebra::common Vc::Vc )
+   INTERFACE algebra::common algebra::common_storage Vc::Vc )
diff --git a/storage/vc/include/algebra/storage/impl/vc_array4.hpp b/storage/vc/include/algebra/storage/impl/vc_array4.hpp
index 014c2524..99de0421 100644
--- a/storage/vc/include/algebra/storage/impl/vc_array4.hpp
+++ b/storage/vc/include/algebra/storage/impl/vc_array4.hpp
@@ -25,6 +25,7 @@ namespace algebra::vc {
  * size of the array.
  *
  * @tparam scalar_t scalar type used in the SimdArray
+ * @todo replace by @c algebra::storage::vector
  */
 template <typename scalar_t>
 struct array4 {
diff --git a/storage/vc/include/algebra/storage/vc.hpp b/storage/vc/include/algebra/storage/vc.hpp
index a1dc84f1..1247ea20 100644
--- a/storage/vc/include/algebra/storage/vc.hpp
+++ b/storage/vc/include/algebra/storage/vc.hpp
@@ -1,6 +1,6 @@
 /** Algebra plugins library, part of the ACTS project
  *
- * (c) 2020-2022 CERN for the benefit of the ACTS project
+ * (c) 2020-2023 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -8,7 +8,8 @@
 #pragma once
 
 // Project include(s).
-#include "algebra/storage/impl/vc_array4.hpp"
+//#include "algebra/storage/impl/vc_array4.hpp"
+#include "algebra/storage/vector.hpp"
 
 // System include(s).
 #include <array>
@@ -25,26 +26,41 @@
 
 namespace algebra::vc {
 
+namespace detail {
+
+/// Hide some template aruments of @c Vc::SimdArray behind type alias
+template <typename T, std::size_t N>
+using simd_array = Vc::SimdArray<T, N>;
+
+}  // namespace detail
+
 /// size type for Vc storage model
 using size_type = std::size_t;
 /// Array type used in the Vc storage model
 template <typename T, size_type N>
-using storage_type = Vc::SimdArray<T, N>;
+using storage_type = detail::simd_array<T, N>;
+/// value type in a linear algebra vector: AoS layout (use scalar directly)
+template <typename T>
+using value_type = T;
+/// Vector type used in the Vc AoS storage model
+template <typename T, std::size_t N>
+using vector_type = storage::vector<N, value_type<T>, detail::simd_array>;
 /// Matrix type used in the Vc storage model
 template <typename T, size_type ROWS, size_type COLS>
 using matrix_type = Vc::array<Vc::array<T, ROWS>, COLS>;
 
-/// 3-element "vector" type, using @c algebra::vc::array4
+/// @todo: Switch back to vectorizing types
+/// 3-element "vector" type, using @c std::array with gather/scatter capability
 template <typename T>
 using vector3 = Vc::array<T, 3>;
-/// Point in 3D space, using @c algebra::vc::array4
+/// Point in 3D space, using @c std::array with gather/scatter capability
 template <typename T>
 using point3 = vector3<T>;
-/// 2-element "vector" type, using @c std::array
+/// 2-element "vector" type, using @c std::array with gather/scatter capability
 template <typename T>
 using vector2 = Vc::array<T, 2>;
-/// Point in 2D space, using @c std::array
+/// Point in 2D space, using @c std::array with gather/scatter capability
 template <typename T>
 using point2 = vector2<T>;
 
-}  // namespace algebra::vc
+}  // namespace algebra::vc
\ No newline at end of file
diff --git a/storage/vc/include/algebra/storage/vc_soa.hpp b/storage/vc/include/algebra/storage/vc_soa.hpp
new file mode 100644
index 00000000..ffc44fc3
--- /dev/null
+++ b/storage/vc/include/algebra/storage/vc_soa.hpp
@@ -0,0 +1,62 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+#include "algebra/storage/vector.hpp"
+
+// System include(s).
+#include <array>
+#include <cstddef>
+
+// Vc include(s).
+#ifdef _MSC_VER
+#pragma warning(push, 0)
+#endif  // MSVC
+#include <Vc/Vc>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif  // MSVC
+
+namespace algebra::vc_soa {
+
+/// 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>
+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>;
+/// Vector type used in the Vc SoA storage model
+template <typename T, std::size_t N>
+using vector_type = storage::vector<N, value_type<T>, storage_type>;
+/// Matrix type used in the Vc SoA storage model
+template <typename T, size_type ROWS, size_type COLS>
+using matrix_type = storage::vector<ROWS * COLS, value_type<T>, storage_type>;
+
+/// 2-element "vector" type, using @c Vc::Vector in every element
+template <typename T>
+using vector2 = vector_type<T, 2>;
+/// Point in 2D space, using @c Vc::Vector in every element
+template <typename T>
+using point2 = vector2<T>;
+/// 3-element "vector" type, using @c Vc::Vector in every element
+template <typename T>
+using vector3 = vector_type<T, 3>;
+/// Point in 3D space, using @c Vc::Vector in every element
+template <typename T>
+using point3 = vector3<T>;
+/// 6-element "vector" type, using @c Vc::Vector in every element
+template <typename T>
+using vector6 = vector_type<T, 6>;
+/// 8-element "vector" type, using @c Vc::Vector in every element
+template <typename T>
+using vector8 = vector_type<T, 8>;
+
+}  // namespace algebra::vc_soa
\ No newline at end of file
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 05d7d405..ba2b7d41 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -39,9 +39,10 @@ endif()
 
 if( ALGEBRA_PLUGINS_INCLUDE_VC )
    algebra_add_test( vc
-      "vc/vc_cmath.cpp" "vc/vc_vc.cpp"
+      "vc/vc_cmath.cpp" "vc/vc_vc.cpp" "vc/vc_soa.cpp"
       LINK_LIBRARIES GTest::gtest_main algebra::tests_common
-                     algebra::vc_cmath algebra::vc_vc )
+                      algebra::vc_storage algebra::vc_cmath algebra::vc_vc
+                      algebra::vc_soa )
 endif()
 
 if( ALGEBRA_PLUGINS_INCLUDE_FASTOR )
diff --git a/tests/vc/vc_soa.cpp b/tests/vc/vc_soa.cpp
new file mode 100644
index 00000000..ac184029
--- /dev/null
+++ b/tests/vc/vc_soa.cpp
@@ -0,0 +1,141 @@
+/** Algebra plugins library, part of the ACTS project
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.f
+ */
+
+// Project include(s).
+#include "algebra/vc_soa.hpp"
+
+// GoogleTest include(s).
+#include <gtest/gtest.h>
+
+#include <array>
+#include <iostream>
+
+using namespace algebra;
+
+using scalar_t = float;
+
+constexpr float tol{1e-6f};
+
+/// 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>;
+  // Value type is Vc::Vector<float>
+  using value_t = typename vector3_v::value_type;
+
+  vector3_v a{value_t::One(), 2.f * value_t::One(), 3.f * value_t::One()};
+  vector3_v b{4.f * value_t::One(), 5.f * value_t::One(), 6.f * value_t::One()};
+
+  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());
+
+  // Masked comparison
+  auto m = a.compare(a);
+  EXPECT_TRUE(m[0].isFull());
+  EXPECT_TRUE(m[1].isFull());
+  EXPECT_TRUE(m[2].isFull());
+
+  m = a.compare(b);
+  EXPECT_FALSE(m[0].isFull());
+  EXPECT_FALSE(m[1].isFull());
+  EXPECT_FALSE(m[2].isFull());
+
+  // Full comparisons
+  EXPECT_TRUE(a == a);
+  EXPECT_FALSE(a == b);
+
+  // 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());
+
+  // 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());
+
+  // 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());
+
+  // 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());
+
+  // 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());
+
+  // 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());
+
+  auto d{vector::dot(a, b)};
+  EXPECT_TRUE((d == value_t(32.f)).isFull());
+
+  value_t norms_a{getter::norm(vector::normalize(a))};
+  value_t norms_b{getter::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());
+
+  static_assert(std::is_convertible_v<decltype(v_expr), vector3_v>,
+                "expression type not convertible");
+}
+
+/// 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>;
+  // Value type is Vc::Vector<float>
+  using value_t = typename vector3_v::value_type;
+
+  vector3_v a{value_t::One(), 2.f * value_t::One(), 3.f * value_t::One()};
+
+  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());
+
+  // All results in the vector are the same, so only check the first one
+
+  // Phi angle
+  auto v_phi = getter::phi(a);
+  EXPECT_NEAR(v_phi[0], static_cast<scalar_t>(std::atan2(2., 1.)), tol);
+
+  // Perpendicular projection
+  auto v_perp = getter::perp(a);
+  EXPECT_NEAR(v_perp[0], std::sqrt(5.), tol);
+
+  // Theta angle
+  auto v_theta = getter::theta(a);
+  EXPECT_NEAR(v_theta[0], static_cast<scalar_t>(std::atan2(std::sqrt(5.), 3.)), tol);
+
+  // Norm of the vector
+  auto v_norm = getter::norm(a);
+  EXPECT_NEAR(v_norm[0], std::sqrt(14.), tol);
+
+  // Eta of the vector
+  auto v_eta = getter::eta(a);
+  EXPECT_NEAR(v_eta[0], static_cast<scalar_t>(std::atanh(1. / std::sqrt(14.) * 3.)), tol);
+}
\ No newline at end of file