前言
参考资料：
高升博客
《CUDA C编程权威指南》
以及 CUDA官方文档
CUDA编程：基础与实践 樊哲勇
我已经更新了我的Github仓库，大家可以前往仓库下载代码
我的CUDA学习仓库
文章、讲解视频同步更新公众《AI知识物语》，B站：出门吃三碗饭
1：共享内存
共享内存是 一种可被程序员直接操控的缓存，主要作用有两个：
（1）一个是减少核函数中对全局内存的访问次数，实现高效的线程块内部的通信；
（2）一个是提高全局内存访问的合并度。
下面是用C++写的一个归约计算
有 N 个元素的数组 x，假如我们需要计算该数组中所有元素的和，
即 sum = x[0] + x[1] + … + x[N - 1]
#include<stdint.h>
#include<cuda.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <math.h>
#include <stdio.h>
#define CHECK(call) \
do \
{ \
const cudaError_t error_code = call; \
if (error_code != cudaSuccess) \
{ \
printf("CUDA Error:\n"); \
printf(" File: %s\n", __FILE__); \
printf(" Line: %d\n", __LINE__); \
printf(" Error code: %d\n", error_code); \
printf(" Error text: %s\n", \
cudaGetErrorString(error_code)); \
exit(1); \
} \
} while (0)
#ifdef USE_DP
typedef double real;
#else
typedef float real;
#endif
const int NUM_REPEATS = 20;
void timing(const real* x, const int N);
real reduce(const real* x, const int N);
int main(void)
{
const int N = 100000000;
const int M = sizeof(real) * N;
real* x = (real*)malloc(M);
for (int n = 0; n < N; ++n)
{
x[n] = 1.23;
}
timing(x, N);
free(x);
return 0;
}
void timing(const real* x, const int N)
{
real sum = 0;
for (int repeat = 0; repeat < NUM_REPEATS; ++repeat)
{
cudaEvent_t start, stop;
CHECK(cudaEventCreate(&start));
CHECK(cudaEventCreate(&stop));
CHECK(cudaEventRecord(start));
cudaEventQuery(start);
sum = reduce(x, N);
CHECK(cudaEventRecord(stop));
CHECK(cudaEventSynchronize(stop));
float elapsed_time;
CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));
printf("Time = %g ms.\n", elapsed_time);
CHECK(cudaEventDestroy(start));
CHECK(cudaEventDestroy(stop));
}
printf("sum = %f.\n", sum);
}
real reduce(const real* x, const int N)
{
real sum = 0.0;
for (int n = 0; n < N; ++n)
{
sum += x[n];
}
return sum;
}
2：线程同步机制
对于多线程的程序， 两个不同线程中指令的执行次序可能和代码中所展现的次序不同。
要保证核函数中语句的执行顺序与出现顺序一致，就必须使用某种同步机制。 在 CUDA 中，提供了一个同步函数 __syncthreads。该函数只能用在核函数中，其最简单的用法是不带任何参数：
__syncthreads();
该函数可保证一个线程块中的所有线程在执行该语句后面的语句之 前都完全执行了该语句前面的语句。然而，该函数只是针对同一个线程块中的线程的，不同线程块中线程的执行次序依然是不确定的。
3：利用线程同步来归约计算
假如数组元素个数是 2 的整数次方（我们稍后会去掉这个假设），我们可以将数 组后半部分的各个元素与前半部分对应的数组元素相加。如果重复此过程，最后得到的第 一个数组元素就是最初的数组中各个元素的和。这就是所谓的折半归约（binary reduction）法。
3.1 全局内存 条件下的归约计算
void __global__ reduce_global(real* d_x, real* d_y)
{
const int tid = threadIdx.x;
//定义指针X，右边表示 d_x 数组第 blockDim.x * blockIdx.x个元素的地址
//该情况下x 在不同线程块，指向全局内存不同的地址---》使用不同的线程块对dx数组不同部分分别进行处理
real* x = d_x + blockDim.x * blockIdx.x;
//blockDim.x >> 1 等价于 blockDim.x /2，核函数中，位操作比 对应的整数操作高效
for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1)
{
if (tid < offset)
{
x[tid] += x[tid + offset];
}
//同步语句，作用：同一个线程块内的线程按照代码先后执行指令（块内同步，块外不用同步）
__syncthreads();
}
if (tid == 0)
{
d_y[blockIdx.x] = x[0];
}
}
3.2 静态共享内存 条件下的归约计算
void __global__ reduce_shared(real* d_x, real* d_y)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
const int n = bid * blockDim.x + tid;
//定义了共享内存数组 s_y[128]，注意关键词 __shared__
__shared__ real s_y[128];
s_y[tid] = (n < N) ? d_x[n] : 0.0;
__syncthreads();
//归约计算用共享内存变量替换了原来的全局内存变量。这里也要记住： 每个线程块都对其中的共享内存变量副本进行操作。在归约过程结束后，每一个线程
//块中的 s_y[0] 副本就保存了若干数组元素的和
for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1)
{
if (tid < offset)
{
s_y[tid] += s_y[tid + offset];
}
__syncthreads();
}
if (tid == 0)
{
d_y[bid] = s_y[0];
}
}
3.3 动态共享内存 条件下的归约计算
在前面的核函数中，我们在定义共享内存数组时指定了一个固定的长度（128） 。 我们的程序假定了这个长度与核函数的执行配置参数 block_size （也就是核函数中 的 blockDim.x）是一样的。如果在定义共享内存变量时不小心把数组长度写错了，就有可能引起错误或者降低核函数性能。
有一种方法可以减少这种错误发生的概率，那就是使用动态的共享内存
- 在调用核函数的执行配置中写下第三个参数：
<<<grid_size, block_size, sizeof(real) * block_size>>>
前面2个参数网格大小和线程块大小，
第三个参数就是核函数中每个线程块需要 定义的动态共享内存的字节数
- 要使用动态共享内存，还需要改变核函数中共享内存变量的声明方式
extern __shared__ real s_y[];这是动态声明
__shared__ real s_y[128]; 这是静态声明
归约计算程序代码
#include<stdint.h>
#include<cuda.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <math.h>
#include <stdio.h>
#define CHECK(call) \
do \
{ \
const cudaError_t error_code = call; \
if (error_code != cudaSuccess) \
{ \
printf("CUDA Error:\n"); \
printf(" File: %s\n", __FILE__); \
printf(" Line: %d\n", __LINE__); \
printf(" Error code: %d\n", error_code); \
printf(" Error text: %s\n", \
cudaGetErrorString(error_code)); \
exit(1); \
} \
} while (0)
#ifdef USE_DP
typedef double real;
#else
typedef float real;
#endif
const int NUM_REPEATS = 100;
const int N = 100000000;
const int M = sizeof(real) * N;
const int BLOCK_SIZE = 128;
void timing(real* h_x, real* d_x, const int method);
int main(void)
{
real* h_x = (real*)malloc(M);
for (int n = 0; n < N; ++n)
{
h_x[n] = 1.23;
}
real* d_x;
CHECK(cudaMalloc(&d_x, M));
printf("\nUsing global memory only:\n");
timing(h_x, d_x, 0);
printf("\nUsing static shared memory:\n");
timing(h_x, d_x, 1);
printf("\nUsing dynamic shared memory:\n");
timing(h_x, d_x, 2);
free(h_x);
CHECK(cudaFree(d_x));
return 0;
}
void __global__ reduce_global(real* d_x, real* d_y)
{
const int tid = threadIdx.x;
//定义指针X，右边表示 d_x 数组第 blockDim.x * blockIdx.x个元素的地址
//该情况下x 在不同线程块，指向全局内存不同的地址---》使用不同的线程块对dx数组不同部分分别进行处理
real* x = d_x + blockDim.x * blockIdx.x;
//blockDim.x >> 1 等价于 blockDim.x /2，核函数中，位操作比 对应的整数操作高效
for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1)
{
if (tid < offset)
{
x[tid] += x[tid + offset];
}
//同步语句，作用：同一个线程块内的线程按照代码先后执行指令（块内同步，块外不用同步）
__syncthreads();
}
if (tid == 0)
{
d_y[blockIdx.x] = x[0];
}
}
void __global__ reduce_shared(real* d_x, real* d_y)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
const int n = bid * blockDim.x + tid;
//定义了共享内存数组 s_y[128]，注意关键词 __shared__
__shared__ real s_y[128];
//将全局内存中的数据复制到共享内存中
//共享内存的特 征：每个线程块都有一个共享内存变量的副本
s_y[tid] = (n < N) ? d_x[n] : 0.0;
//调用函数 __syncthreads 进行线程块内的同步
__syncthreads();
//归约计算用共享内存变量替换了原来的全局内存变量。这里也要记住： 每个线程块都对其中的共享内存变量副本进行操作。在归约过程结束后，每一个线程
//块中的 s_y[0] 副本就保存了若干数组元素的和
for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1)
{
if (tid < offset)
{
s_y[tid] += s_y[tid + offset];
}
__syncthreads();
}
if (tid == 0)
{
d_y[bid] = s_y[0];
}
}
void __global__ reduce_dynamic(real* d_x, real* d_y)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
const int n = bid * blockDim.x + tid;
//声明 动态共享内存 s_y[] 限定词 extern，不能指定数组大小
extern __shared__ real s_y[];
s_y[tid] = (n < N) ? d_x[n] : 0.0;
__syncthreads();
for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1)
{
if (tid < offset)
{
s_y[tid] += s_y[tid + offset];
}
__syncthreads();
}
if (tid == 0)
{//将每一个线程块中归约的结果从共享内存 s_y[0] 复制到全局内 存d_y[bid]
d_y[bid] = s_y[0];
}
}
real reduce(real* d_x, const int method)
{
int grid_size = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
const int ymem = sizeof(real) * grid_size;
const int smem = sizeof(real) * BLOCK_SIZE;
real* d_y;
CHECK(cudaMalloc(&d_y, ymem));
real* h_y = (real*)malloc(ymem);
switch (method)
{
case 0:
reduce_global << <grid_size, BLOCK_SIZE >> > (d_x, d_y);
break;
case 1:
reduce_shared << <grid_size, BLOCK_SIZE >> > (d_x, d_y);
break;
case 2:
reduce_dynamic << <grid_size, BLOCK_SIZE, smem >> > (d_x, d_y);
break;
default:
printf("Error: wrong method\n");
exit(1);
break;
}
CHECK(cudaMemcpy(h_y, d_y, ymem, cudaMemcpyDeviceToHost));
real result = 0.0;
for (int n = 0; n < grid_size; ++n)
{
result += h_y[n];
}
free(h_y);
CHECK(cudaFree(d_y));
return result;
}
void timing(real* h_x, real* d_x, const int method)
{
real sum = 0;
for (int repeat = 0; repeat < NUM_REPEATS; ++repeat)
{
CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));
cudaEvent_t start, stop;
CHECK(cudaEventCreate(&start));
CHECK(cudaEventCreate(&stop));
CHECK(cudaEventRecord(start));
cudaEventQuery(start);
sum = reduce(d_x, method);
CHECK(cudaEventRecord(stop));
CHECK(cudaEventSynchronize(stop));
float elapsed_time;
CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));
printf("Time = %g ms.\n", elapsed_time);
CHECK(cudaEventDestroy(start));
CHECK(cudaEventDestroy(stop));
}
printf("sum = %f.\n", sum);
}
结果比较：
全局内存，耗时25ms，计算结果错误，应该为1.23*10^8，这里后面有很多小数
静态共享内存，耗时28ms，结果也wrong
动态共享内存
耗时29ms
结论：
（1）全局内存的访问速度是所有内存中最低的，应该尽量减少对它的使用。所有设备内存中，寄 存器是最高效的，但在需要线程合作的问题中，用仅对单个线程可见的寄存器是不够的。我们需要使用对整个线程块可见的共享内存。
（2）使用动态共享内存的核函数和使用静态共享内存的核函 数在执行时间上几乎没有差别。所以，使用动态共享内存不会影响程序性能，但有时可提高程序的可维护性。
（3）使用共享内存来改善全局内存的访问方式并不 一定能够提高核函数的性能。所以，在优化CUDA程序时，一般需要对不同的优化方案进
行测试与比较。
（4）关于计算结果SUM出错，这是因为在累加计算中出现了所谓的“大数吃小数”的现象。单精度浮 点数只有 6、7 位精确的有效数字。在上面的函数 reduce 中，将变量 sum 的值累加到 3000 多万后，再将它和 1.23 相加，其值就不再增加了（小数被大数“吃掉了”，但大数并没有
变化）。
目前的有的解决办法比如： Kahan 求和算法
3826
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
