C++23 多维数组切片（std::mdspan）：在 C++ 高性能计算中实现对连续内存的高维度逻辑映射与访问

C++23 多维数组切片（std::mdspan）：在 C++ 高性能计算中实现对连续内存的高维度逻辑映射与访问

各位编程爱好者、C++ 开发者以及高性能计算领域的同仁们，大家好！

今天，我们将深入探讨 C++23 标准库中一个革命性的新特性：std::mdspan。在高性能计算（HPC）领域，处理多维数据是家常便饭，无论是物理模拟、机器学习、图像处理还是科学计算，我们都离不开对矩阵、张量等高维数据的操作。然而，传统 C++ 在处理这些复杂的数据结构时，往往伴随着诸多挑战：内存布局的困扰、切片操作的繁琐、边界检查的缺失以及与底层硬件（特别是缓存）的低效交互。std::mdspan 的出现，正是为了解决这些痛点，它提供了一种类型安全、零开销、极度灵活的多维数据视图，使得 C++ 在 HPC 领域的竞争力得到了显著提升。

本次讲座的目标是帮助大家全面理解 std::mdspan 的设计哲学、核心功能、各种布局策略，以及如何利用它在连续内存上高效地实现高维度逻辑映射与切片操作。我们将通过丰富的代码示例，从基础用法到高级应用，逐步揭示 std::mdspan 的强大威力，并探讨它在实际 HPC 项目中的应用潜力和性能优势。

一、多维数据在 HPC 中的挑战

在深入 std::mdspan 之前，我们首先回顾一下在高性能计算中处理多维数据时，开发者们经常面临的挑战：

1. 内存布局与缓存性能：

   • 行主序 (Row-major) 与列主序 (Column-major)： C/C++ 默认采用行主序（最后维度变化最快），而 Fortran 和许多线性代数库（如 BLAS/LAPACK）则采用列主序（第一维度变化最快）。这种不一致性在跨语言或跨库交互时容易导致混淆和错误。
   • 缓存局部性： 处理器缓存是现代计算机性能的关键。如果数据访问模式不符合内存的物理布局，会导致大量的缓存缺失 (cache miss)，严重拖慢程序执行速度。例如，在行主序存储的矩阵中按列访问元素，会导致不连续的内存跳跃，从而降低缓存效率。

2. 多维切片与子区域访问：

   • 复杂性： 从一个大的多维数组中提取一个子区域（例如，一个子矩阵、一个特定行或列）通常需要复杂的指针算术和维度计算。
   • 安全性： 使用裸指针进行切片操作时，缺乏编译期或运行时的边界检查，极易引入越界访问错误，导致程序崩溃或数据损坏。
   • 灵活性： 传统的 C++ 数组或 std::vector<std::vector<...>> 等结构，难以实现任意维度、任意步长的灵活切片，例如跳跃式切片或降维切片。

3. 泛型编程与维度无关性：

   • 编写能够处理任意维度（2D、3D、ND）多维数据的泛型算法非常困难。通常需要为不同维度编写重载函数或使用复杂的模板元编程技巧。
   • 算法接口需要能够清晰地表达它所操作的数据结构是多维的，并且能够安全地获取其维度信息。

4. 异构计算与互操作性：

   • 在 GPU 等加速器上进行计算时，需要将 CPU 内存中的数据映射到设备内存，并以设备友好的方式访问。mdspan 作为一种视图，其设计使其更容易与 CUDA、OpenCL 等异构计算框架集成。
   • 与现有的高性能库（如 BLAS、LAPACK、TensorFlow C++ API）进行数据交互时，需要一种标准且高效的方式来传递多维数据视图。

5. 零开销抽象：

   • HPC 追求极致的性能，任何额外的开销（如不必要的内存拷贝、虚函数调用、运行时边界检查）都可能成为瓶颈。理想的多维数组抽象应该在提供类型安全和便利性的同时，不引入额外的运行时负担，其性能应与手工优化过的裸指针操作相当。

这些挑战促使 C++ 社区引入 std::mdspan，旨在提供一个标准、高效且类型安全的多维数据视图解决方案。

二、std::mdspan 核心概念与设计哲学

std::mdspan 的核心思想是提供一个 非拥有 (non-owning) 的 多维视图 (multidimensional view)。这意味着 mdspan 不管理底层数据内存的生命周期，它仅仅是提供了一种逻辑上的多维访问方式，指向一块已经存在的连续内存区域。

2.1 核心概念

1. 非拥有视图 (Non-owning View):

   • mdspan 就像一个指针，它指向一块内存，并知道如何将多维索引映射到这块内存中的一维地址。
   • 它的生命周期可以短于底层数据源的生命周期，但必须确保在其活跃期间，底层数据源是有效的。
   • 这与 std::span 的理念相似，std::span 是一维数据的非拥有视图。

2. 维度 (Dimensions) 与范围 (Extents):

   • 一个 mdspan 实例通过其模板参数 Extents 来描述其维度和每个维度的大小。
   • std::extents<IndexType, Extent0, Extent1, ...> 定义了多维数组的形状。Extent_N 可以是编译时常量（静态维度）或 std::dynamic_extent（动态维度）。
   • rank() 方法返回维度数。extent(d) 返回第 d 维的大小。

3. 布局策略 (Layout Policy):

   • 布局策略定义了如何将多维索引 (i, j, k, …) 映射到一维的线性内存地址。这是 mdspan 区别于传统指针算术的关键，也是其灵活性的核心。
   • 标准库提供了三种布局策略：
     • std::layout_right (C-style, 行主序): 默认布局，最后维度变化最快。
     • std::layout_left (Fortran-style, 列主序): 第一维度变化最快。
     • std::layout_stride (Custom, 自定义步长): 最灵活，允许任意步长，是实现高级切片的基础。

