Notes for submdspan LEWG review 2022-10-10

Intro

• submdspan was an important part of the mdspan proposal P0009.
• But: deferred to C++26 due to limits in available LWG review cycles

This opens up an opportunity however: Let us address things we originally deferred for later

• submdspan customization points for user defined layouts
• non-contiguous subrange (e.g. every 3rd element in a dimension)
• static sub-extent with pair/tuple of compile time constants (e.g. tuple of integral_constant)

P2630: Submdspan addresses these concerns

Submdspan Example

// Some data representing 2x3 matrices
double data[6] = { 11, 12, 13,
                   21, 22, 23 };
// An rank-2 mdspan of the data
mdspan matrix(data, 2, 3);


mdspan<double,extents<size_t, dynamic_extent>, layout_stride>
   col_0 = submdspan(matrix, full_extent, 0);
// col_0 == 11, 21
// col_0.stride(0) == 3

mdspan<double,extents<size_t, dynamic_extent>, layout_right>
   row_1 = submdspan(matrix, 1, full_extent);
// row_1 == 21, 22, 23
// row_1.stride(0) == 1

// in practice this is a place where auto is your friend:
auto row_0 = submdspan(matrix, 0, full_extent);

P0009 submdspan

Slice Specifiers

There are a few natural descriptors for describing subranges:

• single integral: reduce the rank by one, sub-index-space consists of all multidimensional indices which match the value in the specified position
• tuple: all values in the begin/end range described in that tuple (half open range like iterator pairs)
• full_extent: all values in that dimension

For example:

// Some data representing 2x3 matrices
double data[6] = { 11, 12, 13, 14,
                   21, 22, 23, 24,
                   31, 32, 33, 34, };
// An rank-2 mdspan of the data
mdspan matrix(data, 3, 4);

// Submatrix excluding first and last column
auto subm = submdspan(matrix, full_extent, pair(1,3));
// subm.rank() ==2 
// subm == 12, 13,
//         22, 23,
//         32, 33

// Column zero
auto col_0 = submdspan(matrix, full_extent, 0);
// col_0 == 11, 21, 31
// col_0.stride(0) == 4

Restrictions in P0009 submdspan

• explicitly constrained on layout being one of layout_left, layout_right, or layout_stride

Return types P0009 submdspan

explicitly specified all the return types - no room for implementation defined behavior

• in all cases where layout_right or layout_left can be returned we mandate such (one could always return layout_stride since every submdspan of the above layouts is representable as such
• maintain static extents where possible (i.e. when using full_extent in a position where the extent is static)

P2630

Slice Specifiers

Maintain slice specifiers from before (note index_type in the following comes from the mdspan argument to submdspan):

• anything convertible to index_type
• anything convertible to tuple<index_type, index_type>
• full_extent_t
• new: strided_index_range

strided_index_range

template<class OffsetType, class ExtentType, class StrideType>
struct strided_index_range {
  using offset_type = OffsetType;
  using extent_type = ExtentType;
  using stride_type = StrideType;

  OffsetType offset;
  ExtentType extent;
  StrideType stride;
};

• enables taking a subrange with a stride in a given dimension
  • matlab: A(start:step:end)
  • numpy A[start:end:step]
  • Fortran A[start:end:step]
• We propose something slightly different:
  • strided_index_range slice{.offset = 2, .extent = 10, .stride = 3} i.e. we don't propose start and end, but start and extent
  • Reason: would enable static extent on the resulting mdspan, with runtime offset (e.g. get a statically sized (4,4) submatrix at runtime offsets)

Somewhat open question: is extent with respect to the original index space, or the resulting sub mdspan?

• Option 1 (current proposal): offset+extent is the end in the other programming models

double data[10] = {0,1,2,3,4,5,6,7,8,9};
mdspan mds(data, 10);
auto sub = submdspan(mds, strided_index_range{1,4,2});
// sub == {1,3}

• Option 2: extent IS the extent of the produced mdspan:

double data[10] = {0,1,2,3,4,5,6,7,8,9};
mdspan mds(data, 10);
auto sub = submdspan(mds, strided_index_range{1,4,2});
// sub == {1,3,5,7}

Benefit of Option 1: all parameters are expressed with respect to the source multidimensional index space
Benefit of Option 2: can get static extent result without static stride

Customization Points

To create the new mdspan we need:

• new data_handle
• new mapping
• new accessor

accessor is copy constructed from source accessor, its type is specified by source accessor
accessor_type::offset_policy
data_handle_type comes from accessor_type::offset_policy::data_handle_type and is computed via src.accessor().offset(src.data_handle(), SOME_OFFSET)
That leaves two things to compute: SOME_OFFSET and the new mapping

We propose "two" customization points for that

• submdspan_mapping
• submdspan_offset

Both take as argument src.mapping() and all the slice specifier arguments.

The customization points are invoked via argument dependent lookup by submdspan

Feedback from Tomasz

Make this a single customization point returning a pair.

I tend to agree: any thoughts from LEWG? If no objection we will make this change for the next revision.

Benefit: one needs to compute some of the same things for offset and the new mapping, could be more efficient to do it in one go.

submdspan_extents

We introduce a function to compute the new extents

• Makes sure we explicitly define preservation of static extents etc.
  • clearly defines new rank, all the extents and static extents.
  • used as a precondition for the customization points: i.e. you better return a mapping with an extents which makes sense ...
  • takes an extents argument and slice specifiers

submdspan

With the above functions we get a submdspan implementation:

template<class T, class E, class L, class A,
         class ... SliceArgs)
