Today’s post is dedicated to a new feature introduced in C++23 called std::mdspan. There are already a few articles online about it, but it caught my personal attention for three reasons:

1. I’ve wanted to get rid of boost::multi_array in my codebase for a long time, as part of the larger effort to reduce dependencies on Boost. Why reduce them? As strange as it may sound, I’m aiming to migrate the codebase toward C++ Modules as soon as possible. Boost is not designed for modularization, and integrating Boost into a modularized project is extremely painful today. boost::multi_array has been one of the most problematic dependencies in this regard.
2. I really wanted to try out the new mdspan feature and the new bracket subscription operator in C++23.
3. I want to prepare the data structures for C++26 std::linalg.

What is std::mdspan?

As you probably know, C++20 introduced std::span, which provides a view over a contiguous array. For example, if you have a std::vector, a raw array, or anything that is laid out contiguously in memory, you can create a span over it and perform operations without copying the data. Suppose you want to sum the elements of an array:

int sum(std::span<const int> input) {
	int output = 0;
	for (auto value : input)
		output += value;
	return output;
}

You can create a span over a part of the array:

std::vector<int> input = { 1, 2, 3, 4, 5 };
auto span = std::span(input.data() + 2, 3);
auto result = sum(span);
std::print("sum of span (3,4,5): {}", result);

Or you can pass a std::array:

std::array array{ 1, 2, 3, 5 };
auto result = sum(array);
std::print("sum of span (1,2,3,5): {}", result);

Notice that adding const to std::span makes the data read-only.
If you want to modify the elements, you must use std::span<int> instead of std::span<const int>.

In essence, std::span is a safe and simple “reference” to an array or a part of it – it doesn’t own the data, it just provides access.

But what about multidimensional arrays?

In many mathematical applications (based on my experience), there is a strong need for something similar, but for multidimensional arrays. For example, my most recent tasks involved working with 3D voxel maps and 3D control grids for free-form deformation. In such cases, it is often necessary to modify a part of the control grid or one of its sides. For this, it is extremely convenient to either extract a slice or create a view over a subregion of the 3D array.

Despite the prevalence of matrices in computational tasks, C++ itself still doesn’t have a built-in matrix type. While mdarray is expected to be added in the future, std::mdspan already provides a clean way to view a multidimensional block of memory as a matrix.

A simple example

Imagine you have a flat array:

std::vector<int> data = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 };

You can interpret it as a 3D array (4 × 3 × 2) using std::mdspan:

std::mdspan<int, std::extents<std::size_t, 4, 3, 2>> mdspan(data.data(), 4, 3, 2);

Now, you can access and modify elements just like in a real 3D array:

mdspan[0, 1, 0] = 42;

Notice how beautiful and natural the indexing looks – no more [][][] proxy chains thanks to the C++23 subscription operator.

Memory layout

In memory, the elements are laid out contiguously, first by Z, then Y, then X (C-style, layout_right). For our 4 × 3 × 2 array:

mdspan[0,0,0] → data[0]
mdspan[0,0,1] → data[1]
mdspan[0,1,0] → data[2]
mdspan[0,1,1] → data[3]
...
mdspan[3,2,1] → data[23]

Linear index formula (for layout_right):

linear_index = i * (extent(1) * extent(2)) + j * extent(2) + k

Useful mdspan properties

std::println("mdspan size: {}", mdspan.size());
std::println("mdspan extents: {} {} {}", mdspan.extent(0), mdspan.extent(1), mdspan.extent(2));
std::println("mdspan rank: {}", mdspan.rank());
std::println("mdspan empty: {}", mdspan.empty());
std::println("mdspan stride 0: {}", mdspan.stride(0));
std::println("mdspan stride 1: {}", mdspan.stride(1));
std::println("mdspan stride 2: {}", mdspan.stride(2));

for (size_t i = 0; i < mdspan.extent(0); ++i) {
	for (size_t j = 0; j < mdspan.extent(1); ++j) {
		for (size_t k = 0; k < mdspan.extent(2); ++k) 
			std::print("[{},{},{}] = {}; ", i,j,k, mdspan[i, j, k]);			
		std::println();
	}
	std::println();
}

I believe the code and properties largely speak for themselves here, but the stride properties are particularly interesting. stride(i) represents the number of memory steps you need to take along axis i to reach the next element in that direction. In other words:

mdspan[i, j, k] == *(data.data() + i * stride(0) + j * stride(1) + k * stride(2))

Where:

Essentially, the main properties have already been covered above, and by now it should be fairly intuitive how to work with std::mdspan. Now let’s take a closer look at the concept of a slice, which is often used in various mathematical tasks. In this context, a slice refers to a 2D array along one of the axes at a fixed level.

For example, along the X-axis, there are 4 slices, as shown earlier.
But how do we write a function that returns the desired slice?

To achieve this, we can utilize std::layout_stride with std::mdspan. The mdspan declaration would then look like this:

using submdspan_t = std::mdspan<
    T,
    std::extents<size_t, std::dynamic_extent, std::dynamic_extent>,
    std::layout_stride
>;

It’s similar to a regular mdspan, but now we explicitly specify std::layout_stride as the third template parameter. Next, we want to extract the i-th layer along the X-axis. First, we specify the shape of such a slice:

std::extents shape{dim1, dim2};

Then, we define how the stride should be calculated:

std::layout_stride::mapping mapping{shape, std::array{dim2, 1uz}};

