From e486e2a8d6601c9d843ab49aee1d694a86adf174 Mon Sep 17 00:00:00 2001 From: Stephen Nicholas Swatman Date: Tue, 19 Nov 2024 12:14:26 +0100 Subject: [PATCH] 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. --- common/include/algebra/concepts.hpp | 35 + .../include/algebra/array_cmath.hpp | 8 + .../include/algebra/eigen_eigen.hpp | 9 + .../include/algebra/eigen_generic.hpp | 8 + .../include/algebra/fastor_fastor.hpp | 9 + .../include/algebra/smatrix_generic.hpp | 8 + .../include/algebra/smatrix_smatrix.hpp | 9 + frontend/vc_aos/include/algebra/vc_aos.hpp | 9 + .../include/algebra/vc_aos_generic.hpp | 8 + frontend/vc_soa/include/algebra/vc_soa.hpp | 9 + .../include/algebra/vecmem_cmath.hpp | 8 + .../algebra/math/impl/generic_matrix.hpp | 190 ++++++ tests/common/test_host_basics.hpp | 637 ++++++++++++++++++ 13 files changed, 947 insertions(+) diff --git a/common/include/algebra/concepts.hpp b/common/include/algebra/concepts.hpp index 251a1f46..c12d72d8 100644 --- a/common/include/algebra/concepts.hpp +++ b/common/include/algebra/concepts.hpp @@ -88,6 +88,41 @@ concept column_matrix = matrix && (algebra::traits::columns == 1); template concept column_matrix3D = column_matrix && (algebra::traits::rows == 3); + +template +concept matrix_multipliable = + matrix&& matrix && + (algebra::traits::columns == + algebra::traits::rows< + MB>)&&std::convertible_to, + algebra::traits::index_t>&& std:: + convertible_to, + algebra::traits::index_t< + MA>>&& requires(algebra::traits::scalar_t sa, + algebra::traits::scalar_t sb) { + {sa * sb}; +}; + +template +concept matrix_multipliable_into = + matrix_multipliable&& matrix && + (algebra::traits::rows == algebra::traits::rows)&&( + algebra::traits::columns == + algebra::traits::columns< + MB>)&&std::convertible_to, + algebra::traits::index_t>&& std:: + convertible_to, + algebra::traits::index_t>&& std:: + convertible_to, + algebra::traits::index_t>&& std:: + convertible_to< + algebra::traits::index_t, + algebra::traits::index_t< + MB>>&& requires(algebra::traits::scalar_t sa, + algebra::traits::scalar_t sb, + algebra::traits::scalar_t& 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 index 0d787dfb..581ce4f3 100644 --- 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 index 328629a4..db285330 100644 --- 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 index 6dfa680d..864e9996 100644 --- 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 index f3621774..34ba3714 100644 --- 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 index 7cc40f83..f9a933fe 100644 --- 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 index d2d0b099..f6946946 100644 --- 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 index cfec4604..66947c85 100644 --- 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 index d1a94a3d..6cc26440 100644 --- 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 index 0a20409e..09b5c0ef 100644 --- 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 index daf11920..c3d5a560 100644 --- 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 diff --git a/math/generic/include/algebra/math/impl/generic_matrix.hpp b/math/generic/include/algebra/math/impl/generic_matrix.hpp index 1bb4ce8d..706c1b14 100644 --- a/math/generic/include/algebra/math/impl/generic_matrix.hpp +++ b/math/generic/include/algebra/math/impl/generic_matrix.hpp @@ -83,6 +83,196 @@ ALGEBRA_HOST_DEVICE inline auto transpose(const M &m) { return ret; } +// Set matrix C to the product AB +template +ALGEBRA_HOST_DEVICE inline void +set_product(MC &C, const MA &A, const MB &B) requires( + algebra::concepts::matrix_multipliable_into) { + using index_t = algebra::traits::index_t; + using value_t = algebra::traits::value_t; + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + for (index_t j = 0; j < algebra::traits::columns; ++j) { + value_t t = 0.f; + + for (index_t k = 0; k < algebra::traits::rows; ++k) { + t += algebra::traits::element_getter_t()(A, i, k) * + algebra::traits::element_getter_t()(B, k, j); + } + + algebra::traits::element_getter_t()(C, i, j) = t; + } + } +} + +// Set matrix C to the product A^TB +template +ALGEBRA_HOST_DEVICE inline void +set_product_left_transpose(MC &C, const MA &A, const MB &B) requires( + algebra::concepts::matrix_multipliable_into< + decltype(transpose(std::declval())), MB, MC>) { + using index_t = algebra::traits::index_t; + using value_t = algebra::traits::value_t; + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + for (index_t j = 0; j < algebra::traits::columns; ++j) { + value_t t = 0.f; + + for (index_t k = 0; k < algebra::traits::rows; ++k) { + t += algebra::traits::element_getter_t()(A, k, i) * + algebra::traits::element_getter_t()(B, k, j); + } + + algebra::traits::element_getter_t()(C, i, j) = t; + } + } +} + +// Set matrix C to the product AB^T +template +ALGEBRA_HOST_DEVICE inline void +set_product_right_transpose(MC &C, const MA &A, const MB &B) requires( + algebra::concepts::matrix_multipliable_into< + MA, decltype(transpose(std::declval())), MC>) { + using index_t = algebra::traits::index_t; + using value_t = algebra::traits::value_t; + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + for (index_t j = 0; j < algebra::traits::columns; ++j) { + value_t t = 0.f; + + for (index_t k = 0; k < algebra::traits::columns; ++k) { + t += algebra::traits::element_getter_t()(A, i, k) * + algebra::traits::element_getter_t()(B, j, k); + } + + algebra::traits::element_getter_t()(C, i, j) = t; + } + } +} + +// Set matrix A to the product AB in place +template +ALGEBRA_HOST_DEVICE inline void +set_inplace_product_right(MA &A, const MB &B) requires( + algebra::concepts::matrix_multipliable_into) { + using index_t = algebra::traits::index_t; + using value_t = algebra::traits::value_t; + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + algebra::traits::get_matrix_t, value_t> + Q; + + for (index_t j = 0; j < algebra::traits::columns; ++j) { + algebra::traits::element_getter_t()(Q, 0, j) = + algebra::traits::element_getter_t()(A, i, j); + } + + for (index_t j = 0; j < algebra::traits::columns; ++j) { + value_t t = 0.f; + + for (index_t k = 0; k < algebra::traits::rows; ++k) { + t += algebra::traits::element_getter_t()(Q, 0, k) * + algebra::traits::element_getter_t()(B, k, j); + } + + algebra::traits::element_getter_t()(A, i, j) = t; + } + } +} + +// Set matrix A to the product BA in place +template +ALGEBRA_HOST_DEVICE inline void +set_inplace_product_left(MA &A, const MB &B) requires( + algebra::concepts::matrix_multipliable_into) { + using index_t = algebra::traits::index_t; + using value_t = algebra::traits::value_t; + + for (index_t j = 0; j < algebra::traits::columns; ++j) { + algebra::traits::get_matrix_t, value_t> + Q; + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + algebra::traits::element_getter_t()(Q, 0, i) = + algebra::traits::element_getter_t()(A, i, j); + } + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + value_t t = 0.f; + + for (index_t k = 0; k < algebra::traits::columns; ++k) { + t += algebra::traits::element_getter_t()(B, i, k) * + algebra::traits::element_getter_t()(Q, 0, k); + } + + algebra::traits::element_getter_t()(A, i, j) = t; + } + } +} + +// Set matrix A to the product AB^T in place +template +ALGEBRA_HOST_DEVICE inline void +set_inplace_product_right_transpose(MA &A, const MB &B) requires( + algebra::concepts::matrix_multipliable_into< + MA, decltype(transpose(std::declval())), MA>) { + using index_t = algebra::traits::index_t; + using value_t = algebra::traits::value_t; + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + algebra::traits::get_matrix_t, value_t> + Q; + + for (index_t j = 0; j < algebra::traits::columns; ++j) { + algebra::traits::element_getter_t()(Q, 0, j) = + algebra::traits::element_getter_t()(A, i, j); + } + + for (index_t j = 0; j < algebra::traits::columns; ++j) { + value_t T = 0.f; + + for (index_t k = 0; k < algebra::traits::columns; ++k) { + T += algebra::traits::element_getter_t()(Q, 0, k) * + algebra::traits::element_getter_t()(B, j, k); + } + + algebra::traits::element_getter_t()(A, i, j) = T; + } + } +} + +// Set matrix A to the product B^TA in place +template +ALGEBRA_HOST_DEVICE inline void +set_inplace_product_left_transpose(MA &A, const MB &B) requires( + algebra::concepts::matrix_multipliable_into< + decltype(transpose(std::declval())), MA, MA>) { + using index_t = algebra::traits::index_t; + using value_t = algebra::traits::value_t; + + for (index_t j = 0; j < algebra::traits::columns; ++j) { + algebra::traits::get_matrix_t, value_t> + Q; + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + algebra::traits::element_getter_t()(Q, 0, i) = + algebra::traits::element_getter_t()(A, i, j); + } + + for (index_t i = 0; i < algebra::traits::rows; ++i) { + value_t T = 0.f; + + for (index_t k = 0; k < algebra::traits::rows; ++k) { + T += algebra::traits::element_getter_t()(B, k, i) * + algebra::traits::element_getter_t()(Q, 0, k); + } + + algebra::traits::element_getter_t()(A, i, j) = T; + } + } +} + /// @returns the determinant of @param m template ALGEBRA_HOST_DEVICE inline algebra::traits::scalar_t determinant( diff --git a/tests/common/test_host_basics.hpp b/tests/common/test_host_basics.hpp index 72282212..24aa6c75 100644 --- a/tests/common/test_host_basics.hpp +++ b/tests/common/test_host_basics.hpp @@ -20,6 +20,7 @@ #include #include #include +#include /// Test case class, to be specialised for the different plugins - vectors template @@ -170,6 +171,8 @@ TYPED_TEST_P(test_host_basics_vector, getter) { } TYPED_TEST_P(test_host_basics_matrix, matrix_2x3) { + static constexpr typename TypeParam::size_type ROWS = 2; + static constexpr typename TypeParam::size_type COLS = 3; using matrix_2x3_t = typename TypeParam::template matrix<2, 3>; @@ -224,12 +227,110 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_2x3) { ASSERT_NEAR(v2[0], 14, this->m_epsilon); ASSERT_NEAR(v2[1], 32, this->m_epsilon); + + // Test the set_product method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } + + // Test the set_product_right_transpose method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < COLS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2; + algebra::matrix::set_product_right_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } + + // Test the set_product_left_transpose method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product_left_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } } TYPED_TEST_P(test_host_basics_matrix, matrix_3x1) { // Print the linear algebra types of this backend using algebra::operator<<; + static constexpr typename TypeParam::size_type ROWS = 3; + static constexpr typename TypeParam::size_type COLS = 1; + // Cross product on vector3 and matrix<3,1> typename TypeParam::template matrix<3, 1> vF; algebra::getter::element(vF, 0, 0) = 5.f; @@ -248,6 +349,101 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_3x1) { // Dot product on vector3 and matrix<3,1> auto dot = algebra::vector::dot(vG, vF); ASSERT_NEAR(dot, 0.f, this->m_epsilon); + + // Test the set_product method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } + + // Test the set_product_right_transpose method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < COLS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2; + algebra::matrix::set_product_right_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } + + // Test the set_product_left_transpose method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product_left_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } } TYPED_TEST_P(test_host_basics_matrix, matrix_6x4) { @@ -366,9 +562,106 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_6x4) { // Test printing std::cout << m << std::endl; + + // Test the set_product method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } + + // Test the set_product_right_transpose method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < COLS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2; + algebra::matrix::set_product_right_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } + + // Test the set_product_left_transpose method. + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < ROWS; ++j) { + algebra::getter::element(m1, i, j) = i * ROWS + j; + } + } + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + algebra::getter::element(m2, i, j) = i * COLS + j; + } + } + + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product_left_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < ROWS; ++i) { + for (std::size_t j = 0; j < COLS; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } } TYPED_TEST_P(test_host_basics_matrix, matrix_3x3) { + static constexpr typename TypeParam::size_type N = 3; + { typename TypeParam::vector3 v = {10.f, 20.f, 30.f}; typename TypeParam::template matrix<3, 3> m33; @@ -425,9 +718,124 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_3x3) { ASSERT_NEAR(algebra::getter::element(m33_inv, 2, 2), -10.f / 20.f, this->m_isclose); } + + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + algebra::getter::element(m1, i, j) = i * N + j; + algebra::getter::element(m2, i, j) = -1 * (i * N + j) + 42; + } + } + + // Test the set_product method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_product_right_transpose method. + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2; + algebra::matrix::set_product_right_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_product_left_transpose method. + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product_left_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_right method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2 = m1; + algebra::matrix::set_inplace_product_right(r2, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_left method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2 = m2; + algebra::matrix::set_inplace_product_left(r2, m1); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_right_transpose method. + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2 = m1; + algebra::matrix::set_inplace_product_right_transpose(r2, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_left_transpose method. + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2 = m2; + algebra::matrix::set_inplace_product_left_transpose(r2, m1); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } } TYPED_TEST_P(test_host_basics_matrix, matrix_2x2) { + static constexpr typename TypeParam::size_type N = 2; typename TypeParam::template matrix<2, 2> m22; algebra::getter::element(m22, 0, 0) = 4.f; @@ -449,9 +857,124 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_2x2) { this->m_isclose); ASSERT_NEAR(algebra::getter::element(m22_inv, 1, 1), 4.f / 16.f, this->m_isclose); + + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + algebra::getter::element(m1, i, j) = i * N + j; + algebra::getter::element(m2, i, j) = -1 * (i * N + j) + 42; + } + } + + // Test the set_product method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_product_right_transpose method. + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2; + algebra::matrix::set_product_right_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_product_left_transpose method. + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product_left_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_right method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2 = m1; + algebra::matrix::set_inplace_product_right(r2, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_left method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2 = m2; + algebra::matrix::set_inplace_product_left(r2, m1); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_right_transpose method. + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2 = m1; + algebra::matrix::set_inplace_product_right_transpose(r2, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_left_transpose method. + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2 = m2; + algebra::matrix::set_inplace_product_left_transpose(r2, m1); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } } TYPED_TEST_P(test_host_basics_matrix, matrix_6x6) { + static constexpr typename TypeParam::size_type N = 6; // Test 6 X 6 big matrix determinant typename TypeParam::template matrix<6, 6> m66_big; @@ -584,6 +1107,120 @@ TYPED_TEST_P(test_host_basics_matrix, matrix_6x6) { auto m66_small_det = algebra::matrix::determinant(m66_small); ASSERT_NEAR((m66_small_det - 4.30636e-11f) / 4.30636e-11f, 0.f, 2.f * this->m_isclose); + + { + typename TypeParam::template matrix m1; + typename TypeParam::template matrix m2; + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + algebra::getter::element(m1, i, j) = i * N + j; + algebra::getter::element(m2, i, j) = -1 * (i * N + j) + 42; + } + } + + // Test the set_product method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_product_right_transpose method. + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2; + algebra::matrix::set_product_right_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_product_left_transpose method. + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2; + algebra::matrix::set_product_left_transpose(r2, m1, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_right method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2 = m1; + algebra::matrix::set_inplace_product_right(r2, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_left method. + { + typename TypeParam::template matrix r1 = m1 * m2; + typename TypeParam::template matrix r2 = m2; + algebra::matrix::set_inplace_product_left(r2, m1); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_right_transpose method. + { + typename TypeParam::template matrix r1 = + m1 * algebra::matrix::transpose(m2); + typename TypeParam::template matrix r2 = m1; + algebra::matrix::set_inplace_product_right_transpose(r2, m2); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + + // Test the set_inplace_product_left_transpose method. + { + typename TypeParam::template matrix r1 = + algebra::matrix::transpose(m1) * m2; + typename TypeParam::template matrix r2 = m2; + algebra::matrix::set_inplace_product_left_transpose(r2, m1); + + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + ASSERT_NEAR(algebra::getter::element(r1, i, j), + algebra::getter::element(r2, i, j), this->m_epsilon); + } + } + } + } } TYPED_TEST_P(test_host_basics_matrix, matrix_small_mixed) {