Make scale work for std::execution::par

CMakeLists.txt

@@ -121,6 +121,14 @@ option(LINALG_ENABLE_BLAS
   "Assume that we are linking with a BLAS library."
   ${BLAS_FOUND})
 
+find_package(TBB)
+option(LINALG_ENABLE_TBB
+  "Enable TBB. Default: autodetect TBB installation."
+  ${TBB_FOUND})
+if(LINALG_ENABLE_TBB)
+  find_package(TBB REQUIRED)
+endif()
+
 find_package(KokkosKernels)
 option(LINALG_ENABLE_KOKKOS
   "Enable Kokkos-based implementation. Default: autodetect Kokkos installation."
@@ -144,6 +152,10 @@ add_library(std::linalg ALIAS linalg)
 
 target_link_libraries(linalg INTERFACE std::mdspan)
 
+if(LINALG_ENABLE_TBB)
+  target_link_libraries(linalg INTERFACE TBB::tbb)
+endif()
+
 if(LINALG_ENABLE_KOKKOS)
   target_link_libraries(linalg INTERFACE Kokkos::kokkos)
   target_link_libraries(linalg INTERFACE Kokkos::kokkoskernels)

README.md

@@ -28,6 +28,11 @@ Other compilers, including MSVC 2019, have been tested in the past.
    - If you want to build examples, set LINALG_ENABLE_EXAMPLES=ON
    - If you have a BLAS installation, set LINALG_ENABLE_BLAS=ON.
      BLAS support is currently experimental.
+   - If you have a TBB (Threading Building Blocks) installation
+     and want to use TBB, set LINALG_ENABLE_TBB=ON (and optionally
+     set TBB_DIR to the lib/cmake/TBB subdirectory of your TBB installation,
+     or wherever the TBBConfig.cmake file happens to live).
+     TBB support is currently experimental.
 4. Build and install as usual
 5. If you enabled tests, use "ctest" to run them
 

examples/01_scale.cpp

@@ -8,11 +8,7 @@
 
 #include <iostream>
 
-#if (! defined(__GNUC__)) || (__GNUC__ > 9)
-#  define MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES 1
-#endif
-
-#ifdef MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES
+#ifdef LINALG_HAS_EXECUTION
 #  include <execution>
 #endif
 
@@ -43,7 +39,7 @@ int main(int argc, char* argv[]) {
 
     // Call linalg::scale x = 2.0*x;
     std::experimental::linalg::scale(2.0, x);
-#ifdef MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES
+#ifdef LINALG_HAS_EXECUTION
     std::experimental::linalg::scale(std::execution::par, 2.0, x);
 #else
     std::experimental::linalg::scale(2.0, x);

examples/02_matrix_vector_product_basic.cpp

@@ -8,11 +8,7 @@
 
 #include <iostream>
 
-#if (! defined(__GNUC__)) || (__GNUC__ > 9)
-#  define MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES 1
-#endif
-
-#ifdef MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES
+#ifdef LINALG_HAS_EXECUTION
 #  include <execution>
 #endif
 
@@ -51,7 +47,7 @@ int main(int argc, char* argv[]) {
     std::experimental::linalg::matrix_vector_product(A, x, y);
 
     // y = 0.5 * y + 2 * A * x
-#ifdef MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES
+#ifdef LINALG_HAS_EXECUTION
     std::experimental::linalg::matrix_vector_product(std::execution::par,
        std::experimental::linalg::scaled(2.0, A), x,
        std::experimental::linalg::scaled(0.5, y), y);

include/experimental/__p1673_bits/blas1_dot.hpp

@@ -129,11 +129,11 @@ Scalar dot(
                 v1.static_extent(0) == v2.static_extent(0));
 
   constexpr bool use_custom = is_custom_dot_avail<
-    decltype(execpolicy_mapper(exec)), decltype(v1), decltype(v2), Scalar
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(v1), decltype(v2), Scalar
     >::value;
 
   if constexpr(use_custom){
-    return dot(execpolicy_mapper(exec), v1, v2, init);
+    return dot(detail::map_execpolicy_with_check(exec), v1, v2, init);
   }
   else{
     return dot(std::experimental::linalg::impl::inline_exec_t(), v1, v2, init);

include/experimental/__p1673_bits/blas1_givens.hpp

@@ -427,13 +427,12 @@ void givens_rotation_apply(
   const Real c,
   const Real s)
 {
-
   constexpr bool use_custom = is_custom_givens_rotation_apply_avail<
-    decltype(execpolicy_mapper(exec)), decltype(x), decltype(y), Real, Real
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(x), decltype(y), Real, Real
     >::value;
 
-  if constexpr(use_custom){
-    givens_rotation_apply(execpolicy_mapper(exec), x, y, c, s);
+  if constexpr (use_custom) {
+    givens_rotation_apply(detail::map_execpolicy_with_check(exec), x, y, c, s);
   }
   else
   {
@@ -524,16 +523,14 @@ void givens_rotation_apply(
   const Real c,
   const complex<Real> s)
 {
-
   constexpr bool use_custom = is_custom_givens_rotation_apply_avail<
-    decltype(execpolicy_mapper(exec)), decltype(x), decltype(y), Real, complex<Real>
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(x), decltype(y), Real, complex<Real>
     >::value;
 
-  if constexpr(use_custom){
-    givens_rotation_apply(execpolicy_mapper(exec), x, y, c, s);
+  if constexpr (use_custom) {
+    givens_rotation_apply(detail::map_execpolicy_with_check(exec), x, y, c, s);
   }
-  else
-  {
+  else {
     givens_rotation_apply(std::experimental::linalg::impl::inline_exec_t(), x, y, c, s);
   }
 }

include/experimental/__p1673_bits/blas1_linalg_add.hpp

@@ -222,17 +222,15 @@ void add(
   std::experimental::mdspan<ElementType_y, std::experimental::extents<SizeType_y, ext_y ...>, Layout_y, Accessor_y> y,
   std::experimental::mdspan<ElementType_z, std::experimental::extents<SizeType_z, ext_z ...>, Layout_z, Accessor_z> z)
 {
-
   constexpr bool use_custom = is_custom_add_avail<
-    decltype(execpolicy_mapper(exec)), decltype(x), decltype(y), decltype(z)
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(x), decltype(y), decltype(z)
     >::value;
 
-  if constexpr(use_custom){
+  if constexpr (use_custom) {
     // for the customization point, it is up to impl to check requirements
-    add(execpolicy_mapper(exec), x, y, z);
+    add(detail::map_execpolicy_with_check(exec), x, y, z);
   }
-  else
-  {
+  else {
     add(std::experimental::linalg::impl::inline_exec_t(), x, y, z);
   }
 }

include/experimental/__p1673_bits/blas1_linalg_copy.hpp

@@ -171,16 +171,14 @@ void copy(
   std::experimental::mdspan<ElementType_x, std::experimental::extents<SizeType_x, ext_x ...>, Layout_x, Accessor_x> x,
   std::experimental::mdspan<ElementType_y, std::experimental::extents<SizeType_y, ext_y ...>, Layout_y, Accessor_y> y)
 {
-
   constexpr bool use_custom = is_custom_copy_avail<
-    decltype(execpolicy_mapper(exec)), decltype(x), decltype(y)
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(x), decltype(y)
     >::value;
 
-  if constexpr(use_custom){
-    copy(execpolicy_mapper(exec), x, y);
+  if constexpr (use_custom) {
+    copy(detail::map_execpolicy_with_check(exec), x, y);
   }
-  else
-  {
+  else {
     copy(std::experimental::linalg::impl::inline_exec_t(), x, y);
   }
 }

include/experimental/__p1673_bits/blas1_linalg_swap.hpp

@@ -179,14 +179,13 @@ void swap_elements(
   std::experimental::mdspan<ElementType_y, std::experimental::extents<SizeType_y, ext_y ...>, Layout_y, Accessor_y> y)
 {
   constexpr bool use_custom = is_custom_vector_swap_elements_avail<
-    decltype(execpolicy_mapper(exec)), decltype(x), decltype(y)
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(x), decltype(y)
     >::value;
 
-  if constexpr(use_custom){
-    return swap_elements(execpolicy_mapper(exec), x, y);
+  if constexpr (use_custom) {
+    return swap_elements(detail::map_execpolicy_with_check(exec), x, y);
   }
-  else
-  {
+  else {
     return swap_elements(std::experimental::linalg::impl::inline_exec_t(), x, y);
   }
 }

include/experimental/__p1673_bits/blas1_matrix_frob_norm.hpp

@@ -138,15 +138,14 @@ Scalar matrix_frob_norm(
   std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, numRows, numCols>, Layout, Accessor> A,
   Scalar init)
 {
-
   constexpr bool use_custom = is_custom_matrix_frob_norm_avail<
-    decltype(execpolicy_mapper(exec)), decltype(A), Scalar
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(A), Scalar
     >::value;
 
-  if constexpr(use_custom){
-    return matrix_frob_norm(execpolicy_mapper(exec), A, init);
+  if constexpr (use_custom) {
+    return matrix_frob_norm(detail::map_execpolicy_with_check(exec), A, init);
   }
-  else{
+  else {
     return matrix_frob_norm(std::experimental::linalg::impl::inline_exec_t(), A, init);
   }
 }

include/experimental/__p1673_bits/blas1_matrix_inf_norm.hpp

@@ -128,15 +128,14 @@ Scalar matrix_inf_norm(
   std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, numRows, numCols>, Layout, Accessor> A,
   Scalar init)
 {
-
   constexpr bool use_custom = is_custom_matrix_inf_norm_avail<
-    decltype(execpolicy_mapper(exec)), decltype(A), Scalar
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(A), Scalar
     >::value;
 
-  if constexpr(use_custom){
-    return matrix_inf_norm(execpolicy_mapper(exec), A, init);
+  if constexpr (use_custom) {
+    return matrix_inf_norm(detail::map_execpolicy_with_check(exec), A, init);
   }
-  else{
+  else {
     return matrix_inf_norm(std::experimental::linalg::impl::inline_exec_t(), A, init);
   }
 }

include/experimental/__p1673_bits/blas1_matrix_one_norm.hpp

@@ -133,13 +133,13 @@ Scalar matrix_one_norm(
 {
 
   constexpr bool use_custom = is_custom_matrix_one_norm_avail<
-    decltype(execpolicy_mapper(exec)), decltype(A), Scalar
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(A), Scalar
     >::value;
 
-  if constexpr(use_custom){
-    return matrix_one_norm(execpolicy_mapper(exec), A, init);
+  if constexpr (use_custom) {
+    return matrix_one_norm(detail::map_execpolicy_with_check(exec), A, init);
   }
-  else{
+  else {
     return matrix_one_norm(std::experimental::linalg::impl::inline_exec_t(), A, init);
   }
 }

include/experimental/__p1673_bits/blas1_scale.hpp

@@ -43,6 +43,12 @@
 #ifndef LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS1_SCALE_HPP_
 #define LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS1_SCALE_HPP_
 
+#include <algorithm> // for __cpp_lib_ranges, if it exists
+#if defined(LINALG_HAS_EXECUTION) && defined(__cpp_lib_ranges)
+#  include <execution>
+#  include <ranges>
+#endif
+
 namespace std {
 namespace experimental {
 inline namespace __p1673_version_0 {
@@ -123,6 +129,40 @@ void scale(
   }
 }
 
+#if defined(LINALG_HAS_EXECUTION) && defined(__cpp_lib_ranges)
+template<class Scalar,
+         class ElementType,
+	 class SizeType,
+         ::std::size_t ... ext,
+         class Layout,
+         class Accessor>
+void scale(
+  std::experimental::linalg::impl::parallel_exec_t&& /* exec */,
+  const Scalar alpha,
+  std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext ...>, Layout, Accessor> x)
+{
+  static_assert(x.rank() <= 2);
+
+  if constexpr (x.rank() == 1) {
+    auto range = std::views::iota(SizeType(0), x.extent(0));
+    std::for_each(std::execution::par, range.begin(), range.end(), [=] (SizeType k) {
+      x[k] *= alpha;
+    });
+  }
+  else if constexpr (x.rank() == 2) {
+#if defined(__cpp_lib_ranges_cartesian_product)
+    auto range = std::views::cartesian_product(std::views::iota(SizeType(0), x.extent(0)), std::views::iota(SizeType(0), x.extent(1)));
+    std::for_each(std::execution::par, range.begin(), range.end(), [=] (auto&& inds) {
+      auto [i, j] = inds;
+      x[i, j] *= alpha;
+    });
+#else
+    linalg_scale_rank_2(alpha, x);
+#endif
+  }
+}
+#endif // LINALG_HAS_EXECUTION && __cpp_lib_ranges
+
 template<class ExecutionPolicy,
          class Scalar,
          class ElementType,
@@ -138,11 +178,11 @@ void scale(
   // Call custom overload if available else call std implementation
 
   constexpr bool use_custom = is_custom_scale_avail<
-    decltype(execpolicy_mapper(exec)), decltype(alpha), decltype(x)
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(alpha), decltype(x)
     >::value;
 
-  if constexpr(use_custom) {
-    scale(execpolicy_mapper(exec), alpha, x);
+  if constexpr (use_custom) {
+    scale(detail::map_execpolicy_with_check(exec), alpha, x);
   } else {
     scale(std::experimental::linalg::impl::inline_exec_t(), alpha, x);
   }

include/experimental/__p1673_bits/blas1_vector_abs_sum.hpp

@@ -104,11 +104,11 @@ Scalar vector_abs_sum(
   Scalar init)
 {
   constexpr bool use_custom = is_custom_vector_abs_sum_avail<
-    decltype(execpolicy_mapper(exec)), decltype(v), Scalar
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(v), Scalar
     >::value;
 
   if constexpr(use_custom){
-    return vector_abs_sum(execpolicy_mapper(exec), v, init);
+    return vector_abs_sum(detail::map_execpolicy_with_check(exec), v, init);
   }
   else
   {

include/experimental/__p1673_bits/blas1_vector_idx_abs_max.hpp

@@ -120,13 +120,13 @@ SizeType idx_abs_max(
   }
 
   constexpr bool use_custom = is_custom_idx_abs_max_avail<
-    decltype(execpolicy_mapper(exec)), decltype(v)
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(v)
     >::value;
 
-  if constexpr(use_custom){
-    return idx_abs_max(execpolicy_mapper(exec), v);
+  if constexpr (use_custom) {
+    return idx_abs_max(detail::map_execpolicy_with_check(exec), v);
   }
-  else{
+  else {
     return idx_abs_max(std::experimental::linalg::impl::inline_exec_t(), v);
   }
 }

include/experimental/__p1673_bits/blas1_vector_norm2.hpp

@@ -111,11 +111,11 @@ Scalar vector_norm2(
   Scalar init)
 {
   constexpr bool use_custom = is_custom_vector_norm2_avail<
-    decltype(execpolicy_mapper(exec)), decltype(x), Scalar
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(x), Scalar
     >::value;
 
   if constexpr(use_custom){
-    return vector_norm2(execpolicy_mapper(exec), x, init);
+    return vector_norm2(detail::map_execpolicy_with_check(exec), x, init);
   }
   else
   {

include/experimental/__p1673_bits/blas1_vector_sum_of_squares.hpp

@@ -136,16 +136,14 @@ sum_of_squares_result<Scalar> vector_sum_of_squares(
   std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> v,
   sum_of_squares_result<Scalar> init)
 {
-
   constexpr bool use_custom = is_custom_vector_sum_of_squares_avail<
-    decltype(execpolicy_mapper(exec)), decltype(v), Scalar
+    decltype(detail::map_execpolicy_with_check(exec)), decltype(v), Scalar
     >::value;
 
-  if constexpr(use_custom){
-    return vector_sum_of_squares(execpolicy_mapper(exec), v, init);
+  if constexpr (use_custom) {
+    return vector_sum_of_squares(detail::map_execpolicy_with_check(exec), v, init);
   }
-  else
-  {
+  else {
     return vector_sum_of_squares(std::experimental::linalg::impl::inline_exec_t(), v, init);
   }
 }