CUDA编程实战04：GPU 内存访问与合并——为什么相同的加法会有不同性能

CUDA编程实战04：GPU 内存访问与合并——为什么相同的加法会有不同性能

系列定位：从能运行的 CUDA 程序出发，逐步建立性能分析与工程优化能力。
本篇适合：已经能编译并运行 Kernel，做过向量加法，但还不理解“为什么只改下标就可能慢很多”的读者。
本篇成果：3 个可直接构建的 CUDA 基准程序、4 张精确原理图、1 张主题插图、完整验证方法、CUDA Event 计时方法和 Nsight Compute 分析路线。
代码目录：blogs/code/04。
本文讨论的 32 字节 Sector 模型主要面向计算能力 6.0 及以上设备；不同架构与工具版本的细节应以对应版本的 NVIDIA 官方文档为准。

---

开篇：代码没有多做一次加法，为什么会慢七倍

先看两段几乎一样的 CUDA 代码。

第一段：

output[i] = input[i];

第二段：

output[i] = input[(i * stride) & mask];

两段代码都只读取一个 float，都只写出一个 float，都由同样多的线程执行，也没有增加复杂数学运算。初学者很容易作出一个看似合理的判断：既然每个线程的工作量一样，运行时间就应该差不多。

然而，本篇附带的 stride_access_benchmark.cu 在一张 NVIDIA GeForce RTX 2060 上得到下面这组 Release 实测结果：

Stride	相邻 Lane 输入地址间隔	Kernel 平均时间	逻辑有效带宽
1	4 B	0.187500 ms	178.96 GB/s
3	12 B	0.331858 ms	101.11 GB/s
5	20 B	0.487751 ms	68.79 GB/s
9	36 B	0.806673 ms	41.60 GB/s
17	68 B	1.316442 ms	25.49 GB/s
33	132 B	1.323668 ms	25.35 GB/s
65	260 B	1.335706 ms	25.12 GB/s

从 Stride 为 1 到 Stride 为 65，数学工作没有变化，输出长度没有变化，访问的输入元素总数也没有变化，平均时间却从约 0.19 毫秒增加到约 1.34 毫秒。两者相差约七倍。

这正是本篇要解决的实际问题：

当 Kernel 已经正确、线程数量也足够多时，怎样判断性能问题是不是来自 Global Memory 的地址模式？怎样用可重复实验，而不是凭感觉，证明相邻线程访问相邻数据更高效？

读完本文，你不应该只会背诵“合并访问更快”。你应该能够拿到一个图像、粒子、矩阵或业务数组 Kernel，完成下面的闭环：

1. 写出一个 Warp 内 32 个 Lane 的实际地址；
2. 判断这些地址覆盖多少个 32 字节 Sector；
3. 区分有用字节、缓存流量与 DRAM 流量；
4. 设计只改变一个变量的基准；
5. 使用 CUDA Event 测 Kernel 时间；
6. 用 CPU 参考或等价结果验证正确性；
7. 使用 Nsight Compute 检查请求数与 Sector 数；
8. 根据访问模式选择 AoS、SoA、转置、分块或共享内存；
9. 知道什么时候不能仅凭一组毫秒数下结论。

GPU 内存系统搬运数据时不会为每个 float 单独派出一辆只装四字节的车。一个 Warp 的请求会按地址落点组织成若干次内存访问。请求越集中，搬来的数据越容易被真正使用；请求越分散，就越可能搬回大量当前线程并不需要的字节。

---

一、先建立本篇的实验纪律

性能文章最容易出现的问题，是先看到一个结果，再为它编故事。例如，版本 A 比版本 B 慢，就直接宣布“版本 A 没有合并”；或者看到理论上多一个 Sector，就断言时间一定增加百分之二十五。这样得到的不是分析能力，而是碰巧说对。

本篇采用三层证据：

第一层：地址推理
    一个 Warp 的 32 个 Lane 到底访问哪些地址？

第二层：时间实验
    在相同数据量、相同运算和相同输出下，CUDA Event 时间是否变化？

第三层：硬件计数
    Nsight Compute 看到的请求数、Sector 数和利用率是否支持推理？

这三层分别回答不同问题。

• 地址推理回答“按照程序语义，可能需要覆盖多少内存区域”；
• Event 计时回答“这个完整 Kernel 在当前机器上实际用了多久”；
• Profiler 计数回答“硬件执行期间，内存层级实际处理了怎样的请求”。

只看地址，不知道缓存和调度怎样改变最终时间；只看时间，不知道变慢究竟来自地址、指令还是频率；只看 Profiler，也可能因为选错 Kernel、采样范围或指标而误读。三层证据结合，才是可以迁移到真实项目的方法。

为了让实验公平，本文所有基准都遵守这些规则：

• 使用 Release 构建；
• 第一次正式计时前先预热；
• 用 CUDA Event 记录 Device 时间；
• 同一配置重复执行多轮；
• 正式计时不包含 cudaMalloc 和 Host/Device 复制；
• 每一组结果都做全量校验；
• 对比时保持 N、线程配置、输出写入量和数学表达式不变；
• 明确“有效带宽”的字节口径；
• 不把一次运行结果当成所有 GPU 的固定结论；
• 不把逻辑请求字节误写成物理 DRAM 字节。

这套纪律比任何一条“优化口诀”都更重要。因为真实项目里的性能变化往往同时受到缓存、寄存器、占用率、指令数、数据布局和工作频率影响。只有控制变量，才能知道自己究竟优化了什么。

---

二、从第三篇的向量加法继续前进

第三篇使用了经典向量加法：

__global__ void vector_add(const float* a,
                           const float* b,
                           float* c,
                           size_t n) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;

    if (i < n) {
        c[i] = a[i] + b[i];
    }
}

第 i 个线程读取 a[i]、读取 b[i]、写入 c[i]。假设一个 Warp 里的 Lane 编号是 0 到 31，并且这个 Warp 从数组某个自然对齐的位置开始，那么它们产生的访问近似如下：

Lane 0  -> a[k + 0]
Lane 1  -> a[k + 1]
Lane 2  -> a[k + 2]
...
Lane 31 -> a[k + 31]

float 占 4 字节，相邻 Lane 的地址相差 4 字节。32 个 Lane 总共需要：

$$B_{useful}=32\times 4=128$$

也就是说，这个 Warp 对数组 A 的一次读取需要 128 字节有用数据。数组 B 的读取方式相同，数组 C 的写入方式也相同。第三篇中的向量加法天然具备适合合并的地址模式，这也是为什么它常被用作最小 CUDA 带宽示例。

现在只改一处下标：

c[i] = a[i * stride] + b[i * stride];

当 stride 大于 1，相邻 Lane 不再读取相邻元素。算术仍是一加法，但内存系统看到的地址集合已经完全不同。CUDA 性能优化的一个重要转折点就在这里：

GPU 执行的不是抽象的“数组元素”，而是指令、地址与数据事务。源代码看起来只差一个下标，硬件工作量却可能差很多。

---

三、Global Memory 不是一次只搬一个 float

为了建立直觉，可以把 Global Memory 想象成距离计算单元较远的大型仓库。线程提出请求时，仓库不会把每个四字节元素装进独立信封逐个送来，而会按一定粒度组织传输。对计算能力 6.0 及以上架构，理解普通 Global Memory 合并访问时，可以从 32 字节事务或 Sector 入手。

需要先区分五个经常混在一起的词。

3.1 元素

元素是程序层概念。例如 float 元素占 4 字节，double 占 8 字节，一个含四个 float 的粒子结构体通常占 16 字节。Kernel 中写的 input[i] 表示某个元素，但硬件最终处理的是它对应的字节地址。

3.2 Lane 请求

Warp 由 32 个 Lane 组成。执行一条 Global Memory 读指令时，当前活跃的 Lane 各自产生地址。某些 Lane 可能因为分支或边界判断不活跃，因此不是每次都恰好有 32 个有效请求。

3.3 Request

Profiler 中的 Request 通常表示某个内存单元收到的一次请求。具体指标命名和聚合层级会随架构与 Nsight Compute 版本变化，不能脱离指标说明机械理解。

3.4 Sector

在 Nsight Compute 的内存分析中，Sector 常以 32 字节为粒度。一个请求可能需要一个或多个 Sector。如果 32 个 Lane 的地址集中在连续 128 字节范围中，就可能由四个 32 字节 Sector 覆盖。

3.5 DRAM 流量

DRAM 流量是最终到显存层级的物理数据移动。它不必等于 Lane 请求字节，也不必等于某一级缓存看到的 Sector 总量。L1 或 L2 命中、写策略、重放以及相邻 Warp 的复用都会改变它。

因此，下面三个数不能随意画等号：

Kernel 语义上需要的有用字节
某一级缓存处理的 Sector 字节
最终到达 DRAM 的物理字节

本文使用“地址级利用率”描述单 Warp 地址覆盖的理想化比例，使用“逻辑有效带宽”描述程序有用字节除以 Event 时间。它们都不是直接测得的 DRAM 总线占用。

---

四、一个 Warp 怎样合并成四个 Sector

看图时按照下面的顺序，不要一开始就背数字。

第一步，确认数据类型。这里每个 Lane 读取一个 float，每个请求是 4 字节。

第二步，列出地址。Lane 0 到 Lane 31 读取连续的 float 0 到 float 31。

第三步，把地址落入 32 字节边界。一个 Sector 可以容纳八个连续 float：

$$N_{float\_per\_sector}=\frac{32}{4}=8$$

第四步，统计覆盖的 Sector。32 个连续 float 恰好覆盖四个 Sector：

$$N_{sector}=\frac{32\times 4}{32}=4$$

第五步，计算地址级利用率。定义：

$${\eta }_{address}=\frac{B_{useful}}{B_{covered}}$$

对自然对齐的连续访问：

$${\eta }_{address}=\frac{128}{4\times 32}=1$$

这个结果表示，在该简化模型中，被覆盖的 128 字节全是当前 Warp 需要的数据。它不表示整个 GPU 已达到理论显存带宽，也不表示一定只有四条机器指令，更不表示缓存层级只有一种请求。它只是说明地址分布非常紧凑。

4.1 合并访问的准确表述

可以把合并访问理解为：

同一个 Warp 在执行某条 Global Memory 指令时，各活跃 Lane 的地址能够由较少且利用率较高的内存事务覆盖。

这个定义有三个关键词。

“同一个 Warp”：不能把不同 Warp 的线程混成一组做最初判断。

“某条内存指令”：同一个 Kernel 可能有多次读取与写入，每一条都要分别分析。

“地址集合”：线程编号连续不等于地址连续。最终地址取决于索引公式、数据布局、类型大小和基址偏移。

4.2 合并访问不等于必须按 Lane 编号递增

“Lane 0 读元素 0，Lane 1 读元素 1”是最容易理解的合并模式，但不能把它误写成唯一合法模式。某些地址重排、广播或缓存访问也可能以较少事务完成。入门阶段最可靠的工程习惯仍然是：

• 让相邻 Lane 尽量访问相邻地址；
• 让数组基址和每行起点保持合理对齐；
• 优先使用连续的标量或向量化数据布局；
• 遇到复杂索引时列出一个 Warp 的地址表；
• 最后用 Profiler 验证，而不是只凭代码外观。

---

五、先学会“地址审计”，再谈优化

拿到一个 Kernel 时，可以使用下面的地址审计模板。

假设 Warp 起始全局线程编号是 warp_base，Lane 编号是 lane：

$$i=warp\_base+lane$$

然后把 Kernel 中每次 Global Memory 访问写成 lane 的函数。例如：

input[i]

对应：

$$index(lane)=warp\_base+lane$$

如果访问：

input[i * stride]

对应：

$$index(lane)=(warp\_base+lane)\times stride$$

相邻 Lane 的元素间隔为：

$$\Delta index=stride$$

相邻 Lane 的字节间隔为：

$$\Delta address=stride\times sizeof(T)$$

如果 T 是 float，Stride 为 8：

$$\Delta address=8\times 4=32$$

这意味着相邻 Lane 的请求刚好落到相邻的 32 字节 Sector 中。一个 Warp 的 32 个 Lane 可能覆盖 32 个不同 Sector，而每个 Sector 当前只使用 4 字节。

