《CUDA编程》8.共享内存的合理使用

共享内存是一种可被程序员直接操控的缓存，主要作用有两个：
①减少核函数中对全局内存的访问次数，实现高效的线程块内部的通信
②提高全局内存访问的合并度
将通过两个具体的例子阐明共享内存的合理使用，一个数组归约的例子和讨矩阵转置的例子

1 例子：数组归约计算

一个有 $N$ 个元素的数组 $x$ 。假如我们需要计算该数组中所有元素的和，即 $s u m = x [0] + x [1] + ... + x [N - 1]$ ，这里先给出C++的代码：

float cumulative_sum(const float* x, int N) {float sum = 0.0;for (int i = 0; i < N; i++) {sum += x[i];}return sum;
}

在这个例子中，我们考虑一个长度为 $10^{8}$ 的一维数组，在主函数中，我们将每个数组元素初始化为 1.23，调用函数 reduce 并计时。

在使用双精度浮点数时，输出：sum = 123000000.110771，该结果前 9 位有效数字都正确，从第 10位开始有错误，运行速度为315ms。
在使用单精度浮点数时，输出：sum = 33554432.000000.，该结果完全错误，运行速度为315ms。

这是因为，在累加计算中出现了所谓的"大数吃小数"的现象。单精度浮点数只有 6、7 位精确的有效数字。在上面的函数中，将变量 sum 的值累加到 3000多万后，再将它和1.23相加，其值就不再增加了（小数被大数"吃掉了"，但大数并没有变化）。
当然现在已经有了其他安全的算法，但我们在CUDA 实现要比上述 C++ 实现稳健（robust）得多，使用单精度浮点数时结果也相当准确。

1.1 只使用全局内存

1.1.1 运行代码

#include <cuda_runtime.h>
#include <iostream>
#include <iomanip> 
#include "error_check.cuh"#define TILE_DIM 32  // 定义每个block的线程块维度// 核函数 
__global__ void Array_sum(float* d_x, float* d_re)
{const int tid = threadIdx.x;float* x = d_x + blockIdx.x * blockDim.x; // 使用 blockIdx.x 来索引数据// 归约求和for (int stride = blockDim.x >> 1; stride > 0; stride >>= 1){if (tid < stride){x[tid] += x[tid + stride]; // 将数据进行归约}__syncthreads(); // 同步线程，确保所有线程完成了归约计算}// 只在第一个线程中写入结果if (tid == 0){d_re[blockIdx.x] = x[0];}
}int main() {// 定义一维数组大小const int N = 100000000;const int size = N * sizeof(float);// 主机上分配内存float* h_A = (float*)malloc(size);const int gridSize = (N + TILE_DIM - 1) / TILE_DIM; // 计算网格数量float* h_re = (float*)malloc(gridSize * sizeof(float)); // 结果数组的大小// 初始化数组数据都为1.23for (int i = 0; i < N; i++){h_A[i] = 1.23f;}// 在设备上分配内存float* d_A, * d_re;CHECK(cudaMalloc((void**)&d_A, size));CHECK(cudaMalloc((void**)&d_re, gridSize * sizeof(float))); // 修正结果数组的大小// 将主机数组数据拷贝到设备CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));// 创建线程块和线程块网格const int blockSize = TILE_DIM; // 每个线程块的线程数量dim3 threads(blockSize); // 一维线程块dim3 gridSize2(gridSize); // 一维网格// 实现计时cudaEvent_t start, stop;CHECK(cudaEventCreate(&start));CHECK(cudaEventCreate(&stop));CHECK(cudaEventRecord(start));// 调用核函数，进行数组归约Array_sum << <gridSize2, threads >> > (d_A, d_re);CHECK(cudaEventRecord(stop));CHECK(cudaEventSynchronize(stop));float milliseconds = 0;CHECK(cudaEventElapsedTime(&milliseconds, start, stop));// 输出运行时间，单位是msstd::cout << "运行时间：" << milliseconds << "ms" << std::endl;// 将结果 re 从设备拷贝到主机CHECK(cudaMemcpy(h_re, d_re, gridSize * sizeof(float), cudaMemcpyDeviceToHost));// 计算最终结果double total_sum = 0.0;for (int i = 0; i < gridSize; i++) {total_sum += h_re[i];}// 输出结果,精度为小数点后10位std::cout << std::fixed << std::setprecision(6); // 你可以根据需要调整精度std::cout << "最终结果：" << total_sum << std::endl;// 释放主机和设备内存free(h_A);free(h_re);CHECK(cudaFree(d_A));CHECK(cudaFree(d_re));return 0;
}