auto submdspan(const mdspan<T,E,L,A>& src, SliceArgs ... args) {
  size_t sub_offset = submdspan_offset(src.mapping(), args...);
  auto sub_map = submdspan_mapping(src.mapping(), args...);
  typename A::offset_policy sub_acc(src.accessor());
  typename A::offset_policy::data_handle_type 
    sub_handle = src.accessor().offset(src.data_handle(), sub_offset);
  return mdspan(sub_handle, sub_map, sub_acc);
}

• Reminder: submdspan_offset(src.mapping(), args...) and submdspan_mapping(src.mapping(), args...) are called using argument dependent lookup
• User defined layouts would simply provide overloads for these functions in the same namespace as the layout.

Preserving compile time information

We extended on the preservation of compile time information

• tuple<integral_constant<T0,val0>, integral_constant<T1,val1> results in static extent of val1-val0
• strided_index_range<OT,integral_constant<ET,ext>,integral_constant<ST,stride>> results in static extent of ext/stride

More feedback from Tomasz

• Some wording improvements (e.g. specify slice specifier somewhere once, instead of redefining legal arguments)
• Don't have submdspan_[extents/mapping/offset] take any slice specifier but just:
  • integral/integral_constant
  • strided_index_range
    • Note strided_index_range covers tuple and full_extent use case - just with more complex spelling
      • tuple<T0, T1> == strided_index_range<T0, T1, integral_constant<size_t,1>>
      • full_extent_t == strided_index_range<integral_constant<size_t,0>, int, integral_constant<size_t,1>> for dynamic extent
      • full_extent_t == strided_index_range<integral_constant<size_t,0>, integral_constant<size_t, Extent>, integral_constant<size_t,1>> for static extent
  • Benefit: makes it easier to implement sumdspan_[extents/mapping/offset]
  • Drawback: the functions are less convenient on their own, but not sure if that is important?

i.e. we would have:

template<class IndexT, size_t ... Extents, class ... SliceSpecifiers>
requires(sizeof...(Extents) == sizeof...(SliceSpecifiers) &&
        ((is-strided-index-range-v<SliceSpecifiers> ||
          is_integral_v<SliceSpecifiers> || 
          is-integral-constant-v<SliceSpecifiers>)...))
auto submdspan_extents(const extents<IndexT, Extents...>& ext, 
                       SliceSpecifiers ... slices);

Polls

1. submdspan should have a single customization point submdspan_mapping which returns a pair<sub_mapping_type, size_t>
2. submdspan customization points should only take integral, integral_constant and strided_index_range slice specifiers.
3. LEWG prefers strided_index_range::extent to be the extent of the created sub mdspan.

No implementation files or tests are identified. Start by reading the P0009 and P2630 discussion in this issue, then determine whether the three open poll questions have received a decision; done would be an agreed direction for a subsequent proposal revision.