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 4, 2024 · 05fd7a7 · 05fd7a7
1 parent 0fb56ea
commit 05fd7a7
Show file tree

Hide file tree

Showing 19 changed files with 565 additions and 7 deletions.
diff --git a/cmake/algebra-plugins-compiler-options-cpp.cmake b/cmake/algebra-plugins-compiler-options-cpp.cmake
@@ -28,7 +28,7 @@ elseif( "${CMAKE_CXX_COMPILER_ID}" MATCHES "MSVC" )
    # Basic flags for all build modes.
    string( REGEX REPLACE "/W[0-9]" "" CMAKE_CXX_FLAGS
       "${CMAKE_CXX_FLAGS}" )
-   algebra_add_flag( CMAKE_CXX_FLAGS "/W4" )
+   algebra_add_flag( CMAKE_CXX_FLAGS "/W4 /bigobj" )
 
    # Fail on warnings, if asked for that behaviour.
    if( ALGEBRA_PLUGINS_FAIL_ON_WARNINGS )

diff --git a/common/include/algebra/concepts.hpp b/common/include/algebra/concepts.hpp
@@ -88,6 +88,34 @@ 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_compatible = matrix<MA>&& matrix<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>>;
+
+template <typename MA, typename MB>
+concept matrix_multipliable =
+    matrix_compatible<MA, MB> &&
+    (algebra::traits::columns<MA> ==
+     algebra::traits::rows<MB>)&&requires(algebra::traits::scalar_t<MA> sa,
+                                          algebra::traits::scalar_t<MB> sb) {
+  {(sa * sb) + (sa * sb)};
+};
+
+template <typename MA, typename MB, typename MC>
+concept matrix_multipliable_into =
+    matrix_multipliable<MA, MB>&& matrix_compatible<MA, MC>&&
+        matrix_compatible<MB, MC> &&
+    (algebra::traits::rows<MC> == algebra::traits::rows<MA>)&&(
+        algebra::traits::columns<MC> ==
+        algebra::traits::columns<
+            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/CMakeLists.txt b/frontend/eigen_eigen/CMakeLists.txt
@@ -8,7 +8,7 @@
 algebra_add_library( algebra_eigen_eigen eigen_eigen
    "include/algebra/eigen_eigen.hpp" )
 target_link_libraries( algebra_eigen_eigen
-   INTERFACE algebra::common algebra::eigen_storage algebra::eigen_math
+   INTERFACE algebra::common algebra::eigen_storage algebra::eigen_math algebra::generic_math
              Eigen3::Eigen )
 algebra_test_public_headers( algebra_eigen_eigen
    "algebra/eigen_eigen.hpp" )
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/CMakeLists.txt b/frontend/fastor_fastor/CMakeLists.txt
@@ -8,6 +8,6 @@
 algebra_add_library( algebra_fastor_fastor fastor_fastor
    "include/algebra/fastor_fastor.hpp" )
 target_link_libraries( algebra_fastor_fastor
-   INTERFACE algebra::common algebra::fastor_storage algebra::fastor_math )
+   INTERFACE algebra::common algebra::fastor_storage algebra::fastor_math algebra::generic_math )
 algebra_test_public_headers( algebra_fastor_fastor
    "algebra/fastor_fastor.hpp" )
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/CMakeLists.txt b/frontend/smatrix_smatrix/CMakeLists.txt
@@ -8,6 +8,6 @@
 algebra_add_library( algebra_smatrix_smatrix smatrix_smatrix
    "include/algebra/smatrix_smatrix.hpp" )
 target_link_libraries( algebra_smatrix_smatrix
-   INTERFACE algebra::common algebra::smatrix_storage algebra::smatrix_math )
+   INTERFACE algebra::common algebra::smatrix_storage algebra::smatrix_math algebra::generic_math )
 algebra_test_public_headers( algebra_smatrix_smatrix
    "algebra/smatrix_smatrix.hpp" )
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/CMakeLists.txt b/frontend/vc_aos/CMakeLists.txt
@@ -8,6 +8,6 @@
 algebra_add_library( algebra_vc_aos vc_aos
    "include/algebra/vc_aos.hpp" )
 target_link_libraries( algebra_vc_aos
-   INTERFACE algebra::common algebra::vc_aos_storage algebra::vc_aos_math )
+   INTERFACE algebra::common algebra::vc_aos_storage algebra::vc_aos_math algebra::generic_math )
 algebra_test_public_headers( algebra_vc_aos
    "algebra/vc_aos.hpp" )
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
@@ -70,6 +70,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/vc_soa/CMakeLists.txt b/frontend/vc_soa/CMakeLists.txt
@@ -8,7 +8,7 @@
 algebra_add_library( algebra_vc_soa vc_soa
    "include/algebra/vc_soa.hpp" )
 target_link_libraries( algebra_vc_soa
-   INTERFACE algebra::common algebra::vc_soa_storage algebra::vc_soa_math
+   INTERFACE algebra::common algebra::vc_soa_storage algebra::vc_soa_math algebra::generic_math
              algebra::vc_aos_math )
 algebra_test_public_headers( algebra_vc_soa
    "algebra/vc_soa.hpp" )
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