地址审计还要检查下面的问题：

1. Warp 是否跨越数组行边界？
2. 行首的 Pitch 是否是元素宽度还是字节宽度？
3. 结构体是否含有填充与对齐？
4. 指针是否做过 + offset？
5. 每个 Lane 读取一个标量还是多个字段？
6. 分支后有哪些 Lane 仍然活跃？
7. 索引是否发生取模、哈希或随机置换？
8. 输入读取与输出写入是否具有相同模式？
9. 数据是否会在缓存中被后续 Warp 重用？
10. 类型转换后真实加载宽度是否变化？

当你能把“看上去复杂的业务代码”压缩成一个 Warp 的地址表时，内存优化就不再神秘。

---

六、实验一：Offset 只偏移一个 float

Offset 实验不改变相邻 Lane 之间的间隔，只改变整个 Warp 的起始位置。

连续对齐版本：

output[i] = input[i];

偏移一个元素：

output[i] = input[i + 1];

两者都让相邻 Lane 访问相邻 float。差别在于：如果 input 的基址按 32 字节对齐，第一种访问可以从 Sector 边界开始；第二种访问会从边界之后 4 字节开始。

对一个完整 Warp，Offset 为 0 时读取 float 0 到 float 31，覆盖四个 Sector。Offset 为 1 时读取 float 1 到 float 32，覆盖五个 Sector：

$$B_{covered,0}=4\times 32=128$$

$$B_{covered,1}=5\times 32=160$$

有用数据都为 128 字节，因此：

$${\eta }_0=\frac{128}{128}=100\%$$

$${\eta }_1=\frac{128}{160}=80\%$$

注意，这里只能得出“单个 Warp 的地址覆盖从四个 Sector 变成五个 Sector”。不能直接得出“整个 Kernel 必慢百分之二十五”。为什么？

因为前一个 Warp 多覆盖的末尾 Sector，可能正好是后一个 Warp 需要的开头 Sector。如果该数据仍在缓存里，后一个 Warp 可以复用。实际硬件还会受请求调度、L1/L2 行为、显存分区映射、动态频率与并发 Warp 影响。

这正是 Offset 实验的价值：它迫使我们区分“地址模型”与“端到端时间”。

---

七、Offset 基准程序逐段解析

完整程序位于：

blogs/code/04/offset_access_benchmark.cu

7.1 Kernel

__global__ void copy_with_offset(const float* input,
                                 float* output,
                                 size_t n,
                                 size_t offset) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;

    if (i < n) {
        output[i] = input[i + offset];
    }
}

这个 Kernel 故意不做加法。原因不是复制比加法高级，而是我们希望把变量减少到最少：只有输入读取的起始地址变化。输出始终是连续写入，因此每组的写入模式相同。

7.2 为什么输入要多分配七个元素

程序测试 Offset 0 到 7。最大访问位置为：

$$i_{max}=N-1+7$$

因此输入长度必须至少为：

$$N_{input}=N+7$$

如果仍只分配 N 个元素，Offset 大于 0 时末尾线程会越界。越界程序可能暂时没有崩溃，但结果和性能都失去意义。性能基准必须先正确，再计时。

7.3 为什么每个 Offset 都预热

第一次 Kernel 启动可能包含上下文初始化影响，GPU 频率也可能尚未稳定。程序对每个 Offset 先执行十次：

for (int warmup = 0; warmup < 10; ++warmup) {
    copy_with_offset<<<blocks, threads>>>(
        d_input, d_output, n, offset);
}

预热不是万能的。消费级 GPU 仍可能因为温度、功耗和桌面图形负载产生波动。因此严谨报告应该保留多次运行，必要时取中位数，并记录 GPU 名称、驱动、Toolkit、功耗模式与温度。

7.4 Event 计时范围

正式计时：

cudaEventRecord(start);

for (int iteration = 0; iteration < iterations; ++iteration) {
    copy_with_offset<<<blocks, threads>>>(
        d_input, d_output, n, offset);
}

cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&total_ms, start, stop);

平均时间：

$$t_{avg}=\frac{t_{total}}{R}$$

其中 R 是迭代次数。循环内没有 cudaMemcpy、Host 验证和内存分配，因此结果主要描述重复 Kernel 的 Device 执行时间。

7.5 有效带宽怎样计算

每个输出元素对应：

• 读取一个 float，4 字节；
• 写入一个 float，4 字节。

逻辑有用字节：

$$B_{logical}=N\times 2\times 4$$

若平均时间以毫秒表示，有效带宽为：

$$B{W}_{effective}=\frac{B_{logical}}{t_{avg}\times 10^6}$$

结果单位是 GB/s，因为分母中的毫秒乘以一百万，等价于把字节每毫秒换算成十进制 GB/s。

这里特意叫“逻辑有效带宽”。Offset 为 1 时，硬件可能为有用的 128 字节覆盖更多 Sector，所以物理流量可能大于逻辑字节。把逻辑字节除以时间，可以公平比较“同样业务数据量被处理得多快”，但不能直接声称得到真实 DRAM 流量。

7.6 全量验证

每个 Offset 完成后，程序把输出复制回 Host，并检查：

for (size_t i = 0; i < n; ++i) {
    if (h_output[i] != h_input[i + offset]) {
        ++mismatches;
    }
}

这里使用精确比较是安全的，因为 Device 只复制已经存在的 float，没有浮点运算舍入。如果 Kernel 做了加法、乘法或函数计算，就应使用绝对误差与相对误差策略，而不是盲目使用 ==。

---

八、Offset 完整源码

下面的文件可以单独保存为 offset_access_benchmark.cu。

#include <cuda_runtime.h>

#include <cstdio>
#include <cstdlib>
#include <vector>

void cuda_check_impl(cudaError_t error, const char* expression,
                     const char* file, int line) {
    if (error != cudaSuccess) {
        std::fprintf(stderr,
                     "CUDA error at %s:%d\n"
                     "  expression: %s\n"
                     "  reason: %s\n",
                     file, line, expression,
                     cudaGetErrorString(error));
        std::exit(EXIT_FAILURE);
    }
}

#define CUDA_CHECK(call) \
    cuda_check_impl((call), #call, __FILE__, __LINE__)

__global__ void copy_with_offset(const float* input, float* output,
                                 size_t n, size_t offset) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (i < n) {
        output[i] = input[i + offset];
    }
}

int main(int argc, char** argv) {
    size_t n = 1ULL << 24;
    int iterations = 100;

    if (argc >= 2) {
        const unsigned long long parsed =
            std::strtoull(argv[1], nullptr, 10);
        if (parsed == 0 || parsed > 50000000ULL) {
            std::fprintf(stderr, "N must be in [1, 50000000].\n");
            return EXIT_FAILURE;
        }
        n = static_cast<size_t>(parsed);
    }

    if (argc >= 3) {
        const long parsed = std::strtol(argv[2], nullptr, 10);
        if (parsed <= 0 || parsed > 10000) {
            std::fprintf(stderr,
                         "iterations must be in [1, 10000].\n");
            return EXIT_FAILURE;
        }
        iterations = static_cast<int>(parsed);
    }

    constexpr size_t max_offset = 7;
    constexpr int threads = 256;
    const int blocks =
        static_cast<int>((n + threads - 1) / threads);
    const size_t input_bytes =
        (n + max_offset) * sizeof(float);
    const size_t output_bytes = n * sizeof(float);

    std::vector<float> h_input(n + max_offset);
    std::vector<float> h_output(n);

    for (size_t i = 0; i < h_input.size(); ++i) {
        h_input[i] =
            static_cast<float>((i * 17 + 3) % 8191) * 0.001f;
    }

    float* d_input = nullptr;
    float* d_output = nullptr;

    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_input), input_bytes));
    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_output), output_bytes));
    CUDA_CHECK(cudaMemcpy(
        d_input, h_input.data(), input_bytes,
        cudaMemcpyHostToDevice));

    cudaEvent_t start{};
    cudaEvent_t stop{};
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));

    int device = 0;
    cudaDeviceProp prop{};
    CUDA_CHECK(cudaGetDevice(&device));
    CUDA_CHECK(cudaGetDeviceProperties(&prop, device));

    std::printf("GPU: %s\n", prop.name);
    std::printf(
        "offset,average_ms,useful_bandwidth_GBps,mismatches\n");

    bool all_passed = true;

    for (size_t offset = 0; offset <= max_offset; ++offset) {
        for (int warmup = 0; warmup < 10; ++warmup) {
            copy_with_offset<<<blocks, threads>>>(
                d_input, d_output, n, offset);
        }
        CUDA_CHECK(cudaGetLastError());
        CUDA_CHECK(cudaDeviceSynchronize());

        CUDA_CHECK(cudaEventRecord(start));
        for (int iteration = 0;
             iteration < iterations;
             ++iteration) {
            copy_with_offset<<<blocks, threads>>>(
                d_input, d_output, n, offset);
        }
        CUDA_CHECK(cudaGetLastError());
        CUDA_CHECK(cudaEventRecord(stop));
        CUDA_CHECK(cudaEventSynchronize(stop));

        float total_ms = 0.0f;
        CUDA_CHECK(cudaEventElapsedTime(
            &total_ms, start, stop));

        const double average_ms =
            static_cast<double>(total_ms) / iterations;
        const double useful_bytes =
            static_cast<double>(n) * 2.0 * sizeof(float);
        const double useful_bandwidth_gb_s =
            useful_bytes / (average_ms * 1.0e6);

        CUDA_CHECK(cudaMemcpy(
            h_output.data(), d_output, output_bytes,
            cudaMemcpyDeviceToHost));

        size_t mismatches = 0;
        for (size_t i = 0; i < n; ++i) {
            if (h_output[i] != h_input[i + offset]) {
                ++mismatches;
            }
        }

        all_passed =
            all_passed && mismatches == 0;

        std::printf("%zu,%.6f,%.2f,%zu\n",
                    offset, average_ms,
                    useful_bandwidth_gb_s, mismatches);
    }

    CUDA_CHECK(cudaEventDestroy(stop));
    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaFree(d_output));
    CUDA_CHECK(cudaFree(d_input));

    std::printf("\nVerification: %s\n",
                all_passed ? "PASSED" : "FAILED");
    return all_passed ? EXIT_SUCCESS : EXIT_FAILURE;
}

---

九、Offset 本机实测：为什么对齐版本没有稳定拿第一

测试环境：

GPU                 : NVIDIA GeForce RTX 2060
Compute Capability  : 7.5
CUDA Toolkit        : 11.6
Compiler            : nvcc 11.6 + Visual Studio 2022
Build               : Release
N                   : 16,777,216 floats
Iterations          : 100
Threads per block   : 256

一次完整运行得到：

offset,average_ms,useful_bandwidth_GBps,mismatches
0,0.719412,186.57,0
1,0.722773,185.70,0
2,0.629693,213.15,0
3,0.647573,207.26,0
4,0.644170,208.36,0
5,0.644321,208.31,0
6,0.643604,208.54,0
7,0.643204,208.67,0

Verification: PASSED

如果只背过“对齐一定最快”，这组数据会让人困惑：Offset 0 并没有稳定胜过 Offset 2 到 7。

正确处理方式不是删掉不符合预期的数据，也不是立刻宣布合并理论错误，而是按证据层级解释。

第一，地址推理仍然成立。在假设基址按 32 字节对齐时，单 Warp 的 Offset 0 覆盖四个 Sector，Offset 1 到 7 往往覆盖五个 Sector。

第二，相邻 Warp 会共享边界 Sector。Offset 版本额外覆盖的末尾部分可能被下一个 Warp 使用，因此缓存能够减轻理论浪费。

第三，Event 时间是整个 Kernel 的综合结果。它同时反映输入读、连续输出写、指令调度、缓存、显存分区、工作频率和桌面系统噪声。

第四，这些时间低于一毫秒。虽然重复一百次已降低启动噪声，但消费级 GPU 的动态频率仍会产生可见差异。更严格实验应多次启动程序、改变测试顺序并报告中位数与离散程度。

第五，要证明“某一级缓存的 Sector 数确实不同”，需要 Nsight Compute 对具体 Kernel 启动采集计数，而不能用 Event 时间倒推。

