Draft: Attempt to specialize matrix_vector_product for parallel_policy

.github/workflows/cmake.yml

@@ -102,8 +102,6 @@ jobs:
 
   test-stdBLAS:
     runs-on: ubuntu-latest
-    container:
-      image: amklinv/mdspan-dependencies:latest
     needs: build-stdblas
 
     steps:

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,8 @@ 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 installation, set LINALG_ENABLE_TBB=ON.
+     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

@@ -2,18 +2,10 @@
 
 #include <iostream>
 
-#if (! defined(__GNUC__)) || (__GNUC__ > 9)
-#  define MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES 1
-#endif
-
-#ifdef MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES
-#  include <execution>
-#endif
-
 // Make mdspan less verbose
-using std::experimental::mdspan;
-using std::experimental::extents;
-using std::experimental::dynamic_extent;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::extents;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent;
 
 int main(int argc, char* argv[]) {
   std::cout << "Scale" << std::endl;
@@ -30,7 +22,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

@@ -2,18 +2,10 @@
 
 #include <iostream>
 
-#if (! defined(__GNUC__)) || (__GNUC__ > 9)
-#  define MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES 1
-#endif
-
-#ifdef MDSPAN_EXAMPLES_USE_EXECUTION_POLICIES
-#  include <execution>
-#endif
-
 // Make mdspan less verbose
-using std::experimental::mdspan;
-using std::experimental::extents;
-using std::experimental::dynamic_extent;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::extents;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent;
 
 int main(int argc, char* argv[]) {
   std::cout << "Matrix Vector Product Basic" << std::endl;
@@ -41,7 +33,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);

examples/03_matrix_vector_product_mixedprec.cpp

@@ -3,11 +3,11 @@
 #include <iostream>
 
 // Make mdspan less verbose
-using std::experimental::mdspan;
-using std::experimental::extents;
-using std::experimental::dynamic_extent;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::extents;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent;
 using std::experimental::submdspan;
-using std::experimental::full_extent;
+using MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent;
 
 int main(int argc, char* argv[]) {
   std::cout << "Matrix Vector Product MixedPrec" << std::endl;

include/experimental/__p1673_bits/blas2_matrix_vector_product.hpp

@@ -44,6 +44,7 @@
 #define LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS2_MATRIX_VECTOR_PRODUCT_HPP_
 
 #include <complex>
+#include <iostream>
 #include <type_traits>
 
 namespace std {
@@ -271,6 +272,7 @@ void matrix_vector_product(
   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)
 {
+  std::cerr << "matrix_vector_product with inline_exec_t\n";
   using size_type = std::common_type_t<
     std::common_type_t<
       std::common_type_t<SizeType_A, SizeType_x>,
@@ -281,8 +283,47 @@ void matrix_vector_product(
       y(i) += A(i,j) * x(j);
     }
   }
+  std::cerr << "exiting matrix_vector_product with inline_exec_t\n";
 }
 
+#ifdef LINALG_HAS_EXECUTION
+MDSPAN_TEMPLATE_REQUIRES(
+         class ElementType_A,
+         class SizeType_A, ::std::size_t numRows_A,
+         ::std::size_t numCols_A,
+         class Layout_A,
+         class Accessor_A,
+         class ElementType_x,
+         class SizeType_x, ::std::size_t ext_x,
+         class Layout_x,
+         class Accessor_x,
+         class ElementType_y,
+         class SizeType_y, ::std::size_t ext_y,
+         class Layout_y,
+         class Accessor_y,
+         /* requires */ (Layout_A::template mapping<extents<SizeType_A, numRows_A, numCols_A> >::is_always_unique())
+)
+void matrix_vector_product(
+  std::execution::parallel_policy&& /* exec */,
+  std::experimental::mdspan<ElementType_A, std::experimental::extents<SizeType_A, numRows_A, numCols_A>, Layout_A, Accessor_A> A,
+  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)
+{
+  std::cerr << "matrix_vector_product with par\n";
+  using size_type = std::common_type_t<
+    std::common_type_t<
+      std::common_type_t<SizeType_A, SizeType_x>,
+      SizeType_y>>;
+  for (size_type i = 0; i < A.extent(0); ++i) {
+    y(i) = ElementType_y{};
+    for (size_type j = 0; j < A.extent(1); ++j) {
+      y(i) += A(i,j) * x(j);
+    }
+  }
+  std::cerr << "exiting matrix_vector_product with par\n";
+}
+#endif
+
 MDSPAN_TEMPLATE_REQUIRES(
          class ExecutionPolicy,
          class ElementType_A,
@@ -310,7 +351,15 @@ void matrix_vector_product(
   constexpr bool use_custom = is_custom_mat_vec_product_avail<
     decltype(execpolicy_mapper(exec)), decltype(A), decltype(x), decltype(y)>::value;
 
+  static unsigned int recursion_count = 0;
+  ++recursion_count;
+  assert(recursion_count == 1u);
+
+  std::cerr << "use_custom: " << use_custom << std::endl;
   if constexpr(use_custom) {
+    static_assert(not std::is_same_v<
+      std::remove_cvref_t<decltype(execpolicy_mapper(exec))>,
+      std::remove_cvref_t<ExecutionPolicy>>);
     matrix_vector_product(execpolicy_mapper(exec), A, x, y);
   } else {
     matrix_vector_product(std::experimental::linalg::impl::inline_exec_t(), A, x, y);
@@ -335,7 +384,9 @@ void matrix_vector_product(
   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)
 {
+  std::cerr << "matrix_vector_product with default_exec_t\n";
   matrix_vector_product(std::experimental::linalg::impl::default_exec_t(), A, x, y);
+  std::cerr << "exiting matrix_vector_product with default_exec_t\n";
 }
 
 namespace impl {
@@ -389,6 +440,7 @@ void matrix_vector_product(
   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)
 {
+  std::cerr << "matrix_vector_product with inline_exec_t\n";
   using size_type = std::common_type_t<
     std::common_type_t<
       std::common_type_t<typename Extents_A::size_type /* SizeType_A */, SizeType_x>,
@@ -400,6 +452,7 @@ void matrix_vector_product(
       z(i) += A(i,j) * x(j);
     }
   }
+  std::cerr << "exiting matrix_vector_product with inline_exec_t\n";
 }
 
 // FIXME (mfh 2022/06/19) Some work-around here for GCC 9 and/or macro insufficiencies.
@@ -434,15 +487,24 @@ void matrix_vector_product(
   std::experimental::mdspan<ElementType_y, std::experimental::extents<SizeType_y, ext_y>, Layout_y, Accessor_y> y,
   std::experimental::mdspan<ElementType_z, Extents_z /* std::experimental::extents<SizeType_z, ext_z> */ , Layout_z, Accessor_z> z)
 {
+  static size_t recursion_count = 0;
+  ++recursion_count;
+
+  assert(recursion_count == 1);
 
   constexpr bool use_custom = is_custom_mat_vec_product_with_update_avail<
     decltype(execpolicy_mapper(exec)), decltype(A), decltype(x), decltype(y), decltype(z)>::value;
 
+  std::cerr << "use_custom: " << use_custom << std::endl;
   if constexpr(use_custom) {
+    static_assert(not std::is_same_v<
+      std::remove_cvref_t<decltype(execpolicy_mapper(exec))>,
+      std::remove_cvref_t<ExecutionPolicy>>);
     matrix_vector_product(execpolicy_mapper(exec), A, x, y, z);
   } else {
     matrix_vector_product(std::experimental::linalg::impl::inline_exec_t(), A, x, y, z);
   }
+  std::cerr << "exiting matrix_vector_product\n";
 }
 
 // FIXME (mfh 2022/06/19) Some work-around here for GCC 9 and/or macro insufficiencies.
@@ -474,7 +536,9 @@ void matrix_vector_product(
   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)
 {
+  std::cerr << "matrix_vector_product with default_exec_t\n";
   matrix_vector_product(std::experimental::linalg::impl::default_exec_t(), A, x, y, z);
+  std::cerr << "exiting matrix_vector_product with default_exec_t\n";
 }
 
 

include/experimental/__p1673_bits/linalg_config.h.in

@@ -5,3 +5,4 @@
 #cmakedefine LINALG_ENABLE_CONCEPTS
 #cmakedefine LINALG_ENABLE_KOKKOS
 #cmakedefine LINALG_ENABLE_KOKKOS_DEFAULT
+#cmakedefine LINALG_ENABLE_TBB

include/experimental/__p1673_bits/linalg_execpolicy_mapper.hpp

@@ -1,9 +1,10 @@
 // For GCC 9 (< 10), <execution> unconditionally includes a TBB header file.
 // If GCC < 10 was not built with TBB support, this causes a build error.
-#if (! defined(__GNUC__)) || (__GNUC__ > 9)
+#ifdef LINALG_HAS_EXECUTION
 #include <execution>
 #endif
 
+#include <iostream>
 #include <type_traits>
 
 namespace std {
@@ -39,7 +40,19 @@ namespace experimental {
 inline namespace __p1673_version_0 {
 namespace linalg {
 template<class T>
-auto execpolicy_mapper(T) { return std::experimental::linalg::impl::inline_exec_t(); }
+auto execpolicy_mapper(const T) {
+  std::cerr << "execpolicy_mapper returning inline_exec_t\n";
+  return std::experimental::linalg::impl::inline_exec_t(); 
+}
+
+#ifdef LINALG_HAS_EXECUTION
+template<>
+auto execpolicy_mapper(const std::execution::parallel_policy) {
+  std::cerr << "execpolicy_mapper returning par\n";
+  return std::execution::par;
+}
+#endif
+
 }
 }
 }

include/experimental/__p1673_bits/macros.hpp

@@ -43,6 +43,30 @@
 #ifndef LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_MACROS_HPP_
 #define LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_MACROS_HPP_
 
+#ifdef _MSVC_LANG
+#define _LINALG_CPLUSPLUS _MSVC_LANG
+#else
+#define _LINALG_CPLUSPLUS __cplusplus
+#endif
+
+#define LINALG_CXX_STD_14 201402L
+#define LINALG_CXX_STD_17 201703L
+#define LINALG_CXX_STD_20 202002L
+
+#define LINALG_HAS_CXX_14 (_LINALG_CPLUSPLUS >= LINALG_CXX_STD_14)
+#define LINALG_HAS_CXX_17 (_LINALG_CPLUSPLUS >= LINALG_CXX_STD_17)
+#define LINALG_HAS_CXX_20 (_LINALG_CPLUSPLUS >= LINALG_CXX_STD_20)
+
+static_assert(_LINALG_CPLUSPLUS >= LINALG_CXX_STD_17, "stdBLAS requires C++17 or later.");
+
+#if ! defined(__clang__) && ! defined(_MSC_VER) && defined(__GNUC__)
+#  if defined(LINALG_ENABLE_TBB)
+#    define LINALG_HAS_EXECUTION 1
+#  endif
+#else
+#  define LINALG_HAS_EXECUTION 1
+#endif
+
 #define P1673_MATRIX_EXTENTS_TEMPLATE_PARAMETERS( MATRIX_NAME ) \
   class SizeType_ ## MATRIX_NAME , \
   ::std::size_t numRows_ ## MATRIX_NAME , \

tests/native/abs_sum.cpp

@@ -1,3 +1,5 @@
+#define MDSPAN_USE_PAREN_OPERATOR 1
+
 #include "gtest/gtest.h"
 #include "gtest_fixtures.hpp"
 

tests/native/add.cpp

@@ -1,3 +1,5 @@
+#define MDSPAN_USE_PAREN_OPERATOR 1
+
 #include "gtest/gtest.h"
 
 #include <experimental/linalg>
@@ -6,9 +8,9 @@
 #include <vector>
 
 namespace {
-  using std::experimental::dynamic_extent;
-  using std::experimental::extents;
-  using std::experimental::mdspan;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::extents;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan;
   using std::experimental::linalg::add;
 
   TEST(BLAS1_add, vector_double)

tests/native/conjugate_transposed.cpp

@@ -1,3 +1,5 @@
+#define MDSPAN_USE_PAREN_OPERATOR 1
+
 #include "gtest/gtest.h"
 
 #include <experimental/linalg>
@@ -6,9 +8,9 @@
 #include <vector>
 
 namespace {
-  using std::experimental::dynamic_extent;
-  using std::experimental::extents;
-  using std::experimental::mdspan;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::extents;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan;
   using std::experimental::linalg::conjugate_transposed;
 
   TEST(conjugate_transposed, mdspan_complex_double)

tests/native/conjugated.cpp

@@ -1,3 +1,6 @@
+
+#define MDSPAN_USE_PAREN_OPERATOR 1
+
 #include "gtest/gtest.h"
 
 #include <experimental/linalg>
@@ -6,10 +9,10 @@
 #include <vector>
 
 namespace {
-  using std::experimental::dextents;
-  using std::experimental::dynamic_extent;
-  using std::experimental::extents;
-  using std::experimental::mdspan;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::dextents;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::extents;
+  using MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan;
   using std::experimental::linalg::conjugated;
 
   template<class ValueType>