Added language standard parallelism to dot product

examples/01_scale.cpp

@@ -26,7 +26,7 @@ int main(int argc, char* argv[]) {
     // With CTAD working we could do, GCC 11.1 works but some others are buggy
     // mdspan x(x_vec.data(), N);
     mdspan<double, extents<std::size_t, dynamic_extent>> x(x_vec.data(),N);
-    for(int i=0; i<x.extent(0); i++) x(i) = i;
+    for(int i=0; i<x.extent(0); i++) x[i] = i;
 
     // Call linalg::scale x = 2.0*x;
     std::experimental::linalg::scale(2.0, x);
@@ -36,6 +36,6 @@ int main(int argc, char* argv[]) {
     std::experimental::linalg::scale(2.0, x);
 #endif
 
-    for(int i=0; i<x.extent(0); i+=5) std::cout << i << " " << x(i) << std::endl;
+    for(int i=0; i<x.extent(0); i+=5) std::cout << i << " " << x[i] << std::endl;
   }
 }

examples/02_matrix_vector_product_basic.cpp

@@ -31,11 +31,11 @@ int main(int argc, char* argv[]) {
     mdspan<double, extents<std::size_t, dynamic_extent>> y(y_vec.data(),N);
     for(int i=0; i<A.extent(0); i++)
       for(int j=0; j<A.extent(1); j++)
-        A(i,j) = 100.0*i+j;
+        A[i,j] = 100.0*i+j;
     for(int i=0; i<x.extent(0); i++)
-      x(i) = 1. * i;
+      x[i] = 1. * i;
     for(int i=0; i<y.extent(0); i++)
-      y(i) = -1. * i;
+      y[i] = -1. * i;
 
     // y = A * x
     std::experimental::linalg::matrix_vector_product(A, x, y);
@@ -50,6 +50,6 @@ int main(int argc, char* argv[]) {
        std::experimental::linalg::scaled(2.0, A), x,
        std::experimental::linalg::scaled(0.5, y), y);
 #endif
-    for(int i=0; i<y.extent(0); i+=5) std::cout << i << " " << y(i) << std::endl;
+    for(int i=0; i<y.extent(0); i+=5) std::cout << i << " " << y[i] << std::endl;
   }
 }

examples/03_matrix_vector_product_mixedprec.cpp

@@ -25,13 +25,13 @@ int main(int argc, char* argv[]) {
     for(int m=0; m<A.extent(0); m++)
       for(int i=0; i<A.extent(1); i++)
         for(int j=0; j<A.extent(2); j++)
-        A(m,i,j) = 1000.0 * m + 100.0 * i + j;
+        A[m,i,j] = 1000.0 * m + 100.0 * i + j;
     for(int i=0; i<x.extent(0); i++)
       for(int m=0; m<x.extent(1); m++)
-        x(i,m) = 33. * i + 0.33 * m;
+        x[i,m] = 33. * i + 0.33 * m;
     for(int m=0; m<y.extent(0); m++)
       for(int i=0; i<y.extent(1); i++)
-        y(m,i) = 33. * m + 0.33 * i;
+        y[m,i] = 33. * m + 0.33 * i;
 
     for(int m = 0; m < M; m++) {
       auto A_m = submdspan(A, m, full_extent, full_extent);
@@ -41,7 +41,7 @@ int main(int argc, char* argv[]) {
       std::experimental::linalg::matrix_vector_product(A_m, x_m, y_m);
     }
 
-    for(int i=0; i<y.extent(0); i+=5) std::cout << i << " " << y(i,1) << std::endl;
+    for(int i=0; i<y.extent(0); i+=5) std::cout << i << " " << y[i,1] << std::endl;
   }
 }
 

include/experimental/__p1673_bits/blas1_dot.hpp

@@ -226,7 +226,7 @@ namespace dot_detail {
   auto dot_return_type_deducer(
     std::experimental::mdspan<ElementType1, std::experimental::extents<SizeType1, ext1>, Layout1, Accessor1> x,
     std::experimental::mdspan<ElementType2, std::experimental::extents<SizeType2, ext2>, Layout2, Accessor2> y)
-  -> decltype(x(0) * y(0));
+  -> decltype(x[0] * y[0]);
 } // namespace dot_detail
 
 

include/experimental/__p1673_bits/blas1_givens.hpp

@@ -399,9 +399,9 @@ void givens_rotation_apply(
   using index_type = ::std::common_type_t<SizeType1, SizeType2>;
   const auto x_extent_0 = static_cast<index_type>(x.extent(0));
   for (index_type i = 0; i < x_extent_0; ++i) {
-    const auto dtemp = c * x(i) + s * y(i);
-    y(i) = c * y(i) - s * x(i);
-    x(i) = dtemp;
+    const auto dtemp = c * x[i] + s * y[i];
+    y[i] = c * y[i] - s * x[i];
+    x[i] = dtemp;
   }
 }
 