4. 访问器 (Accessor Policy):

   • 访问器定义了如何从底层数据指针获取或设置元素。它是一个轻量级的策略对象，封装了对原始指针的解引用和索引操作。
   • std::default_accessor<ElementType> 是默认访问器，它只是简单地对指针进行解引用。
   • 用户可以自定义访问器，例如，实现对原子类型的线程安全访问，或者将数据存储在 GPU 内存中，提供跨设备访问的语义。

2.2 设计哲学

std::mdspan 的设计严格遵循 C++ 的“零开销抽象”原则，其目标是提供与裸指针操作相同甚至更好的性能，同时大大提升代码的安全性、可读性和灵活性。

• 编译期检查与类型安全： 尽可能多的维度信息在编译期确定，允许编译器进行更严格的类型检查和优化。例如，如果尝试访问超出静态维度的索引，编译器会报错。对于动态维度，运行时会进行边界检查（如果开启了调试模式或特定配置）。
• 运行时零开销： mdspan 对象本身非常小，通常只包含一个数据指针和布局映射信息。元素的访问操作通过内联的布局策略方法完成，计算出内存地址，然后通过访问器解引用，不涉及虚函数调用或其他运行时开销。
• 泛型与互操作性： mdspan 的设计使其能够轻松地与泛型算法、模板元编程以及其他库（如 CUDA、BLAS 等）结合，提供统一的多维数据接口。
• 切片支持： 引入 std::submdspan 函数，以安全、高效且符合人体工程学的方式生成原始 mdspan 的子视图，无需内存拷贝。

三、std::mdspan 的基本用法

std::mdspan 的模板参数列表如下：

template<
    class ElementType,
    class Extents,
    class LayoutPolicy = std::layout_right,
    class AccessorPolicy = std::default_accessor<ElementType>
>
class mdspan;

• ElementType: 存储在多维数组中的元素类型（例如 int, double）。
• Extents: 描述多维数组的维度和每个维度的大小。可以使用 std::extents<IndexType, Extent0, Extent1, ...>。Extent_N 可以是 std::dynamic_extent 表示动态维度。
• LayoutPolicy: 布局策略，默认为 std::layout_right。
• AccessorPolicy: 访问器策略，默认为 std::default_accessor<ElementType>。

3.1 构造 mdspan

mdspan 可以从原始指针、std::array、std::vector 等连续内存容器以及其他 mdspan 构造。

示例 1：静态二维 mdspan (行主序)

#include <iostream>
#include <vector>
#include <numeric> // For std::iota
#include <mdspan>  // C++23 header

// 为了方便，使用命名空间别名
namespace stdex = std::experimental; // 在某些编译器上可能仍然是 std::experimental

int main() {
    // 1. 准备底层数据：一个一维的 std::vector 模拟二维数据
    // 3行 x 4列 的矩阵
    std::vector<int> data(3 * 4);
    std::iota(data.begin(), data.end(), 0); // 填充 0, 1, 2, ..., 11

    // 2. 定义 Extents：3行，4列
    // std::extents<std::size_t, 3, 4> 表示一个编译期确定的 3x4 矩阵
    using MyExtents = std::extents<std::size_t, 3, 4>;

    // 3. 构造 mdspan
    // ElementType: int
    // Extents: MyExtents (3x4)
    // LayoutPolicy: std::layout_right (默认，行主序)
    std::mdspan<int, MyExtents> matrix(data.data());

    // 4. 访问元素并打印
    std::cout << "Static 2D mdspan (Row-major):" << std::endl;
    for (std::size_t i = 0; i < matrix.extent(0); ++i) { // 遍历行
        for (std::size_t j = 0; j < matrix.extent(1); ++j) { // 遍历列
            std::cout << matrix(i, j) << "t"; // 使用 operator() 访问元素
        }
        std::cout << std::endl;
    }
    /* 输出:
    0       1       2       3
    4       5       6       7
    8       9       10      11
    */

    // 5. 修改元素
    matrix(0, 0) = 99;
    std::cout << "Modified matrix(0,0): " << matrix(0, 0) << std::endl;
    std::cout << "Underlying data[0]: " << data[0] << std::endl; // 验证底层数据也被修改

    return 0;
}

示例 2：动态二维 mdspan

#include <iostream>
#include <vector>
#include <numeric>
#include <mdspan>

// namespace stdex = std::experimental;

int main() {
    int rows = 3;
    int cols = 4;
    std::vector<double> data(rows * cols);
    std::iota(data.begin(), data.end(), 0.0);

    // std::extents<std::size_t, std::dynamic_extent, std::dynamic_extent>
    // 表示一个运行时确定大小的二维矩阵
    using DynamicExtents = std::extents<std::size_t, std::dynamic_extent, std::dynamic_extent>;

    // 构造 mdspan 时需要提供每个维度的大小
    std::mdspan<double, DynamicExtents> dynamic_matrix(data.data(), rows, cols);

    std::cout << "nDynamic 2D mdspan (Row-major):" << std::endl;
    for (std::size_t i = 0; i < dynamic_matrix.extent(0); ++i) {
        for (std::size_t j = 0; j < dynamic_matrix.extent(1); ++j) {
            std::cout << dynamic_matrix(i, j) << "t";
        }
        std::cout << std::endl;
    }
    /* 输出:
    0       1       2       3
    4       5       6       7
    8       9       10      11
    */

    // 动态维度信息
    std::cout << "Rank: " << dynamic_matrix.rank() << std::endl; // 2
    std::cout << "Extent of dim 0: " << dynamic_matrix.extent(0) << std::endl; // 3
    std::cout << "Extent of dim 1: " << dynamic_matrix.extent(1) << std::endl; // 4
    std::cout << "Total elements: " << dynamic_matrix.size() << std::endl; // 12

    return 0;
}

