Skip to content

Commit

Permalink
Fix symmetric_matrix_rank_k_update, hermitian_matrix_rank_k_update, a…
Browse files Browse the repository at this point in the history 
…nd symmetric_matrix_rank_1_update (#263)

* Fix tests build for submdspan

Header and namespace seem to have changed since last time.

* Add Issue #261 regression test; it fails to build

* Fix #261

Constrain symmetric_matrix_rank_k_update ExecutionPolicy overloads
so that the compiler can distinguish them from alpha overloads.

* Fix symmetric_matrix_rank_k_update algorithm

Algorithm was giving incorrect results.
This commit fixes that and
improves the associated regression test.

* Fix hermitian_matrix_rank_k_update ambiguous overloads

Add a regression test that didn't build before this change,
and now builds correctly.

* Fix hermitian_matrix_rank_k_update

1. Fix ambiguous overload for the no-scalar case
2. Fix the algorithm (so it gives the right answer,
   for both the scalar and no-scalar cases)

* Refactor hermitian_matrix_rank_k_update test

* Fix symmetric_matrix_rank_k_update

1. Fix ambiguous overload for the no-scalar case
2. Fix the algorithm (so it gives the right answer
   for the no-scalar case; the with-scalar case
   was already fixed in a previous commit)

* Add trait for custom execution policies

Add a new trait for testing whether an ExecutionPolicy is a valid
execution policy (either std::is_execution_policy_v or custom).
Use it in symmetric_matrix_rank_k_update and
hermitian_matrix_rank_k_update.

* Fix symmetric_matrix_rank_1_update ambiguous overloads

The reference implementation needs to implement all constraints of
symmetric_matrix_rank_1_update in order to disambiguate overloads.

Add a regression test as well.

* Respond to review feedback
  • Loading branch information
@mhoemmen
mhoemmen committed Jan 16, 2024
1 parent d1a1a11 commit ea57efe
Show file tree
Hide file tree
Showing 10 changed files with 834 additions and 227 deletions.
6 changes: 4 additions & 2 deletions examples/03_matrix_vector_product_mixedprec.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,9 @@
// This must be defined before including any mdspan headers.
#define MDSPAN_USE_PAREN_OPERATOR 1

#include <experimental/mdspan>
#include "experimental/__p2630_bits/submdspan.hpp"
#include <experimental/linalg>

#include <iostream>

#if defined(__cpp_lib_span)
Expand All @@ -17,7 +18,8 @@
using std::experimental::extents;
using std::full_extent; // not in experimental namespace
using std::experimental::mdspan;
using std::experimental::submdspan;

using MDSPAN_IMPL_STANDARD_NAMESPACE :: submdspan;

int main(int argc, char* argv[]) {
std::cout << "Matrix Vector Product MixedPrec" << std::endl;
Expand Down
237 changes: 193 additions & 44 deletions include/experimental/__p1673_bits/blas2_matrix_rank_1_update.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,12 +71,13 @@ struct is_custom_matrix_rank_1_update_avail<
>
: std::true_type{};

template <class Exec, class x_t, class A_t, class Tr_t, class = void>
struct is_custom_symmetric_matrix_rank_1_update_avail : std::false_type {};
template <class Exec, class ScaleFactorType, class x_t, class A_t, class Tr_t, class = void>
struct is_custom_symmetric_matrix_rank_1_update_avail : std::false_type
{};

template <class Exec, class x_t, class A_t, class Tr_t>
struct is_custom_symmetric_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(symmetric_matrix_rank_1_update
Expand All @@ -90,7 +91,28 @@ struct is_custom_symmetric_matrix_rank_1_update_avail<
&& !linalg::impl::is_inline_exec_v<Exec>
>
>
: std::true_type{};
: std::true_type
{};

template <class Exec, class ScaleFactorType, class x_t, class A_t, class Tr_t>
struct is_custom_symmetric_matrix_rank_1_update_avail<
Exec, ScaleFactorType, x_t, A_t, Tr_t,
std::enable_if_t<
std::is_void_v<
decltype(symmetric_matrix_rank_1_update
(std::declval<Exec>(),
std::declval<ScaleFactorType>(),
std::declval<x_t>(),
std::declval<A_t>(),
std::declval<Tr_t>()
)
)
>
&& !linalg::impl::is_inline_exec_v<Exec>
>
>
: std::true_type
{};

template <class Exec, class x_t, class A_t, class Tr_t, class = void>
struct is_custom_hermitian_matrix_rank_1_update_avail : std::false_type {};
Expand Down Expand Up @@ -247,26 +269,142 @@ void matrix_rank_1_update_c(
matrix_rank_1_update(exec, x, conjugated(y), A);
}

// Rank-1 update of a Symmetric matrix

template<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>
// Rank-1 update of a symmetric matrix

// Rank-1 update of a symmetric 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<Triangle, lower_triangle_t> ||
std::is_same_v<Triangle, upper_triangle_t>
)
)
void symmetric_matrix_rank_1_update(
std::experimental::linalg::impl::inline_exec_t&& /* exec */,
ScaleFactorType alpha,
std::experimental::mdspan<ElementType_x, std::experimental::extents<SizeType_x, ext_x>, Layout_x, Accessor_x> x,
std::experimental::mdspan<ElementType_A, std::experimental::extents<SizeType_A, numRows_A, numCols_A>, Layout_A, Accessor_A> A,
Triangle /* t */)
{
using size_type = ::std::common_type_t<SizeType_x, SizeType_A>;

if constexpr (std::is_same_v<Triangle, lower_triangle_t>) {
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) * 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) * 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<ExecutionPolicy> &&
(std::is_same_v<Triangle, lower_triangle_t> ||
std::is_same_v<Triangle, upper_triangle_t>)
)
)
void symmetric_matrix_rank_1_update(
ExecutionPolicy&& exec,
ScaleFactorType alpha,
std::experimental::mdspan<ElementType_x, std::experimental::extents<SizeType_x, ext_x>, Layout_x, Accessor_x> x,
std::experimental::mdspan<ElementType_A, std::experimental::extents<SizeType_A, numRows_A, numCols_A>, Layout_A, Accessor_A> A,
Triangle t)
{
constexpr bool use_custom = is_custom_symmetric_matrix_rank_1_update_avail<
decltype(execpolicy_mapper(exec)), ScaleFactorType, decltype(x), decltype(A), Triangle
>::value;

if constexpr (use_custom) {
symmetric_matrix_rank_1_update(execpolicy_mapper(exec), alpha, x, A, t);
} else {
symmetric_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<ScaleFactorType>) &&
(std::is_same_v<Triangle, lower_triangle_t> ||
std::is_same_v<Triangle, upper_triangle_t>)
)
)
void symmetric_matrix_rank_1_update(
ScaleFactorType alpha,
std::experimental::mdspan<ElementType_x, std::experimental::extents<SizeType_x, ext_x>, Layout_x, Accessor_x> x,
std::experimental::mdspan<ElementType_A, std::experimental::extents<SizeType_A, numRows_A, numCols_A>, Layout_A, Accessor_A> A,
Triangle t)
{
symmetric_matrix_rank_1_update(std::experimental::linalg::impl::default_exec_t(), alpha, x, A, t);
}