@@ -496,9 +496,9 @@ void givens_rotation_apply(
   using index_type = ::std::common_type_t<SizeType1, SizeType2>;
   const auto x_extent_0 = static_cast<index_type>(x.extent(0));
   for (index_type i = 0; i < x_extent_0; ++i) {
-    const auto dtemp = c * x(i) + s * y(i);
-    y(i) = c * y(i) - conj(s) * x(i);
-    x(i) = dtemp;
+    const auto dtemp = c * x[i] + s * y[i];
+    y[i] = c * y[i] - conj(s) * x[i];
+    x[i] = dtemp;
   }
 }
 

include/experimental/__p1673_bits/blas1_linalg_add.hpp

@@ -82,7 +82,7 @@ void add_rank_1(
 
   using size_type = std::common_type_t<SizeType_x, SizeType_y, SizeType_z>;
   for (size_type i = 0; i < z.extent(0); ++i) {
-    z(i) = x(i) + y(i);
+    z[i] = x[i] + y[i];
   }
 }
 
@@ -132,7 +132,7 @@ void add_rank_2(
   using size_type = std::common_type_t<SizeType_x, SizeType_y, SizeType_z>;
   for (size_type j = 0; j < x.extent(1); ++j) {
     for (size_type i = 0; i < x.extent(0); ++i) {
-      z(i,j) = x(i,j) + y(i,j);
+      z[i,j] = x[i,j] + y[i,j];
     }
   }
 }

include/experimental/__p1673_bits/blas1_linalg_copy.hpp

@@ -69,7 +69,7 @@ void copy_rank_1(
                 x.static_extent(0) == y.static_extent(0));
   using size_type = std::common_type_t<SizeType_x, SizeType_y>;
   for (size_type i = 0; i < y.extent(0); ++i) {
-    y(i) = x(i);
+    y[i] = x[i];
   }
 }
 
@@ -98,7 +98,7 @@ void copy_rank_2(
   using size_type = std::common_type_t<SizeType_x, SizeType_y>;
   for (size_type j = 0; j < y.extent(1); ++j) {
     for (size_type i = 0; i < y.extent(0); ++i) {
-      y(i,j) = x(i,j);
+      y[i,j] = x[i,j];
     }
   }
 }

include/experimental/__p1673_bits/blas1_linalg_swap.hpp

@@ -74,7 +74,7 @@ void swap_rank_1(
   using size_type = std::common_type_t<SizeType_x, SizeType_y>;
 
   for (size_type i = 0; i < y.extent(0); ++i) {
-    swap(x(i), y(i));
+    swap(x[i], y[i]);
   }
 }
 
@@ -106,7 +106,7 @@ void swap_rank_2(
 
   for (size_type j = 0; j < y.extent(1); ++j) {
     for (size_type i = 0; i < y.extent(0); ++i) {
-      swap(x(i,j), y(i,j));
+      swap(x[i,j], y[i,j]);
     }
   }
 }

include/experimental/__p1673_bits/blas1_matrix_inf_norm.hpp

@@ -100,14 +100,14 @@ Scalar matrix_inf_norm(
     return result;
   }
   else if(A.extent(0) == size_type(1) && A.extent(1) == size_type(1)) {
-    result += abs(A(0, 0));
+    result += abs(A[0, 0]);
     return result;
   }
 
   for (size_type i = 0; i < A.extent(0); ++i) {
     auto row_sum = init;
     for (size_type j = 0; j < A.extent(1); ++j) {
-      row_sum += abs(A(i,j));
+      row_sum += abs(A[i,j]);
     }
     result = max(row_sum, result);
   }
@@ -170,7 +170,7 @@ namespace matrix_inf_norm_detail {
     class Layout,
     class Accessor>
   auto matrix_inf_norm_return_type_deducer(
-    std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, numRows, numCols>, Layout, Accessor> A) -> decltype(abs(A(0,0)));
+    std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, numRows, numCols>, Layout, Accessor> A) -> decltype(abs(A[0,0]));
 
 } // namespace matrix_inf_norm_detail
 

include/experimental/__p1673_bits/blas1_matrix_one_norm.hpp

@@ -101,7 +101,7 @@ Scalar matrix_one_norm(
     return result;
   }
   else if(A.extent(0) == size_type(1) && A.extent(1) == size_type(1)) {
-    result += abs(A(0, 0));
+    result += abs(A[0, 0]);
     return result;
   }
 
@@ -110,7 +110,7 @@ Scalar matrix_one_norm(
   for (size_type j = 0; j < A.extent(1); ++j) {
     auto col_sum = init;
     for (size_type i = 0; i < A.extent(0); ++i) {
-      col_sum += abs(A(i,j));
+      col_sum += abs(A[i,j]);
     }
     result = max(col_sum, result);
   }
@@ -171,7 +171,7 @@ namespace matrix_one_norm_detail {
     class Layout,
     class Accessor>
   auto matrix_one_norm_return_type_deducer(
-    std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, numRows, numCols>, Layout, Accessor> A) -> decltype(abs(A(0,0)));
+    std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, numRows, numCols>, Layout, Accessor> A) -> decltype(abs(A[0,0]));
 
 } // namespace matrix_one_norm_detail
 

