Remove definition MDSPAN_USE_PAREN_OPERATOR=1 from DDC target (#677)

* Remove definition MDSPAN_USE_PAREN_OPERATOR=1 from DDC target

* Alternative implementation to workaround llvm <17

* Exclude hip oldest with C++ 23

* Add a comment on the logic

.github/workflows/tests.yml

@@ -122,6 +122,9 @@ jobs:
         - image: 'oldest' # nvcc 11 only supports C++-17
           backend: 'cuda'
           cxx_version: '23'
+        - image: 'oldest' # clang bug with multidimensional bracket operator and parameter packs
+          backend: 'hip'
+          cxx_version: '23'
         - image: 'latest' # nvcc 12 only supports C++-20
           backend: 'cuda'
           cxx_version: '23'

CMakeLists.txt

@@ -131,7 +131,6 @@ target_include_directories(
         "$<INSTALL_INTERFACE:include>"
 )
 target_link_libraries(DDC INTERFACE Kokkos::kokkos)
-target_compile_definitions(DDC INTERFACE MDSPAN_USE_PAREN_OPERATOR=1)
 if("${DDC_BUILD_DOUBLE_PRECISION}")
     target_compile_definitions(DDC INTERFACE DDC_BUILD_DOUBLE_PRECISION)
 endif()

include/ddc/chunk.hpp

@@ -13,6 +13,7 @@
 #include "ddc/chunk_common.hpp"
 #include "ddc/chunk_span.hpp"
 #include "ddc/chunk_traits.hpp"
+#include "ddc/detail/kokkos.hpp"
 #include "ddc/kokkos_allocator.hpp"
 
 namespace ddc {
@@ -192,7 +193,7 @@ class Chunk<ElementType, DiscreteDomain<DDims...>, Allocator>
         static_assert((is_discrete_element_v<DElems> && ...), "Expected DiscreteElements");
         assert(((select<DDims>(take<DDims>(delems...)) >= front<DDims>(this->m_domain)) && ...));
         assert(((select<DDims>(take<DDims>(delems...)) <= back<DDims>(this->m_domain)) && ...));
-        return this->m_internal_mdspan(uid<DDims>(take<DDims>(delems...))...);
+        return DDC_MDSPAN_ACCESS_OP(this->m_internal_mdspan, uid<DDims>(take<DDims>(delems...))...);
     }
 
     /** Element access using a list of DiscreteElement
@@ -208,7 +209,7 @@ class Chunk<ElementType, DiscreteDomain<DDims...>, Allocator>
         static_assert((is_discrete_element_v<DElems> && ...), "Expected DiscreteElements");
         assert(((select<DDims>(take<DDims>(delems...)) >= front<DDims>(this->m_domain)) && ...));
         assert(((select<DDims>(take<DDims>(delems...)) <= back<DDims>(this->m_domain)) && ...));
-        return this->m_internal_mdspan(uid<DDims>(take<DDims>(delems...))...);
+        return DDC_MDSPAN_ACCESS_OP(this->m_internal_mdspan, uid<DDims>(take<DDims>(delems...))...);
     }
 
     /** Returns the label of the Chunk

include/ddc/chunk_span.hpp

@@ -330,7 +330,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
         static_assert((is_discrete_element_v<DElems> && ...), "Expected DiscreteElements");
         assert(((select<DDims>(take<DDims>(delems...)) >= front<DDims>(this->m_domain)) && ...));
         assert(((select<DDims>(take<DDims>(delems...)) <= back<DDims>(this->m_domain)) && ...));
-        return this->m_internal_mdspan(uid<DDims>(take<DDims>(delems...))...);
+        return DDC_MDSPAN_ACCESS_OP(this->m_internal_mdspan, uid<DDims>(take<DDims>(delems...))...);
     }
 
     /** Access to the underlying allocation pointer

include/ddc/detail/kokkos.hpp

@@ -13,6 +13,20 @@
 
 #include "macros.hpp"
 
+// if MDSPAN_USE_BRACKET_OPERATOR is defined then we are likely
+// using the Kokkos implementation of mdspan
+#if defined(MDSPAN_USE_BRACKET_OPERATOR)
+#if MDSPAN_USE_BRACKET_OPERATOR
+#define DDC_MDSPAN_ACCESS_OP(mds, ...) mds[__VA_ARGS__]
+#else
+#define DDC_MDSPAN_ACCESS_OP(mds, ...) mds(__VA_ARGS__)
+#endif
+// else we are likely using an other implementation
+// then we can only use the bracket operator
+#else
+#define DDC_MDSPAN_ACCESS_OP(mds, ...) mds[__VA_ARGS__]
+#endif
+
 namespace ddc::detail {
 
 template <class T>
@@ -133,8 +147,8 @@ KOKKOS_FUNCTION mdspan_to_kokkos_layout_t<typename MP::layout_type> build_kokkos
                 std::size_t,
                 Kokkos::extents<std::size_t, sizeof...(Is), 2>,
                 Kokkos::layout_right> const interleaved_extents_strides(storage.data());
-        ((interleaved_extents_strides(Is, 0) = ep.extent(Is),
-          interleaved_extents_strides(Is, 1) = mapping.stride(Is)),
+        ((DDC_MDSPAN_ACCESS_OP(interleaved_extents_strides, Is, 0) = ep.extent(Is),
+          DDC_MDSPAN_ACCESS_OP(interleaved_extents_strides, Is, 1) = mapping.stride(Is)),
          ...);
         return make_layout_stride(storage, std::make_index_sequence<sizeof...(Is) * 2> {});
     } else {

include/ddc/kernels/splines/bsplines_non_uniform.hpp

@@ -600,59 +600,63 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDim> NonUniformBSplines<CDim, D>::
 
     double saved;
     double temp;
-    ndu(0, 0) = 1.0;
+    DDC_MDSPAN_ACCESS_OP(ndu, 0, 0) = 1.0;
     for (std::size_t j = 0; j < degree(); ++j) {
         left[j] = x - ddc::coordinate(icell - j);
         right[j] = ddc::coordinate(icell + j + 1) - x;
         saved = 0.0;
         for (std::size_t r = 0; r < j + 1; ++r) {
             // compute inverse of knot differences and save them into lower
             // triangular part of ndu
-            ndu(r, j + 1) = 1.0 / (right[r] + left[j - r]);
+            DDC_MDSPAN_ACCESS_OP(ndu, r, j + 1) = 1.0 / (right[r] + left[j - r]);
             // compute basis functions and save them into upper triangular part
             // of ndu
-            temp = ndu(j, r) * ndu(r, j + 1);
-            ndu(j + 1, r) = saved + right[r] * temp;
+            temp = DDC_MDSPAN_ACCESS_OP(ndu, j, r) * DDC_MDSPAN_ACCESS_OP(ndu, r, j + 1);
+            DDC_MDSPAN_ACCESS_OP(ndu, j + 1, r) = saved + right[r] * temp;
             saved = left[j - r] * temp;
         }
-        ndu(j + 1, j + 1) = saved;
+        DDC_MDSPAN_ACCESS_OP(ndu, j + 1, j + 1) = saved;
     }
     // Save 0-th derivative
     for (std::size_t j = 0; j < degree() + 1; ++j) {
-        derivs(j, 0) = ndu(degree(), j);
+        DDC_MDSPAN_ACCESS_OP(derivs, j, 0) = DDC_MDSPAN_ACCESS_OP(ndu, degree(), j);
     }
 
     for (int r = 0; r < int(degree() + 1); ++r) {
         int s1 = 0;
         int s2 = 1;
-        a(0, 0) = 1.0;
+        DDC_MDSPAN_ACCESS_OP(a, 0, 0) = 1.0;
         for (int k = 1; k < int(n + 1); ++k) {
             double d = 0.0;
             int const rk = r - k;
             int const pk = degree() - k;
             if (r >= k) {
-                a(0, s2) = a(0, s1) * ndu(rk, pk + 1);
-                d = a(0, s2) * ndu(pk, rk);
+                DDC_MDSPAN_ACCESS_OP(a, 0, s2)
+                        = DDC_MDSPAN_ACCESS_OP(a, 0, s1) * DDC_MDSPAN_ACCESS_OP(ndu, rk, pk + 1);
+                d = DDC_MDSPAN_ACCESS_OP(a, 0, s2) * DDC_MDSPAN_ACCESS_OP(ndu, pk, rk);
             }
             int const j1 = rk > -1 ? 1 : (-rk);
             int const j2 = (r - 1) <= pk ? k : (degree() - r + 1);
             for (int j = j1; j < j2; ++j) {
-                a(j, s2) = (a(j, s1) - a(j - 1, s1)) * ndu(rk + j, pk + 1);
-                d += a(j, s2) * ndu(pk, rk + j);
+                DDC_MDSPAN_ACCESS_OP(a, j, s2)
+                        = (DDC_MDSPAN_ACCESS_OP(a, j, s1) - DDC_MDSPAN_ACCESS_OP(a, j - 1, s1))
+                          * DDC_MDSPAN_ACCESS_OP(ndu, rk + j, pk + 1);
+                d += DDC_MDSPAN_ACCESS_OP(a, j, s2) * DDC_MDSPAN_ACCESS_OP(ndu, pk, rk + j);
             }
             if (r <= pk) {
-                a(k, s2) = -a(k - 1, s1) * ndu(r, pk + 1);
-                d += a(k, s2) * ndu(pk, r);
+                DDC_MDSPAN_ACCESS_OP(a, k, s2) = -DDC_MDSPAN_ACCESS_OP(a, k - 1, s1)
+                                                 * DDC_MDSPAN_ACCESS_OP(ndu, r, pk + 1);
+                d += DDC_MDSPAN_ACCESS_OP(a, k, s2) * DDC_MDSPAN_ACCESS_OP(ndu, pk, r);
             }
-            derivs(r, k) = d;
+            DDC_MDSPAN_ACCESS_OP(derivs, r, k) = d;
             Kokkos::kokkos_swap(s1, s2);
         }
     }
 
     int r = degree();
     for (int k = 1; k < int(n + 1); ++k) {
         for (std::size_t i = 0; i < derivs.extent(0); ++i) {
-            derivs(i, k) *= r;
+            DDC_MDSPAN_ACCESS_OP(derivs, i, k) *= r;
         }
         r *= degree() - k;
     }

include/ddc/kernels/splines/bsplines_uniform.hpp

@@ -413,17 +413,17 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDim> UniformBSplines<CDim, D>::
     double xx;
     double temp;
     double saved;
-    values(0) = 1.0;
+    DDC_MDSPAN_ACCESS_OP(values, 0) = 1.0;
     for (std::size_t j = 1; j < values.size(); ++j) {
         xx = -offset;
         saved = 0.0;
         for (std::size_t r = 0; r < j; ++r) {
             xx += 1;
-            temp = values(r) / j;
-            values(r) = saved + xx * temp;
+            temp = DDC_MDSPAN_ACCESS_OP(values, r) / j;
+            DDC_MDSPAN_ACCESS_OP(values, r) = saved + xx * temp;
             saved = (j - xx) * temp;
         }
-        values(j) = saved;
+        DDC_MDSPAN_ACCESS_OP(values, j) = saved;
     }
 
     return discrete_element_type(jmin);
@@ -447,29 +447,29 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDim> UniformBSplines<CDim, D>::
     double xx;
     double temp;
     double saved;
-    derivs(0) = 1.0 / ddc::step<knot_discrete_dimension_type>();
+    DDC_MDSPAN_ACCESS_OP(derivs, 0) = 1.0 / ddc::step<knot_discrete_dimension_type>();
     for (std::size_t j = 1; j < degree(); ++j) {
         xx = -offset;
         saved = 0.0;
         for (std::size_t r = 0; r < j; ++r) {
             xx += 1.0;
-            temp = derivs(r) / j;
-            derivs(r) = saved + xx * temp;
+            temp = DDC_MDSPAN_ACCESS_OP(derivs, r) / j;
+            DDC_MDSPAN_ACCESS_OP(derivs, r) = saved + xx * temp;
             saved = (j - xx) * temp;
         }
-        derivs(j) = saved;
+        DDC_MDSPAN_ACCESS_OP(derivs, j) = saved;
     }
 
     // Compute derivatives
     double bjm1 = derivs[0];
     double bj = bjm1;
-    derivs(0) = -bjm1;
+    DDC_MDSPAN_ACCESS_OP(derivs, 0) = -bjm1;
     for (std::size_t j = 1; j < degree(); ++j) {
-        bj = derivs(j);
-        derivs(j) = bjm1 - bj;
+        bj = DDC_MDSPAN_ACCESS_OP(derivs, j);
+        DDC_MDSPAN_ACCESS_OP(derivs, j) = bjm1 - bj;
         bjm1 = bj;
     }
-    derivs(degree()) = bj;
+    DDC_MDSPAN_ACCESS_OP(derivs, degree()) = bj;
 
     return discrete_element_type(jmin);
 }
@@ -507,45 +507,47 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDim> UniformBSplines<CDim, D>::
     double xx;
     double temp;
     double saved;
-    ndu(0, 0) = 1.0;
+    DDC_MDSPAN_ACCESS_OP(ndu, 0, 0) = 1.0;
     for (std::size_t j = 1; j < degree() + 1; ++j) {
         xx = -offset;
         saved = 0.0;
         for (std::size_t r = 0; r < j; ++r) {
             xx += 1.0;
-            temp = ndu(j - 1, r) / j;
-            ndu(j, r) = saved + xx * temp;
+            temp = DDC_MDSPAN_ACCESS_OP(ndu, j - 1, r) / j;
+            DDC_MDSPAN_ACCESS_OP(ndu, j, r) = saved + xx * temp;
             saved = (j - xx) * temp;
         }
-        ndu(j, j) = saved;
+        DDC_MDSPAN_ACCESS_OP(ndu, j, j) = saved;
     }
     for (std::size_t i = 0; i < ndu.extent(1); ++i) {
-        derivs(i, 0) = ndu(degree(), i);
+        DDC_MDSPAN_ACCESS_OP(derivs, i, 0) = DDC_MDSPAN_ACCESS_OP(ndu, degree(), i);
     }
 
     for (int r = 0; r < int(degree() + 1); ++r) {
         int s1 = 0;
         int s2 = 1;
-        a(0, 0) = 1.0;
+        DDC_MDSPAN_ACCESS_OP(a, 0, 0) = 1.0;
         for (int k = 1; k < int(n + 1); ++k) {
             double d = 0.0;
             int const rk = r - k;
             int const pk = degree() - k;
             if (r >= k) {
-                a(0, s2) = a(0, s1) / (pk + 1);
-                d = a(0, s2) * ndu(pk, rk);
+                DDC_MDSPAN_ACCESS_OP(a, 0, s2) = DDC_MDSPAN_ACCESS_OP(a, 0, s1) / (pk + 1);
+                d = DDC_MDSPAN_ACCESS_OP(a, 0, s2) * DDC_MDSPAN_ACCESS_OP(ndu, pk, rk);
             }
             int const j1 = rk > -1 ? 1 : (-rk);
             int const j2 = (r - 1) <= pk ? k : (degree() - r + 1);
             for (int j = j1; j < j2; ++j) {
-                a(j, s2) = (a(j, s1) - a(j - 1, s1)) / (pk + 1);
-                d += a(j, s2) * ndu(pk, rk + j);
+                DDC_MDSPAN_ACCESS_OP(a, j, s2)
+                        = (DDC_MDSPAN_ACCESS_OP(a, j, s1) - DDC_MDSPAN_ACCESS_OP(a, j - 1, s1))
+                          / (pk + 1);
+                d += DDC_MDSPAN_ACCESS_OP(a, j, s2) * DDC_MDSPAN_ACCESS_OP(ndu, pk, rk + j);
             }
             if (r <= pk) {
-                a(k, s2) = -a(k - 1, s1) / (pk + 1);
-                d += a(k, s2) * ndu(pk, r);
+                DDC_MDSPAN_ACCESS_OP(a, k, s2) = -DDC_MDSPAN_ACCESS_OP(a, k - 1, s1) / (pk + 1);
+                d += DDC_MDSPAN_ACCESS_OP(a, k, s2) * DDC_MDSPAN_ACCESS_OP(ndu, pk, r);
             }
-            derivs(r, k) = d;
+            DDC_MDSPAN_ACCESS_OP(derivs, r, k) = d;
             Kokkos::kokkos_swap(s1, s2);
         }
     }
@@ -557,7 +559,7 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDim> UniformBSplines<CDim, D>::
     double d = degree() * inv_dx;
     for (int k = 1; k < int(n + 1); ++k) {
         for (std::size_t i = 0; i < derivs.extent(0); ++i) {
-            derivs(i, k) *= d;
+            DDC_MDSPAN_ACCESS_OP(derivs, i, k) *= d;
         }
         d *= (degree() - k) * inv_dx;
     }

include/ddc/kernels/splines/spline_builder.hpp

@@ -500,7 +500,10 @@ int SplineBuilder<
             offset = jmin.uid() - start.uid() + bsplines_type::degree() / 2 - BSplines::degree();
         } else {
             int const mid = bsplines_type::degree() / 2;
-            offset = jmin.uid() - start.uid() + (values(mid) > values(mid + 1) ? mid : mid + 1)
+            offset = jmin.uid() - start.uid()
+                     + (DDC_MDSPAN_ACCESS_OP(values, mid) > DDC_MDSPAN_ACCESS_OP(values, mid + 1)
+                                ? mid
+                                : mid + 1)
                      - BSplines::degree();
         }
     } else {
@@ -677,14 +680,17 @@ void SplineBuilder<
         // all derivatives by multiplying the i-th derivative by dx^i
         for (std::size_t i = 0; i < bsplines_type::degree() + 1; ++i) {
             for (std::size_t j = 1; j < bsplines_type::degree() / 2 + 1; ++j) {
-                derivs(i, j) *= ddc::detail::ipow(m_dx, j);
+                DDC_MDSPAN_ACCESS_OP(derivs, i, j) *= ddc::detail::ipow(m_dx, j);
             }
         }
 
         // iterate only to deg as last bspline is 0
         for (std::size_t i = 0; i < s_nbc_xmin; ++i) {
             for (std::size_t j = 0; j < bsplines_type::degree(); ++j) {
-                matrix->set_element(i, j, derivs(j, s_nbc_xmin - i - 1 + s_odd));
+                matrix->set_element(
+                        i,
+                        j,
+                        DDC_MDSPAN_ACCESS_OP(derivs, j, s_nbc_xmin - i - 1 + s_odd));
             }
         }
     }
@@ -703,7 +709,7 @@ void SplineBuilder<
             int const j = ddc::detail::
                     modulo(int(jmin.uid() - m_offset + s),
                            static_cast<int>(ddc::discrete_space<BSplines>().nbasis()));
-            matrix->set_element(ix.uid() - start + s_nbc_xmin, j, values(s));
+            matrix->set_element(ix.uid() - start + s_nbc_xmin, j, DDC_MDSPAN_ACCESS_OP(values, s));
         }
     });
 
@@ -727,15 +733,15 @@ void SplineBuilder<
         // all derivatives by multiplying the i-th derivative by dx^i
         for (std::size_t i = 0; i < bsplines_type::degree() + 1; ++i) {
             for (std::size_t j = 1; j < bsplines_type::degree() / 2 + 1; ++j) {
-                derivs(i, j) *= ddc::detail::ipow(m_dx, j);
+                DDC_MDSPAN_ACCESS_OP(derivs, i, j) *= ddc::detail::ipow(m_dx, j);
             }
         }
 
         int const i0 = ddc::discrete_space<BSplines>().nbasis() - s_nbc_xmax;
         int const j0 = ddc::discrete_space<BSplines>().nbasis() - bsplines_type::degree();
         for (std::size_t j = 0; j < bsplines_type::degree(); ++j) {
             for (std::size_t i = 0; i < s_nbc_xmax; ++i) {
-                matrix->set_element(i0 + i, j0 + j, derivs(j + 1, i + s_odd));
+                matrix->set_element(i0 + i, j0 + j, DDC_MDSPAN_ACCESS_OP(derivs, j + 1, i + s_odd));
             }
         }
     }

tests/splines/bsplines.cpp

@@ -75,7 +75,7 @@ TYPED_TEST(BSplinesFixture, PartitionOfUnityUniform)
         ddc::discrete_space<BSplinesX>().eval_basis(values, test_point);
         double sum = 0.0;
         for (std::size_t j(0); j < degree + 1; ++j) {
-            sum += values(j);
+            sum += DDC_MDSPAN_ACCESS_OP(values, j);
         }
         EXPECT_LE(fabs(sum - 1.0), 1.0e-15);
     }
@@ -110,7 +110,7 @@ TYPED_TEST(BSplinesFixture, PartitionOfUnityNonUniform)
         ddc::discrete_space<BSplinesX>().eval_basis(values, test_point);
         double sum = 0.0;
         for (std::size_t j(0); j < degree + 1; ++j) {
-            sum += values(j);
+            sum += DDC_MDSPAN_ACCESS_OP(values, j);
         }
         EXPECT_LE(fabs(sum - 1.0), 1.0e-15);
     }