Add new generic matrix operators

Our current approach for e.g. matrix multiplication reads very naturally, e.g. `A = B*C` models $A = BC$, but this has a problem on GPU code. Indeed, if these matrices have size $N \times N$, then we first concretize an $N \times N$ matrix which we then have to copy element-by-element into $A$. This means that we need to keep that many registers live, which is quite large for e.g. $7 \times 7$ free matrices (which thus occupy 49 registers). This problem can be alleviated by implementing optimized operators. More precisely, this PR implements the following new methods: * `set_product(A, B, C)` computes $A = BC$ without intermediate values. * `set_product_left_transpose(A, B, C)` computes $A = B^TC$ without intermediate values. * `set_product_right_transpose(A, B, C)` computes $A = BC^T` without intermediate values. * `set_inplace_product_right(A, B)` computes $A = AB$ in place. * `set_inplace_product_left(A, B)` computes $A = BA$ in place. * `set_inplace_product_right_transpose(A, B)` computes $A = AB^T$ in place. * `set_inplace_product_left_transpose(A, B)` computes $A = B^TA$ in place.
acts-project · Dec 3, 2024 · e486e2a · e486e2a
1 parent 0fb56ea
commit e486e2a
Show file tree

Hide file tree

Showing 13 changed files with 947 additions and 0 deletions.
diff --git a/common/include/algebra/concepts.hpp b/common/include/algebra/concepts.hpp
@@ -88,6 +88,41 @@ concept column_matrix = matrix<M> && (algebra::traits::columns<M> == 1);
 
 template <typename M>
 concept column_matrix3D = column_matrix<M> && (algebra::traits::rows<M> == 3);
+
+template <typename MA, typename MB>
+concept matrix_multipliable =
+    matrix<MA>&& matrix<MB> &&
+    (algebra::traits::columns<MA> ==
+     algebra::traits::rows<
+         MB>)&&std::convertible_to<algebra::traits::index_t<MA>,
+                                   algebra::traits::index_t<MB>>&& std::
+        convertible_to<algebra::traits::index_t<MB>,
+                       algebra::traits::index_t<
+                           MA>>&& requires(algebra::traits::scalar_t<MA> sa,
+                                           algebra::traits::scalar_t<MB> sb) {
+  {sa * sb};
+};
+
+template <typename MA, typename MB, typename MC>
+concept matrix_multipliable_into =
+    matrix_multipliable<MA, MB>&& matrix<MC> &&
+    (algebra::traits::rows<MC> == algebra::traits::rows<MA>)&&(
+        algebra::traits::columns<MC> ==
+        algebra::traits::columns<
+            MB>)&&std::convertible_to<algebra::traits::index_t<MA>,
+                                      algebra::traits::index_t<MC>>&& std::
+        convertible_to<algebra::traits::index_t<MC>,
+                       algebra::traits::index_t<MA>>&& std::
+            convertible_to<algebra::traits::index_t<MB>,
+                           algebra::traits::index_t<MC>>&& std::
+                convertible_to<
+                    algebra::traits::index_t<MC>,
+                    algebra::traits::index_t<
+                        MB>>&& requires(algebra::traits::scalar_t<MA> sa,
+                                        algebra::traits::scalar_t<MB> sb,
+                                        algebra::traits::scalar_t<MC>& sc) {
+  {sc += (sa * sb)};
+};
 /// @}
 
 /// Transform concept

diff --git a/frontend/array_cmath/include/algebra/array_cmath.hpp b/frontend/array_cmath/include/algebra/array_cmath.hpp
@@ -93,6 +93,14 @@ using cmath::determinant;
 using cmath::inverse;
 using cmath::transpose;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix

diff --git a/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp b/frontend/eigen_eigen/include/algebra/eigen_eigen.hpp
@@ -9,6 +9,7 @@
 
 // Project include(s).
 #include "algebra/math/eigen.hpp"
+#include "algebra/math/generic.hpp"
 #include "algebra/storage/eigen.hpp"
 
 // Eigen include(s).
@@ -65,6 +66,14 @@ using eigen::math::set_zero;
 using eigen::math::transpose;
 using eigen::math::zero;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix

diff --git a/frontend/eigen_generic/include/algebra/eigen_generic.hpp b/frontend/eigen_generic/include/algebra/eigen_generic.hpp
@@ -71,6 +71,14 @@ using generic::math::set_zero;
 using generic::math::transpose;
 using generic::math::zero;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix

diff --git a/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp b/frontend/fastor_fastor/include/algebra/fastor_fastor.hpp
@@ -9,6 +9,7 @@
 
 // Project include(s).
 #include "algebra/math/fastor.hpp"
+#include "algebra/math/generic.hpp"
 #include "algebra/storage/fastor.hpp"
 
 // Fastor include(s).
@@ -67,6 +68,14 @@ using fastor::math::set_zero;
 using fastor::math::transpose;
 using fastor::math::zero;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix

diff --git a/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp b/frontend/smatrix_generic/include/algebra/smatrix_generic.hpp
@@ -62,6 +62,14 @@ using generic::math::set_zero;
 using generic::math::transpose;
 using generic::math::zero;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix

diff --git a/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp b/frontend/smatrix_smatrix/include/algebra/smatrix_smatrix.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/math/generic.hpp"
 #include "algebra/math/smatrix.hpp"
 #include "algebra/storage/smatrix.hpp"
 
@@ -58,6 +59,14 @@ using smatrix::math::set_zero;
 using smatrix::math::transpose;
 using smatrix::math::zero;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix

diff --git a/frontend/vc_aos/include/algebra/vc_aos.hpp b/frontend/vc_aos/include/algebra/vc_aos.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/math/generic.hpp"
 #include "algebra/math/vc_aos.hpp"
 #include "algebra/storage/vc_aos.hpp"
 
@@ -63,6 +64,14 @@ using vc_aos::math::set_zero;
 using vc_aos::math::transpose;
 using vc_aos::math::zero;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix

diff --git a/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp b/frontend/vc_aos_generic/include/algebra/vc_aos_generic.hpp
@@ -53,6 +53,14 @@ using generic::math::perp;
 using generic::math::phi;
 using generic::math::theta;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace vector

diff --git a/frontend/vc_soa/include/algebra/vc_soa.hpp b/frontend/vc_soa/include/algebra/vc_soa.hpp
@@ -8,6 +8,7 @@
 #pragma once
 
 // Project include(s).
+#include "algebra/math/generic.hpp"
 #include "algebra/math/impl/vc_aos_transform3.hpp"
 #include "algebra/math/vc_soa.hpp"
 #include "algebra/storage/vc_soa.hpp"
@@ -71,6 +72,14 @@ using vc_soa::math::set_zero;
 using vc_soa::math::transpose;
 using vc_soa::math::zero;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 }  // namespace matrix
 
 namespace vc_soa {

diff --git a/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp b/frontend/vecmem_cmath/include/algebra/vecmem_cmath.hpp
@@ -92,6 +92,14 @@ using generic::math::determinant;
 using generic::math::inverse;
 using generic::math::transpose;
 
+using generic::math::set_inplace_product_left;
+using generic::math::set_inplace_product_left_transpose;
+using generic::math::set_inplace_product_right;
+using generic::math::set_inplace_product_right_transpose;
+using generic::math::set_product;
+using generic::math::set_product_left_transpose;
+using generic::math::set_product_right_transpose;
+
 /// @}
 
 }  // namespace matrix