示例 3：三维 mdspan (动态维度)

#include <iostream>
#include <vector>
#include <numeric>
#include <mdspan>

// namespace stdex = std::experimental;

int main() {
    int dim0 = 2;
    int dim1 = 3;
    int dim2 = 4;
    std::vector<float> data(dim0 * dim1 * dim2);
    std::iota(data.begin(), data.end(), 0.0f);

    using Dynamic3DExtents = std::extents<std::size_t, std::dynamic_extent, std::dynamic_extent, std::dynamic_extent>;

    // 构造 mdspan 时提供所有动态维度的大小
    std::mdspan<float, Dynamic3DExtents> volume(data.data(), dim0, dim1, dim2);

    std::cout << "nDynamic 3D mdspan (Row-major):" << std::endl;
    for (std::size_t i = 0; i < volume.extent(0); ++i) {
        for (std::size_t j = 0; j < volume.extent(1); ++j) {
            for (std::size_t k = 0; k < volume.extent(2); ++k) {
                std::cout << "volume(" << i << "," << j << "," << k << ") = " << volume(i, j, k) << std::endl;
            }
        }
    }
    /* 部分输出:
    volume(0,0,0) = 0
    volume(0,0,1) = 1
    ...
    volume(1,2,3) = 23
    */

    return 0;
}

3.2 访问元素

mdspan 提供了 operator() 用于访问元素，传入与维度数相等的索引参数。这比 operator[] 链式访问（如 matrix[i][j]）更清晰和高效，因为它避免了多层解引用。

// 访问二维 mdspan 的元素
int val = matrix(row_index, col_index);

// 访问三维 mdspan 的元素
float val3d = volume(idx0, idx1, idx2);

四、深入理解布局策略

布局策略是 mdspan 最强大的特性之一，它决定了多维索引如何映射到一维内存地址。选择合适的布局策略对于缓存性能至关重要。

4.1 std::layout_right (行主序)

• 特点： 这是 C/C++ 多维数组的默认布局。在内存中，最后维度（最右边的索引）变化最快，相邻元素在内存中也是相邻的。
• 内存地址计算： 对于一个 N 维 mdspan，索引为 (i_0, i_1, ..., i_{N-1})，其在内存中的偏移量计算公式为：
  offset = i_0 * S_1 * S_2 * ... * S_{N-1} + i_1 * S_2 * ... * S_{N-1} + ... + i_{N-2} * S_{N-1} + i_{N-1}
  其中 S_k 是第 k 维的大小。
• 适用场景： 当算法主要按行（或最后维度）顺序访问数据时，例如遍历矩阵的每一行，或者进行图像处理中的像素扫描。

示例：

#include <iostream>
#include <vector>
#include <numeric>
#include <mdspan>

// namespace stdex = std::experimental;

int main() {
    std::vector<int> data(2 * 3); // 2行3列
    std::iota(data.begin(), data.end(), 0); // 0, 1, 2, 3, 4, 5

    // 默认就是 layout_right
    std::mdspan<int, std::extents<std::size_t, 2, 3>, std::layout_right> matrix_right(data.data());

    std::cout << "mdspan with std::layout_right (Row-major):" << std::endl;
    // 打印内存中的实际值
    for (std::size_t i = 0; i < data.size(); ++i) {
        std::cout << data[i] << " ";
    }
    std::cout << std::endl; // Output: 0 1 2 3 4 5

    // 逻辑访问
    std::cout << "Logical access (i,j) -> value:" << std::endl;
    for (std::size_t i = 0; i < matrix_right.extent(0); ++i) {
        for (std::size_t j = 0; j < matrix_right.extent(1); ++j) {
            std::cout << "(" << i << "," << j << ")=" << matrix_right(i, j) << "t";
        }
        std::cout << std::endl;
    }
    /* Output:
    (0,0)=0       (0,1)=1       (0,2)=2
    (1,0)=3       (1,1)=4       (1,2)=5
    */

    // 验证物理地址和逻辑访问的对应关系
    // matrix_right(0,0) -> data[0]
    // matrix_right(0,1) -> data[1]
    // matrix_right(0,2) -> data[2]
    // matrix_right(1,0) -> data[3]
    // matrix_right(1,1) -> data[4]
    // matrix_right(1,2) -> data[5]

    return 0;
}

4.2 std::layout_left (列主序)

• 特点： Fortran 语言和许多线性代数库（如 BLAS、LAPACK）采用的布局。在内存中，第一维度（最左边的索引）变化最快，相邻元素在内存中也是相邻的。
• 内存地址计算： 对于一个 N 维 mdspan，索引为 (i_0, i_1, ..., i_{N-1})，其在内存中的偏移量计算公式为：
  offset = i_0 + i_1 * S_0 + i_2 * S_0 * S_1 + ... + i_{N-1} * S_0 * S_1 * ... * S_{N-2}
  其中 S_k 是第 k 维的大小。
• 适用场景： 当算法主要按列（或第一维度）顺序访问数据时，例如进行矩阵乘法时，如果内层循环遍历列，使用列主序可以获得更好的缓存局部性。与 BLAS/LAPACK 库进行数据交互时，通常需要列主序。

示例：

#include <iostream>
#include <vector>
#include <numeric>
#include <mdspan>

// namespace stdex = std::experimental;

