注：本markdown代码部分格式高亮基于c
代码详见(https://github.com/H-Freax/CUDA_optimization)
试试看优化CUDA上的加法吧(cuda_sum2.cu)
文章目录
先来试试并行化（优化1.0）
在cuda_sum1.cu中，并没有进行并行化，整个程序只有一个thread，效果并不是很好。
这主要是因为GPU的架构导致的，在CUDA中，一般的内容复制到显示记忆体中的部分，即global memory，这些部分是没有cache的，而且存取global memory的时间较长，通常是数百个cycle。由于程序只有一个thread，每次读取global memory 的内容需要等到读取到内容、累加以后才能进行下一步。
由于global memory没有cache，如果想避开存取的巨量时间，就需要利用大量threads。
要如何把计算平方和的程序并行化呢？我们可以把数字分成若干组，分别计算平方和以后进行相加。
一开始，我们可以把最后的和的相加放在CPU上执行。
首先，基于first_cuda.cu，添加以下代码
#define DATA_SIZE 1048576
#define THREAD_NUM 256
修改sumOfSquares为以下内容，即__优化1.0__
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
const int size = DATA_SIZE / THREAD_NUM;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid * size; i < (tid + 1) * size; i++) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
其中，threadIdx.x是CUDA中表示目前的thread是第几个thread的变量，该变量从0开始计算，由于我们设置有256个thread，所以计数为0~255，利用这个变量，我们可以计算各个组的平方和，另外当threadIdx.x=0的时候可以进行计算时间的start。
因为会有256个计算结果，所以result也需要扩大占用的位置。需要修改main函数为
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM);
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time);
int sum[THREAD_NUM];
clock_t time_used;
cudaMemcpy(&sum, result, sizeof(int) * THREAD_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t),
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
最后，在CPU端计算和
int final_sum = 0;
for(int i = 0; i < THREAD_NUM; i++) {
final_sum += sum[i];
}
printf("sum: %d time: %d\n", final_sum, time_used);
final_sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("sum (CPU): %d\n", final_sum);
编译后可以发现结果相同的前提下，速度快了77倍！
来试试基于记忆体的存取模式的优化吧（优化2.0）
显示卡上的记忆体是DRAM（即动态随机存取存储器，主要的作用原理是利用电容内存储电荷的多寡来代表一个二进制比特是1还是0。）因此最有效率的存储方式，是连续进行存储。
上面编写的程序看是连续存储记忆体的位置，但是考虑到thread的执行方式，当一个thread在等待内容的时候，GPU会切换到下一个thread。也就是说实际上执行的顺序是类似
thread 0 -> thread 1 -> thread 2 -> ……
因此，在同一个thread中连续存储记忆体，实际执行的时候并不连续，要让实际执行时连续，应该要让thread 0读取第一个数，thread 1读取第二个数，以此类推。因此需要修改sumOfSquares如下
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid; i < DATA_SIZE; i += THREAD_NUM) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
编译后执行结果相同，又比上一版快了三倍！
如果增加thread数目的数目，就可以看到更好的效率，例如512个，主要取决于GPU的block中最多能有几个thread，同时如果thread数目增加太多，CPU端相加的工作也会变多。
还能有更多的并行化吗？（优化3.0）
上面提到了block，接下来我们来介绍一下block。
在CUDA中，thread可以进行分组，也就是block，一个block中的thread有一个共用的shared memory，可以进行同步工作。不同的block之间的thread不行。接下来我们试试用多个block来进一步增加thread的数目。
首先在#define的位置修改代码
#define DATA_SIZE 1048576
#define BLOCK_NUM 32
#define THREAD_NUM 256
表示接下来我们会用到32个block，每个block有256个threads，一共有32*256=8192个threads
接下来修改sumOfSquares部分
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int sum = 0;
int i;
if(tid == 0) time[bid] = clock();
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
sum += num[i] * num[i];
}
result[bid * THREAD_NUM + tid] = sum;
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
blockIdx.x的用法跟threadIdx.x相同，表示的是block的编号，在这个版本中我们记录每个block的开始以及结束时间。
最后修改main函数部分
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result,
sizeof(int) * THREAD_NUM * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<BLOCK_NUM, THREAD_NUM, 0>>>(gpudata, result,
time);
int sum[THREAD_NUM * BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(&sum, result, sizeof(int) * THREAD_NUM * BLOCK_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
int final_sum = 0;
for(int i = 0; i < THREAD_NUM * BLOCK_NUM; i++) {
final_sum += sum[i];
}
clock_t min_start, max_end;
min_start = time_used[0];
max_end = time_used[BLOCK_NUM];
for(int i = 1; i < BLOCK_NUM; i++) {
if(min_start > time_used[i])
min_start = time_used[i];
if(max_end < time_used[i + BLOCK_NUM])
max_end = time_used[i + BLOCK_NUM];
}
printf("sum: %d time: %d\n", final_sum, max_end - min_start);
基本上就是增加result的大小，修改计算时间的方法，把每个block最早的开始时间减去最晚的结束时间，得到最终的时间。相较上版本又快了很多，但是在CPU上的部分时间增加了，因为CPU上需要加的数字更多了，为了避免这个问题，我们可以让每个block都计算自己的threads的计算结果的和。
来试试Thread的同步（优化4.0）
这个版本中，我们让每个block都计算自己的threads的计算结果的和，把sumOfSquares修改如下
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
__syncthreads();
if(tid == 0) {
for(i = 1; i < THREAD_NUM; i++) {
shared[0] += shared[i];
}
result[bid] = shared[0];
}
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
利用__shared__表示这个变量存在于shared memory，是一个block中每个thread都公用的记忆体，会使用GPU上的记忆体，可以不用担心存取时间的问题。
__syncthreads()是一个CUDA内置的函数，表示把block中的所有thread都同步到这个点再执行。
接下来把main函数部分改成
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<BLOCK_NUM, THREAD_NUM,
THREAD_NUM * sizeof(int)>>>(gpudata, result, time);
int sum[BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(&sum, result, sizeof(int) * BLOCK_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
int final_sum = 0;
for(int i = 0; i < BLOCK_NUM; i++) {
final_sum += sum[i];
}
可以发现，现在CPU上只需要加32个数字就可以了，又有了优化。
但是还是有优化的空间，最终的相加的工作，被分配给每个block的thread 0来进行，并不是最有效的方法，这个相加的动作是可以进行并行化的。
来试试用树状加法优化（优化5.0）
树状加法即透过树型结构的启发
![[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ahpZpI6R-1650191714396)(C:\Users\QianYaoyao\AppData\Roaming\Typora\typora-user-images\image-20220417182758993.png)]](https://code84.com/wp-content/uploads/2022/09/b4c8b6b32e244e5fa628f3726e4eaa99.png)
把sumOfSquares修改如下
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
int offset = 1, mask = 1;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
__syncthreads();
while(offset < THREAD_NUM) {
if((tid & mask) == 0) {
shared[tid] += shared[tid + offset];
}
offset += offset;
mask = offset + mask;
__syncthreads();
}
if(tid == 0) {
result[bid] = shared[0];
time[bid + BLOCK_NUM] = clock();
}
}
还有什么改进空间吗？
上个版本的树状加法在GPU执行的时候可能会存在share memory的bank conflict的问题。
bank conflict是什么？：
在CUDA装置中，shared memory被分成数个bank，如果同时每个thread存取不同的bank，不会存在问题，当两个或多个threads存同一个bank的时候，就会产生bank conflict
可以进行如下改写
offset = THREAD_NUM / 2;
while(offset > 0) {
if(tid < offset) {
shared[tid] += shared[tid + offset];
}
offset >>= 1;
__syncthreads();
}
这时候省去了mask，也有进一步提升，当然，再优化的话，可以展开树状加法
if(tid < 128) { shared[tid] += shared[tid + 128]; }
__syncthreads();
if(tid < 64) { shared[tid] += shared[tid + 64]; }
__syncthreads();
if(tid < 32) { shared[tid] += shared[tid + 32]; }
__syncthreads();
if(tid < 16) { shared[tid] += shared[tid + 16]; }
__syncthreads();
if(tid < 8) { shared[tid] += shared[tid + 8]; }
__syncthreads();
if(tid < 4) { shared[tid] += shared[tid + 4]; }
__syncthreads();
if(tid < 2) { shared[tid] += shared[tid + 2]; }
__syncthreads();
if(tid < 1) { shared[tid] += shared[tid + 1]; }
__syncthreads();