这组结果带来的真正收获是：

理论模型用于提出可检验假设，不负责保证某台机器上的毫秒数严格按简单比例排列。

如果博客只展示一组漂亮而整齐的数据，读者到了真实项目中遇到缓存与频率影响就会无所适从。保留反直觉结果，反而更接近工程实践。

---

十、实验二：Stride 把相邻 Lane 主动拆散

Offset 只改变起点，Stride 则改变相邻 Lane 之间的地址距离。设：

source_index = i * stride;

对一个 Warp，Lane 0 到 Lane 31 的元素下标是：

0 × stride
1 × stride
2 × stride
...
31 × stride

相邻 Lane 的字节间隔：

$$\Delta address=stride\times 4$$

当 Stride 为 1，间隔为 4 字节；当 Stride 为 8，间隔为 32 字节。

在简化的单 Warp 模型中，Stride 为 8 时，每个 Lane 的一个 float 可能位于不同 Sector。32 个 Lane 需要的有用数据仍为：

$$B_{useful}=32\times 4=128$$

覆盖的 Sector 数据可能达到：

$$B_{covered}=32\times 32=1024$$

地址级利用率：

$${\eta }_{address}=\frac{128}{1024}=12.5\%$$

这不是说 Kernel 的有效带宽一定只剩峰值的八分之一。它表示从当前 Warp 的地址集合看，每个 32 字节区域只消费其中 4 字节，是一种极不紧凑的访问型态。缓存复用和其他流量仍会影响最终时间。

---

十一、怎样让 Stride 实验保持公平

很多网上的 Stride 示例写成：

output[i] = input[i * stride];

然后只分配 N 个输入元素。这样会越界。有人用：

output[i] = input[(i * stride) % n];

这样虽然不越界，但如果 N 与 Stride 有公因数，就只会反复访问部分输入。例如 N 为 16、Stride 为 4 时，访问序列只有：

0, 4, 8, 12, 0, 4, 8, 12, ...

这会让数据集突然缩小，缓存命中率大幅上升，实验不再只比较合并访问。

本篇使用两个条件：

1. N 是 2 的幂；
2. Stride 是奇数。

因为任何奇数与 2 的幂互质：

$$\gcd (stride,N)=1$$

映射：

$$j=(i\times stride)\mathrm{mod}N$$

就会形成 0 到 N 减 1 的一个排列。也就是说，每组 Stride 都会恰好访问全部 N 个输入元素一次，没有因为重复小子集而获得额外缓存优势。

当 N 是 2 的幂时，取模还可以写成位与：

source_index = (i * stride) & (n - 1);

这不是为了炫技，而是明确实验的数学结构。程序启动时会检查 N 确实是 2 的幂，测试的 Stride 列表全部为奇数。

此外，输出始终写到 output[i]，确保写入连续。这样性能变化主要来自输入读取地址，而不是读写两边同时改变。

---

十二、Stride 基准的关键代码

Kernel：

__global__ void gather_with_stride(const float* input,
                                   float* output,
                                   size_t n,
                                   size_t mask,
                                   size_t stride) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;

    if (i < n) {
        const size_t source_index = (i * stride) & mask;
        output[i] = input[source_index];
    }
}

测试步长：

const std::vector<size_t> strides{
    1, 3, 5, 9, 17, 33, 65
};

为什么不使用 2、4、8、16？因为 N 是 2 的幂，偶数 Stride 与 N 不互质，会重复访问子集。这里使用 9、17、33、65 等奇数，让相邻 Lane 的间隔跨过一个或多个 32 字节区域，同时仍保持全排列。

验证逻辑与 Device 完全相同：

for (size_t i = 0; i < n; ++i) {
    const size_t source_index = (i * stride) & mask;
    if (h_output[i] != h_input[source_index]) {
        ++mismatches;
    }
}

每个 Stride 都必须输出零个 Mismatch。否则没有资格比较性能。

---

十三、Stride 完整源码

完整文件位于 blogs/code/04/stride_access_benchmark.cu。下面的版本可以直接复制。

#include <cuda_runtime.h>

#include <cstdio>
#include <cstdlib>
#include <vector>

void cuda_check_impl(cudaError_t error,
                     const char* expression,
                     const char* file,
                     int line) {
    if (error != cudaSuccess) {
        std::fprintf(stderr,
                     "CUDA error at %s:%d\n"
                     "  expression: %s\n"
                     "  reason: %s\n",
                     file, line, expression,
                     cudaGetErrorString(error));
        std::exit(EXIT_FAILURE);
    }
}

#define CUDA_CHECK(call) \
    cuda_check_impl((call), #call, __FILE__, __LINE__)

__global__ void gather_with_stride(const float* input,
                                   float* output,
                                   size_t n,
                                   size_t mask,
                                   size_t stride) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;

    if (i < n) {
        const size_t source_index =
            (i * stride) & mask;
        output[i] = input[source_index];
    }
}

bool is_power_of_two(size_t value) {
    return value != 0 &&
           (value & (value - 1)) == 0;
}

int main(int argc, char** argv) {
    size_t n = 1ULL << 22;
    int iterations = 100;

    if (argc >= 2) {
        const unsigned long long parsed =
            std::strtoull(argv[1], nullptr, 10);

        if (parsed == 0 ||
            parsed > 67108864ULL ||
            !is_power_of_two(
                static_cast<size_t>(parsed))) {
            std::fprintf(
                stderr,
                "N must be a power of two "
                "in [1, 67108864].\n");
            return EXIT_FAILURE;
        }

        n = static_cast<size_t>(parsed);
    }

    if (argc >= 3) {
        const long parsed =
            std::strtol(argv[2], nullptr, 10);

        if (parsed <= 0 || parsed > 10000) {
            std::fprintf(
                stderr,
                "iterations must be in [1, 10000].\n");
            return EXIT_FAILURE;
        }

        iterations = static_cast<int>(parsed);
    }

    const std::vector<size_t> strides{
        1, 3, 5, 9, 17, 33, 65
    };

    constexpr int threads = 256;
    const int blocks =
        static_cast<int>(
            (n + threads - 1) / threads);
    const size_t bytes = n * sizeof(float);
    const size_t mask = n - 1;

    std::vector<float> h_input(n);
    std::vector<float> h_output(n);

    for (size_t i = 0; i < n; ++i) {
        h_input[i] =
            static_cast<float>(
                (i * 13 + 5) % 16381) * 0.0001f;
    }

    float* d_input = nullptr;
    float* d_output = nullptr;

    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_input), bytes));
    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_output), bytes));
    CUDA_CHECK(cudaMemcpy(
        d_input, h_input.data(), bytes,
        cudaMemcpyHostToDevice));

    cudaEvent_t start{};
    cudaEvent_t stop{};
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));

    int device = 0;
    cudaDeviceProp prop{};
    CUDA_CHECK(cudaGetDevice(&device));
    CUDA_CHECK(cudaGetDeviceProperties(&prop, device));

    std::printf("GPU: %s\n", prop.name);
    std::printf(
        "stride,byte_gap_per_lane,"
        "average_ms,useful_bandwidth_GBps,"
        "mismatches\n");

    bool all_passed = true;

    for (const size_t stride : strides) {
        for (int warmup = 0; warmup < 10; ++warmup) {
            gather_with_stride<<<blocks, threads>>>(
                d_input, d_output,
                n, mask, stride);
        }
        CUDA_CHECK(cudaGetLastError());
        CUDA_CHECK(cudaDeviceSynchronize());

        CUDA_CHECK(cudaEventRecord(start));

        for (int iteration = 0;
             iteration < iterations;
             ++iteration) {
            gather_with_stride<<<blocks, threads>>>(
                d_input, d_output,
                n, mask, stride);
        }

        CUDA_CHECK(cudaGetLastError());
        CUDA_CHECK(cudaEventRecord(stop));
        CUDA_CHECK(cudaEventSynchronize(stop));

        float total_ms = 0.0f;
        CUDA_CHECK(cudaEventElapsedTime(
            &total_ms, start, stop));

        const double average_ms =
            static_cast<double>(total_ms) /
            iterations;
        const double useful_bytes =
            static_cast<double>(n) *
            2.0 * sizeof(float);
        const double useful_bandwidth_gb_s =
            useful_bytes /
            (average_ms * 1.0e6);

        CUDA_CHECK(cudaMemcpy(
            h_output.data(), d_output, bytes,
            cudaMemcpyDeviceToHost));

        size_t mismatches = 0;

        for (size_t i = 0; i < n; ++i) {
            const size_t source_index =
                (i * stride) & mask;

            if (h_output[i] !=
                h_input[source_index]) {
                ++mismatches;
            }
        }

        all_passed =
            all_passed && mismatches == 0;

        std::printf(
            "%zu,%zu,%.6f,%.2f,%zu\n",
            stride,
            stride * sizeof(float),
            average_ms,
            useful_bandwidth_gb_s,
            mismatches);
    }

    CUDA_CHECK(cudaEventDestroy(stop));
    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaFree(d_output));
    CUDA_CHECK(cudaFree(d_input));

    std::printf(
        "\nVerification: %s\n",
        all_passed ? "PASSED" : "FAILED");

    return all_passed
        ? EXIT_SUCCESS
        : EXIT_FAILURE;
}

13.1 为什么参数上限必须存在

程序允许从命令行传入 N 和迭代次数，但没有接受任意大值。原因很实际：

• std::vector 需要 Host 内存；
• cudaMalloc 需要 Device 内存；
• N 参与字节数和 Grid 大小计算；
• 迭代次数过大会让程序长时间占用 GPU；
• 一个错误的输入不应该把机器拖入内存交换或显示驱动超时。

教程代码也应该有输入边界。能运行不等于可以忽略资源管理。

13.2 为什么索引使用 size_t

数组长度和字节数使用 size_t，避免在较大 N 下由 32 位有符号整数产生溢出。全局索引先把 blockIdx.x 转为 size_t：

const size_t i =
    static_cast<size_t>(blockIdx.x) *
    blockDim.x +
    threadIdx.x;

如果先用 32 位整数完成乘法，再转换为 64 位，溢出已经发生，事后转换无法补救。这个细节在小数组中看不出来，在大规模图像、体数据或稀疏索引中很重要。

13.3 为什么不在计时循环中验证

Device 到 Host 的复制和 CPU 验证放在计时循环之后。否则测到的是：

Kernel + 同步 + D2H + CPU 循环

这不能回答访问模式怎样影响 Kernel。另一方面，完全不验证也不行，因为编写错误的索引可能碰巧跑得更快。正确做法是把“性能范围”和“正确性检查”分开执行，但两者都保留。

---

十四、Stride 实测结果逐行读懂

本机参数：

GPU                 : NVIDIA GeForce RTX 2060
N                   : 4,194,304 floats
Input size          : 16 MiB
Output size         : 16 MiB
Iterations          : 100
Warmup per stride   : 10
Threads per block   : 256
Build               : Release

程序输出：

stride,byte_gap_per_lane,average_ms,useful_bandwidth_GBps,mismatches
1,4,0.187500,178.96,0
3,12,0.331858,101.11,0
5,20,0.487751,68.79,0
9,36,0.806673,41.60,0
17,68,1.316442,25.49,0
33,132,1.323668,25.35,0
65,260,1.335706,25.12,0

Verification: PASSED

14.1 Stride 为 1

相邻 Lane 地址间隔 4 字节，输入读取连续；输出也连续。逻辑有效带宽约 179 GB/s，是本组基线。

14.2 Stride 为 3

相邻 Lane 相隔 12 字节。Warp 覆盖范围扩大，其中夹着不属于当前请求的元素。有效带宽降到约 101 GB/s。此时并非所有 Sector 都只使用一个 float，但利用率已经明显低于连续访问。

14.3 Stride 为 5

地址间隔变为 20 字节，覆盖区域继续扩大，有效带宽约 69 GB/s。算术不变，下降主要与内存请求分散相关。

14.4 Stride 为 9

相邻地址相隔 36 字节，已经超过一个 32 字节 Sector。多数相邻 Lane 会落入不同 Sector，逻辑有效带宽下降到约 42 GB/s。

14.5 Stride 为 17、33、65