int main() {
    std::vector<int> data(2 * 3); // 2行3列
    // 注意：这里的data仍然是按行主序填充的，但mdspan会将其解释为列主序
    // 所以逻辑上的(0,0) (1,0) (0,1) (1,1) ... 会对应到物理内存的 0 1 2 3 ...
    std::iota(data.begin(), data.end(), 0); // 0, 1, 2, 3, 4, 5

    // 使用 layout_left 策略
    std::mdspan<int, std::extents<std::size_t, 2, 3>, std::layout_left> matrix_left(data.data());

    std::cout << "nmdspan with std::layout_left (Column-major):" << std::endl;
    // 打印内存中的实际值
    for (std::size_t i = 0; i < data.size(); ++i) {
        std::cout << data[i] << " ";
    }
    std::cout << std::endl; // Output: 0 1 2 3 4 5

    // 逻辑访问
    std::cout << "Logical access (i,j) -> value:" << std::endl;
    for (std::size_t i = 0; i < matrix_left.extent(0); ++i) {
        for (std::size_t j = 0; j < matrix_left.extent(1); ++j) {
            std::cout << "(" << i << "," << j << ")=" << matrix_left(i, j) << "t";
        }
        std::cout << std::endl;
    }
    /* Output:
    (0,0)=0       (0,1)=2       (0,2)=4
    (1,0)=1       (1,1)=3       (1,2)=5
    */

    // 验证物理地址和逻辑访问的对应关系
    // matrix_left(0,0) -> data[0]
    // matrix_left(1,0) -> data[1]
    // matrix_left(0,1) -> data[2]
    // matrix_left(1,1) -> data[3]
    // matrix_left(0,2) -> data[4]
    // matrix_left(1,2) -> data[5]

    return 0;
}

4.3 std::layout_stride (自定义步长)

• 特点： 这是最灵活的布局策略，它允许用户为每个维度指定一个步长 (stride)。步长表示当该维度的索引增加 1 时，底层数据指针需要跳过的元素数量。
• 内存地址计算： 对于一个 N 维 mdspan，索引为 (i_0, i_1, ..., i_{N-1})，如果步长数组为 (s_0, s_1, ..., s_{N-1})，其在内存中的偏移量计算公式为：
  offset = i_0 * s_0 + i_1 * s_1 + ... + i_{N-1} * s_{N-1}
• 适用场景：
  • 任意切片： 它是实现 submdspan 的基础，可以灵活地创建子视图、跳跃视图、降维视图等。
  • 不规则布局： 当数据在内存中不是严格的行主序或列主序时，例如稀疏矩阵的部分视图。
  • 转置视图： 通过交换步长和维度大小，可以轻松实现矩阵的转置视图，而无需复制数据。

构造 layout_stride 的 mdspan 需要一个 layout_stride::mapping 对象，该对象在构造时接收 extents 和 strides。

示例：自定义步长实现转置视图

#include <iostream>
#include <vector>
#include <numeric>
#include <mdspan>

// namespace stdex = std::experimental;

int main() {
    std::vector<int> data(3 * 4); // 3行4列
    std::iota(data.begin(), data.end(), 0); // 0, 1, ..., 11

    // 原始矩阵 (行主序)
    std::mdspan<int, std::extents<std::size_t, 3, 4>> matrix(data.data());

    std::cout << "Original Matrix (3x4, Row-major):" << std::endl;
    for (std::size_t i = 0; i < matrix.extent(0); ++i) {
        for (std::size_t j = 0; j < matrix.extent(1); ++j) {
            std::cout << matrix(i, j) << "t";
        }
        std::cout << std::endl;
    }
    /* Output:
    0       1       2       3
    4       5       6       7
    8       9       10      11
    */

    // 构造转置视图 (4x3)
    // 新的 Extents 是 4x3
    using TransposedExtents = std::extents<std::size_t, 4, 3>;

    // 原始矩阵的行步长是 4 (跳过4个元素到下一行)
    // 原始矩阵的列步长是 1 (跳过1个元素到下一列)
    // 对于转置矩阵 (4x3)，它的 "行" 对应原矩阵的列，"列" 对应原矩阵的行。
    // 所以转置矩阵的第一个维度 (新行) 步长应该与原矩阵的列步长相同 (1)。
    // 转置矩阵的第二个维度 (新列) 步长应该与原矩阵的行步长相同 (4)。
    std::array<std::size_t, 2> strides = {1, 4}; // 注意顺序，对应新维度0和维度1

    // 构造 layout_stride::mapping
    std::layout_stride::mapping<TransposedExtents> transposed_mapping(
        TransposedExtents(4, 3), // 新的维度大小
        strides
    );

    // 构造 mdspan
    std::mdspan<int, TransposedExtents, std::layout_stride> transposed_matrix(
        data.data(), transposed_mapping
    );

    std::cout << "nTransposed Matrix View (4x3, using layout_stride):" << std::endl;
    for (std::size_t i = 0; i < transposed_matrix.extent(0); ++i) {
        for (std::size_t j = 0; j < transposed_matrix.extent(1); ++j) {
            std::cout << transposed_matrix(i, j) << "t";
        }
        std::cout << std::endl;
    }
    /* Output:
    0       4       8
    1       5       9
    2       6       10
    3       7       11
    */

    return 0;
}

通过这个例子，我们可以看到 layout_stride 的强大之处：它允许我们在不复制数据的情况下，以完全不同的逻辑布局来解释同一块内存。

五、std::mdspan 的多维切片

多维切片是 std::mdspan 的核心应用场景之一。std::submdspan 函数模板提供了一种安全、高效且直观的方式来创建现有 mdspan 的子视图。submdspan 返回的仍然是一个 mdspan，它共享原始数据，但具有新的 Extents 和 LayoutPolicy (通常是 layout_stride)。

std::submdspan 的函数签名类似于：

template<class ElementType, class Extents, class LayoutPolicy, class AccessorPolicy, class... SliceSpecs>
constexpr auto submdspan(
    mdspan<ElementType, Extents, LayoutPolicy, AccessorPolicy> m,
    SliceSpecs... slices
);

SliceSpecs 可以是以下几种类型：

