数据内存布局#

虽然线程数是可以用2D 3D的（逻辑上），但是显存是1D的（物理意义上）。

在前面我们提到过这个cudaMalloc函数，用于GPU分配空间。

1
cudaError_t cudaMalloc (void** devPtr, size_t size)

这个是在GPU显存上，申请一个size字节的连续内存，返回的是是否分配成功，然后分配的GPU显存的指针，会被存在devPtr指向的地址中。

这里分配的显存，是线性的，也就是1D的。

函数限定符#

前面我们用__global__来定义了和函数。CUDA C中还有更多的函数限定符。

__device__ 和 __host__ 可以同时使用，让同一个函数在CPU和GPU上都能用。

SM：流式处理器#

前面说了，当CPU调用一次核函数，比如

1
train_kernel<<<grid_size, block_size>>>(...)

这时候，我们要调用非常非常多的线程。而在调用时，<<<grid, block>>> 是一个软件代码里划分的层次，这里我们可以用二维或者三维来组织这些线程。

不过，在扔进显卡后，这些线程实际上还都是一维的，那么这时候是怎么对应的呢？

实际上，我们用一个Grid（线程格）来组织这么多的线程。而其中，一个Block（线程块）会被分配到一个指定的SM。

注意：是一个Block会被分配给一个SM，但一个SM可以容纳不止一个Block。

一个SM有这样的一些东西：

1. CUDA核心
2. 共享内存、一级缓存（这里共享内存，是每个block占用一份，不是SM上所有Block共享的。所谓共享，是同一个Block内的线程共享。）
3. 寄存器文件
4. 加载/存储单元
5. 特殊功能单元
6. 线程束调度器

一个家用级别的显卡，比如RTX3060，大概有28个SM。而例如A100这种专业的显卡，也只有108个SM左右。

那么，一个任务，会分配给几个SM呢？

先来看硬件上限。一个Block最多占用一个SM，不会跨SM，所以Grid有多少个Block，就最多能用多少个SM（当然也不能超过GPU的SM总数）。

这是上限，那实际呢？每个SM资源有限，主要是三样东西：寄存器数量、共享内存大小、以及 SM 支持的最大线程数。一个 Block 进来会占用这些资源，如果一个 SM 的资源还够，就可以再塞第二个 Block 进来同时跑，这叫做占用率（Occupancy）。

所以实际用几个 SM，是由“Block 数量”和“每个 SM 能同时容纳几个 Block” 共同决定的，硬件调度器自动分配，人是控制不了具体哪个 Block 去哪个 SM。

如果都分配完了呢？

假设一个SM有 48KB 共享内存，每个 Block 申请 20KB，那这个 SM 只能同时容纳 2 个 Block（40KB），第 3 个 Block 塞不进来，就要等前面的 Block 执行完释放资源，才能进来。

如果所有 SM 都满了，剩下的 Block 就排队等待。

而且共享内存只是卡住的原因之一，寄存器也会卡。（比如超级无敌多的局部变量！）

SIMT架构、warp#

CUDA采用的是SIMT，也就是单指令多线程架构。和这个相对的是SIMD，单指令多数据架构。

这看上去比较难懂，对比一下可能比较好。

SIMD（Single Instruction, Multiple Data） 是 CPU 上的概念。一条指令，同时对多个数据做同样的操作。比如一条 AVX 指令可以同时加 8 个 float。但这 8 个数据必须打包在一起，用专门的寄存器，程序员（或编译器）要显式地去安排这件事。

举个例子，

1
for (int i = 0; i < 8; i++) {
2
    c[i] = a[i] + b[i];
3
}

这看上去是循环加了 8 次，但编译器开了 -O2 优化之后，可能直接生成一条 AVX 指令，把 8 个 float 一次加完。没写任何特殊代码，但 SIMD 已经在用了。

SIMT（Single Instruction, Multiple Threads） 是 GPU 的概念。同样是一条指令操作多份数据，但包装方式不同——它伪装成”多个独立线程”，每个线程有自己独立的寄存器、自己的 PC 指针，你写代码的时候感觉像在写普通的单线程逻辑，硬件偷偷把 32 个线程捆成一个 Warp 一起跑。

还是举个例子，要对32个数各自+1，核函数是：

1
__global__ void add_one(int *a) {
2
    int i = threadIdx.x;
3
    a[i] = a[i] + 1;
4
}

这里写的是 第 i 个线程处理第 i 个数，每个线程是独立的逻辑，也就是单线程视角的代码。

但硬件实际干的是，把 Thread 0 到 Thread 31 捆成一个 Warp，同时对 32 个数据执行同一条加法指令。本质上和 SIMD 做的事一模一样。