Note: uz (introduced in C++23) is used as a suffix for size_t literals.

This means:

• Moving along the Y-axis requires a jump of dim2 elements (i.e., one row of Z elements).
• Moving along the Z-axis requires a jump of 1 element.

The slice is then constructed as:

subspan_t slice{ptr, mapping};

where ptr is a pointer to span[sliceIndex, 0, 0].

For the X-axis, slice number 1 corresponds to the following elements (shown both in 2D slice coordinates, 3D mdspan coordinates, and linear index):

slice[0,0] → mdspan[1,0,0] → data[6]
slice[0,1] → mdspan[1,0,1] → data[7]
slice[1,0] → mdspan[1,1,0] → data[8]
slice[1,1] → mdspan[1,1,1] → data[9]
slice[2,0] → mdspan[1,2,0] → data[10]
slice[2,1] → mdspan[1,2,1] → data[11]

It would be even more convenient if C++ eventually adopted a slicing mechanism similar to what NumPy offers, where slices and subarrays can be selected using concise and expressive syntax. For example, in NumPy, extracting a slice along the X-axis at index 1 would look like this:

slice = data[1, :, :]

In fact, there is a proposal – submdspan. However, it has not yet been officially included in the standard.

The complete function looks like this (I intentionally made the function a bit verbose to make the stride logic more visible):

template <typename T, typename Extents>
auto extract_slice(std::mdspan<T, Extents> mdspan_view, size_t slice_index, size_t dimension_index)
-> std::mdspan<T, std::dextents<size_t, 2>, std::layout_stride>
{
	static_assert(Extents::rank() == 3, "extract_slice supports only 3D mdspans");

	if (slice_index >= mdspan_view.extent(dimension_index))
		throw std::out_of_range("slice_index out of bounds");

	size_t extent0 = mdspan_view.extent(0);
	size_t extent1 = mdspan_view.extent(1);
	size_t extent2 = mdspan_view.extent(2);

	switch (dimension_index) {
	case 0: {
		std::extents shape{ extent1, extent2 };
		std::array strides{ extent2, 1uz };
		return { &mdspan_view[slice_index, 0, 0], {shape, strides} };
	}
	case 1: {
		std::extents shape{ extent0, extent2 };
		std::array strides{ extent1 * extent2, 1uz };
		return { &mdspan_view[0, slice_index, 0], {shape, strides} };
	}
	case 2: {
		std::extents shape{ extent0, extent1 };
		std::array strides{ extent1 * extent2, extent2 };
		return { &mdspan_view[0, 0, slice_index], {shape, strides} };
	}
	default:
		throw std::invalid_argument("Invalid dimension_index");
	}
}

When using extract_slice, it might seem convenient to copy data directly via std::copy(slice.data_handle(), slice.data_handle() + slice.size(), dst.begin()), and in simple cases this works because the slice memory is dense and contiguous. However, since extract_slice creates slices using std::layout_stride, there is no general guarantee that the elements are physically continuous in memory.

Replacing boost::multi_array with std::mdspan

Simply replacing boost::multi_array with std::mdspan is not enough, because mdspan is just a view. You need to manage the actual storage separately, typically with std::vector.

Here, I will describe the operations that I had to rewrite when replacing boost::multi_array with std::mdspan and std::vector.
The main operations include retrieving the size of the 3D array, reading and writing an element by a 3D index, extracting the i-th slice along a specified axis, and resizing the 3D array.

Resizing a boost::multi_array is fairly straightforward:

data.resize(boost::extents[newShape[0]][newShape[1]][newShape[2]]);

Whereas with the std::vector + std::mdspan solution, you need to manually manage the size of the vector and rebind the std::mdspan to the resized storage after each resize:

storage.resize(newShape[0] * newShape[1] * newShape[2]);
mdspan = mdspan(storage.data(), newShape);

And finally, here’s the equivalent of extract_slice for boost::multi_array (analogous to the version presented above for mdspan):

auto extract_slice(size_t slice_index, size_t dimension_index) {
	typename boost::multi_array_types::index_gen::gen_type<2, 3>::type slice;

	if (dimension_index == 0) {
		slice.ranges_[0] = boost::multi_array_types::index_range(slice_index);
		slice.ranges_[1] = boost::multi_array_types::index_range();
		slice.ranges_[2] = boost::multi_array_types::index_range();
	}
	else if (dimension_index == 1) {
		slice.ranges_[0] = boost::multi_array_types::index_range();
		slice.ranges_[1] = boost::multi_array_types::index_range(slice_index);
		slice.ranges_[2] = boost::multi_array_types::index_range();
	}
	else if (dimension_index == 2) {
		slice.ranges_[0] = boost::multi_array_types::index_range();
		slice.ranges_[1] = boost::multi_array_types::index_range();
		slice.ranges_[2] = boost::multi_array_types::index_range(slice_index);
	}
	else
		throw std::invalid_argument("Invalid dimension_index");	

	return slice;
}

And its usage:

typename boost::multi_array<T, 3>::template array_view<2>::type slice = data[extract_slice(index, dir)];

Migration turned out to be relatively simple and, in many places, significantly cleaned up the legacy code. It would also be great if submdspan gets added, as it would simplify the code even further and eliminate the need for the monstrous extract_slice method.

Bonus: eliminating boost::multi_array helps avoid internal linkage conflicts that are problematic for C++ Modules.