1. std::integral_constant<IndexType, Value>： 选择一个固定索引，将该维度降维。
   • 例如：std::integral_constant<std::size_t, 0>() 选择第0个索引。
2. std::full_extent_t： 选择整个维度，不改变该维度。
   • 例如：std::full_extent 或 std::full_extent_t{}。
3. std::pair<IndexType, IndexType>： 选择一个范围 [start, end)。
   • 例如：std::pair{1, 3} 选择从索引 1 到 2 的元素。
4. std::tuple<...>： 允许在 layout_stride 的基础上，对切片维度进一步指定步长。这在 submdspan 的 SliceSpecs 参数中较少直接使用，更多是 layout_stride 本身的能力。

重要提示： submdspan 的返回值类型会根据 SliceSpecs 自动推导。如果切片导致所有维度都是静态的，则返回的 mdspan 也会是静态的；如果包含动态切片，则返回动态 mdspan。如果切片中包含 std::full_extent 或范围，则新的 mdspan 的 LayoutPolicy 会变成 std::layout_stride，因为其步长可能不再是标准行主序或列主序。

6.1 切片示例

我们使用一个 4×5 的矩阵作为基础数据。

#include <iostream>
#include <vector>
#include <numeric>
#include <mdspan> // C++23 header
#include <array>

// namespace stdex = std::experimental;

// 辅助函数：打印 mdspan
template <class MDSpan>
void print_mdspan(const std::string& title, const MDSpan& m) {
    std::cout << title << " (extents: ";
    for (std::size_t d = 0; d < m.rank(); ++d) {
        std::cout << m.extent(d) << (d == m.rank() - 1 ? "" : "x");
    }
    std::cout << "):" << std::endl;

    if (m.rank() == 2) {
        for (std::size_t i = 0; i < m.extent(0); ++i) {
            for (std::size_t j = 0; j < m.extent(1); ++j) {
                std::cout << m(i, j) << "t";
            }
            std::cout << std::endl;
        }
    } else if (m.rank() == 1) {
        for (std::size_t i = 0; i < m.extent(0); ++i) {
            std::cout << m(i) << "t";
        }
        std::cout << std::endl;
    } else {
        std::cout << "Cannot print mdspan with rank " << m.rank() << std::endl;
    }
    std::cout << std::endl;
}

int main() {
    // 原始数据：4行5列的矩阵
    std::vector<int> data(4 * 5);
    std::iota(data.begin(), data.end(), 10); // 10, 11, ..., 29

    using MatrixExtents = std::extents<std::size_t, 4, 5>;
    std::mdspan<int, MatrixExtents> original_matrix(data.data());
    print_mdspan("Original Matrix", original_matrix);
    /* Output:
    Original Matrix (extents: 4x5):
    10      11      12      13      14
    15      16      17      18      19
    20      21      22      23      24
    25      26      27      28      29
    */

    // 1. 选择特定行 (降维切片)
    // 选择第 1 行 (索引为 1 的行)
    auto row_1 = std::submdspan(original_matrix, std::integral_constant<std::size_t, 1>(), std::full_extent);
    print_mdspan("Slice: Row 1", row_1);
    /* Output:
    Slice: Row 1 (extents: 5):
    15      16      17      18      19
    */

    // 2. 选择特定列 (降维切片)
    // 选择第 2 列 (索引为 2 的列)
    auto col_2 = std::submdspan(original_matrix, std::full_extent, std::integral_constant<std::size_t, 2>());
    print_mdspan("Slice: Column 2", col_2);
    /* Output:
    Slice: Column 2 (extents: 4):
    12      17      22      27
    */

    // 3. 获取子矩阵 (范围切片)
    // 从 (1, 1) 开始，取 2行3列 的子矩阵
    // 行范围 [1, 3) (索引 1, 2)
    // 列范围 [1, 4) (索引 1, 2, 3)
    auto sub_matrix = std::submdspan(original_matrix, std::pair{1, 3}, std::pair{1, 4});
    print_mdspan("Slice: Sub-matrix (rows 1-2, cols 1-3)", sub_matrix);
    /* Output:
    Slice: Sub-matrix (rows 1-2, cols 1-3) (extents: 2x3):
    16      17      18
    21      22      23
    */

    // 4. 跳跃切片 (通过 layout_stride 实现)
    // 获取所有偶数行和偶数列的元素
    // 这需要先创建一个 layout_stride::mapping，然后用它来构造 mdspan
    // 或者，利用 submdspan 的高级切片能力 (C++23)

    // C++23 引入了带有步长的切片范围，这使得跳跃切片更加直接
    // std::tuple<Start, End, Step>
    // std::submdspan(original_matrix, std::tuple{0, 4, 2}, std::tuple{0, 5, 2}); // 并非所有编译器都已完全支持此语法
    // 另一种更通用的方式是先用 layout_stride 构造新的 mdspan

    // 创建一个从 original_matrix 中跳过行和列的视图
    // 新的 Extents: 2x3 (原来4行，每隔一行取，就是2行；原来5列，每隔一列取，就是3列)
    std::array<std::size_t, 2> new_extents = {2, 3};
    // 原始矩阵的行步长是 5 (下一个元素的内存地址 = 当前地址 + 5 * sizeof(int))
    // 原始矩阵的列步长是 1
    // 如果要跳过一行，步长就是 2 * 5 = 10
    // 如果要跳过一列，步长就是 2 * 1 = 2
    std::array<std::size_t, 2> strides_skip = {original_matrix.stride(0) * 2, original_matrix.stride(1) * 2};

    std::layout_stride::mapping<std::extents<std::size_t, 2, 3>> skip_mapping(
        std::extents<std::size_t, 2, 3>{new_extents[0], new_extents[1]},
        strides_skip
    );

    // 基指针仍然是原始数据的起始指针
    std::mdspan<int, std::extents<std::size_t, 2, 3>, std::layout_stride> skip_matrix(
        data.data(), skip_mapping
    );
    print_mdspan("Slice: Skip Rows/Cols (0,0), (0,2), (0,4); (2,0), (2,2), (2,4)", skip_matrix);
    /* Output:
    Slice: Skip Rows/Cols (0,0), (0,2), (0,4); (2,0), (2,2), (2,4) (extents: 2x3):
    10      12      14
    20      22      24
    */

    // 5. 降维切片并选择特定元素
    // 获取 original_matrix(1, 2) 这个单个元素
    // std::submdspan 返回的 mdspan 的 rank 会减少到 0
    auto element_view = std::submdspan(original_matrix, std::integral_constant<std::size_t, 1>(), std::integral_constant<std::size_t, 2>());
    std::cout << "Slice: Single Element (1,2) = " << element_view() << std::endl; // rank 0 的 mdspan 用 operator()() 访问
    std::cout << "Element view rank: " << element_view.rank() << std::endl; // 0
    std::cout << std::endl;
    /* Output:
    Slice: Single Element (1,2) = 17
    Element view rank: 0
    */

    // 6. 混合切片：例如，选择所有行，但只选择第1列到第3列的子集
    auto mixed_slice = std::submdspan(original_matrix, std::full_extent, std::pair{1, 4});
    print_mdspan("Slice: All Rows, Cols 1-3", mixed_slice);
    /* Output:
    Slice: All Rows, Cols 1-3 (extents: 4x3):
    11      12      13
    16      17      18
    21      22      23
    26      27      28
    */

    return 0;
}