// Rank-1 update of a symmetric matrix without scaling factor alpha

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<Triangle, lower_triangle_t> ||
std::is_same_v<Triangle, upper_triangle_t>
)
)
void symmetric_matrix_rank_1_update(
std::experimental::linalg::impl::inline_exec_t&& /* exec */,
std::experimental::mdspan<ElementType_x, std::experimental::extents<SizeType_x, ext_x>, Layout_x, Accessor_x> x,
std::experimental::mdspan<ElementType_A, std::experimental::extents<SizeType_A, numRows_A, numCols_A>, Layout_A, Accessor_A> A,
Triangle /* t */)
{
using size_type = std::common_type_t<SizeType_x, SizeType_A>;

if constexpr (std::is_same_v<Triangle, lower_triangle_t>) {
for (size_type j = 0; j < A.extent(1); ++j) {
for (size_type i = j; i < A.extent(0); ++i) {
Expand All @@ -283,47 +421,57 @@ void symmetric_matrix_rank_1_update(
}
}

template<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>
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<ExecutionPolicy> &&
(std::is_same_v<Triangle, lower_triangle_t> ||
std::is_same_v<Triangle, upper_triangle_t>)
)
)
void symmetric_matrix_rank_1_update(
ExecutionPolicy&& exec,
std::experimental::mdspan<ElementType_x, std::experimental::extents<SizeType_x, ext_x>, Layout_x, Accessor_x> x,
std::experimental::mdspan<ElementType_A, std::experimental::extents<SizeType_A, numRows_A, numCols_A>, Layout_A, Accessor_A> A,
Triangle t)
{

constexpr bool use_custom = is_custom_symmetric_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) {
symmetric_matrix_rank_1_update(execpolicy_mapper(exec), x, A, t);
}
else
{
} else {
symmetric_matrix_rank_1_update(std::experimental::linalg::impl::inline_exec_t(), x, A, t);
}
}

template<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>
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<Triangle, lower_triangle_t> ||
std::is_same_v<Triangle, upper_triangle_t>
)
)
void symmetric_matrix_rank_1_update(
std::experimental::mdspan<ElementType_x, std::experimental::extents<SizeType_x, ext_x>, Layout_x, Accessor_x> x,
std::experimental::mdspan<ElementType_A, std::experimental::extents<SizeType_A, numRows_A, numCols_A>, Layout_A, Accessor_A> A,
Expand All @@ -332,8 +480,9 @@ void symmetric_matrix_rank_1_update(
symmetric_matrix_rank_1_update(std::experimental::linalg::impl::default_exec_t(), x, A, t);
}

// Rank-k update of a Hermitian matrix

// Rank-1 update of a Hermitian matrix
// Rank-1 update of a Hermitian matrix without scaling factor alpha

template<class ElementType_x,
class SizeType_x, ::std::size_t ext_x,
Expand Down
Loading

0 comments on commit ea57efe

Please sign in to comment.