Fix getters

acts-project · Nov 21, 2024 · 3d95d67 · 3d95d67
1 parent b9286d3
commit 3d95d67
Show file tree

Hide file tree

Showing 12 changed files with 204 additions and 213 deletions.
diff --git a/benchmarks/vc_soa/vc_soa_vector.cpp b/benchmarks/vc_soa/vc_soa_vector.cpp
@@ -70,17 +70,6 @@ int main(int argc, char** argv) {
   //
   // Register all benchmarks
   //
-  algebra::register_benchmark<phi_f_t>(cfg_s, "_single");
-  algebra::register_benchmark<phi_d_t>(cfg_d, "_double");
-  algebra::register_benchmark<theta_f_t>(cfg_s, "_single");
-  algebra::register_benchmark<theta_d_t>(cfg_d, "_double");
-  algebra::register_benchmark<perp_f_t>(cfg_s, "_single");
-  algebra::register_benchmark<perp_d_t>(cfg_d, "_double");
-  algebra::register_benchmark<norm_f_t>(cfg_s, "_single");
-  algebra::register_benchmark<norm_d_t>(cfg_d, "_double");
-  algebra::register_benchmark<eta_f_t>(cfg_s, "_single");
-  algebra::register_benchmark<eta_d_t>(cfg_d, "_double");
-
   algebra::register_benchmark<add_f_t>(cfg_s, "_single");
   algebra::register_benchmark<add_d_t>(cfg_d, "_double");
   algebra::register_benchmark<sub_f_t>(cfg_s, "_single");
@@ -92,6 +81,17 @@ int main(int argc, char** argv) {
   algebra::register_benchmark<normlz_f_t>(cfg_s, "_single");
   algebra::register_benchmark<normlz_d_t>(cfg_d, "_double");
 
+  algebra::register_benchmark<phi_f_t>(cfg_s, "_single");
+  algebra::register_benchmark<phi_d_t>(cfg_d, "_double");
+  algebra::register_benchmark<theta_f_t>(cfg_s, "_single");
+  algebra::register_benchmark<theta_d_t>(cfg_d, "_double");
+  algebra::register_benchmark<perp_f_t>(cfg_s, "_single");
+  algebra::register_benchmark<perp_d_t>(cfg_d, "_double");
+  algebra::register_benchmark<norm_f_t>(cfg_s, "_single");
+  algebra::register_benchmark<norm_d_t>(cfg_d, "_double");
+  algebra::register_benchmark<eta_f_t>(cfg_s, "_single");
+  algebra::register_benchmark<eta_d_t>(cfg_d, "_double");
+
   ::benchmark::Initialize(&argc, argv);
   ::benchmark::RunSpecifiedBenchmarks();
   ::benchmark::Shutdown();

diff --git a/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp b/math/vc_aos/include/algebra/math/impl/vc_aos_vector.hpp
@@ -11,6 +11,7 @@
 #include "algebra/concepts.hpp"
 #include "algebra/math/common.hpp"
 #include "algebra/qualifiers.hpp"
+#include "algebra/storage/vc_aos.hpp"
 #include "algebra/storage/vector.hpp"
 
 // Vc include(s).
@@ -30,44 +31,42 @@ namespace algebra::vc_aos::math {
 
 /// This method retrieves phi from a vector @param v
 template <typename vector_t>
-requires(Vc::is_simd_vector<vector_t>::value ||
-         algebra::detail::is_storage_vector_v<vector_t>) ALGEBRA_HOST_DEVICE
-    inline auto phi(const vector_t &v) {
+  requires(Vc::is_simd_vector<vector_t>::value ||
+           algebra::detail::is_storage_vector_v<vector_t>)
+ALGEBRA_HOST_DEVICE inline auto phi(const vector_t &v) {
   return algebra::math::atan2(v[1], v[0]);
 }
 
 /// This method retrieves the perpendicular magnitude of a vector @param v
 template <typename vector_t>
-requires(Vc::is_simd_vector<vector_t>::value ||
-         algebra::detail::is_storage_vector_v<vector_t>) ALGEBRA_HOST_DEVICE
-    inline auto perp(const vector_t &v) {
+  requires(Vc::is_simd_vector<vector_t>::value ||
+           algebra::detail::is_storage_vector_v<vector_t>)
+ALGEBRA_HOST_DEVICE inline auto perp(const vector_t &v) {
   return algebra::math::sqrt(algebra::math::fma(v[0], v[0], v[1] * v[1]));
 }
 
 /// This method retrieves theta from a vector @param v
 template <typename vector_t>
-requires(Vc::is_simd_vector<vector_t>::value ||
-         algebra::detail::is_storage_vector_v<vector_t>) ALGEBRA_HOST_DEVICE
-    inline auto theta(const vector_t &v) {
+  requires(Vc::is_simd_vector<vector_t>::value ||
+           algebra::detail::is_storage_vector_v<vector_t>)
+ALGEBRA_HOST_DEVICE inline auto theta(const vector_t &v) {
   return algebra::math::atan2(perp(v), v[2]);
 }
 
 /// Dot product between two input vectors
 ///
-/// @tparam vector_type generic input vector type
+/// @tparam vector_t generic input vector type
 ///
 /// @param a the first input vector
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
-template <typename vector_type1, typename vector_type2>
-requires(
-    (Vc::is_simd_vector<vector_type1>::value ||
-     algebra::detail::is_storage_vector_v<
-         vector_type1>)&&(Vc::is_simd_vector<vector_type2>::value ||
-                          algebra::detail::is_storage_vector_v<vector_type2>))
-    ALGEBRA_HOST_DEVICE
-    inline auto dot(const vector_type1 &a, const vector_type2 &b) {
+template <typename vector_t1, typename vector_t2>
+  requires((Vc::is_simd_vector<vector_t1>::value ||
+            algebra::detail::is_storage_vector_v<vector_t1>) &&
+           (Vc::is_simd_vector<vector_t2>::value ||
+            algebra::detail::is_storage_vector_v<vector_t2>))
+ALGEBRA_HOST_DEVICE inline auto dot(const vector_t1 &a, const vector_t2 &b) {
 
   return (a * b).sum();
 }
@@ -76,22 +75,22 @@ requires(
 ///
 /// @param v the input vector
 template <typename vector_t>
-requires(Vc::is_simd_vector<vector_t>::value ||
-         algebra::detail::is_storage_vector_v<vector_t>) ALGEBRA_HOST_DEVICE
-    inline auto norm(const vector_t &v) {
+  requires(Vc::is_simd_vector<vector_t>::value ||
+           algebra::detail::is_storage_vector_v<vector_t>)
+ALGEBRA_HOST_DEVICE inline auto norm(const vector_t &v) {
 
   return algebra::math::sqrt(dot(v, v));
 }
 
 /// Get a normalized version of the input vector
 ///
-/// @tparam vector_type generic input vector type
+/// @tparam vector_t generic input vector type
 ///
 /// @param v the input vector
-template <typename vector_type>
-requires(Vc::is_simd_vector<vector_type>::value ||
-         algebra::detail::is_storage_vector_v<vector_type>) ALGEBRA_HOST_DEVICE
-    inline auto normalize(const vector_type &v) {
+template <typename vector_t>
+  requires(Vc::is_simd_vector<vector_t>::value ||
+           algebra::detail::is_storage_vector_v<vector_t>)
+ALGEBRA_HOST_DEVICE inline vector_t normalize(const vector_t &v) {
 
   return v / norm(v);
 }
@@ -101,47 +100,47 @@ requires(Vc::is_simd_vector<vector_type>::value ||
 ///
 /// @param v the input vector
 template <typename vector_t>
-requires(Vc::is_simd_vector<vector_t>::value ||
-         algebra::detail::is_storage_vector_v<vector_t>) ALGEBRA_HOST_DEVICE
-    inline auto eta(const vector_t &v) noexcept {
+  requires(Vc::is_simd_vector<vector_t>::value ||
+           algebra::detail::is_storage_vector_v<vector_t>)
+ALGEBRA_HOST_DEVICE inline auto eta(const vector_t &v) noexcept {
 
   return algebra::math::atanh(v[2] / norm(v));
 }
 
 /// Cross product between two input vectors - 3 Dim
 ///
-/// @tparam vector_type generic input vector type
+/// @tparam vector_t generic input vector type
 ///
 /// @param a the first input vector
 /// @param b the second input vector
 ///
 /// @return a vector representing the cross product
-template <typename vector_type1, typename vector_type2>
-requires(
-    (Vc::is_simd_vector<vector_type1>::value ||
-     algebra::detail::is_storage_vector_v<
-         vector_type1>)&&(Vc::is_simd_vector<vector_type2>::value ||
-                          algebra::detail::is_storage_vector_v<vector_type2>))
-    ALGEBRA_HOST_DEVICE
-    inline auto cross(const vector_type1 &a, const vector_type2 &b)
-        -> decltype(a * b - a * b) {
-
-  return {algebra::math::fma(a[1], b[2], -b[1] * a[2]),
-          algebra::math::fma(a[2], b[0], -b[2] * a[0]),
-          algebra::math::fma(a[0], b[1], -b[0] * a[1]), 0.f};
+template <typename vector_t1, typename vector_t2>
+  requires((Vc::is_simd_vector<vector_t1>::value ||
+            algebra::detail::is_storage_vector_v<vector_t1>) &&
+           (Vc::is_simd_vector<vector_t2>::value ||
+            algebra::detail::is_storage_vector_v<vector_t2>))
+ALGEBRA_HOST_DEVICE inline auto cross(const vector_t1 &a, const vector_t2 &b) {
+
+  using scalar_t = decltype(algebra::math::fma(a[1], b[2], -b[1] * a[2]));
+
+  return vc_aos::vector_type<scalar_t, 3>{
+      algebra::math::fma(a[1], b[2], -b[1] * a[2]),
+      algebra::math::fma(a[2], b[0], -b[2] * a[0]),
+      algebra::math::fma(a[0], b[1], -b[0] * a[1])};
 }
 
 /// Elementwise sum
 ///
-/// @tparam vector_type generic input vector type
+/// @tparam vector_t generic input vector type
 ///
 /// @param v the vector whose elements should be summed
 ///
 /// @return the sum of the elements
-template <typename vector_type>
-requires(Vc::is_simd_vector<vector_type>::value ||
-         algebra::detail::is_storage_vector_v<vector_type>) ALGEBRA_HOST_DEVICE
-    inline auto sum(const vector_type &v) {
+template <typename vector_t>
+  requires(Vc::is_simd_vector<vector_t>::value ||
+           algebra::detail::is_storage_vector_v<vector_t>)
+ALGEBRA_HOST_DEVICE inline auto sum(const vector_t &v) {
   return v.get().sum();
 }
 

diff --git a/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp b/storage/cmath/include/algebra/storage/impl/cmath_getter.hpp
@@ -169,7 +169,7 @@ struct block_getter {
   template <std::size_t SIZE, std::size_t ROWS, std::size_t COLS,
             concepts::scalar scalar_t,
             template <typename, std::size_t> class array_t>
-  ALGEBRA_HOST_DEVICE inline array_t<scalar_t, SIZE> operator()(
+  ALGEBRA_HOST_DEVICE inline array_t<scalar_t, SIZE> vector(
       const array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
       std::size_t col) {
 
@@ -204,7 +204,7 @@ ALGEBRA_HOST_DEVICE inline array_t<scalar_t, SIZE> vector(
     const array_t<array_t<scalar_t, ROWS>, COLS> &m, std::size_t row,
     std::size_t col) {
 
-  return block_getter().template operator()<SIZE>(m, row, col);
+  return block_getter().template vector<SIZE>(m, row, col);
 }
 
 /// Sets a matrix of dimension @tparam ROW and @tparam COL as submatrix of

diff --git a/storage/common/include/algebra/storage/matrix.hpp b/storage/common/include/algebra/storage/matrix.hpp
@@ -32,9 +32,9 @@ struct alignas(alignof(storage::vector<ROW, scalar_t, array_t>)) matrix {
 
   /// Construct from given column vectors @param v
   template <concepts::vector... vector_t>
-  ALGEBRA_HOST_DEVICE requires(sizeof...(vector_t) ==
-                               COL) explicit matrix(vector_t &&... v)
-      : m_storage{std::forward<vector_t>(v)...} {}
+  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, concepts::scalar S,
@@ -76,17 +76,17 @@ struct alignas(alignof(storage::vector<ROW, scalar_t, array_t>)) matrix {
   /// @{
   // AoS
   template <std::size_t... I>
-  ALGEBRA_HOST_DEVICE requires(
-      !std::is_scalar_v<scalar_t>) constexpr bool equal(const matrix &rhs,
-                                                        std::index_sequence<
-                                                            I...>) const {
+  ALGEBRA_HOST_DEVICE
+    requires(!std::is_scalar_v<scalar_t>)
+  constexpr bool equal(const matrix &rhs, std::index_sequence<I...>) const {
     return (... && (m_storage[I] == rhs[I]));
   }
 
   // SoA
   template <std::size_t... I>
-  ALGEBRA_HOST requires(std::is_scalar_v<scalar_t>) constexpr bool equal(
-      const matrix &rhs, std::index_sequence<I...>) const {
+  ALGEBRA_HOST
+    requires(std::is_scalar_v<scalar_t>)
+  constexpr bool equal(const matrix &rhs, std::index_sequence<I...>) const {
     return (... && ((m_storage[I].get() == rhs[I].get()).isFull()));
   }
   /// @}
@@ -105,15 +105,15 @@ struct alignas(alignof(storage::vector<ROW, scalar_t, array_t>)) matrix {
 
   template <std::size_t R, std::size_t C, typename S1, typename S2,
             template <typename, std::size_t> class A>
-  requires(std::is_scalar_v<S2> || std::is_same_v<S1, S2>) ALGEBRA_HOST_DEVICE
-      friend constexpr decltype(auto)
-      operator*(S2 a, const matrix<A, S1, R, C> &rhs) noexcept;
+  /*requires(std::is_scalar_v<S2> || std::is_same_v<S1, S2>)*/
+  ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator*(
+      const S2 a, const matrix<A, S1, R, C> &rhs) noexcept;
 
-  template <std::size_t R, std::size_t C, typename S1, typename S2,
-            template <typename, std::size_t> class A>
-  requires(std::is_scalar_v<S2> || std::is_same_v<S1, S2>) ALGEBRA_HOST_DEVICE
-      friend constexpr decltype(auto)
-      operator*(const matrix<A, S1, R, C> &lhs, S2 a) noexcept;
+  template <std::size_t R, std::size_t C, concepts::scalar S1,
+            concepts::scalar S2, template <typename, std::size_t> class A>
+  /*requires(std::is_scalar_v<S2> || std::is_same_v<S1, S2>)*/
+  ALGEBRA_HOST_DEVICE friend constexpr decltype(auto) operator*(
+      const matrix<A, S1, R, C> &lhs, const S2 a) noexcept;
 
   /// Matrix-vector multiplication
   template <std::size_t R, std::size_t C, typename S,
@@ -212,8 +212,9 @@ template <concepts::matrix matrix_t, concepts::scalar scalar_t,
           std::size_t... J>
 ALGEBRA_HOST_DEVICE constexpr matrix_t matrix_scalar_mul(
     scalar_t a, const matrix_t &rhs, std::index_sequence<J...>) noexcept {
+  using mat_scalar_t = algebra::traits::value_t<matrix_t>;
 
-  return matrix_t{(a * rhs[J])...};
+  return matrix_t{(static_cast<mat_scalar_t>(a) * rhs[J])...};
 }
 
 /// Matrix addition
@@ -257,30 +258,22 @@ ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator-(
   return matrix_sub(lhs, rhs, std::make_index_sequence<matrix_t::columns()>());
 }
 
-template <std::size_t ROW, std::size_t COL, concepts::scalar scalar1_t,
-          concepts::scalar scalar2_t,
-          template <typename, std::size_t> class array_t>
-requires(std::is_scalar_v<scalar1_t> ||
-         std::is_same_v<scalar1_t, scalar2_t>) ALGEBRA_HOST_DEVICE
-    constexpr decltype(auto)
-    operator*(scalar1_t a,
-              const matrix<array_t, scalar2_t, ROW, COL> &rhs) noexcept {
+template <std::size_t R, std::size_t C, typename S1, typename S2,
+          template <typename, std::size_t> class A>
+ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
+    const S2 a, const matrix<A, S1, R, C> &rhs) noexcept {
 
-  using matrix_t = matrix<array_t, scalar2_t, ROW, COL>;
+  using matrix_t = matrix<A, S2, R, C>;
 
-  return matrix_scalar_mul(a, rhs,
+  return matrix_scalar_mul(static_cast<S1>(a), rhs,
                            std::make_index_sequence<matrix_t::columns()>());
 }
 
-template <std::size_t ROW, std::size_t COL, concepts::scalar scalar1_t,
-          concepts::scalar scalar2_t,
-          template <typename, std::size_t> class array_t>
-requires(std::is_scalar_v<scalar1_t> ||
-         std::is_same_v<scalar1_t, scalar2_t>) ALGEBRA_HOST_DEVICE
-    constexpr decltype(auto)
-    operator*(const matrix<array_t, scalar1_t, ROW, COL> &lhs,
-              scalar2_t a) noexcept {
-  return a * lhs;
+template <std::size_t R, std::size_t C, concepts::scalar S1,
+          concepts::scalar S2, template <typename, std::size_t> class A>
+ALGEBRA_HOST_DEVICE constexpr decltype(auto) operator*(
+    const matrix<A, S1, R, C> &lhs, const S2 a) noexcept {
+  return static_cast<S1>(a) * lhs;
 }
 
 /// Matrix-vector multiplication