这些示例展示了 submdspan 的强大和灵活性。通过组合不同的切片参数，我们可以创建出满足各种需求的多维数据视图，而这一切都无需进行任何数据拷贝，极大地提升了效率。

六、std::mdspan 在 HPC 中的高级应用与性能优化

std::mdspan 的价值不仅仅在于提供一个方便的语法，更在于它在高性能计算领域的深远影响。

6.1 与泛型算法结合

mdspan 使得编写泛型、维度无关的算法成为可能。函数可以接收一个 mdspan 模板参数，并通过其 rank() 和 extent() 方法来适应不同维度的输入。

template <class ElementType, class Extents, class LayoutPolicy, class AccessorPolicy>
void scale_mdspan(std::mdspan<ElementType, Extents, LayoutPolicy, AccessorPolicy> m, ElementType factor) {
    // 遍历所有元素，无论维度如何
    // m.size() 返回总元素数
    // m.rank() 返回维度数
    // m.extent(d) 返回第 d 维的大小
    // m(idx0, idx1, ...) 访问元素

    // 对于简单的线性遍历，可以利用 mdspan 的迭代器或直接使用布局映射的线性索引
    // 但更常见的做法是嵌套循环，以保持缓存局部性
    if constexpr (Extents::rank() == 2) {
        for (std::size_t i = 0; i < m.extent(0); ++i) {
            for (std::size_t j = 0; j < m.extent(1); ++j) {
                m(i, j) *= factor;
            }
        }
    } else if constexpr (Extents::rank() == 3) {
        for (std::size_t i = 0; i < m.extent(0); ++i) {
            for (std::size_t j = 0; j < m.extent(1); ++j) {
                for (std::size_t k = 0; k < m.extent(2); ++k) {
                    m(i, j, k) *= factor;
                }
            }
        }
    } else {
        // 对于更高维度或未知维度，可以递归或使用线性迭代器（如果布局允许）
        // 对于 layout_right 和 layout_left，线性遍历是可行的
        // for (std::size_t i = 0; i < m.size(); ++i) { /* ... */ }
        // 但为了缓存局部性，通常还是推荐嵌套循环
        std::cerr << "Scaling for rank " << m.rank() << " not explicitly optimized." << std::endl;
        // Fallback or error
    }
}

// 调用示例
// scale_mdspan(original_matrix, 2.0);
// scale_mdspan(sub_matrix, 0.5);

6.2 异构计算接口 (CUDA/OpenCL)

mdspan 非常适合作为 GPU 核函数的输入。其非拥有视图的特性意味着我们可以将 mdspan 指向 GPU 内存中的数据。

1. 数据传输： 首先，将数据从 CPU 传输到 GPU 设备内存。
2. 创建设备 mdspan： 使用 mdspan 包装设备内存指针。如果需要，可以定义自定义 AccessorPolicy 来处理设备内存访问的特殊性（例如，在某些框架中可能需要特定的指针类型或修饰符）。
3. 传递给核函数： 将设备 mdspan 对象（或其底层指针和维度信息）传递给 CUDA 核函数。核函数可以直接使用 mdspan 的多维索引访问语义。