include/experimental/__p1673_bits/blas1_scale.hpp

@@ -61,7 +61,7 @@ void linalg_scale_rank_1(
   std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
 {
   for (SizeType i = 0; i < x.extent(0); ++i) {
-    x(i) *= alpha;
+    x[i] *= alpha;
   }
 }
 
@@ -78,7 +78,7 @@ void linalg_scale_rank_2(
 {
   for (SizeType j = 0; j < A.extent(1); ++j) {
     for (SizeType i = 0; i < A.extent(0); ++i) {
-      A(i,j) *= alpha;
+      A[i,j] *= alpha;
     }
   }
 }

include/experimental/__p1673_bits/blas1_vector_abs_sum.hpp

@@ -87,7 +87,7 @@ Scalar vector_abs_sum(
   const SizeType numElt = v.extent(0);
   for (SizeType i = 0; i < numElt; ++i) {
     using std::abs;
-    init += abs(v(i));
+    init += abs(v[i]);
   }
   return init;
 }
@@ -140,7 +140,7 @@ namespace vector_abs_detail {
     class Accessor>
   auto vector_abs_return_type_deducer(
     std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
-  -> decltype(abs(x(0)));
+  -> decltype(abs(x[0]));
 } // namespace vector_abs_detail
 
 

include/experimental/__p1673_bits/blas1_vector_idx_abs_max.hpp

@@ -76,17 +76,17 @@ SizeType idx_abs_max_default_impl(
   std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> v)
 {
   using std::abs;
-  using magnitude_type = decltype(abs(v(0)));
+  using magnitude_type = decltype(abs(v[0]));
 
   if (v.extent(0) == 0) {
     return std::numeric_limits<SizeType>::max();
   }
 
   SizeType maxInd = 0;
-  magnitude_type maxVal = abs(v(0));
+  magnitude_type maxVal = abs(v[0]);
   for (SizeType i = 1; i < v.extent(0); ++i) {
-    if (maxVal < abs(v(i))) {
-      maxVal = abs(v(i));
+    if (maxVal < abs(v[i])) {
+      maxVal = abs(v[i]);
       maxInd = i;
     }
   }

include/experimental/__p1673_bits/blas1_vector_norm2.hpp

@@ -82,7 +82,7 @@ template<class ElementType,
          class Scalar>
 Scalar vector_norm2(
   std::experimental::linalg::impl::inline_exec_t&& exec,
-  std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x,
+  mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x,
   Scalar init)
 {
   // Initialize the sum of squares result
@@ -107,7 +107,7 @@ template<class ExecutionPolicy,
          class Scalar>
 Scalar vector_norm2(
   ExecutionPolicy&& exec,
-  std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x,
+  mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x,
   Scalar init)
 {
   constexpr bool use_custom = is_custom_vector_norm2_avail<
@@ -129,7 +129,7 @@ template<class ElementType,
          class Accessor,
          class Scalar>
 Scalar vector_norm2(
-  std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x,
+  mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x,
   Scalar init)
 {
   return vector_norm2(std::experimental::linalg::impl::default_exec_t(), x, init);
@@ -147,16 +147,16 @@ namespace vector_norm2_detail {
     class Layout,
     class Accessor>
   auto vector_norm2_return_type_deducer(
-    std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
-  -> decltype(abs(x(0)) * abs(x(0)));
+    mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
+  -> decltype(abs(x[0]) * abs(x[0]));
 } // namespace vector_norm2_detail
 
 template<class ElementType,
          class SizeType, ::std::size_t ext0,
          class Layout,
          class Accessor>
 auto vector_norm2(
-  std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
+  mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
 -> decltype(vector_norm2_detail::vector_norm2_return_type_deducer(x))
 {
   using return_t = decltype(vector_norm2_detail::vector_norm2_return_type_deducer(x));
@@ -170,7 +170,7 @@ template<class ExecutionPolicy,
          class Accessor>
 auto vector_norm2(
   ExecutionPolicy&& exec,
-  std::experimental::mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
+  mdspan<ElementType, std::experimental::extents<SizeType, ext0>, Layout, Accessor> x)
 -> decltype(vector_norm2_detail::vector_norm2_return_type_deducer(x))
 {
   using return_t = decltype(vector_norm2_detail::vector_norm2_return_type_deducer(x));

include/experimental/__p1673_bits/blas1_vector_sum_of_squares.hpp

@@ -105,8 +105,8 @@ sum_of_squares_result<Scalar> vector_sum_of_squares(
   Scalar scale = init.scaling_factor;
   Scalar ssq = init.scaled_sum_of_squares;
   for (SizeType i = 0; i < x.extent(0); ++i) {
-    if (abs(x(i)) != 0.0) {
-      const auto absxi = abs(x(i));
+    if (abs(x[i]) != 0.0) {
+      const auto absxi = abs(x[i]);
       const auto quotient = scale / absxi;
       if (scale < absxi) {
           ssq = Scalar(1.0) + ssq * quotient * quotient;