间隔继续扩大，但时间没有无限增长，而是在约 1.32 到 1.34 毫秒附近形成平台，有效带宽约 25 GB/s。

平台并不表示更大 Stride 被自动优化回来了。更合理的解释是：当访问已经足够分散、每个 Lane 基本需要独立 Sector 后，再继续扩大地址间隔不会让单 Warp 的 Sector 数超过活跃 Lane 数。系统已经进入另一种瓶颈区间，延迟隐藏、并发请求数、缓存与显存吞吐共同决定上限。

14.6 性能倍率

Stride 65 相对 Stride 1 的时间倍率：

$$R_{time}=\frac{1.335706}{0.187500}\approx 7.12$$

对应的逻辑有效带宽比例：

$$R_{bw}=\frac{178.96}{25.12}\approx 7.12$$

因为每组的逻辑字节完全相同，时间倍率与有效带宽倍率互为对应。这里才能公平地说，相同有用数据量在分散访问下处理得慢约七倍。

14.7 不要把 7.12 倍推广成固定常数

换成其他 GPU、数据规模、缓存配置、功耗状态或 Stride 集合，倍率会变化。教学结论不是“Stride 65 必慢 7.12 倍”，而是：

• 地址分散会增加覆盖 Sector；
• 增加的流量与请求降低有效利用率；
• 变化可由控制变量实验观察；
• 绝对倍率必须在目标设备和真实规模上测量。

如果在你的 GPU 上连续访问只有 90 GB/s，或者平台出现在 35 GB/s，不代表程序错误。先确认 Release、N、迭代次数、校验结果和后台负载，再比较趋势。

---

十五、实验三：AoS 与 SoA 不是编码风格之争

粒子、点云、物体、传感器记录和数据库行经常使用结构体：

struct Particle {
    float x;
    float y;
    float z;
    float w;
};

一种布局叫 Array of Structures，简称 AoS：

Particle particles[N];

内存排列：

x0 y0 z0 w0 | x1 y1 z1 w1 | x2 y2 z2 w2 | ...

另一种布局叫 Structure of Arrays，简称 SoA：

float x[N];
float y[N];
float z[N];
float w[N];

内存排列：

x0 x1 x2 x3 ... | y0 y1 y2 y3 ... | z0 z1 ... | w0 w1 ...

如果 Kernel 同时使用同一个粒子的 x、y、z、w，AoS 可能很自然；如果 Kernel 当前只使用所有粒子的 x，SoA 可以让一个 Warp 的 x 请求连续排列。

本篇结构体明确对齐为 16 字节：

struct __align__(16) Particle {
    float x;
    float y;
    float z;
    float w;
};

AoS 中，相邻粒子的 x 相隔：

$$\Delta addres{s}_{AoS}=sizeof(Particle)=16$$

一个 Warp 的 32 个 x 位于大约 512 字节结构体足迹中：

$$B_{footprint,AoS}=32\times 16=512$$

SoA 中，相邻 x 相隔 4 字节：

$$\Delta addres{s}_{SoA}=sizeof(float)=4$$

一个 Warp 的 32 个 x 连续占用：

$$B_{footprint,SoA}=32\times 4=128$$

注意“足迹”不等于编译器必定加载全部 y、z、w。AoS Kernel 源代码只读取 x，生成的加载指令通常也只需要 x 的标量。但这些 x 分散在更宽的地址区域，所覆盖 Sector 中的其他字节对当前 Kernel 没有贡献。

15.1 业务需求决定布局

不要把本节总结为“SoA 永远比 AoS 快”。正确问题是：

当前热 Kernel 在同一时刻需要跨很多对象读取同一个字段，还是需要对单个对象读取大多数字段？

典型倾向如下：

场景	可能更合适的起点	原因
更新所有粒子的 x 坐标	SoA	相邻 Lane 读取连续 x
只统计所有对象的类别	SoA	避免为其他字段扩大地址足迹
每线程完整处理一个小结构体	AoS 或 AoSoA	同一线程字段局部性可能更好
SIMD/SIMT 分批处理固定宽度对象	AoSoA	在字段连续和对象分组间折中
CPU 与 GPU 共享序列化格式	需要实测	转换成本可能抵消 Kernel 收益
多个 Kernel 各自消费不同字段	SoA	每个阶段只读取需要的数组

AoSoA 是 Array of Structures of Arrays。它把对象按固定小块分组，每块内部按字段连续存放。大型粒子系统、向量化 CPU 与 GPU 共用数据时，经常用它在多种访问需求之间折中。本篇先掌握 AoS 与 SoA，后续再讨论更复杂布局。

---

十六、AoS/SoA 基准怎样保证数学一致

两个 Kernel 都执行：

output[i] = x * scale + bias;

AoS 版本：

__global__ void transform_aos(const Particle* particles,
                              float* output,
                              size_t n,
                              float scale,
                              float bias) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;

    if (i < n) {
        output[i] =
            particles[i].x * scale + bias;
    }
}

SoA 版本：

__global__ void transform_soa(const float* x,
                              float* output,
                              size_t n,
                              float scale,
                              float bias) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;

    if (i < n) {
        output[i] =
            x[i] * scale + bias;
    }
}

两个输入中的 x 值完全相同，Scale 与 Bias 相同，Block 数和线程数相同，输出都是连续数组。区别只有 x 在内存中的布局。

为了避免“两个版本都以同一种方式写错”，验证包含三层：

1. AoS 输出与 CPU 参考比较；
2. SoA 输出与 CPU 参考比较；
3. AoS 输出与 SoA 输出逐元素比较。

允许误差：

$$ϵ=10^{-6}$$

因为 CPU 与 GPU 都执行同样的简单乘加，最大绝对误差很小。但程序仍然不使用“只看前十个”的抽样验证，而是检查全部 N 个元素。

16.1 为什么有效带宽仍按 8N 字节计算

从算法语义看，每个元素只消费一个 x，写出一个结果：

$$B_{logical}=N\times (4+4)$$

AoS 数组整体占 16N 字节，但 Kernel 并没有在语义上消费 y、z、w。因此用 8N 字节计算“逻辑有效带宽”，可以表达同样业务工作在两种布局中完成得多快。

如果把 AoS 的输入字节直接写成 16N，再与 SoA 的 4N 比较，就会人为给 AoS 增加“完成了更多有用工作”的分子，产生误导。物理流量应该由 Profiler 测，不应该用结构体总大小猜。

---

十七、AoS/SoA 完整源码

完整文件位于 blogs/code/04/aos_soa_benchmark.cu。

#include <cuda_runtime.h>

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>

void cuda_check_impl(cudaError_t error,
                     const char* expression,
                     const char* file,
                     int line) {
    if (error != cudaSuccess) {
        std::fprintf(stderr,
                     "CUDA error at %s:%d\n"
                     "  expression: %s\n"
                     "  reason: %s\n",
                     file, line, expression,
                     cudaGetErrorString(error));
        std::exit(EXIT_FAILURE);
    }
}

#define CUDA_CHECK(call) \
    cuda_check_impl((call), #call, __FILE__, __LINE__)

struct __align__(16) Particle {
    float x;
    float y;
    float z;
    float w;
};

__global__ void transform_aos(
    const Particle* particles,
    float* output,
    size_t n,
    float scale,
    float bias) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) *
        blockDim.x +
        threadIdx.x;

    if (i < n) {
        output[i] =
            particles[i].x * scale + bias;
    }
}

__global__ void transform_soa(
    const float* x,
    float* output,
    size_t n,
    float scale,
    float bias) {
    const size_t i =
        static_cast<size_t>(blockIdx.x) *
        blockDim.x +
        threadIdx.x;

    if (i < n) {
        output[i] =
            x[i] * scale + bias;
    }
}

int main(int argc, char** argv) {
    size_t n = 1ULL << 23;
    int iterations = 100;

    if (argc >= 2) {
        const unsigned long long parsed =
            std::strtoull(argv[1], nullptr, 10);

        if (parsed == 0 ||
            parsed > 30000000ULL) {
            std::fprintf(
                stderr,
                "N must be in [1, 30000000].\n");
            return EXIT_FAILURE;
        }

        n = static_cast<size_t>(parsed);
    }

    if (argc >= 3) {
        const long parsed =
            std::strtol(argv[2], nullptr, 10);

        if (parsed <= 0 || parsed > 10000) {
            std::fprintf(
                stderr,
                "iterations must be in [1, 10000].\n");
            return EXIT_FAILURE;
        }

        iterations = static_cast<int>(parsed);
    }

    constexpr float scale = 1.25f;
    constexpr float bias = -0.5f;
    constexpr int threads = 256;

    const int blocks =
        static_cast<int>(
            (n + threads - 1) / threads);
    const size_t aos_bytes =
        n * sizeof(Particle);
    const size_t field_bytes =
        n * sizeof(float);

    std::vector<Particle> h_particles(n);
    std::vector<float> h_x(n);
    std::vector<float> h_aos_output(n);
    std::vector<float> h_soa_output(n);

    for (size_t i = 0; i < n; ++i) {
        const float x =
            static_cast<float>(
                (i * 29 + 7) % 32749) *
            0.0001f;

        h_particles[i] = Particle{
            x,
            static_cast<float>(i % 101),
            static_cast<float>(i % 211),
            1.0f
        };

        h_x[i] = x;
    }

    Particle* d_particles = nullptr;
    float* d_x = nullptr;
    float* d_aos_output = nullptr;
    float* d_soa_output = nullptr;

    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_particles),
        aos_bytes));
    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_x),
        field_bytes));
    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_aos_output),
        field_bytes));
    CUDA_CHECK(cudaMalloc(
        reinterpret_cast<void**>(&d_soa_output),
        field_bytes));

    CUDA_CHECK(cudaMemcpy(
        d_particles,
        h_particles.data(),
        aos_bytes,
        cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(
        d_x,
        h_x.data(),
        field_bytes,
        cudaMemcpyHostToDevice));

    cudaEvent_t start{};
    cudaEvent_t stop{};
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));

    for (int warmup = 0;
         warmup < 10;
         ++warmup) {
        transform_aos<<<blocks, threads>>>(
            d_particles,
            d_aos_output,
            n, scale, bias);
        transform_soa<<<blocks, threads>>>(
            d_x,
            d_soa_output,
            n, scale, bias);
    }

    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    CUDA_CHECK(cudaEventRecord(start));

    for (int iteration = 0;
         iteration < iterations;
         ++iteration) {
        transform_aos<<<blocks, threads>>>(
            d_particles,
            d_aos_output,
            n, scale, bias);
    }

    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));

    float aos_total_ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(
        &aos_total_ms, start, stop));

    CUDA_CHECK(cudaEventRecord(start));

    for (int iteration = 0;
         iteration < iterations;
         ++iteration) {
        transform_soa<<<blocks, threads>>>(
            d_x,
            d_soa_output,
            n, scale, bias);
    }

    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));

    float soa_total_ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(
        &soa_total_ms, start, stop));

    CUDA_CHECK(cudaMemcpy(
        h_aos_output.data(),
        d_aos_output,
        field_bytes,
        cudaMemcpyDeviceToHost));
    CUDA_CHECK(cudaMemcpy(
        h_soa_output.data(),
        d_soa_output,
        field_bytes,
        cudaMemcpyDeviceToHost));

    size_t aos_mismatches = 0;
    size_t soa_mismatches = 0;
    size_t cross_mismatches = 0;
    float max_abs_error = 0.0f;

    for (size_t i = 0; i < n; ++i) {
        const float reference =
            h_x[i] * scale + bias;
        const float aos_error =
            std::fabs(
                h_aos_output[i] -
                reference);
        const float soa_error =
            std::fabs(
                h_soa_output[i] -
                reference);

        max_abs_error =
            std::max(
                max_abs_error,
                std::max(aos_error, soa_error));

        if (aos_error > 1.0e-6f) {
            ++aos_mismatches;
        }
        if (soa_error > 1.0e-6f) {
            ++soa_mismatches;
        }
        if (h_aos_output[i] !=
            h_soa_output[i]) {
            ++cross_mismatches;
        }
    }

    const double aos_average_ms =
        static_cast<double>(
            aos_total_ms) /
        iterations;
    const double soa_average_ms =
        static_cast<double>(
            soa_total_ms) /
        iterations;
    const double logical_useful_bytes =
        static_cast<double>(n) *
        2.0 * sizeof(float);
    const double aos_useful_bandwidth =
        logical_useful_bytes /
        (aos_average_ms * 1.0e6);
    const double soa_useful_bandwidth =
        logical_useful_bytes /
        (soa_average_ms * 1.0e6);

    int device = 0;
    cudaDeviceProp prop{};
    CUDA_CHECK(cudaGetDevice(&device));
    CUDA_CHECK(cudaGetDeviceProperties(
        &prop, device));

    std::printf("GPU: %s\n", prop.name);
    std::printf(
        "layout,average_ms,"
        "logical_useful_bandwidth_GBps,"
        "mismatches\n");
    std::printf(
        "AoS,%.6f,%.2f,%zu\n",
        aos_average_ms,
        aos_useful_bandwidth,
        aos_mismatches);
    std::printf(
        "SoA,%.6f,%.2f,%zu\n",
        soa_average_ms,
        soa_useful_bandwidth,
        soa_mismatches);
    std::printf(
        "\nSoA speedup over AoS: %.2fx\n",
        aos_average_ms /
        soa_average_ms);
    std::printf(
        "Cross-layout mismatches: %zu\n",
        cross_mismatches);
    std::printf(
        "Maximum absolute error: %.9e\n",
        max_abs_error);

    CUDA_CHECK(cudaEventDestroy(stop));
    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaFree(d_soa_output));
    CUDA_CHECK(cudaFree(d_aos_output));
    CUDA_CHECK(cudaFree(d_x));
    CUDA_CHECK(cudaFree(d_particles));

    const bool passed =
        aos_mismatches == 0 &&
        soa_mismatches == 0 &&
        cross_mismatches == 0;

    std::printf(
        "\nVerification: %s\n",
        passed ? "PASSED" : "FAILED");

    return passed
        ? EXIT_SUCCESS
        : EXIT_FAILURE;
}

