Make mdspan size_type controllable

Document #: P2553R1
Date: 2022-03-15
Project: Programming Language C++
LEWG
Reply-to: Christian Trott
<>
Damien Lebrun-Grandie
<>
Mark Hoemmen
<>
Daniel Sunderland
<>

1 Revision History

1.1 R1

1.2 Initial Version 2022-02 Mailing

1.2.1 Review in LEWG 2022-03-01

2 Description

P0009 explicitly sets the size_type of extents to size_t, which is then used by layout mappings and mdspan. While this matches span whose extent function returns size_t, this behavior has significant performance impact on various architectures where 64-bit integer throughput is significantly lower than 32-bit integer computation throughput.

While that problem is present for span it is a much bigger issue for mdspan, since every data access into a multi-dimensional array is accompanied by an offset calculation which typically costs two integer operations per dimension.

On current GPUs, which are an essential hardware component in machines used for high performance computing, machine learning and graphics (arguably the core constituency for mdspan), the ratio between 32-bit and 64-bit integer throughput is often 2x-4x.

2.1 Example

To gauge the impact, we investigated a simple stencil code benchmark, which is hosted in the mdspan reference implementation. That benchmark is using CUDA and compares a variant with raw pointers and explicit index calculation against a version which uses mdspan.

The mdspan variant does in each CUDA thread the following code:

for(size_t i = blockIdx.x+d; i < s.extent(0)-d; i += gridDim.x) {
  for(size_t j = threadIdx.z+d; j < s.extent(1)-d; j += blockDim.z) {
    for(size_t k = threadIdx.y+d; k < s.extent(2)-d; k += blockDim.y) {
      value_type sum_local = 0;
      for(size_t di = i-d; di < i+d+1; di++) {
      for(size_t dj = j-d; dj < j+d+1; dj++) {
      for(size_t dk = k-d; dk < k+d+1; dk++) {
        sum_local += s(di, dj, dk);
      }}}
      o(i,j,k) = sum_local;
    }
  }
}

The raw pointer variant looks like this:

for(size_t i = blockIdx.x+d; i < x-d; i += gridDim.x) {
  for(size_t j = threadIdx.z+d; j < y-d; j += blockDim.z) {
    for(size_t k = threadIdx.y+d; k < z-d; k += blockDim.y) {
      value_type sum_local = 0;
      for(size_t di = i-d; di < i+d+1; di++) {
      for(size_t dj = j-d; dj < j+d+1; dj++) {
      for(size_t dk = k-d; dk < k+d+1; dk++) {
        sum_local += data[dk + dj*z + di*z*y];
      }}}
      data_o[k + j*z + i*z*y] = sum_local;
    }
  }
}

Running the raw pointer variant with unsigned vs size_t as the loop indices, results in a timing of 31ms vs 56ms. The same is observed for the mdspan variant when switching in the mdspan implementation the size_type from size_t to unsigned. The 31ms result can also be obtained when leaving size_type as size_t but casting extents.extent(r) to the user provided index type inside the layout mappings index calculation operator while using unsigned as the loop index type in the algorithm.

2.2 Possible Ways To Address The Issue

2.2.1 Mappings Doing Offset Calculation With Argument Type

One way to address this issue would be for mappings to do all their internal calculations with the common_type of the user-provided indices. That includes casting extents.extent(i). However the drawback of this approach is that it is hard to identify overflows, which depend on layout as well.

// The following is ok, extents converts to size_t, required_span_size returns size_t too
mdspan<double,dextents<3>,layout_right> r(ptr,2000,2000,2000); 
mdspan<double,dextents<3>,layout_left>  l(ptr,2000,2000,2000);
...
r(1,1,1000) = 5; // ok
r(1000,1,1) = 5; // overflow
l(1,1,1000) = 5; // overflow
l(1000,1,1) = 5; // ok

In particular, in situations where allocations and mdspan creation happens in another code, location this could be an issue.

2.2.2 Template extents on size_type

In order to make overflow a better controllable artifact, and avoid accidental overflows we can make the index type part of the type. The natural place for this is extents. Every index calculation related piece of the mdspan proposal gets its size_type from extents, specifically both layout mappings, and mdspan itself is required to get its public size_type type member from extents_type::size_type. Furthermore, extents defines the complete iteration space for mdspan. Note, that a mapping might require a larger integer range that the product of extents (e.g. layout_stride::required_span_size can return a number larger than the product of its extents).

2.2.3 Second Extents Type Templated on size_type

Instead of modifying extents, we could introduce a new type basic_extents which is templated on the size type and the extents, but otherwise is identical to extents: When we can make anything in the mdspan proposal which accepts extents also accept basic_extents.

Potentially, extents could be just a template alias to basic_extents:

template<size_t ... Extents>
using extents = basic_extents<size_t, Extents...>;

Unfortunately that means that the following type of code would not work:

