CUDA共享内存 SM(Streaming Multiprocessors)与板块(block) 在先前我们利用网格跨步循环实现了数组的加法,但是其中利用了atomicAdd()这样的原子化操作导致实际上每一个线程都是串行化进行了,并不能发挥GPU的并行计算能力。
如果不使用原子计算的方式,我们就要将sum变为数组,下边的这段代码每一个线程计算一段地址0-1023, 1024-2045 …,将求和的大小缩减到CPU可接受的范围,最后在CPU上进行求和。
1 2 3 4 5 6 7 8 9 __global void parallel_sum (int * sum, int const * arr, int n) { for (int i = blockIdx.x * blockDim.x + threadIdx.x;i < n / 1024 ;i += blockDim.x * gridDim.x){ int local_sum = 0 ; for (int j = i * 1024 ; j < i * 1024 + 1024 ;j++){ local_sum += arr[j]; } sum[i] = local_sum; } }
这种求和的内存访问有很大的缺点,每个线程内部仍旧是串行化的:
线程0依次访问arr[0], arr[1]…
线程1依次访问arr[1024], arr[1025]…
在上述的方式里,local_sum对前一时刻是有依赖关系的,去除依赖关系是提升效率的关键。
利用线程局部数组来并行化处理
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 __global__ void parallel_sum (int * sum, int const * arr, int n) { for (int i = threadIdx.x + blockDim.x * blockIdx.x; i < n / 1024 ; i += blockDim.x * gridDim.x){ int local_sum[1024 ]; for (int j = 0 ;j < 1024 ;j++){ local_sum[j] = arr[i * 1024 + j]; } for (int j = 0 ;j < 512 ;j++){ local_sum[j] += local_sum[j + 512 ]; } for (int j = 0 ;j < 256 ;j++){ local_sum[j] += local_sum[j + 256 ]; } for (int j = 0 ;j < 128 ;j++){ local_sum[j] += local_sum[j + 128 ]; } for (int j = 0 ;j < 64 ;j++){ local_sum[j] += local_sum[j + 64 ]; } for (int j = 0 ;j < 32 ;j++){ local_sum[j] += local_sum[j + 32 ]; } for (int j = 0 ;j < 16 ;j++){ local_sum[j] += local_sum[j + 16 ]; } for (int j = 0 ;j < 8 ;j++){ local_sum[j] += local_sum[j + 8 ]; } for (int j = 0 ;j < 4 ;j++){ local_sum[j] += local_sum[j + 4 ]; } for (int j = 0 ;j < 2 ;j++){ local_sum[j] += local_sum[j + 2 ]; } for (int j = 0 ;j < 1 ;j++){ local_sum[j] += local_sum[j + 1 ]; } sum[i] = local_sum[0 ]; } }
Shared Memory 上面已经实现了无数据依赖可以并行的for循环, 那么如何将他变成真正的并行?对于上面的无数据依赖的for,我们可以将线程升级为板块(block),线程局部数组升级为板块局部数组。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 __global__void parallel_sum (int *sum, int const * arr, int n) { __shared__ int local_sum[1024 ]; int i = threadIdx.x; int j = blockIdx.x; local_sum[j] = arr[i * 1024 + j]; __syncthreads(); if (j < 512 ){ local_sum[j] += local_sum[j + 512 ]; } __syncthreads(); if (j < 256 ){ local_sum[j] += local_sum[j + 256 ]; } __syncthreads(); if (j < 128 ){ local_sum[j] += local_sum[j + 128 ]; } __syncthreads(); if (j < 64 ){ local_sum[j] += local_sum[j + 64 ]; } __syncthreads(); if (j < 32 ){ local_sum[j] += local_sum[j + 32 ]; } __syncthreads(); if (j < 16 ){ local_sum[j] += local_sum[j + 16 ]; } __syncthreads(); if (j < 8 ){ local_sum[j] += local_sum[j + 8 ]; } __syncthreads(); if (j < 4 ){ local_sum[j] += local_sum[j + 4 ]; } __syncthreads(); if (j < 2 ){ local_sum[j] += local_sum[j + 2 ]; } __syncthreads(); if (j == 0 ){ sum[i] = local_sum[0 ] + local_sum[1 ]; } }
这里代码是线程并行的, 如果不加__syncthreads(), 有的代码运行到了j < 32,而有的才运行到 j < 64,这就导致了归约顺序出现了问题。
我们可以进一步优化这个代码,因为线程组是以32个为单位进行了, 一个线程组内部本身就是同步的。在j < 32的时候,不需要加入sync..
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 __global__void parallel_sum (int *sum, int const * arr, int n) { __shared__ int local_sum[1024 ]; int i = threadIdx.x; int j = blockIdx.x; local_sum[j] = arr[i * 1024 + j]; __syncthreads(); if (j < 512 ){ local_sum[j] += local_sum[j + 512 ]; } __syncthreads(); if (j < 256 ){ local_sum[j] += local_sum[j + 256 ]; } __syncthreads(); if (j < 128 ){ local_sum[j] += local_sum[j + 128 ]; } __syncthreads(); if (j < 64 ){ local_sum[j] += local_sum[j + 64 ]; } __syncthreads(); if (j < 32 ){ local_sum[j] += local_sum[j + 32 ]; } if (j < 16 ){ local_sum[j] += local_sum[j + 16 ]; } if (j < 8 ){ local_sum[j] += local_sum[j + 8 ]; } if (j < 4 ){ local_sum[j] += local_sum[j + 4 ]; } if (j < 2 ){ local_sum[j] += local_sum[j + 2 ]; } if (j == 0 ){ sum[i] = local_sum[0 ] + local_sum[1 ]; } }
线程组分歧 在 NVIDIA GPU 中,基本的调度单位是 Warp(线程束/线程组),一个 Warp 包含 32 个线程。 线程组分歧是指:当同一个 Warp 中的 32 个线程在执行代码时,遇到了条件分支(如 if-else、switch 或 while),且一部分线程走向了 A 分支,另一部分线程走向了 B 分支。由于一个 Warp 在同一时刻只能执行同一条指令,硬件无法让线程 0 跑 if 的同时让线程 1 跑 else。
当分歧发生时,GPU 硬件会采取“牺牲时间”的策略:
执行分支 A:此时,所有走向分支 B 的线程会被屏蔽(Masked),处于空闲等待状态。
执行分支 B:此时,所有走向分支 A 的线程被屏蔽,轮到分支 B 的线程工作。
合并:两个分支都跑完后,32 个线程重新汇合,继续执行后面的指令。
结果: 原本应该并行完成的任务,因为分歧变成了串行执行。如果 if 和 else 各占一半,那么该 Warp 的理论计算性能直接损失 50%。
所以我们上边的一段代码实际上还可以进行优化,将j < 32的部分都合并为一个,避免了线程分歧:
1 2 3 4 5 6 7 8 9 10 11 12 13 if (j < 32 ){ local_sum[j] += local_sum[j + 32 ]; local_sum[j] += local_sum[j + 16 ]; local_sum[j] += local_sum[j + 8 ]; local_sum[j] += local_sum[j + 4 ]; local_sum[j] += local_sum[j + 2 ]; if (j == 0 ){ sum[i] = local_sum[0 ] + local_sum[1 ]; } }
寄存器打翻(register spill) 板块中线程数量太多
GPU一个板块中的线程,共享一个共同的寄存器仓库(实际上也是一段比较小的内存), 所以当把模块中的线程数量(blockDim)太多的时候,会导致每个线程分配到的寄存器数量急剧缩小。另一方面, 如果程序需要大量的寄存器,就没有办法全部装在高效的寄存器仓库里,导致一部分寄存器“打翻”到一级缓存中,这些寄存器与寄存器仓库中的读写速度就不一致,相对而言就低效
如果线程局部分配了一个数组,并通过动态下标访问,那么这个数组会被分配到一级缓存中,因为寄存器不能动态寻址。
对于Fermi架构, 每个线程最多有63个寄存器。
板块中的线程太少 每个SM一次调度板块中的一个线程组(warp),也就是32个线程。当线程陷入内存等待的时候,可以切换到另外一个线程,这样就算一个warp的内存延迟就被另外一个warp隐藏起来了。因此如果线程数量太少,则无法通过多个warp之间的调度来隐藏内存等待的延迟。
block中的线程数量最好是32的整数倍, 因为如果是33个线程的情况下, 第二个warp中只有一个线程,非常浪费。
综合上面的情况:对于使用寄存器少、访存为主的核函数,使用大的blockDim为宜。反之,使用较小的blockDim。
Example矩阵转置 首先实现一个简单的矩阵转置版本用一维数组的方式存储。[nx, ny] -> [ny, nx]
1 2 3 4 5 6 7 8 9 10 11 12 template <typename T>__global__ void parallel_transpose (T* out, T* in,int nx, int ny) { int linearized = blockIdx.x * blockDim.x + threadIdx.x; if (linearized > nx * ny) return ; int row = linearized / nx; int col = linearized % nx; if (row < ny && col < nx) { out[col * ny + row] = in[row * nx + col]; } }
这段代码对于内存的访问仍不够高效, 因为in/out都存在于GPU的显存当中,并且对于out的访问是跨步的,不符合内存局部性原理。这个问题可以通过板块共享内存解决。下边这段代码通过跨blockSize步数来进行访存(实际和blockDim没什么区别)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 template <typename T, int blockSize = 32 >__global__ void parallel_transpose (T* out, T* in, int nx, int ny) { int x = blockIdx.x * blockSize + threadIdx.x; int y = blockIdx.y * blockSize + threadIdx.y; if (x >= nx || y >= ny) return ; __shared__ T memory[blockSize * blockSize]; memory[threadIdx.y * blockSize + threadIdx.x] = in[y * nx + x]; __syncthreads(); int rx = blockIdx.y * blockSize + threadIdx.x; int ry = blockIdx.x * blockSize + threadIdx.y; if (rx < ny && ry < nx) { out[ry * ny + rx] = memory[threadIdx.x * blockSize + threadIdx.y]; } }
这里比较反直觉的是为什么多了一步的搬运步骤,反而要比直接从in -> out快呢 gemini: GPU 读写全局显存(Global Memory)时,不是一个线程一个线程读的,而是 Warp(32个线程) 一起读。
合并访问(Coalesced):如果 32 个线程访问的是连续的 128 字节,硬件>只需要 1 次内存事务。
跨步访问(Strided):如果 32 个线程访问的地址散落在各地(比如转置时的列访问),硬件可能需要 32 次内存事务。
结论:跨步访问的带宽浪费可能高达 32 倍。
那么这样看起来像是从显存将数据搬移到共享内存中的速度要比从显存到显存快,事实是这样吗?
在指令的条数上,确实是显存 -> 共享内存 -> 显存的比较长(数据需要从显存搬移到寄存器中再进入到共享内存上),不过现在的GPU似乎可以省去搬到寄存器这个步骤(待求证)
真正使得共享内存快的原因是:全局内存访问是以32/128字节为一个整体为单位访问的,如果跨步很大, 但是只读取一个数据(比如int),那么有效的数据只有4/128,损失的大量的带宽。
Bank Conflict 前一节实现了一个简单的矩阵转置的例子,但是其中仍旧有一个问题,GPU的shared_memory是按照一个个bank为单位进行划分的,总共会划分为32个banks,而线程的访存方式为 。我们期望的方式应该是两个不同的线程i与线程j的访存不冲突,即i与j同时访问不同的bank(如果访问同一个bank,bank内部则需要串行化处理,否则会带来data-race)。
$$\text{Bank Index} = (\text{Address} / 4\text{ bytes}) \pmod{32}$$
PS 就是说
arr[0], arr[32], arr[64]…一系列的数据会位于bank_0
arr[1], arr[33],arr[65]…会位于bank_1
那么在我们转置的过程中,一个warp对于全局内存的访问方式恰好就是这种按列 读取,导致一个warp的所有线程都在同时的读取同一个bank。解决这个问题的方式为padding,让共享内存多出一行来,这样就可以恰好的让每一个线程都”错位”。
给出没有bank conflict的代码:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 template <typename T, int blockSize = 32 >__global__ void parallel_transpose (T* out, T* in, int nx, int ny) { __shared__ T memory[blockSize * (blockSize + 1 )]; int x = blockIdx.x * blockSize + threadIdx.x; int y = blockIdx.y * blockSize + threadIdx.y; if (x < nx && y < ny) { memory[threadIdx.y * (blockSize + 1 ) + threadIdx.x] = in[y * ny + x]; } __syncthreads(); int rx = blockIdx.y * blockSize + threadIdx.x; int ry = blockIdx.x * blockSize + threadIdx.y; if (rx < ny && ry < nx) { out[ry * ny + rx] = memory[threadIdx.x * (blockSize + 1 ) + threadIdx.y]; } }