---

十八、AoS/SoA 本机结果

测试参数：

GPU                 : NVIDIA GeForce RTX 2060
N                   : 8,388,608 particles
Particle size       : 16 bytes
Iterations          : 100
Expression          : output[i] = x * 1.25 - 0.5
Build               : Release

实测输出：

layout,average_ms,logical_useful_bandwidth_GBps,mismatches
AoS,0.819582,81.88,0
SoA,0.319387,210.12,0

SoA speedup over AoS: 2.57x
Cross-layout mismatches: 0
Maximum absolute error: 2.384185791e-07

Verification: PASSED

速度比：

$$S_{SoA/AoS}=\frac{0.819582}{0.319387}\approx 2.57$$

为什么不是理论足迹比例的四倍？

第一，AoS 的 512 字节足迹不等于每个 Warp 一定从 DRAM 读取完整 512 字节且毫无复用。

第二，两个 Kernel 都包含连续输出写、索引计算、边界判断和乘加，布局只影响其中的输入读取部分。

第三，GPU 可以并发处理多个 Warp，通过延迟隐藏减轻单次读取等待。

第四，缓存层级与内存请求调度会合并、复用或重排部分流量。

第五，时间还受到频率与测量噪声影响。

因此，2.57 倍是这台机器、这组规模和这个只读取 x 的 Kernel 的实际结果，不是所有 AoS/SoA 场景的固定比例。

不过，趋势已经足够清楚：当业务只消费一个字段时，SoA 让有用字段连续排列，能够显著提高逻辑有效带宽。

18.1 什么时候改布局可能不值得

如果上游数据以 AoS 到达，GPU 只执行一次很短的 Kernel，然后立刻把数据发走，为转换成 SoA 额外做一次全量重排可能得不偿失。需要比较：

$$t_{saved}>t_{convert}$$

如果数据在 GPU 上会被十几个 Kernel 反复使用，转换成本可以被多次收益摊薄：

$$R\times (t_{AoS}-t_{SoA})>t_{convert}$$

其中 R 是受益 Kernel 的执行次数。工程优化不能只看单个 Kernel 的最快布局，还要看完整流水线。

18.2 布局也是接口设计

AoS 改 SoA 会影响：

• 内存分配数量；
• Host 到 Device 复制方式；
• API 参数；
• 序列化格式；
• CPU 侧处理；
• 多个 Kernel 之间的数据契约；
• 对齐与内存池；
• 调试和可读性。

因此在大型项目里，最好的优化方案可能不是一次性推翻结构，而是先为热点字段建立连续镜像、只在阶段边界转换，或使用 AoSoA 分块。性能收益必须与维护成本一起评估。

---

十九、完整构建与运行指南

本篇目录结构：

blogs/
└── code/
    └── 04/
        ├── CMakeLists.txt
        ├── offset_access_benchmark.cu
        ├── stride_access_benchmark.cu
        └── aos_soa_benchmark.cu

CMake 文件：

cmake_minimum_required(VERSION 3.18)

project(cuda_blog_04 LANGUAGES CXX CUDA)

set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CUDA_STANDARD 14)
set(CMAKE_CUDA_STANDARD_REQUIRED ON)

add_executable(
    offset_access_benchmark
    offset_access_benchmark.cu)
add_executable(
    stride_access_benchmark
    stride_access_benchmark.cu)
add_executable(
    aos_soa_benchmark
    aos_soa_benchmark.cu)

foreach(target IN ITEMS
        offset_access_benchmark
        stride_access_benchmark
        aos_soa_benchmark)
    if(MSVC)
        target_compile_options(
            ${target}
            PRIVATE
            $<$<COMPILE_LANGUAGE:CUDA>:-Xcompiler=/utf-8>
            $<$<COMPILE_LANGUAGE:CXX>:/utf-8>)
    endif()
endforeach()

19.1 Windows + Visual Studio 2022

在 blogs/code/04 目录打开 PowerShell：

cmake -S . -B build `
  -G "Visual Studio 17 2022" `
  -A x64 `
  -DCMAKE_CUDA_ARCHITECTURES=75

cmake --build build --config Release

75 对应本文实测 RTX 2060 的计算能力。你的设备如果不同，应改成自己的目标架构。例如不要因为复制教程命令就把所有 GPU 都编译为 75。可以先运行第二篇的设备查询程序，或查询设备 Compute Capability。

运行：

.\build\Release\offset_access_benchmark.exe
.\build\Release\stride_access_benchmark.exe
.\build\Release\aos_soa_benchmark.exe

自定义 N 和迭代次数：

.\build\Release\offset_access_benchmark.exe 16777216 200
.\build\Release\stride_access_benchmark.exe 4194304 200
.\build\Release\aos_soa_benchmark.exe 8388608 200

Stride 程序要求 N 为 2 的幂。下面是合法示例：

1048576
2097152
4194304
8388608

如果传入 3,000,000，程序会拒绝，而不是悄悄破坏全排列实验。

19.2 Linux

cmake -S . -B build \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_CUDA_ARCHITECTURES=75

cmake --build build -j

./build/offset_access_benchmark
./build/stride_access_benchmark
./build/aos_soa_benchmark

多配置生成器和单配置生成器对 Release 的指定方式不同。Visual Studio 在构建阶段使用 --config Release，常见 Linux Makefile 或 Ninja 在配置阶段使用 -DCMAKE_BUILD_TYPE=Release。

19.3 不使用 CMake

Windows 开发者命令行：

nvcc -O3 -std=c++14 `
  -arch=sm_75 `
  offset_access_benchmark.cu `
  -o offset_access_benchmark.exe

nvcc -O3 -std=c++14 `
  -arch=sm_75 `
  stride_access_benchmark.cu `
  -o stride_access_benchmark.exe

nvcc -O3 -std=c++14 `
  -arch=sm_75 `
  aos_soa_benchmark.cu `
  -o aos_soa_benchmark.exe

Linux：

nvcc -O3 -std=c++14 -arch=sm_75 \
  offset_access_benchmark.cu \
  -o offset_access_benchmark

nvcc -O3 -std=c++14 -arch=sm_75 \
  stride_access_benchmark.cu \
  -o stride_access_benchmark

nvcc -O3 -std=c++14 -arch=sm_75 \
  aos_soa_benchmark.cu \
  -o aos_soa_benchmark

19.4 成功验收标准

不要只看程序打印了表格。每个可用结果至少满足：

• 进程退出码为 0；
• 最后一行是 Verification: PASSED；
• 每一组 mismatches 都为 0；
• N 和 Iterations 与报告一致；
• 使用 Release；
• GPU 名称正确；
• 没有 CUDA API 错误；
• 测量期间没有同时运行其他重负载；
• 连续运行多次趋势大致一致。

如果 mismatches 不为 0，先修正确性。不要把错误版本的高带宽当成优化成果。

---

二十、有效带宽到底“有效”在哪里

本篇三个程序都输出有效带宽。这个指标非常有用，也非常容易被滥用。

20.1 算法有用字节

复制 Kernel 每个元素逻辑上读 4 字节、写 4 字节：

$$B_{logical}=8N$$

向量加法读两个输入、写一个输出：

$$B_{logical,add}=12N$$

SAXPY 读取 X、读取并写回 Y，按实现可能按 12N 字节作为基础口径。关键不是背某个数字，而是明确列出每个数组的逻辑读写。

20.2 有效带宽

$$B{W}_{effective}=\frac{B_{logical}}{t}$$

它回答：

每秒完成了多少“对算法有意义的数据工作”？

当两个版本完成同一任务时，这个指标很适合比较。访问浪费越多，完成相同逻辑字节所需时间通常越长，有效带宽就越低。

20.3 它不是显存芯片实际流量

假设 Stride 访问让每个 4 字节请求覆盖一个 32 字节 Sector，硬件某一级处理的字节可能远大于 8N。缓存命中又可能让部分请求不去 DRAM。只用源码无法得到精确物理流量。

要研究物理层，应查看 Nsight Compute 对应架构支持的：

• Global Load Requests；
• Global Load Sectors；
• L1/TEX 请求和 Sector；
• L2 请求和 Sector；
• DRAM Read Bytes；
• DRAM Write Bytes；
• 各层吞吐率；
• Excessive Sectors 或类似派生指标。

指标名称会随版本变化。应在当前工具中阅读指标说明，不要从旧博客复制一串名字后默认含义永远不变。

20.4 为什么带宽可能超过你以为的数

如果数据大部分命中缓存，有效带宽使用逻辑字节计算时可能很高，甚至高于设备标称 DRAM 带宽。这不违反物理规律，因为分母中的数据没有全部经过 DRAM。它说明当前工作集从更靠近 SM 的层级得到服务。

反过来，有效带宽很低也不一定只因合并差。可能原因包括：

• N 太小，启动开销占比高；
• 分支严重；
• 指令依赖长；
• Occupancy 不足；
• 寄存器溢出到 Local Memory；
• 原子操作串行化；
• 数据传输被误计入；
• GPU 降频；
• Debug 构建；
• 桌面图形或其他进程争用 GPU。

所以有效带宽是起点，不是诊断终点。

---

二十一、用 Nsight Compute 把推理变成硬件证据

CUDA Event 已经告诉我们 Stride 版本变慢，但它没有直接显示“一个请求用了多少 Sector”。这一步交给 Nsight Compute。

Nsight Compute 是面向单个 CUDA Kernel 的性能分析器。它有图形界面 ncu-ui，也有命令行工具 ncu。它与 Nsight Systems 的侧重点不同：

• Nsight Systems 擅长查看整个应用时间线、CPU/GPU 并发、传输、同步与空闲；
• Nsight Compute 擅长深入某个 Kernel，查看内存、指令、Warp、Occupancy 和管线利用率。

实际工作常用顺序是：

先用 Nsight Systems 找到耗时热点
    ↓
再用 Nsight Compute 分析那个 Kernel
    ↓
提出修改并重新测量

本篇程序很小，已经知道热点就是 gather_with_stride，可以直接用 Nsight Compute。

21.1 先确认工具存在

Windows PowerShell：

ncu --version
ncu-ui

Linux：

ncu --version
ncu-ui

如果命令找不到，检查 Nsight Compute 是否随 CUDA Toolkit 安装，以及它的目录是否加入 PATH。Windows 常见安装位置在 NVIDIA Corporation 的 Nsight Compute 目录中，具体路径取决于版本。