// 伪代码示例：CUDA集成
/*
#include <cuda_runtime.h>
#include <mdspan>

// 假设有一个自定义访问器，用于标记数据在设备内存上
template <class ElementType>
struct device_accessor {
    using element_type = ElementType;
    using reference = ElementType&; // 注意：在设备代码中这可能是 __device__ 引用
    using pointer = ElementType*;   // 在设备代码中这可能是 __device__ 指针

    constexpr reference access(pointer p, size_t i) const noexcept {
        return p[i];
    }
    constexpr pointer offset(pointer p, size_t i) const noexcept {
        return p + i;
    }
};

// CUDA 核函数
__global__ void add_vectors_kernel(
    std::mdspan<float, std::extents<std::size_t, std::dynamic_extent>, std::layout_right, device_accessor<float>> a,
    std::mdspan<float, std::extents<std::size_t, std::dynamic_extent>, std::layout_right, device_accessor<float>> b,
    std::mdspan<float, std::extents<std::size_t, std::dynamic_extent>, std::layout_right, device_accessor<float>> c
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < a.extent(0)) {
        c(idx) = a(idx) + b(idx);
    }
}

int main() {
    size_t N = 1024;
    std::vector<float> h_a(N), h_b(N), h_c(N); // Host data
    // ... initialize h_a, h_b ...

    float *d_a, *d_b, *d_c; // Device data
    cudaMalloc(&d_a, N * sizeof(float));
    cudaMalloc(&d_b, N * sizeof(float));
    cudaMalloc(&d_c, N * sizeof(float));

    cudaMemcpy(d_a, h_a.data(), N * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b.data(), N * sizeof(float), cudaMemcpyHostToDevice);

    // 创建设备 mdspan
    std::mdspan<float, std::extents<std::size_t, std::dynamic_extent>, std::layout_right, device_accessor<float>>
        d_mdspan_a(d_a, N), d_mdspan_b(d_b, N), d_mdspan_c(d_c, N);

    // 启动核函数
    int blockSize = 256;
    int numBlocks = (N + blockSize - 1) / blockSize;
    add_vectors_kernel<<<numBlocks, blockSize>>>(d_mdspan_a, d_mdspan_b, d_mdspan_c);
    cudaDeviceSynchronize();

    cudaMemcpy(h_c.data(), d_c, N * sizeof(float), cudaMemcpyDeviceToHost);
    // ... process h_c ...

    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
    return 0;
}
*/

注意： 上述 CUDA 示例中的 device_accessor 和 std::mdspan 在设备代码中直接使用，需要编译器（如 NVCC）的支持。C++23 mdspan 的设计确实考虑了异构计算，但具体的设备支持细节取决于编译器和平台。

6.3 缓存优化

选择正确的布局策略对缓存性能至关重要。

• *矩阵乘法 (`C = A B`)：**
  • 如果 A 是行主序，B 是列主序，那么 C 也是行主序，这可能是一个较好的组合，因为 A 的行和 B 的列都可以顺序访问。
  • 如果所有矩阵都是行主序，那么在计算 C(i,j) = sum(A(i,k) * B(k,j)) 时，B(k,j) 的访问会跳跃式地访问内存（因为 j 变化时，k 保持不变，在行主序中 B(k,j) 和 B(k,j+1) 是相邻的，但 B(k,j) 和 B(k+1,j) 则不是）。
  • 通过将 B 视为列主序 (layout_left) 或对其进行转置视图，可以显著改善缓存局部性。

示例：矩阵乘法与布局策略

假设我们计算 C = A * B，其中 A 是 M x K 矩阵，B 是 K x N 矩阵，C 是 M x N 矩阵。
标准的朴素矩阵乘法算法：

for (int i = 0; i < M; ++i) {
    for (int j = 0; j < N; ++j) {
        C(i, j) = 0;
        for (int k = 0; k < K; ++k) {
            C(i, j) += A(i, k) * B(k, j);
        }
    }
}

如果 A, B, C 都是 layout_right (行主序)：

• A(i,k) 访问良好（k 变化，最后维度）。
• B(k,j) 访问糟糕（k 变化，第一维度）。
• C(i,j) 写入良好（j 变化，最后维度）。

优化方案：

1. 将 B 存储为 layout_left (列主序)： 这样 B(k,j) 访问时，k 变化是按列访问，可以利用缓存局部性。

#include <iostream>
#include <vector>
#include <numeric>
#include <chrono>
#include <mdspan>

// namespace stdex = std::experimental;

using namespace std::chrono;

// 辅助函数：填充 mdspan
template <class MDSpan>
void fill_mdspan(MDSpan& m, int start_val) {
    // 假设是 2D mdspan
    for (std::size_t i = 0; i < m.extent(0); ++i) {
        for (std::size_t j = 0; j < m.extent(1); ++j) {
            m(i, j) = start_val++;
        }
    }
}

// 朴素矩阵乘法函数
template <class MatA, class MatB, class MatC>
void multiply_matrices(const MatA& A, const MatB& B, MatC& C) {
    for (std::size_t i = 0; i < A.extent(0); ++i) { // M
        for (std::size_t j = 0; j < B.extent(1); ++j) { // N
            C(i, j) = 0;
            for (std::size_t k = 0; k < A.extent(1); ++k) { // K
                C(i, j) += A(i, k) * B(k, j);
            }
        }
    }
}

int main() {
    const std::size_t M = 500, K = 500, N = 500;

    // 1. All Row-major
    std::vector<double> data_A_rm(M * K);
    std::vector<double> data_B_rm(K * N);
    std::vector<double> data_C_rm(M * N);

    std::mdspan<double, std::extents<std::size_t, M, K>> A_rm(data_A_rm.data());
    std::mdspan<double, std::extents<std::size_t, K, N>> B_rm(data_B_rm.data());
    std::mdspan<double, std::extents<std::size_t, M, N>> C_rm(data_C_rm.data());

    fill_mdspan(A_rm, 1);
    fill_mdspan(B_rm, 100);

    auto start_rm = high_resolution_clock::now();
    multiply_matrices(A_rm, B_rm, C_rm);
    auto end_rm = high_resolution_clock::now();
    std::cout << "All Row-major time: " << duration_cast<milliseconds>(end_rm - start_rm).count() << " ms" << std::endl;
    // Expected output: ~1000-2000 ms for 500x500 matrices (depends on hardware)

    // 2. A and C Row-major, B Column-major
    std::vector<double> data_B_cm(K * N); // Re-initialize for B in column-major
    std::mdspan<double, std::extents<std::size_t, K, N>, std::layout_left> B_cm(data_B_cm.data());
    fill_mdspan(B_cm, 100); // Fill B in column-major logical order

    auto start_cm = high_resolution_clock::now();
    multiply_matrices(A_rm, B_cm, C_rm); // A_rm (Row), B_cm (Col), C_rm (Row)
    auto end_cm = high_resolution_clock::now();
    std::cout << "A(Row), B(Col), C(Row) time: " << duration_cast<milliseconds>(end_cm - start_cm).count() << " ms" << std::endl;
    // Expected output: Significantly faster, ~500-1000 ms

    return 0;
}