template<class T, class L, class A, size_t ... Extents>
void foo(mdspan<T,extents<Extents...>,L,A> a) {...{

However we believe the common use case would be to template on the extents object itself, mitigating this issue:

template<class T, class E, class L, class A>
void foo(mdspan<T,E,L,A> a) {...}

2.2.4 Template LayoutPolicy on size_type

We could also template the layout policy class on size_type, and use that type for the offset calculation, casting extents::extent explicitly on use. However this still means that the size of the object is larger (i.e. we would still store 64bit extents, instead of 32bit) and that additional casts will happen.

2.2.5 What we prefer:

All in all we prefer the option of making extents require the additional argument (2.2.2), with the next best thing being the introduction basic_extents and making extents an alias to basic_extents with size_t as the size_type. If LEWG would prefer the second option, the wording is largely the same with the following changes at the end:

template<size_t ... Extents>
using extents = basic_extents<size_t, Extents...>;

LEWG would need to decide whether to make dextents have a size_type template parameter or not.

2.3 Why we can’t fix this later

In principle we could add a second extents type later, though it may break code such as the one shown before (in the sense that it wouldn’t generally work for every instance of mdspan anymore):

template<class T, class L, class A, size_t ... Extents>
void foo(mdspan<T,extents<Extents...>,L,A> a) {...{

3 Editing Notes

The proposed changes are relative to the working draft of the standard as of P0009 R16.

4 Wording

4.1 In 22.7.X [mdspan.syn]

Replace:

template<size_t... Extents>
  class extents;

template<size_t Rank>
  using dextents = see below;

with:

template<class SizeT, size_t... Extents>
  class extents;

template<class SizeT, size_t Rank>
  using dextents = see below;

4.2 In 22.7.X.1 [mdspan.extents.overview]:

template<size_t... Extents>
class extents {
public:
  using size_type = size_t;
  using rank_type = size_t;

with

template<class SizeT, size_t... Extents>
class extents {
public:
  using size_type = SizeT;
  using rank_type = size_t;
  static constexpr size_type static_extent(rank_type) noexcept;

with

  static constexpr size_t static_extent(rank_type) noexcept;
  template<size_t... OtherExtents>
    explicit(see below)
    constexpr extents(const extents<OtherExtents...>&) noexcept;

with:

  template<class OtherSizeT, size_t... OtherExtents>
    explicit(see below)
    constexpr extents(const extents<OtherSizeT, OtherExtents...>&) noexcept;
  // [mdspan.extents.cmp], extents comparison operators
  template<size_t... OtherExtents>
    friend constexpr bool operator==(const extents&, const extents<OtherExtents...>&) noexcept;

with:

  // [mdspan.extents.cmp], extents comparison operators
  template<class OtherSizeT, size_t... OtherExtents>
    friend constexpr bool operator==(const extents&, const extents<OtherSizeT, OtherExtents...>&) noexcept;
  template<class... Integrals>
  static constexpr bool are-valid-extents(Integrals... values) noexcept; // exposition only

  template<class OtherSizeT, size_t... OtherExtents>
  static constexpr bool are-valid-extents(extents<OtherSizeT, OtherExtents>) noexcept; // exposition only

  template<class OtherSizeT, size_t N>
  static constexpr bool are-valid-extents(array<OtherSizeT, N>) noexcept; // exposition only

If SizeT is not an integral type other than bool, then the program is ill-formed.

4.3 In 22.7.X.2 [mdspan.extents.helpers]:

template<class... Integrals>
static constexpr bool are-valid-extents(Integrals... values) noexcept; // exposition only

Returns:

template<class OtherSizeT, size_t... OtherExtents>
static constexpr bool are-valid-extents(extents<OtherSizeT, OtherExtents> e) noexcept; // exposition only

Returns:

template<class OtherSizeT, size_t N>
static constexpr bool are-valid-extents(array<OtherSizeT, N> arr) noexcept; // exposition only

Returns:

4.4 In subsection 22.7.X.3 [mdspan.extents.cons]

template<size_t... OtherExtents>
  explicit((((Extents!=dynamic_extent) && (OtherExtents==dynamic_extent)) || ... ))
  constexpr extents(const extents<OtherExtents...>& other) noexcept;

to:

template<class OtherSizeT, size_t... OtherExtents>
  explicit(see below)
  constexpr extents(const extents<OtherSizeT, OtherExtents...>& other) noexcept;

4.5 In subsection 22.7.X.4 [mdspan.extents.obs]

with:

  static constexpr size_t static_extent(rank_type i) const noexcept;

4.6 In subsection 22.7.X.5 [mdspan.extents.compare]

with:

  template<class OtherSizeT, size_t... OtherExtents>
    friend constexpr bool operator==(const extents& lhs,
                                     const extents<OtherSizeT, OtherExtents...>& rhs) noexcept;

Returns: true if lhs.rank() == rhs.rank() is true and lhs.extents(r) equals rhs.extents(r) for every rank index r of rhs, otherwise false.

4.7 In subsection 22.7.X.6 [mdspan.extents.dextents]

4.8 In subsection 22.7.X.3.1 [mdspan.layout.left.ctors]

Preconditions: other.required_span_size() is a representable value of size_type ([basic.fundamental]).

Preconditions: other.required_span_size() is a representable value of size_type ([basic.fundamental]).

Preconditions:

4.9 In subsection 22.7.X.4.1 [mdspna.layout.right.ctors]

Preconditions: other.required_span_size() is a representable value of size_type ([basic.fundamental]).

Preconditions: other.required_span_size() is a representable value of size_type ([basic.fundamental]).

Preconditions:

4.10 In Subsection 22.7.X.5.2 [mdspan.layout.stride.cotrs]

4.11 In subsection 22.7.X.1 [mdspan.mdspan.overview]

with:

  static constexpr size_t static_extent(size_t r) { return Extents::static_extent(r); }

with

  constexpr size_t size() const;

with:

  template <class ElementType, class SizeType, size_t N>
  mdspan(ElementType*, const array<SizeType, N>&)
    -> mdspan<ElementType, dextents<size_t, N>>;

  template <class ElementType, class SizeT, size_t... ExtentsPack>
    mdspan(ElementType*, const extents<SizeT, ExtentsPack...>&)
    -> mdspan<ElementType, extents<SizeT, ExtentsPack...>>;

4.12 In subsection 22.7.X.1 [mdspan.mdspan.ctors]

4.13 In Subsection 22.7.X [mdspan.submdspan]

5 Implementation

There is an mdspan implementation available at https://github.com/kokkos/mdspan/.

6 Related Work