21.2 为什么只分析一个 Kernel 启动

Profiler 采集完整指标时，可能对同一个 Kernel 进行多次 Replay。若一次分析几百个启动，运行时间会非常长，报告也难以阅读。正确方法是使用 Kernel 过滤、跳过数量和采集数量，只选一个代表性启动。

本篇 Stride 程序在迭代次数为 1 时，每个 Stride 包含：

10 次预热 + 1 次计时 = 11 次 Kernel 启动

要分析 Stride 1 的计时 Kernel，可以跳过前十个预热：

ncu --set full `
  --launch-skip 10 `
  --launch-count 1 `
  .\build\Release\stride_access_benchmark.exe `
  4194304 1

要分析 Stride 65，前面有六组 Stride，每组十一条启动，再加当前组十次预热：

$$N_{skip}=6\times 11+10=76$$

命令：

ncu --set full `
  --launch-skip 76 `
  --launch-count 1 `
  .\build\Release\stride_access_benchmark.exe `
  4194304 1

Linux 只需改可执行文件路径和换行符：

ncu --set full \
  --launch-skip 10 \
  --launch-count 1 \
  ./build/stride_access_benchmark \
  4194304 1

如果你的代码修改了预热次数、Stride 列表或迭代结构，跳过数量必须重新计算。最稳妥的方法是在 Kernel 名称中区分版本、使用 NVTX 标记范围，或让程序接受单独 Stride 参数。教学程序保持接口简单，所以这里演示手动计算。

21.3 生成报告文件

ncu --set full `
  --launch-skip 10 `
  --launch-count 1 `
  --export stride1_report `
  .\build\Release\stride_access_benchmark.exe `
  4194304 1

再为 Stride 65 生成一份报告：

ncu --set full `
  --launch-skip 76 `
  --launch-count 1 `
  --export stride65_report `
  .\build\Release\stride_access_benchmark.exe `
  4194304 1

用 ncu-ui 打开两份报告，可以把 Stride 1 设为 Baseline，再比较 Stride 65。报告比较比手抄几个数字可靠，因为它能同时显示变化比例和其他关联指标。

21.4 优先看什么

进入 Memory Workload Analysis 或对应版本的内存分析页面，重点寻找：

1. Global Load 的 Request；
2. Global Load 对应的 Sector；
3. 每个 Request 平均使用多少 Sector；
4. 理想 Sector 与过量 Sector；
5. L1/TEX 命中情况；
6. L2 命中情况；
7. DRAM Read/Write Bytes；
8. 内存吞吐占峰值比例；
9. 是否有 Local Memory 流量；
10. Source 页面中是哪一行产生相关访问。

指标名称可能包含很长的硬件前缀。例如某些版本会显示 L1/TEX Global Load Requests 和 Global Load Sectors，某些版本还提供 Ideal 或 Excessive Sectors。不要因为名字长就只看最后一个百分比。先点开指标说明，确认：

• 统计的是 Load 还是 Store；
• 统计层级是 L1/TEX、L2 还是 DRAM；
• 单位是 Count、Byte、Percent 还是每请求平均值；
• 是全部 Kernel 还是某一条指令；
• 是否受到 Replay 影响。

21.5 用 Sectors per Request 建立直觉

可以用一个派生量帮助比较：

$$R_{sector/request}=\frac{N_{sector}}{N_{request}}$$

它表示平均每个请求需要多少 Sector。对一个读取连续 float 的完整 Warp，理想地址覆盖可由四个 32 字节 Sector 完成；分散访问可能需要更多。

但不要把这个派生量与“每线程事务数”混为一谈。Request 的定义由所看硬件单元决定，某些层级会对请求重新组织。最可靠的判断方式是同时看工具给出的 Ideal、Actual 与 Excessive 指标，并阅读当前版本说明。

21.6 Profiler 下的时间不能代替正常运行时间

Nsight Compute 采集多组指标时可能重复执行 Kernel，插入同步并改变缓存状态。报告里的 Duration 可以用于同一次分析的参考，但发布性能结果时仍应使用无 Profiler 的 Release 程序和 CUDA Event。

正确分工是：

无 Profiler 的 Event 时间：报告用户会感受到的 Kernel 性能
Profiler 计数：解释为什么性能会这样

不要拿 ncu --set full 运行后的总进程时间当优化前后速度。

---

二十二、遇到 ERR_NVGPUCTRPERM 怎么办

在本文写作机器上，ncu 可以启动，但采集硬件计数时返回：

ERR_NVGPUCTRPERM
The user does not have permission to access
NVIDIA GPU Performance Counters

这不是 CUDA Kernel 错误，也不是 Nsight Compute 没安装，而是当前用户没有读取 GPU 性能计数器的权限。

Windows 上有两种常见解决思路。

第一种是以管理员身份启动 Nsight Compute，或从管理员 PowerShell/命令提示符运行完整 ncu 命令。

第二种是在当前 NVIDIA App 中进入：

System
  → Advanced
    → Developer
      → Manage GPU Performance Counters

然后按组织安全策略授予合适权限。修改这一设置需要管理员权限。

旧版驱动或控制界面中，也可能在 NVIDIA Control Panel 的 Developer 设置下管理 GPU Performance Counters。新版官方说明已经把 NVIDIA App 作为当前方式，并提示 Windows Control Panel 方式会逐步弃用。

Linux 上通常需要管理员能力或由系统管理员开放性能计数器权限。企业服务器、容器和共享集群的策略可能不同，不要为了跑教程擅自修改公共机器的安全设置。

本机不能采集计数并不影响三个 CUDA 示例的正确运行与 Event 计时。文章没有伪造 Profiler 数字，而是明确报告权限限制。你获得权限后，可以按照上一节命令补齐第三层证据。

这也是工程报告应有的态度：

没有测到的数据就标记为未测，不用理论估算冒充硬件计数。

---

二十三、缓存与合并是两个维度

初学者常问：“如果缓存命中，合并访问还重要吗？”答案是重要，但两者解决的问题不同。

合并关注同一 Warp 的地址如何组织成请求；缓存关注请求的数据是否已经位于更近的存储层级。

可以构造四种组合：

合并情况	缓存情况	典型表现
好	命中高	请求紧凑且大多由近端缓存服务，通常理想
好	命中低	请求紧凑，但需要更多访问 L2 或 DRAM
差	命中高	请求分散，但部分代价被复用掩盖
差	命中低	请求分散又缺少复用，通常最糟

Offset 实验可能属于“地址覆盖略差，但相邻 Warp 有缓存复用”，因此 Event 时间未严格按四 Sector和五 Sector的理论比例变化。

Stride 大且形成全排列时，一个 Warp 内请求分散，跨 Warp 的近期空间局部性也被打乱，更容易暴露带宽下降。

23.1 不要用小数组证明 DRAM 性能

如果输入只有几十 KB，它可能完全驻留在缓存里。你测到的主要是缓存带宽与 Kernel 启动，而不是 Global Memory 到 DRAM 的行为。

本篇 Stride 默认输入和输出各 16 MiB，总工作集超过 RTX 2060 的常见 L2 容量，目的是降低“整个数组常驻 L2”的可能性。但即使如此，单次运行中的局部复用仍存在。

调整 N 时可以观察三个区间：

1. 很小：启动和缓存占主导；
2. 中等：跨越缓存容量，变化明显；
3. 很大：DRAM 与系统稳定性占主导。

把 N 从 2 的 16 次方逐步增加到 2 的 24 次方，记录每个 Stride 的曲线，比只报告一个 N 更有分析价值。

23.2 冷缓存和热缓存回答不同问题

重复执行同一 Kernel 会使一部分数据或元数据处于热状态。如果你的真实业务会反复处理同一数据，热缓存结果具有意义；如果每次处理全新数据，应该设计轮换缓冲区或扩大工作集。

不要把“清空缓存”当成必须遵守的仪式。先写清楚要模拟哪种业务：

• 单次冷启动；
• 同一帧多次处理；
• 多帧轮换；
• 多 Kernel 共享数据；
• 数据流只读一次。

基准环境要代表目标场景。

---

二十四、二维图像：Thread x 应该对应连续列

假设一张行优先存储的灰度图像，宽为 W，高为 H：

pixel = image[y * width + x];

最自然的二维线程映射：

const int x =
    blockIdx.x * blockDim.x + threadIdx.x;
const int y =
    blockIdx.y * blockDim.y + threadIdx.y;

在 CUDA 的线性线程组织中，threadIdx.x 是最快变化维度。一个 Warp 通常先沿 x 方向覆盖线程。因此相邻 Lane 倾向于拥有相邻 x，并读取相邻像素。

如果错误地交换映射：

const int y =
    blockIdx.x * blockDim.x + threadIdx.x;
const int x =
    blockIdx.y * blockDim.y + threadIdx.y;

相邻 Lane 可能沿 y 变化。行优先数组中，相邻行同一列的地址间隔约为：

$$\Delta address=W\times sizeof(T)$$

对宽 1920 的 float 图像：

$$\Delta address=1920\times 4=7680$$

这会让一个 Warp 的请求极度分散。

24.1 行宽不是 Pitch 时要小心

使用 cudaMallocPitch 时，每行起点之间的距离是 Pitch 字节，不一定等于 width * sizeof(T)：

char* row =
    reinterpret_cast<char*>(base) +
    y * pitch;
float* row_values =
    reinterpret_cast<float*>(row);
float value = row_values[x];

常见错误是把 Pitch 当元素数再次乘 sizeof(float)，或者把元素宽度当字节 Pitch。两者都会产生错误地址。

Pitch 的存在通常帮助每行满足合适对齐，但它不能自动修复错误的线程映射。仍要确保相邻 Lane 沿连续列访问。

24.2 图像邻域访问

卷积或滤波中，每个线程读取周围多个像素。某一条读指令，例如所有 Lane 同时读取自己的中心像素，通常仍可合并；读取左邻居时地址整体偏移一个元素，可能出现类似 Offset 的边界 Sector；读取上邻居时整个 Warp 访问上一行的连续区域，也可能保持合并。

问题不是“用了邻域就一定不合并”，而是逐条列出：

center : image[y * pitch + x]
left   : image[y * pitch + x - 1]
right  : image[y * pitch + x + 1]
up     : image[(y - 1) * pitch + x]
down   : image[(y + 1) * pitch + x]

只要同一条指令中相邻 Lane 的 x 仍相邻，多数访问仍有良好的空间聚集。重复读取邻域会带来冗余，这时共享内存分块才可能进一步减少 Global Memory 流量。

---

二十五、矩阵：行访问和列访问为什么差别巨大

C/C++ 常见二维矩阵按行优先排列：

index = row * leading_dimension + column;

如果相邻 Lane 的 column 递增：

$$index(lane)=row\times L+colum{n}_0+lane$$

地址连续。

如果相邻 Lane 的 row 递增：

$$index(lane)=(ro{w}_0+lane)\times L+column$$

相邻地址间隔：

$$\Delta address=L\times sizeof(T)$$

当 L 很大，这与 Stride 实验完全同构。

矩阵转置是最典型例子。朴素转置很难让输入读取和输出写入同时连续：

output[x * height + y] =
    input[y * width + x];

相邻线程读取输入行时连续，但写入输出列时跨步；若交换线程映射，写入连续但读取跨步。解决思路不是让一个线程做更多随机访问，而是使用共享内存 Tile：

第一阶段：按行从 Global Memory 合并读取到 Shared Memory
同步
第二阶段：交换 Tile 内索引
按行向 Global Memory 合并写出

这将是下一篇的核心实战。你会看到共享内存不仅是“更快的小内存”，还是地址重排的中转站。

25.1 Leading Dimension 不能忽略

真实矩阵可能有 Padding，行跨度称为 Leading Dimension。即使列数是 1000，行跨度也可能是 1024。地址审计要使用真实跨度：

matrix[row * ld + column]

不要在 Kernel 里把 ld 写死为逻辑列数。错误既可能导致结果错误，也可能改变对齐和访问效率。

---

二十六、结构体对齐、填充与向量类型

26.1 sizeof 才是实际步长

看见：

struct Record {
    char flag;
    float value;
};