注意： 实际性能提升会受到编译器优化、CPU 架构、缓存大小、以及其他系统负载的影响。但通常情况下，B 采用列主序（或对 B 进行转置）会带来显著的性能提升。

6.4 与现有库的互操作性

mdspan 可以轻松地与 BLAS (Basic Linear Algebra Subprograms) 等库进行互操作。许多 BLAS 函数期望接收指向数据起始位置的指针以及每个维度的步长。mdspan 可以直接提供这些信息：

• m.data() 返回底层数据指针。
• m.extent(d) 返回维度 d 的大小。
• m.stride(d) 返回维度 d 的步长（以元素为单位）。

例如，一个 GEMM (Generalized Matrix Multiplication) 函数可能需要 lda, ldb, ldc (leading dimensions)。如果 mdspan 是 layout_right，其 stride(0) 就是 lda。如果 mdspan 是 layout_left，其 stride(1) 就是 lda。mdspan 提供了 mapping().stride(d) 和 mapping().required_span_size() 等方法来获取这些底层布局信息。

七、std::mdspan 与其他 C++ 特性的结合

1. std::span： std::mdspan 可以看作是 std::span 的多维泛化。std::span 提供了对一维连续内存的非拥有视图。当 mdspan 的 rank() 为 1 时，它在很多方面与 std::span 行为相似，甚至可以从 std::span 构造 mdspan。

2. std::vector / std::array： mdspan 通常以 std::vector 或 std::array 作为其底层数据源，利用它们的内存管理能力。

3. Concepts (C++20)： mdspan 大量使用了 Concepts 来约束其模板参数，确保类型安全并提供更好的编译错误信息。开发者在编写泛型算法时，也可以利用 Concepts 来约束 mdspan 参数，例如要求其 rank() 或 extent() 满足特定条件。

   template <typename T>
   concept IsRank2MDSPAN = std::is_mdspan_v<T> && (T::rank() == 2);
   
   template <IsRank2MDSPAN Mat>
   void process_2d_matrix(Mat m) {
       // 保证 m 是一个 2 维 mdspan
       // ...
   }

4. Ranges (C++20)： 虽然 mdspan 本身不直接是 Range，但未来可能会有适配器将 mdspan 转化为 Range，从而可以利用 Range-V3 或 C++20 Ranges 库的强大功能进行链式操作和转换。

八、潜在的陷阱与注意事项

1. 生命周期管理： mdspan 是一个视图，它不拥有底层数据。因此，必须确保 mdspan 对象在底层数据（如 std::vector）的生命周期内保持有效。如果底层数据被销毁，mdspan 将指向无效内存，导致未定义行为。

   // 错误示例：悬空 mdspan
   std::mdspan<int, std::extents<std::size_t, 2, 2>> dangling_mdspan;
   {
       std::vector<int> data(4);
       dangling_mdspan = std::mdspan<int, std::extents<std::size_t, 2, 2>>(data.data());
   } // data 在这里被销毁，dangling_mdspan 变得无效
   // dangling_mdspan(0,0); // 未定义行为

2. 线程安全： mdspan 本身不提供任何线程安全保证。如果多个线程同时读写同一个 mdspan 的共享数据，需要外部同步机制（如互斥锁、原子操作）来避免数据竞争。

3. 维度不匹配与索引越界： mdspan 在编译期对静态维度进行严格检查。对于动态维度，默认情况下不强制进行运行时边界检查，以保持零开销。在调试模式下或通过特定构建选项，可能会启用运行时检查。始终确保索引操作在有效范围内。

4. 性能分析： 尽管 mdspan 被设计为零开销，但在复杂的场景下，性能瓶颈可能仍然存在于算法本身或缓存访问模式。始终进行基准测试和性能分析，以验证优化效果。

九、展望与未来发展

std::mdspan 的引入是 C++ 在高性能计算领域迈出的重要一步。它不仅为多维数据处理提供了标准化的、类型安全的、零开销的解决方案，也为 C++ 语言在科学计算、机器学习等数据密集型应用中提供了更强大的工具。

未来的发展方向包括：

• std::mdarray： 正在提案中的 std::mdarray 将是 mdspan 的完美搭档。mdarray 将提供拥有内存的多维数组，与 mdspan 共同构成完整的、标准的多维数据管理方案。mdarray 负责内存管理，mdspan 负责提供视图和切片。
• 与并行编程模型更紧密的集成： mdspan 的设计使其易于与 OpenMP、MPI、SYCL 等并行编程模型结合，简化在异构和分布式系统上处理多维数据的复杂性。
• 更高级的切片功能： 可能会有更丰富的切片语法或辅助函数，例如直接支持负索引、更复杂的跳跃模式等。
• 与 P2300 (std::execution) 结合： 随着 C++ 标准库中并行执行模型的发展，mdspan 有望与此结合，为多维数据上的并行算法提供更自然的表达方式。

结语

std::mdspan 是 C++23 献给高性能计算领域的一份厚礼。它以其卓越的类型安全、零开销抽象和无与伦比的灵活性，彻底改变了 C++ 中处理多维数据的方式。掌握并善用 std::mdspan，将使我们能够编写出更安全、更高效、更易于维护的 C++ HPC 应用程序，更好地驾驭现代计算的复杂性。现在，我们有了更强大的工具，去探索更高维度的数据世界！