输出结果
在这里插入图片描述
观察结果发现，在CUDA中使用单精度进行计算，不仅运算结果正确，而且速度也比c++代码快了5倍（我的设备是GTX 1650）

1.1.2 分析代码

①__syncthreads()函数
在核函数中：

for (int stride = blockDim.x >> 1; stride > 0; stride >>= 1){if (tid < stride){x[tid] += x[tid + stride]; // 将数据进行归约}__syncthreads(); // 同步线程，确保所有线程完成了归约计算}

在归约操作后面使用了__syncthreads()函数，是因为核函数操作是多线程计算的，所以可能上一个归约操作还没有完成，下一个归约操作就开启了，可能导致计算错误，为了保证顺序进行，所以使用该函数。

②位运算

for (int stride = blockDim.x >> 1; stride > 0; stride >>= 1){}

这里使用“右移一位”来替代“除以2”的操作，因为位运算的速度更快

1.2 使用共享内存

全局内存的访问速度是所有内存中最低的，应该尽量减少对它的使用。寄存器是最高效的，但在需要线程合作的问题中，用仅对单个线程可见的寄存器是不够的。

所以共享内存成为最佳选择，因为它提供了一个全局可见、快速且高效的存储空间，供同一个线程块内的所有线程使用。

在核函数中，要将一个变量定义为共享内存变量，就要在定义语句中加上一个限定符__shared__。一般情况下，我们需要的是一个长度等于线程块大小的数组。

1.2.1 修改代码

修改了核函数代码和TILE_DIM 的值，保证block的大小和共享内存的大小一致，其余不变

#include <cuda_runtime.h>
#include <iostream>
#include <iomanip> 
#include "error_check.cuh"#define TILE_DIM 128  // 定义每个block的线程块维度// 核函数 
__global__ void Array_sum(float* d_x, float* d_re)
{const int tid = threadIdx.x;const int bid = blockIdx.x;const int n   = bid * blockDim.x + tid;__shared__ float s_y[TILE_DIM];s_y[tid] = (n < 100000000) ? d_x[n] : 0.0;__syncthreads();// 归约求和for (int stride = blockDim.x >> 1; stride > 0; stride >>= 1){if (tid < stride){s_y[tid] += s_y[tid + stride]; // 将数据进行归约}__syncthreads(); // 同步线程，确保所有线程完成了归约计算}// 只在第一个线程中写入结果if (tid == 0){d_re[blockIdx.x] = s_y[0];}
}int main() {// 定义一维数组大小const int N = 100000000;const int size = N * sizeof(float);// 主机上分配内存float* h_A = (float*)malloc(size);const int gridSize = (N + TILE_DIM - 1) / TILE_DIM; // 计算网格数量float* h_re = (float*)malloc(gridSize * sizeof(float)); // 结果数组的大小// 初始化数组数据都为1.23for (int i = 0; i < N; i++){h_A[i] = 1.23f;}// 在设备上分配内存float* d_A, * d_re;CHECK(cudaMalloc((void**)&d_A, size));CHECK(cudaMalloc((void**)&d_re, gridSize * sizeof(float))); // 修正结果数组的大小// 将主机数组数据拷贝到设备CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));// 创建线程块和线程块网格const int blockSize = TILE_DIM; // 每个线程块的线程数量dim3 threads(blockSize); // 一维线程块dim3 gridSize2(gridSize); // 一维网格// 实现计时cudaEvent_t start, stop;CHECK(cudaEventCreate(&start));CHECK(cudaEventCreate(&stop));CHECK(cudaEventRecord(start));// 调用核函数，进行数组归约Array_sum << <gridSize2, threads >> > (d_A, d_re);CHECK(cudaEventRecord(stop));CHECK(cudaEventSynchronize(stop));float milliseconds = 0;CHECK(cudaEventElapsedTime(&milliseconds, start, stop));// 输出运行时间，单位是msstd::cout << "运行时间：" << milliseconds << "ms" << std::endl;// 将结果 re 从设备拷贝到主机CHECK(cudaMemcpy(h_re, d_re, gridSize * sizeof(float), cudaMemcpyDeviceToHost));// 计算最终结果double total_sum = 0.0;for (int i = 0; i < gridSize; i++) {total_sum += h_re[i];}// 输出结果,精度为小数点后10位std::cout << std::fixed << std::setprecision(6); // 你可以根据需要调整精度std::cout << "最终结果：" << total_sum << std::endl;// 释放主机和设备内存free(h_A);free(h_re);CHECK(cudaFree(d_A));CHECK(cudaFree(d_re));return 0;
}

输出结果
在这里插入图片描述
观察结果又比使用全局内存的CUDA程序快了20ms

1.2.2 分析代码

主要看核函数代码：

核函数中定义了一个共享内存数组s_y[TILE_DIM];
通过一个三元运算s_y[tid] = (n < 100000000) ? d_x[n] : 0.0;，来把全局内存的数据复制到共享内存中，数组大小之外的赋值为0，即不对计算产生影响
调用函数 __syncthreads 进行线程块内的同步
归约计算过程中，用共享内存变量替换了原来的全局内存变量
因为共享内存变量的生命周期仅仅在核函数内，所以必须在核函数结束之前将共享内存中的某些结果保存到全局内存，所以if (tid == 0)判断语句会把共享内存中的数据复制给全局内存

以上的共享内存使用方式叫做使用静态共享内存，因为共享内存数组的长度是固定的。

1.3 使用动态共享内存

在上面的核函数中，我们在定义共享内存数组是一个固定的长度，且程序让该长度和block_size是一致的，但是如果在定义共享内存变量时不小心把数组长度写错了，就有可能引起错误或者降低核函数性能。

有一种方法可以减少这种错误发生的概率，那就是使用动态的共享内存，使用方法如下：

①在核函数调用代码中，写下第三个参数：

// 调用核函数，进行数组归约
Array_sum << <gridSize2, threads,sizeof(float)*blockSize>> > (d_A, d_re);

第三个参数就是核函数中每个线程块需要定义的动态共享内存的字节数，没写的时候默认是0

②改变核函数中共享内存的声明方式：
使用extern限定词，且不能指定数组大小

extern __shared__ float s_y[];

输出结果变化不大：
在这里插入图片描述

2 例子：矩阵转置

2.1 运行代码

在矩阵转置问题中，对全局内存的读和写这两个操作，总有一个是合并的，另一个是非合并的，那么利用共享内存可以改善全局内存的访问模式，使得对全局内存的读和写都是合并的，依然使用行索引转列索引，代码如下：

#include <cuda_runtime.h>
#include <iostream>
#include "error_check.cuh"#define TILE_DIM 32  // 定义每个block的线程块维度__global__ void cpy_matrix(const float* A, float* B, const int N) {__shared__ float S[TILE_DIM][TILE_DIM ];  // 动态共享内存不能直接定义二维数组int nx1 = blockIdx.x * TILE_DIM + threadIdx.x;  // 计算当前线程的列索引int ny1 = blockIdx.y * TILE_DIM + threadIdx.y;  // 计算当前线程的行索引if (nx1 < N && ny1 < N) {// 列索引转行索引，实现矩阵转置S[threadIdx.y][threadIdx.x] = A[ny1 * N + nx1];}__syncthreads();// 转置后的线程索引（交换 x 和 y）int nx2 = blockIdx.y * TILE_DIM + threadIdx.x;int ny2 = blockIdx.x * TILE_DIM + threadIdx.y;if (nx2 < N && ny2 < N) {B[ny2 * N + nx2] = S[threadIdx.x][threadIdx.y];  // 从共享内存写入全局内存}
}int main() {// 定义矩阵大小const int N = 1024;const int size = N * N * sizeof(float);// 主机上分配内存float* h_A = (float*)malloc(size);float* h_B = (float*)malloc(size);// 初始化矩阵数据for (int i = 0; i < N * N; i++) {h_A[i] = 1.0f;}// 在设备上分配内存float* d_A, * d_B;CHECK(cudaMalloc((void**)&d_A, size));CHECK(cudaMalloc((void**)&d_B, size));// 将主机矩阵数据拷贝到设备CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));CHECK(cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice));// 创建线程块和线程块网格dim3 threads(TILE_DIM, TILE_DIM);  // 使用32x8的线程块来减少银行冲突dim3 gridSize((N + TILE_DIM - 1) / TILE_DIM, (N + TILE_DIM - 1) / TILE_DIM);// 计算核函数运行时间cudaEvent_t start, stop;CHECK(cudaEventCreate(&start));CHECK(cudaEventCreate(&stop));CHECK(cudaEventRecord(start));// 调用核函数，将矩阵 A 转置并复制到矩阵 Bcpy_matrix << <gridSize, threads >> > (d_A, d_B, N);CHECK(cudaEventRecord(stop));CHECK(cudaEventSynchronize(stop));float milliseconds = 0;CHECK(cudaEventElapsedTime(&milliseconds, start, stop));std::cout << "kernel time: " << milliseconds << "ms" << std::endl;// 将矩阵 B 从设备拷贝到主机CHECK(cudaMemcpy(h_B, d_B, size, cudaMemcpyDeviceToHost));// 释放主机和设备内存free(h_A);free(h_B);cudaFree(d_A);cudaFree(d_B);return 0;
}

输出结果：
在这里插入图片描述
比《CUDA编程》7.全局内存的合理使用中的行索引转列索引的0.35ms要快，但是还是比列索引慢。

2.2 分析代码

if (nx1 < N && ny1 < N) {// 列索引转行索引，实现矩阵转置S[threadIdx.y][threadIdx.x] = A[ny1 * N + nx1];}__syncthreads();// 转置后的线程索引（交换 x 和 y）int nx2 = blockIdx.y * TILE_DIM + threadIdx.x;int ny2 = blockIdx.x * TILE_DIM + threadIdx.y;if (nx2 < N && ny2 < N) {B[ny2 * N + nx2] = S[threadIdx.x][threadIdx.y];  // 从共享内存写入全局内存}

首先按行索引把全局内存的数据复制到共享内存，这一操作是顺序操作，所以是合并访问（第7章已经讨论过），速度较快
共享内存数据按照列索引将数据复制回全局内存中去，这一步不是顺序访问，但是由于共享内存速度快，弥补了非合并访速度慢的缺点，所以最后的运行速度也快上不少

3 避免共享内存的bank冲突

关于共享内存，有一个内存 bank 的概念值得注意，。为了获得高的内存带宽，共享内存在物理上被分为 32 个同样宽度的、能被同时访问的内存 bank。

①bank冲突定义： 多个线程同时访问共享内存的同一个bank，导致这些访问不能被并行处理，从而降低性能，如下示意图：
在这里插入图片描述
②为什么上述代码会产生bank冲突

__shared__ float S[TILE_DIM][TILE_DIM]; // TILE_DIM = 32创建了一个
32x32 的共享内存数组，表示一个总共 1024 个浮点数的数组。

S[threadIdx.y][threadIdx.x] = A[ny1 * N + nx1];这段代码中，每个线程根据其 threadIdx.x 和 threadIdx.y 的值访问共享内存。如果 threadIdx.y 为 0，threadIdx.x 从 0 到 31 的所有线程将依次访问：S[0][0]、S[0][1]、S[0][31]

和上图一比，就发现实际是访问的同步一个bank，所以产生了bank冲突，而且每个bank收到了32个访问，所以是32路bank冲突

③解决方法
通常可以用改变共享内存数组大小的方式来消除或减轻共享内存的 bank 冲突

__shared__ float S[TILE_DIM][TILE_DIM + 1];  // +1 用于避免银行冲突

输出结果：
在这里插入图片描述
解决了bank冲突后，运行速度又提升了，现在只比列索引慢了1ms，所以合理的使用共享内存，可以有效改善非合并访问的性能瓶颈

因为让原本的32列共享内存数组变成了33列，线程的访问模式将会变得更加分散，bank 的访问更加均匀，从而避免了多个线程同时请求同一 bank 的情况。