不能直接认为它占 5 字节。编译器通常会为了 float 对齐插入填充，sizeof(Record) 可能为 8。相邻 value 的地址间隔由 sizeof(Record) 决定。

Host 与 Device 编译布局通常保持一致，但跨编译器、跨语言或手工序列化时必须明确验证：

static_assert(sizeof(Particle) == 16,
              "Unexpected Particle layout");

教程中的 Particle 使用四个 float 并显式 16 字节对齐，使实验布局清晰。

26.2 float4 不会自动解决所有问题

把四个 float 改成 float4 可以表达 16 字节向量，并可能帮助编译器生成合适的向量化加载，但前提包括：

• 指针满足所需对齐；
• 每线程确实需要四个分量；
• 数组长度和尾部处理正确；
• 数据布局与访问模式匹配。

如果每线程只需要 x，读取 float4 反而明确搬入 y、z、w，增加无用数据。不要为了“向量类型看起来更 GPU”就替换标量。

26.3 reinterpret_cast 的对齐风险

下面的代码只有在地址满足 float4 对齐要求时才安全：

const float4* values4 =
    reinterpret_cast<const float4*>(values);

若 values 实际是基址加一个 float 的偏移，地址可能不再按 16 字节对齐。错误的向量化转换可能导致错误结果、未定义行为或更差性能。

优化前先检查：

基址对齐
元素偏移
结构体大小
编译器生成代码
尾部元素

26.4 cudaMalloc 对齐不代表所有子指针对齐

cudaMalloc 返回的基址具有较强对齐保证，但：

float* shifted = d_input + 1;

会把地址移动 4 字节。基址对齐不能自动传递给任意偏移后的子指针。Offset 实验正是在利用这个事实。

---

二十七、十个常见错误与纠正方法

错误一：只看线程编号，不算地址

错误判断：“线程是连续的，所以访问合并。”

纠正：把每条 Global Memory 访问写成 Lane 的函数，计算相邻 Lane 的字节间隔。

错误二：把一个元素等同于一次事务

错误判断：“每线程只读 4 字节，所以带宽浪费很小。”

纠正：内存系统按事务或 Sector 组织，一个 4 字节请求可能使某一级覆盖 32 字节。

错误三：把 Sector 字节等同于 DRAM 字节

纠正：注明所看层级。L1/TEX、L2 与 DRAM 的计数不同，缓存命中会改变下游流量。

错误四：用不同数据量比较 Stride

例如大 Stride 只访问输入的一个小子集，缓存更容易命中。

纠正：使用互质映射形成全排列，或为每个 Stride 分配足够大的独立输入并确保访问元素数相同。

错误五：输出写入也同时改成 Stride

这样读写两侧都变化，无法判断哪一侧造成下降。

纠正：先让输出连续，只改变输入；再做单独的跨步写实验。

错误六：只执行一次极短 Kernel

几十微秒的结果容易受计时分辨率、启动和频率影响。

纠正：预热、重复多轮、计算平均或中位数，并扩大数据规模。

错误七：在 Debug 下比较

Debug 构建可能关闭优化、增加调试信息并改变代码生成。

纠正：正确性调试可用 Debug，性能发布必须用 Release，并记录编译参数。

错误八：只验证第一个元素

错误索引可能让第一个元素恰好正确，其他位置全部错误。

纠正：教学基准做全量验证；大型生产测试至少结合全量校验、哈希、随机抽样和边界样本。

错误九：看到毫秒下降就宣布合并改善

也许新版少算了元素、删除了边界、改变了精度或提高了缓存复用。

纠正：先证明两个版本输出等价、逻辑工作量一致，再用 Profiler 检查 Sector。

错误十：为了合并访问破坏算法语义

例如强行重排数据，却忘记上游和下游使用的是原顺序。

纠正：把索引映射、逆映射和数据契约写清楚；任何布局转换都必须有独立正确性测试。

---

二十八、从本篇方法迁移到真实项目

28.1 点云属性过滤

业务：从千万点中读取强度 intensity，筛选超过阈值的点。

若数据是：

struct Point {
    float x;
    float y;
    float z;
    float intensity;
    int label;
};

而筛选阶段只读 intensity，AoS 会让相邻强度间隔至少为结构体大小。可以尝试单独维护：

float* intensity;

但还要考虑筛选后的索引是否需要回到坐标数组。合理方案可能是：

SoA 保存热点字段
筛选输出索引
后续只对通过索引读取坐标

性能评估必须包含后续 Gather，不能只优化第一阶段。

28.2 图神经网络邻接访问

图算法天然包含间接索引：

feature = node_features[neighbors[e]];

相邻 Lane 的 neighbors[e] 可能完全随机，合并访问很差。不能简单把下标改为连续，因为那会改变图。

可考虑：

• 对节点重新编号提高局部性；
• 按社区或分区重排；
• 对边排序；
• 缓存高频节点特征；
• 使用适合稀疏访问的库；
• 让 Warp 协作处理同一邻接表。

这些优化本质上仍是控制地址集合，只是约束比向量复制更复杂。

28.3 Embedding 查表

推荐系统与语言模型常见：

embedding[id]

ID 随机时，每个 Lane 可能访问不同长向量。可以根据批次排序 ID、合并重复 ID、提高请求局部性，但排序本身有成本，并且必须恢复输出顺序。

收益条件近似为：

$$t_{memory,saved}>t_{sort}+t_{restore}$$

只有完整流水线测量才能判断。

28.4 数据库列式计算

GPU 上做某列过滤、聚合时，SoA 或列式布局非常自然，因为每个 Kernel 只读取需要的列。行式布局更方便取出完整记录。数据库为什么区分行存与列存，与本篇 AoS/SoA 的矛盾高度相似：访问模式决定布局。

28.5 粒子模拟

一个阶段只更新位置，另一个阶段需要位置与速度，碰撞阶段还需要半径与类型。没有一种布局在所有阶段都绝对最好。

可以：

• 热字段分离；
• 冷字段留在 AoS；
• 按粒子块使用 AoSoA；
• 在阶段边界重排；
• 以总帧时间选择方案。

这就是从教程到工程的关键升级：不问“哪个布局最好”，而问“哪个布局最符合热点阶段的访问模式，转换代价能否被收益覆盖”。

---

二十九、如何写一份可信的内存性能报告

建议每次实验记录下面内容：

实验名称：
日期：
GPU：
Compute Capability：
驱动版本：
CUDA Toolkit：
操作系统：
构建类型：
编译架构：
功耗模式：
数据类型：
N：
输入/输出逻辑字节：
工作集总大小：
Block 大小：
Grid 大小：
预热次数：
计时次数：
计时工具：
计时是否包含复制：
访问公式：
相邻 Lane 字节间隔：
理论 Sector 覆盖：
Event 平均/中位时间：
逻辑有效带宽：
Profiler 请求：
Profiler Sector：
缓存命中：
DRAM 字节：
最大误差：
Mismatch 数：
结论：
未解决问题：

29.1 为什么记录未解决问题

一份专业报告不要求解释所有现象，但必须区分：

• 已由数据支持的结论；
• 依据模型作出的推断；
• 尚未测量的因素；
• 因权限或工具限制未完成的项目。

例如本文 Offset 结果可以写：

已确认：所有 Offset 输出正确，Event 时间处于约 0.63–0.72 ms。
地址推理：Offset 0 的单 Warp 覆盖四个 Sector，非零小偏移可覆盖五个。
推断：相邻 Warp 缓存复用和运行频率可能弱化或改变时间排序。
未完成：本机因 ERR_NVGPUCTRPERM 未采集硬件性能计数。

这比“Offset 2 神奇地比对齐快”更可信。

29.2 报告趋势，不隐藏离散

至少运行五次独立进程，记录每组时间。可以计算：

$$t_{mean}=\frac{1}{K}\sum _{k=1}^K{t}_k$$

标准差：

$$\sigma =\sqrt{\frac{1}{K}\sum _{k=1}^K(t_k-t_{mean}{)}^2}$$

性能数据受到偶发干扰时，中位数往往比平均数更稳健。还可以报告最小值作为接近理想稳定状态的参考，但不能只挑最小值而隐瞒其他结果。

29.3 改动一次只回答一个问题

错误做法：

AoS 改 SoA
Block 从 128 改 512
float 改 half
加入共享内存
删除边界判断

然后看到快了五倍，却不知道每项贡献多少。

正确做法：

基线
  → 只改布局
  → 只改 Block
  → 只改精度
  → 只加共享内存

每一步都保存正确性与性能结果。最后可以组合，但分析阶段必须能够归因。

---

三十、进阶实验任务

下面十二个任务可以把本文从“读懂”变成“掌握”。每项都应保留源码、输出和结论。

任务一：扩展 Offset 到 0–31

把最大 Offset 改为 31，输入多分配 31 个元素。运行五次，画出 Offset 与有效带宽曲线。

思考：

• 是否每八个 Offset 出现周期？
• Offset 0、8、16、24 是否具有相似对齐？
• 测试顺序会不会影响前几组？
• 重启程序后趋势是否稳定？

任务二：改变测试顺序

分别使用：

0 → 1 → 2 → ... → 7
7 → 6 → 5 → ... → 0
随机顺序

如果 Offset 0 总在第一组表现不同，可能与预热或频率状态有关；如果无论顺序都稳定不同，再研究分区与缓存。

任务三：扫描 N

Stride 程序测试：

2^16
2^18
2^20
2^22
2^24

观察工作集跨越缓存容量时曲线如何变化。每个 N 都要保持为 2 的幂。

任务四：增加跨步写

新增 Kernel：

output[(i * stride) & mask] = input[i];

输入读取连续，输出写入跨步。与本文“跨步读、连续写”比较。注意某些写入映射虽然是排列，但写合并和缓存策略与读不同，不能直接套用相同结论。

任务五：同时跨步读写

output[(i * stride) & mask] =
    input[(i * stride) & mask];

与前两项组合，判断两侧地址变化的叠加效果。

任务六：AoS 读取两个字段

改为：

output[i] =
    particles[i].x +
    particles[i].y;

SoA 版本读取 x[i] 和 y[i]。比较只读一个字段、两个字段、四个字段时两种布局的差距。你可能发现随着 AoS 中有用字段增加，它的相对浪费减少。

任务七：实现 AoSoA

每块放 32 或 64 个粒子：

struct ParticleBlock {
    float x[32];
    float y[32];
    float z[32];
    float w[32];
};

让一个 Warp 处理一个块。比较 AoS、SoA 和 AoSoA，并考虑尾部不足一块的处理。

任务八：加入 CSV 输出

输出：

gpu,n,pattern,stride,run,ms,gbps,mismatches

不要只保存格式化表格。CSV 便于绘图、回归检测和不同机器汇总。

任务九：报告中位数

每个配置执行七轮独立计时，排序后取中位数。额外报告最小值、最大值和标准差，判断当前环境是否稳定。

任务十：采集 Nsight Compute

获得性能计数权限后，对 Stride 1、9、65 各采一个 Kernel，记录：

Load Requests
Load Sectors
Sectors per Request
L1/TEX Hit Rate
L2 Hit Rate
DRAM Read Bytes
Kernel Duration

验证 Sector 趋势是否支持地址推理。

任务十一：二维矩阵行列复制

实现两个 Kernel：

一个 Warp 沿行读取
一个 Warp 沿列读取

保证读写元素总数相同，验证结果后比较。再为列读取设计转置或共享内存版本。

任务十二：建立性能回归门禁

不要用固定绝对 GB/s 作为所有机器的失败线。可以在同一机器上比较相对基线：

连续访问带宽不能低于历史中位数的某个合理比例
Stride 趋势必须保持
Mismatch 必须为零

性能回归阈值要允许正常波动，但要能发现 Debug 构建、意外索引变化和缓存行为突变。

---

三十一、问题排查手册

31.1 编译器找不到 CUDA

检查：

nvcc --version
cmake --version

Visual Studio 版本必须与 CUDA Toolkit 支持范围匹配。第一篇已经详细介绍环境搭建；如果 CMake 的 CUDA Compiler 检测失败，优先确认 Toolkit 与 Host Compiler，而不是修改 Kernel。

31.2 no kernel image is available

