diff --git a/include/experimental/__p1673_bits/blas2_matrix_rank_1_update.hpp b/include/experimental/__p1673_bits/blas2_matrix_rank_1_update.hpp index 45073e79..196ae650 100644 --- a/include/experimental/__p1673_bits/blas2_matrix_rank_1_update.hpp +++ b/include/experimental/__p1673_bits/blas2_matrix_rank_1_update.hpp @@ -114,12 +114,13 @@ struct is_custom_symmetric_matrix_rank_1_update_avail< : std::true_type {}; -template -struct is_custom_hermitian_matrix_rank_1_update_avail : std::false_type {}; +template +struct is_custom_hermitian_matrix_rank_1_update_avail : std::false_type +{}; template struct is_custom_hermitian_matrix_rank_1_update_avail< - Exec, x_t, A_t, Tr_t, + Exec, void, x_t, A_t, Tr_t, std::enable_if_t< std::is_void_v< decltype(hermitian_matrix_rank_1_update @@ -133,7 +134,28 @@ struct is_custom_hermitian_matrix_rank_1_update_avail< && !linalg::impl::is_inline_exec_v > > - : std::true_type{}; + : std::true_type +{}; + +template +struct is_custom_hermitian_matrix_rank_1_update_avail< + Exec, ScaleFactorType, x_t, A_t, Tr_t, + std::enable_if_t< + std::is_void_v< + decltype(hermitian_matrix_rank_1_update + (std::declval(), + std::declval(), + std::declval(), + std::declval(), + std::declval() + ) + ) + > + && !linalg::impl::is_inline_exec_v + > + > + : std::true_type +{}; } // end anonymous namespace @@ -297,7 +319,7 @@ void symmetric_matrix_rank_1_update( std::experimental::mdspan, Layout_A, Accessor_A> A, Triangle /* t */) { - using size_type = ::std::common_type_t; + using size_type = std::common_type_t; if constexpr (std::is_same_v) { for (size_type j = 0; j < A.extent(1); ++j) { @@ -482,25 +504,139 @@ void symmetric_matrix_rank_1_update( // Rank-k update of a Hermitian matrix +// Rank-1 update of a Hermitian matrix with scaling factor alpha + +MDSPAN_TEMPLATE_REQUIRES( + class ScaleFactorType, + class ElementType_x, + class SizeType_x, ::std::size_t ext_x, + class Layout_x, + class Accessor_x, + class ElementType_A, + class SizeType_A, ::std::size_t numRows_A, + ::std::size_t numCols_A, + class Layout_A, + class Accessor_A, + class Triangle, + /* requires */ ( + std::is_same_v || + std::is_same_v + ) +) +void hermitian_matrix_rank_1_update( + std::experimental::linalg::impl::inline_exec_t&& /* exec */, + ScaleFactorType alpha, + std::experimental::mdspan, Layout_x, Accessor_x> x, + std::experimental::mdspan, Layout_A, Accessor_A> A, + Triangle /* t */) +{ + using size_type = std::common_type_t; + + if constexpr (std::is_same_v) { + for (size_type j = 0; j < A.extent(1); ++j) { + for (size_type i = j; i < A.extent(0); ++i) { + A(i,j) += alpha * x(i) * impl::conj_if_needed(x(j)); + } + } + } + else { + for (size_type j = 0; j < A.extent(1); ++j) { + for (size_type i = 0; i <= j; ++i) { + A(i,j) += alpha * x(i) * impl::conj_if_needed(x(j)); + } + } + } +} + +MDSPAN_TEMPLATE_REQUIRES( + class ExecutionPolicy, + class ScaleFactorType, + class ElementType_x, + class SizeType_x, ::std::size_t ext_x, + class Layout_x, + class Accessor_x, + class ElementType_A, + class SizeType_A, ::std::size_t numRows_A, + ::std::size_t numCols_A, + class Layout_A, + class Accessor_A, + class Triangle, + /* requires */ ( + impl::is_linalg_execution_policy_other_than_inline_v && + (std::is_same_v || + std::is_same_v) + ) +) +void hermitian_matrix_rank_1_update( + ExecutionPolicy&& exec, + ScaleFactorType alpha, + std::experimental::mdspan, Layout_x, Accessor_x> x, + std::experimental::mdspan, Layout_A, Accessor_A> A, + Triangle t) +{ + constexpr bool use_custom = is_custom_hermitian_matrix_rank_1_update_avail< + decltype(execpolicy_mapper(exec)), ScaleFactorType, decltype(x), decltype(A), Triangle + >::value; + + if constexpr (use_custom) { + hermitian_matrix_rank_1_update(execpolicy_mapper(exec), alpha, x, A, t); + } else { + hermitian_matrix_rank_1_update(std::experimental::linalg::impl::inline_exec_t(), alpha, x, A, t); + } +} + +MDSPAN_TEMPLATE_REQUIRES( + class ScaleFactorType, + class ElementType_x, + class SizeType_x, ::std::size_t ext_x, + class Layout_x, + class Accessor_x, + class ElementType_A, + class SizeType_A, ::std::size_t numRows_A, + ::std::size_t numCols_A, + class Layout_A, + class Accessor_A, + class Triangle, + /* requires */ ( + (! impl::is_linalg_execution_policy_other_than_inline_v) && + (std::is_same_v || + std::is_same_v) + ) +) +void hermitian_matrix_rank_1_update( + ScaleFactorType alpha, + std::experimental::mdspan, Layout_x, Accessor_x> x, + std::experimental::mdspan, Layout_A, Accessor_A> A, + Triangle t) +{ + hermitian_matrix_rank_1_update(std::experimental::linalg::impl::default_exec_t(), alpha, x, A, t); +} + // Rank-1 update of a Hermitian matrix without scaling factor alpha -template +MDSPAN_TEMPLATE_REQUIRES( + class ElementType_x, + class SizeType_x, ::std::size_t ext_x, + class Layout_x, + class Accessor_x, + class ElementType_A, + class SizeType_A, ::std::size_t numRows_A, + ::std::size_t numCols_A, + class Layout_A, + class Accessor_A, + class Triangle, + /* requires */ ( + std::is_same_v || + std::is_same_v + ) +) void hermitian_matrix_rank_1_update( std::experimental::linalg::impl::inline_exec_t&& /* exec */, std::experimental::mdspan, Layout_x, Accessor_x> x, std::experimental::mdspan, Layout_A, Accessor_A> A, Triangle /* t */) { - using size_type = ::std::common_type_t; + using size_type = std::common_type_t; if constexpr (std::is_same_v) { for (size_type j = 0; j < A.extent(1); ++j) { @@ -518,47 +654,57 @@ void hermitian_matrix_rank_1_update( } } -template +MDSPAN_TEMPLATE_REQUIRES( + class ExecutionPolicy, + class ElementType_x, + class SizeType_x, ::std::size_t ext_x, + class Layout_x, + class Accessor_x, + class ElementType_A, + class SizeType_A, ::std::size_t numRows_A, + ::std::size_t numCols_A, + class Layout_A, + class Accessor_A, + class Triangle, + /* requires */ ( + impl::is_linalg_execution_policy_other_than_inline_v && + (std::is_same_v || + std::is_same_v) + ) +) void hermitian_matrix_rank_1_update( ExecutionPolicy&& exec, std::experimental::mdspan, Layout_x, Accessor_x> x, std::experimental::mdspan, Layout_A, Accessor_A> A, Triangle t) { - constexpr bool use_custom = is_custom_hermitian_matrix_rank_1_update_avail< - decltype(execpolicy_mapper(exec)), decltype(x), decltype(A), Triangle + decltype(execpolicy_mapper(exec)), void, decltype(x), decltype(A), Triangle >::value; - if constexpr(use_custom){ + if constexpr (use_custom) { hermitian_matrix_rank_1_update(execpolicy_mapper(exec), x, A, t); - } - else - { + } else { hermitian_matrix_rank_1_update(std::experimental::linalg::impl::inline_exec_t(), x, A, t); } } -template +MDSPAN_TEMPLATE_REQUIRES( + class ElementType_x, + class SizeType_x, ::std::size_t ext_x, + class Layout_x, + class Accessor_x, + class ElementType_A, + class SizeType_A, ::std::size_t numRows_A, + ::std::size_t numCols_A, + class Layout_A, + class Accessor_A, + class Triangle, + /* requires */ ( + std::is_same_v || + std::is_same_v + ) +) void hermitian_matrix_rank_1_update( std::experimental::mdspan, Layout_x, Accessor_x> x, std::experimental::mdspan, Layout_A, Accessor_A> A, diff --git a/tests/native/CMakeLists.txt b/tests/native/CMakeLists.txt index 3c14e984..e0f910c5 100644 --- a/tests/native/CMakeLists.txt +++ b/tests/native/CMakeLists.txt @@ -24,6 +24,7 @@ linalg_add_test(gemv_no_ambig) linalg_add_test(givens) linalg_add_test(hemm) linalg_add_test(herk) +linalg_add_test(her) linalg_add_test(idx_abs_max) linalg_add_test(iterator) linalg_add_test(matrix_inf_norm) diff --git a/tests/native/her.cpp b/tests/native/her.cpp new file mode 100644 index 00000000..b0cf8e51 --- /dev/null +++ b/tests/native/her.cpp @@ -0,0 +1,111 @@ +#include "./gtest_fixtures.hpp" + +#include +#include +#include + +namespace { + using std::experimental::linalg::lower_triangle; + using std::experimental::linalg::hermitian_matrix_rank_1_update; + using std::experimental::linalg::upper_triangle; + using std::extents; + using std::layout_right; + using std::mdspan; + + // Regression test for ambiguous overloads of + // hermitian_matrix_rank_1_update (related to + // https://github.com/kokkos/stdBLAS/issues/261 ). + // + // The reference implementation needs to implement all constraints + // of hermitian_matrix_rank_1_update in order to disambiguate + // overloads. + TEST(BLAS3_her, AmbiguousOverloads) + { + constexpr auto map_A = layout_right::mapping{extents{}}; + constexpr auto map_expected = map_A; + constexpr auto map_x = layout_right::mapping{extents{}}; + using V = std::complex; + + // A = [-1.0 -2.0 -4.0] + // [-2.0 -3.0 -5.0] + // [-4.0 -5.0 -6.0] + // + // x = [ 2.0 + 7.0i] + // [ 5.0] + // [ 11.0] + // + // x x^H = [53.0 10.0 + 35.0i 22.0 + 77.0i] + // [10.0 - 35.0i 25.0 55.0 ] + // [22.0 - 77.0i 55.0 121.0 ] + // + // A + x x^H = [52.0 8.0 + 35.0i 18.0 + 77.0i] + // [ 8.0 - 35.0i 22.0 50.0] + // [18.0 - 77.0i 50.0 115.0] + constexpr std::array A_storage_original{ + V(-1.0, 0.0), V(-2.0, 0.0), V(-4.0, 0.0), + V(-2.0, 0.0), V(-3.0, 0.0), V(-5.0, 0.0), + V(-4.0, 0.0), V(-5.0, 0.0), V(-6.0, 0.0) + }; + constexpr std::array x_storage_original{ + V( 2.0, 7.0), + V( 5.0, 0.0), + V(11.0, 0.0) + }; + constexpr std::array expected_storage_original{ + V(52.0, 0.0), V( 8.0, 35.0), V( 18.0, 77.0), + V( 8.0, -35.0), V(22.0, 0.0), V( 50.0, 0.0), + V(18.0, -77.0), V(50.0, 0.0), V(115.0, 0.0) + }; + + auto A_storage = A_storage_original; + mdspan A{A_storage.data(), map_A}; + + auto expected_storage = expected_storage_original; + mdspan expected{expected_storage.data(), map_expected}; + + auto x_storage = x_storage_original; + mdspan x{x_storage.data(), map_x}; + + auto check_upper_triangle = [&] () { + for (std::size_t row = 0; row < A.extent(0); ++row) { + for (std::size_t col = row; col < A.extent(1); ++col) { + const auto expected_rc = expected(row, col); + const auto A_rc = A(row, col); + EXPECT_EQ(expected_rc, A_rc) << "at (" << row << ", " << col << ")"; + } + } + }; + auto check_lower_triangle = [&] () { + for (std::size_t row = 0; row < A.extent(0); ++row) { + for (std::size_t col = 0; col <= row; ++col) { + const auto expected_rc = expected(row, col); + const auto A_rc = A(row, col); + EXPECT_EQ(expected_rc, A_rc) << "at (" << row << ", " << col << ")"; + } + } + }; + + hermitian_matrix_rank_1_update(1.0, x, A, upper_triangle); + check_upper_triangle(); + + // Reset values, just in case some bug might have overwritten them. + A_storage = A_storage_original; + expected_storage = expected_storage_original; + x_storage = x_storage_original; + hermitian_matrix_rank_1_update(1.0, x, A, lower_triangle); + check_lower_triangle(); + + A_storage = A_storage_original; + expected_storage = expected_storage_original; + x_storage = x_storage_original; + hermitian_matrix_rank_1_update(x, A, upper_triangle); + check_upper_triangle(); + + A_storage = A_storage_original; + expected_storage = expected_storage_original; + x_storage = x_storage_original; + hermitian_matrix_rank_1_update(x, A, lower_triangle); + check_lower_triangle(); + } + +} // end anonymous namespace diff --git a/tests/native/herk.cpp b/tests/native/herk.cpp index 289fb2c4..838b7409 100644 --- a/tests/native/herk.cpp +++ b/tests/native/herk.cpp @@ -17,7 +17,7 @@ namespace { // The reference implementation needs to implement all constraints // of hermitian_matrix_rank_k_update in order to disambiguate // overloads. - TEST(BLAS3_herk, Issue261_FollowOn) + TEST(BLAS3_herk, AmbiguousOverloads) { constexpr auto map_C = layout_left::mapping{extents{}}; constexpr auto map_expected = map_C; diff --git a/tests/native/syr.cpp b/tests/native/syr.cpp index 2bc0b079..ee9c8547 100644 --- a/tests/native/syr.cpp +++ b/tests/native/syr.cpp @@ -13,12 +13,12 @@ namespace { // Regression test for ambiguous overloads of // symmetric_matrix_rank_1_update (related to - // https://github.com/kokkos/stdBLAS/issues/261 ). + // https://github.com/kokkos/stdBLAS/issues/261). // // The reference implementation needs to implement all constraints // of symmetric_matrix_rank_1_update in order to disambiguate // overloads. - TEST(BLAS3_syrk, Issue261) + TEST(BLAS3_syr, AmbiguousOverloads) { constexpr auto map_A = layout_right::mapping{extents{}}; constexpr auto map_expected = map_A; diff --git a/tests/native/syrk.cpp b/tests/native/syrk.cpp index d8864954..40a6641c 100644 --- a/tests/native/syrk.cpp +++ b/tests/native/syrk.cpp @@ -16,7 +16,7 @@ namespace { // The reference implementation needs to implement all constraints // of symmetric_matrix_rank_k_update in order to disambiguate // overloads. - TEST(BLAS3_syrk, Issue261) + TEST(BLAS3_syrk, AmbiguousOverloads_Issue261) { constexpr auto map_C = layout_left::mapping{extents{}}; constexpr auto map_expected = map_C;