但是这样就很显然，有一个问题——如果这32个线程，我写的不是顺序执行呢？比如中间有if，那他们有可能同时无法执行同样的指令。

这就是线程束分化（Warp Divergence），这是 SIMT 的核心痛点。

例如：

1
__global__ void kernel(int *a) {
2
    int i = threadIdx.x;
3
    if (i % 2 == 0) {
4
        a[i] = a[i] + 1;  // 偶数线程走这里
5
    } else {
6
        a[i] = a[i] * 2;  // 奇数线程走这里
7
    }
8
}

同一个 Warp 里，Thread 0、2、4… 要执行加法，Thread 1、3、5… 要执行乘法。但硬件要求 32 个线程同时执行同一条指令，怎么办？

硬件的解决方式是串行化：

第一步：让偶数线程执行加法，奇数线程强制闲置（屏蔽掉）

第二步：让奇数线程执行乘法，偶数线程强制闲置（屏蔽掉）

如果我们再极端一点，假设一共有 $n$ 个分叉，总共有 $2^n$ 种可能的指令，就在最坏情况下，需要 $O(2^n)$ 个时间周期。不过这是理论，因为一个warp只有32个线程，不可能有超过 32 条不同的路径同时出现。所以最坏情况不是 $O(2^n)$ ，而是 $O(\min(2^n, 32))$ ，也就是最坏情况是 32 步串行执行。

NOTE

当 $n$ 比较小的时候（比如 $n=1$ 就是 2 条路径，$n=2$ 就是 4 条），$2^n$ 是瓶颈。但 n 一旦超过 5，$2^n$ 就超过 32 了，后面再多分支也不会更坏，因为线程数就这么多。

所以现实中的最坏情况是：Warp 里 32 个线程每个人走的路径都不一样，硬件要串行执行 32 次，每次只有 1 个线程在干活，其余 31 个闲置。效率变成 1/32，退化成完全串行了。

这也是为什么 GPU 并不适合分支密集、逻辑复杂的任务——这类任务交给 CPU 反而更合适。

在并行线程中，可能多个线程都要访问同一个数据，或者写入同一个数据。假设这是共享内存，那怎么办？

比如：

1
__global__ void kernel(int *a) {
2
    a[0] = threadIdx.x;  // 32个线程都往a[0]写
3
}

32 个线程同时写同一个地址，最终 a[0] 是多少？

不确定，取决于谁最后写进去，结果是随机的。

如果确实需要这么做的话，可以用原子操作。比如求和，可以这样写：

1
atomicAdd(&result, value);

原子操作保证每次只有一个线程在操作这个地址，操作完了下一个线程才能进来。这样结果是正确的，但代价是这些操作被串行化了，性能会下降。

算法优化#

1
for (int offset = 4; offset > 0; offset /= 2) {
2
    partial += offset;
3
}

1
partial += 4;
2
partial += 2;
3
partial += 1;

1
#pragma unroll
2
for (int offset = 4; offset > 0; offset /= 2) {
3
    partial += offset;
4
}

1
// 雇佣了 512 个线程
2
int tid = threadIdx.x;
3

4
// 第一轮循环：步长是 256 (512的一半)
5
for (int stride = 256; stride > 0; stride /= 2) {
6
    if (tid < stride) {  // 只有编号小于 256 的线程才准干活！
7
        idata[tid] += idata[tid + stride];
8
    }
9
    __syncthreads(); // 所有人必须等干活的人做完
10
}

1
// 我们现在只雇佣 256 个线程了！网格大小直接砍半。
2
int tid = threadIdx.x;
3

4
// 在进入循环前，让 0 号线程去把 0 号和 256 号加起来；
5
// 让 1 号线程去把 1 号和 257 号加起来...
6
g_idata[idx] += g_idata[idx + blockDim.x];
7

8
__syncthreads();
9

10
// 此时，512 个数字已经被折叠成 256 个了。
11
// 接下来再舒舒服服地进入原先的循环。
12
for (int stride = 128; stride > 0; stride /= 2) {
13
    // ... 原来的逻辑 ...
14
}

1
// 只让前 32 个线程进来
2
if (tid < 32) {
3
    volatile int *vmem = idata;
4

5
    // 以下 6 行代码，没有 for 循环，没有 __syncthreads() 互相等待
6
    // 32 个线程靠着天然的硬件同步，像机关枪一样一口气打完！
7
    vmem[tid] += vmem[tid + 32]; // 这是应对如果刚好剩 64 个数据的情况
8
    vmem[tid] += vmem[tid + 16];
9
    vmem[tid] += vmem[tid + 8];
10
    vmem[tid] += vmem[tid + 4];
11
    vmem[tid] += vmem[tid + 2];
12
    vmem[tid] += vmem[tid + 1];
13
}