可执行文件没有包含当前 GPU 可执行的代码。重新设置 CMAKE_CUDA_ARCHITECTURES 或 -arch。不要把本文 RTX 2060 的 75 当成通用值。

31.3 cudaErrorMemoryAllocation

减小 N。AoS/SoA 程序同时分配结构体数组、字段数组和两个输出，Host 侧也有多个向量。总占用不只是屏幕上看到的单个数组大小。

粗略 Device 分配：

$$B_{device}=16N+4N+4N+4N=28N$$

N 为 8,388,608 时约为 224 MiB，不含上下文和其他进程。

31.4 结果全部正确但 GB/s 很低

检查：

• 是否 Release；
• N 是否太小；
• 是否把迭代次数设为 1；
• 是否正在屏幕录制、训练模型或运行游戏；
• 笔记本是否使用电池或节能模式；
• 是否发生热降频；
• 是否远程桌面环境改变 GPU；
• Event 范围是否包含同步之外的操作；
• 当前 GPU 是否真的是预期设备。

31.5 每次数据波动很大

增加预热和迭代，关闭重负载，运行多轮取中位数。不要为了得到整齐结果无限增加迭代，导致 GPU 温度持续上升并进入另一种频率状态。最好记录温度与功率。

31.6 Stride 程序拒绝 N

N 必须是 2 的幂，目的是保证奇数 Stride 与 N 互质，从而形成全排列。若需要任意 N，可以选择与 N 互质的 Stride，并使用取模；程序要计算最大公约数后拒绝不互质组合。

31.7 Nsight Compute 没有采集到 Kernel

可能原因：

• --launch-skip 太大；
• Kernel 名称过滤不匹配；
• 可执行文件路径错误；
• 进程在采集前退出；
• 目标 Kernel 由子进程启动；
• 性能计数权限不足。

先用较小 N 和 --launch-count 1，暂时不跳过，确认能捕获第一条 Kernel。然后逐步加过滤。

31.8 Nsight Compute 分析非常慢

--set full 会采集大量指标，可能需要多次 Replay。先查看可用 Section：

ncu --list-sections

选择内存相关 Section，或只指定需要的 Metrics。不同版本名称不同，以当前工具输出为准。

31.9 Mismatch 不为零

不要先怀疑显卡。逐项检查：

• 输入分配是否包含最大 Offset；
• N 与字节数是否一致；
• Grid 是否覆盖 N；
• 边界判断是否存在；
• Stride 映射的 Host 和 Device 是否一致；
• mask 是否为 N 减 1；
• N 是否为 2 的幂；
• cudaMemcpy 方向是否正确；
• Kernel 后是否检查异步错误；
• 输出是否从正确指针复制。

用很小 N，例如 32 或 64，打印每个 i、Source Index、Expected 和 Actual，更容易定位。

---

三十二、初学者最容易混淆的问答

问：线程越多，内存合并就越好吗？

不是。合并由一个 Warp 内活跃 Lane 的地址分布决定。增加 Grid 可以提供更多并行度和延迟隐藏，但不会自动把分散地址变连续。

问：Block 大小是 256，所以按 256 个线程合并吗？

不是。分析基本单位是 Warp，通常 32 个线程。一个 256 线程 Block 包含八个 Warp，每个 Warp 分别形成请求。

问：只要数组连续存储就一定合并吗？

不是。数组本身连续，但索引可以随机或跨步。必须看相邻 Lane 实际选择的元素。

问：使用 Unified Memory 会自动合并吗？

不会。Unified Memory 改变内存管理和迁移方式，不会改变 Kernel 中 Lane 生成的地址模式。页面迁移与合并访问是不同层面。

问：Pinned Memory 能改善 Device Global Memory 合并吗？

Pinned Host Memory 主要影响 Host/Device 传输和异步复制条件。Kernel 在 Device 内访问 Global Memory时是否合并，仍由 Device 地址模式决定。

问：共享内存一定比合并 Global Memory 快吗？

共享内存延迟低、带宽高，但需要先把数据装入，而且可能有 Bank Conflict、同步和容量限制。若数据只读一次且 Global Memory 已良好合并，盲目经过共享内存可能没有收益。

问：L1 命中很高就不用改 Stride 吗？

要看总时间和业务规模。缓存也许已经掩盖代价，优化收益有限；也可能在更大数据上命中率下降。用代表性工作集测量，不按单一指标做决定。

问：为什么 Offset 的理论利用率是 80%，实测却没有慢 20%？

地址利用率只描述单 Warp 覆盖；相邻 Warp 可能复用边界 Sector。实际时间还包含连续写和其他开销，并受缓存与频率影响。理论比例不是时间比例。

问：为什么 AoS 只读 x，不是只加载 4 字节吗？

单个 Lane 语义上只需要 4 字节，但不同 Lane 的 x 相隔结构体大小，整个 Warp 的地址覆盖更分散。合并分析关注所有 Lane 的地址集合。

问：SoA 快 2.57 倍，是不是所有项目都应重构？

不是。这个结果针对只读 x 的 Kernel。若后续每个线程读取全部字段，或转换成本很高，整体收益可能不同。先在热点流水线上基准。

问：怎样知道优化完成？

至少满足：结果正确；目标数据规模下时间下降；Profiler 指标支持预期原因；没有把成本转移到更慢的转换或传输阶段；代码仍可维护。

---

三十三、本篇必须带走的十八条结论

1. GPU 内存性能由地址模式决定，不只由算术数量决定。
2. 合并访问的分析单位是同一 Warp 的某条内存指令。
3. 计算能力 6.0 及以上普通 Global Memory 合并可从 32 字节事务或 Sector 模型入手。
4. 32 个 Lane 连续读取 32 个 float，有用数据是 128 字节。
5. 自然对齐的连续 128 字节可由四个 32 字节 Sector 覆盖。
6. Offset 可能让四个 Sector 变成五个，但时间不必严格慢百分之二十五。
7. 相邻 Warp 的缓存复用能够缓解边界 Sector 浪费。
8. Stride 增大会把相邻 Lane 地址拆散，显著降低逻辑有效带宽。
9. 公平 Stride 实验必须保持访问元素总数和工作集规模一致。
10. N 为 2 的幂且 Stride 为奇数时，乘法取模映射形成全排列。
11. 有效带宽使用算法有用字节计算，不等于物理 DRAM 带宽。
12. Cache Hit、Sector 和 DRAM Bytes 属于不同层级，不能混用。
13. AoS 与 SoA 的优劣取决于热点 Kernel 消费哪些字段。
14. 只读一个字段时，SoA 往往让相邻 Lane 地址更连续。
15. 布局转换有成本，必须用完整流水线而不是单 Kernel 决策。
16. CUDA Event 负责测正常执行时间，Nsight Compute 负责解释硬件行为。
17. 没有权限采集的计数应标为未测，不能用估算冒充。
18. 可靠优化流程是地址推理、正确性验证、Event 测量、Profiler 证据和整体复测。

一句话总结：

CUDA 内存优化不是让每个线程少写一行代码，而是让一个 Warp 请求的数据尽量集中，让每次搬运都服务更多真正有用的字节。

---

三十四、学习成果验收

完成本篇后，请不要只说“我知道合并访问”。使用下面模板提交自己的实验记录：

我的 GPU：
Compute Capability：
CUDA Toolkit：
构建类型：

Offset N：
Offset 迭代次数：
Offset 0 时间：
Offset 1 时间：
Offset 趋势是否稳定：
我的解释：

Stride N：
Stride 列表：
Stride 1 时间与带宽：
最大 Stride 时间与带宽：
时间倍率：
为什么每组访问相同数量元素：

AoS 结构体大小：
AoS 时间：
SoA 时间：
SoA 加速比：
两个版本最大误差：
两个版本 Mismatch：

Nsight Compute 是否安装：
是否有性能计数权限：
连续访问 Request/Sector：
跨步访问 Request/Sector：
缓存与 DRAM 观察：

我准备优化的真实 Kernel：
一个 Warp 的地址公式：
相邻 Lane 字节间隔：
我准备控制的变量：
我不会改变的变量：
验收指标：

如果你能独立填写这些内容，并解释为什么逻辑带宽不等于 DRAM 带宽，就已经跨过“只会写 Kernel”的阶段，开始具备 GPU 性能工程能力。

---

三十五、下一篇：用共享内存修复矩阵转置

本篇已经看到一个难题：行优先矩阵转置时，输入读取和输出写入很难同时连续。只改变线程映射，常常只是把跨步问题从读端搬到写端。

下一篇：

CUDA编程实战05：共享内存与矩阵转置——把跨步访问改造成合并访问

将完成：

• __shared__ 的声明、生命周期与作用域；
• 一个 Block 怎样协作加载 Tile；
• 为什么需要 __syncthreads()；
• 朴素转置的跨步写问题；
• 共享内存转置怎样让读写两端都合并；
• Shared Memory Bank 与 Bank Conflict；
• 为什么 Tile 常加入一列 Padding；
• 朴素版、分块版、Padding 版完整程序；
• CUDA Event 性能对比；
• Nsight Compute 中 Global Memory 与 Shared Memory 指标；
• 非整数倍矩阵尺寸的边界保护；
• 从转置迁移到卷积、矩阵乘和图像处理。

请保留本篇的 Stride 程序和实验报告。下一篇会先用它证明“跨步为什么慢”，再用共享内存完成地址重排。你会第一次看到一个非常有代表性的 CUDA 优化过程：

发现跨步写
  → 用 Tile 合并读取
  → 在 Shared Memory 中转置
  → 同步
  → 合并写出
  → 处理 Bank Conflict
  → 正确性与性能同时验收

这条路线也是很多高性能图像和矩阵 Kernel 的基本骨架。

---

三十六、官方参考资料

• CUDA Programming Guide：Writing SIMT Kernels
• CUDA C++ Best Practices Guide：Coalesced Access to Global Memory
• NVIDIA Nsight Compute User Guide
• NVIDIA Nsight Compute Profiling Guide
• ERR_NVGPUCTRPERM 官方解决说明

官方文档会随 CUDA 与 Nsight Compute 版本更新。特别是指标名称、缓存行为说明、权限界面与架构规则，应优先查阅与你实际安装版本匹配的文档。本文的地址模型和本机数据用于建立方法，不替代目标架构文档。

---

结语：从“代码一样”看到“地址完全不同”

第四篇真正希望改变的，是你观察 CUDA 程序的角度。

过去看到：

output[i] = input[i];

你可能只看见一次赋值。现在你应该同时看到：

一个 Warp
32 个活跃 Lane
32 个 4 字节请求
连续 128 字节地址范围
四个 32 字节 Sector
连续输出写
逻辑字节口径
缓存复用可能性
Profiler 验证入口

过去看到：

output[i] = input[(i * stride) & mask];

你可能只看见一个更复杂的下标。现在你应该追问：

• Stride 是多少？
• 相邻 Lane 相隔多少字节？
• 地址覆盖多少 Sector？
• 映射是否访问全部输入？
• 输出是否仍连续？
• 工作集是否跨越缓存？
• Event 时间如何变化？
• Nsight Compute 是否看到更多 Sector？
• 优化是否改变结果顺序？

这就是 CUDA 从入门走向实战的分界线：不再把 GPU 当成自动加速的黑盒，而是能够把源代码、线程、地址、事务、缓存、时间和业务价值连成一条因果链。

本文的三个程序分别提供三种能力：

• Offset 教你尊重对齐，同时警惕简单理论与实际时间之间的差距；
• Stride 教你设计公平实验，并直观看到地址分散造成的巨大损失；
• AoS/SoA 教你把底层事务知识转化为业务数据布局决策。

如果这三组实验都在你的机器上通过，并且你能解释自己的结果，那么你已经拥有一套可以反复使用的内存访问分析框架。下一次看到图像处理慢、矩阵列访问慢、粒子字段慢或查表慢时，不必先搜索某个神秘编译参数。先画出一个 Warp 的地址，再用时间和计数器让假设接受检验。

关注并继续完成第五篇。下一次我们不只观察跨步访问有多慢，而会亲手用 Shared Memory 把它重新组织成高效的合并访问。那会是你第一次从“发现瓶颈”走到“设计数据通路并修复瓶颈”的完整 CUDA 优化闭环。