最近打算进行hipe相关的开发，hipe与CUDA还是比较像的，所以打算入门一下CUDA，希望顺着CUDA可以入门CUDA的相关编程。
我是通过《GPU高性能计算之CUDA》这本书来学习这些内容的，这本书已经比较早了，所以还保留着很多CUDA已经抛弃的旧特性，在这篇笔记中不会涉及这些内容。
目录
1、CUDA基础
CUDA是为英伟达GPU进行通用计算编程的工具（CUDA的WiKi）。
这个章节会简单过一下CUDA的基本概念，会占据近一半的篇幅。
1.1、CUDA编程模型
1.1.1、主机与设备
CUDA编程模型分为主机和设备两个部分，主机部分就是CPU代码，编程的方式和一般的C代码一模一样，主要执行串行代码，而设备部分就是GPU要运行的部分，一般执行有较大并行潜力的代码。在程序中主机代码和设备代码交替执行。
1.1.2、kernal函数的定义与调用
设备部分代码，也就是GPU要执行的内容。是一种经过扩展的C语言语法。这里先引入一个例子，讲述了使用GPU来完成两个大的一维数组input_arr1和input_arr2对应为相加，将按位相加的结果存在output_arr中。也就是假设：
输入：input_arr1 = [1，1，1]，input_arr2 = [2，2，2]
输出：output_arr = [3，3，3]
但是数组其实很大就是了。
下面就是这个例子的原码（原码链接）：
#include
#include
#include
//这个程序展露了CUDA最为基本的编程模型，主要包括将数据在CPU中建好，然后在GPU上申请空间，
//将CPU的数据拷贝到GPU中，然后在GPU中进行计算，最后将数据从GPU拷贝到CPU中。
//https://blog.csdn.net/xiaohu2022/article/details/79599947#cuda编程模型基础
//Tesla K40使用的是15个SM，每个SM有192个SP
//这个程序的主要工作就是将两个数组进行并行叠加。
//传入三个指针作为形参，两个是输入数组，一个是输出数组
__global__ void add_arr(int* input1, int* input2, int* output){
//这个函数处理拷贝进来的两个数组
//首先查看当前线程的编号，让一个线程去处理一个位置的计算
//这里注意几个用来算索引位置的结构体
//threadIdx、blockIdx、blockDim、gridDim，分别记录维度和索引，内建变量还有warpSize。
//因为只是一个维度的网格，所以就不需要考虑网格层面的事情
// printf("进入内核函数\n");
int thread_idx = threadIdx.x + blockIdx.x * blockDim.x;
//执行一次加法操作
output[thread_idx] = input1[thread_idx] + input2[thread_idx];
//查看计算的结果
// if(thread_idx == 1){
// printf("%d", output[thread_idx]);
// }
}
int main(void)
{
// cudaSetDevice(1);
//数组总大小
int N = 15 * 192;
// int N = 1;
//创建三个数组
int* input_arr1;
int* input_arr2;
int* output_arr;
//分配三个数组
input_arr1 = (int *)malloc(N*sizeof(int));
input_arr2 = (int *)malloc(N*sizeof(int));
output_arr = (int *)malloc(N*sizeof(int));
//用来遍历的变量
int i;
//为三个数组初始化空间
for(i = 0; i < N; i++){
input_arr1[i] = 1;
input_arr2[i] = 2;
output_arr[i] = 0;
}
printf("初始化完毕\n");
//申请GPU空间
int* device_input_arr1 = NULL;
int* device_input_arr2 = NULL;
int* device_output_arr = NULL;
//使用cudaMalloc来申请空间
//函数原型__host____device__cudaError_t cudaMalloc (void** devPtr, size_t size)
//注意这里要传入的是指针的地址
cudaMalloc((void **)&device_input_arr1, N*sizeof(int));
cudaMalloc((void **)&device_input_arr2, N*sizeof(int));
cudaMalloc((void **)&device_output_arr, N*sizeof(int));
printf("空间分配完毕\n");
//这里进行数据拷贝
//函数原型 __host__cudaError_t cudaMemcpy (void *dst, const void *src, size_t count, cudaMemcpyKind kind)
//方向分为两个：cudaMemcpyHostToDevice、cudaMemcpyDeviceToHost
cudaMemcpy((void *)device_input_arr1, (void *)input_arr1, N*sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy((void *)device_input_arr2, (void *)input_arr2, N*sizeof(int), cudaMemcpyHostToDevice);
printf("数据拷贝完毕\n");
//这里开始建立网格和block，我们建立一维网格，一个block内也建立一维线程，使用dim3这个结构
dim3 grid_dimension(15);
dim3 block_dimension(192);
//将两个网格结构传入函数中
add_arr<<>>(device_input_arr1, device_input_arr2, device_output_arr);
printf("退出内核函数\n");
//将结果数组从拷出来
cudaMemcpy((void *)output_arr, (void *)device_output_arr, N*sizeof(int), cudaMemcpyDeviceToHost);
//看看结果对不对
for(i = 0; i < N; i++){
if(output_arr[i] != 3){
printf("发现错误,%d\n", i);
break;
}
}
//将三个数组注销
free(input_arr1);
free(input_arr2);
free(output_arr);
cudaFree(device_input_arr1);
cudaFree(device_input_arr2);
cudaFree(device_output_arr);
return 0;
}
先不看运行代码的内容，就看一个函数的定义与调用。第14行__global__ void add_arr(int* input1, int* input2, int* output)所对应的就是kernal函数的定义，这段函数对应的代码会在GPU在运行。第96行的代码就是对于这个kernal函数的调用。在main函数中的其他代码和一般的C语言是无异的。而kernal函数就是在串行代码中调用内核函数，这个时候GPU就会开始运行kernal函数所规定的内容。在<<<>>>>中实际上就规定了要开多少个线程，线程的结构是什么样的。
1.1.3、线程结构
与CPU不同，GPU拥有成百上千个可以并发运行的计算单元，在理想状况下，GPU可以同时高效运行成百上千个线程，这是CPU做不到的。但是GPU也有一个限制，那就是在同时启动的线程中，每个线程都要执行同一套代码，这套代码就是在kernal函数中规定的。
线程结构的第一个层次叫做网格，每个网格中有很多的线程块，可以是一维、二维也可以是三维的分布。可以针对实际我们面对的问题，来建立不同维度的线程网格。每个线程块中又有比较多的线程，线程在线程块中也是可以按照多个维度来分布的。
在上面这个例子中我们就建立了一维的线程网格和一维的线程块。
//这里开始建立网格和block，我们建立一维网格，一个block内也建立一维线程，使用dim3这个结构
dim3 grid_dimension(15);
dim3 block_dimension(192);
//将两个网格结构传入函数中
add_arr<<>>(device_input_arr1, device_input_arr2, device_output_arr);
在这个网格中一共15个线程块依次排开，每个线程块负责线性数组input_arr1、input_arr2、output_arr的一部分连续的内容。而每一个线程又分别应对一个线程块负责内容的一部分更小的内容，在这个例子中就是某一位的相加。
有一系列的内建结构体可以在kernal函数中被使用，来让当前线程知道自己在全局上到底是第几号线程，并且负责哪些数据的计算，这些内建变量就是：threadIdx、blockIdx、blockDim、gridDim。代表线程在块内的位置以及线程所在块在网格中的位置。
1.1.4、硬件映射
线程网格与硬件对应
在GPU中有一系列的结构来负责一个线程块和一个线程要计算的内容。GPU的硬件结构如下：
一个GPU又多个SM构成，每个SM又由一系列的SP构成。一个SM的不同SP是可以很方便地进行通信的。在执行内核函数的时候，一个线程块会被分配给一个SM，而这个线程块中的线程又会分配给一系列的SP。
线程束
CUDA会将相邻的32个线程合并为一个“线程束”，线程束是GPU指令执行的基本单位，在一个线程束中每个线程执行的指令是严格一致，并且是完全同步的。这就使得我们需要谨慎对待分支语句和合并访存的代码。
1.2、CUDA软件体系
CUDA是一个CPU程序，通常来讲我们使用CUDA运行时API来进行开发。
CUDA C语言
CUDA在C语言上做了如下扩展：
1、函数类型限定符：device、host__和__global
2、变量类型限定符：包括：device、shared 和__constant__
3、矢量类型：char4、 ushort3、 double2、 dim3
4、内建变量：blockIdx、threadIdx、gridDim、blockDim、warpSize
5、<<<>>>，第一个参数是网格的结构，第二个是线程块的结构，第三个参数是共享内存的大小，第四个参数是这个内核函数在哪个流上运行的流的编号。前两个参数在前文的例子（链接）中我们用一个一维网格和一维线程块已经举了例子。
6、一些新的函数
1.3、存储器模型
寄存器
每个SP中都有一定数量的寄存器，拥有极快的访问速度。在上面的例子（链接）中，int thread_idx = threadIdx.x + blockIdx.x * blockDim.x;中thread_idx就是使用寄存器来存储的。
局部存储器
在内核函数中，申请的数组、寄存器存不下的数据都被放在SP的局部存储器中，局部存储器本质上就是显存，实际上还是比较慢的。
共享存储器
这是一种高速存储，并且可以被同一个线程块中的线程共享。在数组前面加上__shared__作为限定符，就可以吧数组声明为一个在共享存储器中的数组。
共享存储器一共有两种初始化和声明方式：__shared__ float array[10]还有extern __shared__ float array[]。这两种方式有点像CPU编程中的堆区和栈区数组，不加extern关键字的类似于栈区数组，需要这一上来，而加extern关键字的类似于堆区数组，可以在程序运行的时候根据实际情况来申请空间，而申请空间的大小需要在核函数调用的时候在<<<>>>第三个参数处规定。
全局存储器
全局储存器就是显存，在上文的例子中（链接），我们需要将两个大数组的各自位依次相加，这需要将两个数组传到GPU，然后需要将数组传回来：
//申请GPU空间
int* device_input_arr1 = NULL;
int* device_input_arr2 = NULL;
int* device_output_arr = NULL;
//使用cudaMalloc来申请空间
//函数原型__host____device__cudaError_t cudaMalloc (void** devPtr, size_t size)
//注意这里要传入的是指针的地址
cudaMalloc((void **)&device_input_arr1, N*sizeof(int));
cudaMalloc((void **)&device_input_arr2, N*sizeof(int));
cudaMalloc((void **)&device_output_arr, N*sizeof(int));
printf("空间分配完毕\n");
//这里进行数据拷贝
//函数原型 __host__cudaError_t cudaMemcpy (void *dst, const void *src, size_t count, cudaMemcpyKind kind)
//方向分为两个：cudaMemcpyHostToDevice、cudaMemcpyDeviceToHost
cudaMemcpy((void *)device_input_arr1, (void *)input_arr1, N*sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy((void *)device_input_arr2, (void *)input_arr2, N*sizeof(int), cudaMemcpyHostToDevice);
printf("数据拷贝完毕\n");
//这里开始建立网格和block，我们建立一维网格，一个block内也建立一维线程，使用dim3这个结构
dim3 grid_dimension(15);
dim3 block_dimension(192);
//将两个网格结构传入函数中
add_arr<<>>(device_input_arr1, device_input_arr2, device_output_arr);
printf("退出内核函数\n");
//将结果数组从拷出来
cudaMemcpy((void *)output_arr, (void *)device_output_arr, N*sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(device_input_arr1);
cudaFree(device_input_arr2);
cudaFree(device_output_arr);
这里展现了和全局内存相关的三个函数：cudaMalloc、cudaMemcpy、cudaFree。
除了使用cudaMalloc来分配线性数组之外，我们还可以使用__device__限定符来声明全局存储器，这种声明方式需要在所有的函数需要在所有函数之外进行声明，作为全局变量使用。
当然对于二维和三维的数组，我们需要使用cudaMallocPitch()和cudaMalloc3D()。这些函数在分配上可以满足对齐要求。至于为什么要“对齐”，什么叫“对齐”，后文自有分解。
主机端内存
主机端内存就是CPU的RAM。总体的使用方法和C语言申请内存的方法是一致的。但是CUDA将一种新形式的内存空间暴露给用户，那就是pinned内存。这种内存不会通过置换算法来swap到磁盘中，这段空间永远都保存在内存中。使用pinned内存可以很好地保证内存与显存的数据传输效率，并且pinned内存可以很好地映射到GPU的地址空间，让GPU可以直接使用主机端内存，从而减少大容量的数据拷贝。
pinned内存使用通过cudaHostAlloc()和cudaFreeHost()来分配和释放。通过在cudaHostAlloc()加入不同的参数来获得pinned memory的更多功能。首先是cudaHostAllocPortable参数，这个参数会让申请的pinned内存在不同的CPU线程之间共享，一般pinned内存是一个CPU线程独享的，通过这个参数可以在多个CPU线程之间共享数组。使用cudaHostAllocWriteCombined参数可以申请write-combined内存，这种内存在CPU端没有不会进行缓存，这样子CPU就不需要对这段内存进行一致性扫描，从而进一步提高CPU和GPU之间的传输效率。
GPU还支持零拷贝内存功能，通过cudaHostAlloc()函数，我们可以得到一个分配空间的主机端地址，然后我们使用cudaSetDeviceFlags()来进行内存映射，最后我们使用cudaHostGetDevicePointer()函数来获得这个主机端空间的设备端地址，这样子我们就不一定要将数据拷贝到设备端再进行计算了。
常数存储器
这是一种只读存储器，对于同一个线程束的不同线程访问同一个常数存储器的数据时，只要一个周期，这就是一个读优化的存储器。常数存储器在数组和变量前面加一个__constant__，就可以将这个数组和变量放到常数存储器中了，我们可以直接初始化，也可以使用cudaMemcpyToSymbol函数利用拷贝初始化。
1.4、CUDA通信
线程束
线程束是自带同步的，因为在GPU中，一个线程束中所有线程运行的指令是完全相同的。
__syncthreads()
__syncthreads()函数可以实现一个线程块内的线程同步。经常和共享内存的操作配套使用，来保证共享内存的一致性。
memory fence函数
这个函数主要有两个，用来保证线程间数据通信的可靠性。__threadfence()是保证这个线程在这个语句之前的所有对于全局内存和共享内存的访问已经完成，执行结果对于所有线程都是可见的。__threadfence_block()也是做相同的保证，但是执行结果对于一个block中的所有线程可见。
我们为什么要使用memory fence函数呢？那是因为无论是CPU还是GPU，处理器的指令都是乱序发射的。如果几行代码之间没有明确的先后顺序，可能就会被打乱执行，这在多线程中是非常致命的，因为这很可能会导致明明共享资源没有被更新，但是因为其他的语句已经运行，根据线程之间的逻辑激活了本来不应该在资源没有更新的时候激活的逻辑，这就导致了一定的错误。
关于memory fence，这里有一个讲的比较清楚的博客：Memory Fence。memory fence在CPU和GPU中都是存在的，都是解决在多线程环境下指令乱序发射导致的问题。
对于共享资源的更新是不是成功进行需要通过其他的线程来确认，__threadfence()需要等待整个网格中所有线程的确认，而__threadfence_block()只需要一个线程块的所有线程确认。之所以有这两种形式的确认，我认为这个和GPU的缓存结构有关系，因为一个SM中所有的SP都共享一级缓存，所以有时候我们并不需要将数据写会全局内存，修改的结果就可以在一个线程块中的所有线程可见，但是对于其他块的线程就不是这样了。所以在GPU中有两个层次的memory fence操作。
CPU与GPU之间的同步
有时候我们会调用一些一些异步的接口来操作GPU，这会使得一些操作GPU的函数立即返回，CPU继续执行剩下的代码，这就是需要加入一些机制来解决CPU与GPU之间的同步问题。
cudaThreadSynchronize()是一个类似的函数，CPU计算会在这行代码处卡住，然后等待之前所有的GPU操作全部结束。
cudaStreamSynchronize()也是同样的道理，CPU会在这一行卡住，等待某一个流的操作全部结束。
cudaEventSynchronize()会等待某一个事件的发生。
Volatile
这是一个关键字，我们可以使用这个关键字来声明一个变量是“敏感的”，是其他线程很有可能修改的，防止某个线程因为直接从缓存读取而忽略了其他线程对于某个资源在全局上的修改。
ATOM操作
原子性操作是完全互斥的。防止多个线程对于同一个资源的修改出现问题。
1.5、异步执行
1.5.1、流
流有点像GPU的“多线程”。对于GPU来说，同一个流的所有操作是串行的，但是不同流的操作之间的并行的。我们可以在调用CUDA的异步API来激活GPU的异步操作，并且将这个异步操作分配给某一个流。这样子就可以实现CPU与GPU之间的并行。如果流有多个，那么就可以实现GPU不同任务之间的并行。
流的创建和使用是比较简单的，流的创建使用函数cudaStreamCreate()，流的销毁使用函数cudaStreamDestroy()。
流的使用也非常简单，调用Async为后缀的函数，然后将流作为参数传入就好了。在kernal函数的调用中<<<>>>>的第四个参数就是流，可以将内核函数分配给一个流进行调度。
1.5.2、事件
事件是和流紧密配合的东西，可以用来进行CPU与GPU之间的同步。流使用接口cudaEventCreate()来创建，使用cudaStreamDestroy()来销毁，使用cudaEventRecord()来打时间戳，这个函数可以传入时间和流，这个事件会在当前事件和流的cudaEventRecord()之前的对应流的所有的异步GPU操作全部完成后激活，我们不仅可以记录两个时间戳之间的时间，也可以使用cudaEventSynchronize()来利用事件进行CPU与GPU之间的同步。
我们将会使用一个例子来进行CPU与GPU之间的异步计算，并且进行适当地同步工作。我们将利用GPU计算两个数组的大小，分别使用流来进行并行计算。在每一个block中我们将利用共享内存进行规约求和，不同block之间我们将使用全局内存进行并行串行求和。在这个过程中，不仅会用到共享内存来进行规约求和，也会用到共享变量来记录当前已经完成的求和情况，方便最后对每一个block结果的串行归约。我们需要根据输出数组的情况开辟合适的线程网格来解决问题。并且使用事件来进行时间的记录。这个例子将会利用到几乎所有的CUDA的基本内容。
因为使用规约求和这种算法，这个博客讲述了规约求和的算法思想：理解cuda并行程序的规约思想
当然想要一步到位还是很难的，我们首先先编写了一个单流的GPU规约求和。
1.6、例子：CUDA规约求和
源代码放在了我的码云上：源码链接
/**
我们将会使用一个例子来进行CPU与GPU之间的异步计算，并且进行适当地同步工作。我们将利用GPU计算两个数组的大小，
分别使用流来进行并行计算。在每一个block中我们将利用共享内存进行规约求和，不同block之间我们将使用全局内存进行并行串行求和。
在这个过程中，不仅会用到共享内存来进行规约求和，也会用到共享变量来记录当前已经完成的求和情况，方便最后对每一个block结果的串行归约。
我们需要根据输出数组的情况开辟合适的线程网格来解决问题。并且使用事件来进行时间的记录。这个例子将会利用到几乎所有的CUDA的基本内容。
**/
#include
#include
#include
//需要求和的数组大小
#define N 5000
//用来记录计算完成情况的变量
__device__ int count = 0;
//判断是不是可以将过程量相加的变量
// __device__ int judge = 0;
//首先声明内核函数，这个函数处理的就是进行一个数组中所有元素的相加操作
//传入参数就是当前数组的指针以及当前数组的大小。
//并且传入一个数组result来保存每个block计算的中间结果，将最后的计算结果存在这个中间结果数组的第0号位置。
//中间结果数组的大小为gridDim.x
__global__ void arr_sum(int* input_arr, int size, long* result){
// printf("进入内核函数\n");
//这里进入内核函数
//首先先进行规约求和，首先先申请共享内存，一共192位
//注意关键字是shared而不是share
__shared__ long temp_arr[192];
//计算总的线程编号，来索引全局内存
unsigned int thread_idx = threadIdx.x + blockIdx.x * blockDim.x;
// //进行数据拷贝，因为数组的大小可能不是对齐的，加了判断防止数组越界
if(thread_idx < size){
temp_arr[threadIdx.x] = input_arr[thread_idx];
}else {
//因为输出的数组可能
temp_arr[threadIdx.x] = 0;
}
// printf("共享内存拷贝完毕\n");
// 拷贝结束之后必须进行同步，保证所有warp的操作全部完成
__syncthreads();
//这里进行规约求和，规约求和的特点就是每个线程只执行一个计算
//所以只有数据规模一半的线程是活跃的
int active = 192 / 2;
//这里开始进行规约求和
//我们需要一个循环，让工作线程的边界不断缩减
int i;
//当规约到奇数个项的时候就难以规约了，这个时候使用暴力相加的方式就好
//当active的边界是3的时候，还有6个元素，还可以规约
//所有i的边界包括3
for(i = active; i >= 3; i = i / 2){
//如果当前线程是活跃线程就参与计算
if(threadIdx.x < i){
//开始计算两个点的相加
temp_arr[threadIdx.x] = temp_arr[threadIdx.x] + temp_arr[threadIdx.x + i];
//每次算完都必须同步
//注意不要在分支语句中加同步，因为有些线程进不来
// __syncthreads();
}
//不要在分支中同步
__syncthreads();
}
// // printf("规约求和完毕\n");
int judge = 0;
if(threadIdx.x == 0){
//然后使用每个block的第一个线程来做最后三个元素的相加
// printf("%ld,%ld,%ld\n", temp_arr[0], temp_arr[1], temp_arr[2]);
result[blockIdx.x] = temp_arr[0] + temp_arr[1] + temp_arr[2];
// printf("%ld\n", result[blockIdx.x]);
//这里强制等待写回
__threadfence();
//这里记录一个变量，每个块完成计算之后，就为这个值自增
atomicAdd(&count,1);
if(count == gridDim.x){
judge = 1;
}
}
//这套实现是错的，不能让一个不知道何时开始运行的线程等待那么久
//这会导致程序的崩溃。
// if(thread_idx == 0){
// //这里如果是整个网格的第一个线程，那就完成这个result数组的求和操作
// while(count != gridDim.x){
// }
// for(i = 1; i < gridDim.x; i++){
// printf("%ld\n", result[i]);
// result[0] = result[0] + result[i];
// }
// __threadfence();
// // printf("%ld\n", result[0]);
// }
//这里代表最后一个块，最后一个块进行收尾
if(judge == 1){
if(threadIdx.x == 0){
for(i = 1; i < gridDim.x; i++){
result[0] = result[0] + result[i];
}
}
}
}
//主函数
int main(){
//这里计算网格结构，我们使用的是Tesla K40，15SM，192SP
//我们建立一维网格，计算需要几个block
int block_num = (N / 192);
//如果数组的大小没有办法与192取整，那么就需要一个新的块
if(N % 192 != 0){
block_num++;
}
//创建两个数组为了方便验证结果，分别是正序的和倒序的等差数列。数组的大小为N
//使用pinned memory
int* arr1;
long* result1;
cudaMallocHost((void **)&arr1, N * sizeof(int), cudaHostAllocWriteCombined);
//用来存储结果的数组
cudaMallocHost((void **)&result1, block_num * sizeof(long));
//用来迭代的变量
int i;
//将两个数组分别初始化
for(i = 0; i < N; i++){
arr1[i] = i;
// arr2[i] = N - i - 1;
}
printf("主机数据初始化完毕\n");
//申请GPU空间
int* device_input_arr1 = NULL;
long* device_result1 = NULL;
//使用cudaMalloc来申请空间
//函数原型__host____device__cudaError_t cudaMalloc (void** devPtr, size_t size)
//注意这里要传入的是指针的地址
cudaMalloc((void **)&device_input_arr1, N*sizeof(int));
cudaMalloc((void **)&device_result1, block_num*sizeof(long));
printf("设备空间分配完毕\n");
//完成数据拷贝，从这里开始使用流
cudaMemcpy((void *)device_input_arr1, (void *)arr1, N*sizeof(int), cudaMemcpyHostToDevice);
//初始化中间结果数组
cudaMemset((void *)device_result1, 0, block_num * sizeof(long));
printf("激活初始化\n");
//网格建立，一维网格
dim3 grid_dimension(block_num);
//每个block拥有和sp数量一致的thread，也是一维的
dim3 block_dimension(192);
cudaError_t error_check;
//使用内核函数，不分配动态共享内存
arr_sum<<>>(device_input_arr1, N, device_result1);
error_check = cudaGetLastError();
if( error_check != cudaSuccess ){
printf("%s\n" , cudaGetErrorString( error_check ) );
system("pause") ;
return 0 ;
}
//将数据拷贝出来
cudaMemcpy((void *)result1, (void *)device_result1, block_num * sizeof(long), cudaMemcpyDeviceToHost);
printf("激活数据拷出\n");
//检查数据是不是正确
long real_result = (0 + N - 1) * N / 2;
printf("结果%ld\n", result1[0]);
if(result1[0] == real_result){
printf("结果正确\n");
}else{
printf("结果不正确\n");
}
//析构pinned memory
cudaFreeHost(arr1);
cudaFreeHost(result1);
//析构在设备端分配的函数
cudaFree(device_input_arr1);
cudaFree(device_result1);
return 0;
}
在串行代码中我们创造一个等差数列，然后将数据传入内核函数中。并且我们申请了一个数组result来保存计算的中间结果。每个线程块会产生一个计算结果，归并排序的过程就不赘述了。
我们在这段代码的编写过程中遇到了很多问题，GPU的代码确实非常难以调试，下面我们做一些总结。
__syncthreads()的使用雷区
这个函数在使用的时候要注意，在分支语句中使用要格外注意，这是一个block内部所有线程的同步的函数。但是如果一个block钟的不同线程在不同分支的话，在不同分支的__syncthreads()就极有可能造成程序崩溃。所以__syncthreads()一定要放在所有线程都经过的路径上。
比如这次遇到的错误：
for(i = active; i >= 3; i = i / 2){
//如果当前线程是活跃线程就参与计算
if(threadIdx.x < i){
//开始计算两个点的相加
temp_arr[threadIdx.x] = temp_arr[threadIdx.x] + temp_arr[threadIdx.x + i];
//每次算完都必须同步
//注意不要在分支语句中加同步，因为有些线程进不来
// __syncthreads();
}
//不要在分支中同步
__syncthreads();
}
本来在__syncthreads()是放在if语句中的，结果出错了。
禁止让单个线程进入死循环的等待
这个错误主要来自于这段已经被注释掉的错误代码：
//这套实现是错的，不能让一个不知道何时开始运行的线程等待那么久
//这会导致程序的崩溃。
if(thread_idx == 0){
//这里如果是整个网格的第一个线程，那就完成这个result数组的求和操作
while(count != gridDim.x){
}
for(i = 1; i < gridDim.x; i++){
printf("%ld\n", result[i]);
result[0] = result[0] + result[i];
}
__threadfence();
// printf("%ld\n", result[0]);
}
我们让整个网格的第一个线程进入死循环等待所有块的中间结果全部计算完毕。但是这种方式的错误的，让一个线程陷入这样的等待会导致一个线程块长时间占用一个SM，导致显卡崩溃。正确的方式是我们找出最后一个运行的线程块，然后让最后一个线程块的某一个线程来完成所有中间结果的规约操作。这样子就避免的某一个线程长时间的等待。
1.7、例子：多流的规约求和
/**
我们将会使用一个例子来进行CPU与GPU之间的异步计算，并且进行适当地同步工作。我们将利用GPU计算两个数组的大小，
分别使用流来进行并行计算。在每一个block中我们将利用共享内存进行规约求和，不同block之间我们将使用全局内存进行并行串行求和。
在这个过程中，不仅会用到共享内存来进行规约求和，也会用到共享变量来记录当前已经完成的求和情况，方便最后对每一个block结果的串行归约。
我们需要根据输出数组的情况开辟合适的线程网格来解决问题。并且使用事件来进行时间的记录。这个例子将会利用到几乎所有的CUDA的基本内容。
**/
#include
#include
#include
//需要求和的数组大小
#define N 1000
//用来记录计算完成情况的变量
__device__ int count1 = 0;
__device__ int count2 = 0;
//判断是不是可以将过程量相加的变量
// __device__ int judge = 0;
//首先声明内核函数，这个函数处理的就是进行一个数组中所有元素的相加操作
//传入参数就是当前数组的指针以及当前数组的大小。
//并且传入一个数组result来保存每个block计算的中间结果，将最后的计算结果存在这个中间结果数组的第0号位置。
//中间结果数组的大小为gridDim.x
//第四个参数传入流的编号
__global__ void arr_sum(int* input_arr, int size, long* result, int stream_id){
// printf("进入内核函数\n");
//这里进入内核函数
//首先先进行规约求和，首先先申请共享内存，一共192位
//注意关键字是shared而不是share
__shared__ long temp_arr[192];
//计算总的线程编号，来索引全局内存
unsigned int thread_idx = threadIdx.x + blockIdx.x * blockDim.x;
// //进行数据拷贝，因为数组的大小可能不是对齐的，加了判断防止数组越界
if(thread_idx < size){
temp_arr[threadIdx.x] = input_arr[thread_idx];
}else {
//因为输出的数组可能
temp_arr[threadIdx.x] = 0;
}
// printf("共享内存拷贝完毕\n");
// 拷贝结束之后必须进行同步，保证所有warp的操作全部完成
__syncthreads();
//这里进行规约求和，规约求和的特点就是每个线程只执行一个计算
//所以只有数据规模一半的线程是活跃的
int active = 192 / 2;
//这里开始进行规约求和
//我们需要一个循环，让工作线程的边界不断缩减
int i;
//当规约到奇数个项的时候就难以规约了，这个时候使用暴力相加的方式就好
//当active的边界是3的时候，还有6个元素，还可以规约
//所有i的边界包括3
for(i = active; i >= 3; i = i / 2){
//如果当前线程是活跃线程就参与计算
if(threadIdx.x < i){
//开始计算两个点的相加
temp_arr[threadIdx.x] = temp_arr[threadIdx.x] + temp_arr[threadIdx.x + i];
}
//不要在分支中同步
__syncthreads();
}
// // printf("规约求和完毕\n");
int judge = 0;
if(threadIdx.x == 0){
//然后使用每个block的第一个线程来做最后三个元素的相加
// printf("%ld,%ld,%ld\n", temp_arr[0], temp_arr[1], temp_arr[2]);
result[blockIdx.x] = temp_arr[0] + temp_arr[1] + temp_arr[2];
// printf("%ld\n", result[blockIdx.x]);
//这里强制等待写回
__threadfence();
//这里记录一个变量，每个块完成计算之后，就为这个值自增
if(stream_id == 1){
atomicAdd(&count1,1);
if(count1 == gridDim.x){
judge = 1;
}
}else{
atomicAdd(&count2,1);
if(count2 == gridDim.x){
judge = 1;
}
}
}
//这里代表最后一个块，最后一个块进行收尾
if(judge == 1){
if(threadIdx.x == 0){
for(i = 1; i < gridDim.x; i++){
// printf("%ld\n", result[i]);
result[0] = result[0] + result[i];
}
}
}
}
//主函数
int main(){
//这里计算网格结构，我们使用的是Tesla K40，15SM，192SP
//我们建立一维网格，计算需要几个block
int block_num = (N / 192);
//如果数组的大小没有办法与192取整，那么就需要一个新的块
if(N % 192 != 0){
block_num++;
}
//创建两个数组为了方便验证结果，分别是正序的和倒序的等差数列。数组的大小为N
//使用pinned memory
int* arr1;
int* arr2;
long* result1;
long* result2;
cudaMallocHost((void **)&arr1, N * sizeof(int), cudaHostAllocWriteCombined);
cudaMallocHost((void **)&arr2, N * sizeof(int), cudaHostAllocWriteCombined);
//用来存储结果的数组
cudaMallocHost((void **)&result1, block_num * sizeof(long));
cudaMallocHost((void **)&result2, block_num * sizeof(long));
//用来迭代的变量
int i;
//将两个数组分别初始化
for(i = 0; i < N; i++){
arr1[i] = i;
arr2[i] = N - i - 1;
}
printf("主机数据初始化完毕\n");
//申请GPU空间
int* device_input_arr1 = NULL;
int* device_input_arr2 = NULL;
long* device_result1 = NULL;
long* device_result2 = NULL;
//使用cudaMalloc来申请空间
//函数原型__host____device__cudaError_t cudaMalloc (void** devPtr, size_t size)
//注意这里要传入的是指针的地址
cudaMalloc((void **)&device_input_arr1, N*sizeof(int));
cudaMalloc((void **)&device_input_arr2, N*sizeof(int));
cudaMalloc((void **)&device_result1, block_num*sizeof(long));
cudaMalloc((void **)&device_result2, block_num*sizeof(long));
printf("设备空间分配完毕\n");
//创建两个流
cudaStream_t stream[2];
for(i = 0; i < 2; i++){
cudaStreamCreate(&stream[i]);
}
//完成数据拷贝，从这里开始使用流
cudaMemcpyAsync((void *)device_input_arr1, (void *)arr1, N*sizeof(int), cudaMemcpyHostToDevice, stream[0]);
cudaMemcpyAsync((void *)device_input_arr2, (void *)arr2, N*sizeof(int), cudaMemcpyHostToDevice, stream[1]);
//初始化中间结果数组
cudaMemsetAsync((void *)device_result1, 0, block_num * sizeof(long), stream[0]);
cudaMemsetAsync((void *)device_result2, 0, block_num * sizeof(long), stream[1]);
printf("激活初始化\n");
//网格建立，一维网格
dim3 grid_dimension(block_num);
//每个block拥有和sp数量一致的thread，也是一维的
dim3 block_dimension(192);
cudaError_t error_check;
//使用内核函数，不分配动态共享内存
arr_sum<<>>(device_input_arr1, N, device_result1, 1);
arr_sum<<>>(device_input_arr2, N, device_result2, 2);
error_check = cudaGetLastError();
if( error_check != cudaSuccess ){
printf("%s\n" , cudaGetErrorString( error_check ) );
system("pause") ;
return 0 ;
}
printf("激活内核\n");
//将数据拷贝出来
cudaMemcpyAsync((void *)result1, (void *)device_result1, block_num * sizeof(long), cudaMemcpyDeviceToHost, stream[0]);
cudaMemcpyAsync((void *)result2, (void *)device_result2, block_num * sizeof(long), cudaMemcpyDeviceToHost, stream[1]);
printf("激活数据拷出\n");
printf("进入等待\n");
//这里等待所有的GPU代码执行完毕
cudaThreadSynchronize();
//检查数据是不是正确
long real_result = (0 + N - 1) * N / 2;
printf("结果%ld\n", result2[0]);
if(result1[0] == real_result){
printf("第一个结果正确\n");
}else{
printf("第一个结果不正确\n");
}
if(result2[0] == real_result){
printf("第二个结果正确\n");
}else{
printf("第二个结果不正确\n");
}
//析构两个流
for(i = 0; i < 2; i++){
//注意这个函数的使用，流的创建输入的是
cudaStreamDestroy(stream[i]);
}
//析构pinned memory
cudaFreeHost(arr1);
cudaFreeHost(arr2);
cudaFreeHost(result1);
cudaFreeHost(result2);
//析构在设备端分配的函数
cudaFree(device_input_arr1);
cudaFree(device_input_arr2);
cudaFree(device_result1);
cudaFree(device_result2);
return 0;
}
多流的使用方法也是比较简单的，我们只需要稍微修改一下原有的代码，在内核函数中加入一个变量来标示当前流的编号。并且使用同步机制来等待两个流计算完毕。
2、显卡的硬件结构
2.1、显卡的基本结构
下面的一张图展现了显卡的硬件构造：
在编程当中，我们通常只能看到SM与SP两个硬件的层次，但是，但实际上显卡的硬件层次比我们编程的时候所看见的要复杂。
下面这张图代表了G80架构的显卡构造，虽然已经非常过时，但是我们还是可以看到一些端倪：
我们可以看到在SM之上还有一个叫做TPC的结构。而TPC又包含在一个叫做SPA的结构当中。SPA被称作流处理器阵列，TPC被称作线程处理器群，而SM被称作流多处理器。在GPU中所有的计算单元都在SPA中，而SPA中又有多层的结构：
TPC中有多个SM，一个TPC中的所有SM共享指令与常量的二级缓存：
在一个SM中拥有指令和常量的以及缓存以及共享内存，每个SP都有一个对应的寄存器文件。
2.2、warp（线程束）与共享内存
线程束是SM指令发射的基本单位，在现有的GPU中线程束包含32个线程，这32个线程会保证指令集的同步，也就是说这32个线程是在同一时刻执行的指令是一模一样的。要明白GPU的指令发射方法，我们需要看明白下面这张图：
指令在SM中的各个元器件中是按照swap为单位来传递的。Entry Warp Instruction Buffer，存储着即将进行计算的一系列指令，这些指令会按照一定的优先级重新排序。SM的发射逻辑会找到优先级最高的指令进行计算。
warp的意义除了指令发射之外，还体现在SM内的共享内存的使用上。共享内存里面有一种叫做bank的结构，不同bank之间是可以并发运作的，但是相同bank之间只能串行读写数据。下面的示意图展现了这个过程，这是一个共享内存，每一列就是一个bank，地址根据bank的顺序顺序分布。所以我们要保证一个线程束的不同线程在共享内存中一次读取的位置正好在不同的bank中（0-31正好顺序排在不同的bank中）。这样子一个线程束中所有线程的可以在同时执行一条指令的时候一口气取出所有需要的数据了。
2.3、全局内存（显存）
共享内存是最慢的存储介质，但是实际上它也有可以提升性能的方法。因为GPU的显存很大，所以会由多个内存控制器来分别控制，每个内存控制器管辖的空间是连续的，叫做分区，而不同分区之间也是连续的。
上面这张图就展现了一个数组的分区，我们可以看到分区是取模分布的，分区比较小，但是用户申请的空间可能很大。所以一个空间中不同分区会连续地、取模地分布。而如果我们让我们要同时取的数据正好在不同的分区，这样子就可以极大地促成IO的并行化，从而提高IO效率。这个问题在GPU中被称作“分区冲突”问题。
3、CUDA的优化
一般来讲第一次编写的CUDA程序是难以发挥硬件的完整性能的，我们需要根据硬件的特点来对CUDA程序做优化。
3.1、测量程序的运行时间
在CUDA程序中我们可以在设备端也可以在主机端为设备测试运行时间。因为核函数只能规定一个线程的运行内容，所以在核函数中记录运行时间，我们需要在第一个运行的warp中记录开始时间，然后在最后一个运行的warp中记录结束时间。在内核函数中clock函数来进行计时。
我们运行一个程序来进行计时的演示。这个程序非常简单，我们为每一个线程进行一个等差数列的求和操作，然后通过clock函数来计算时间。使用clock函数来进行计时，这个函数的返回值不是时间，而是每个线程在执行代码用了多少个时钟周期。所以时钟周期还要除以GPU频率才能算出真正的时间。这里有一个博客介绍得比较详细：CUDA之clock()方法详解。
这里我编写了一个程序来演示了一下在CUDA核函数内进行计时工作的过程，在核函数中我们申请了共享内存，并且对全局内存进行存取，来尽可能模仿一个正常的核函数操作。
此外我们还使用了另外一个基于事件的CUDA接口，cudaEventElapsedTime来进行时间的记录。对于事件的利用可以查看这个博客：GPU编程系列之三】cuda stream和event相关内容。
下面是这段计时函数的源码，在我们码云上有最新的版本（https://gitee.com/ZhenDu_ICT/CUDA_practice/blob/master/CUDA程序计时/clock_timer.cu）：
//这个程序我们打算在内核函数中记录时间戳，然后在CPU逻辑中查看第一个线程和最后一个线程运行的时间。
#include
#include
#include
//就传入两个数组，一个数组记录每一个block开始的时间，一个数组记录每一个block结束的时间
//使用clock函数来进行计时，这个函数的返回值不是时间，而是每个线程在执行代码用了多少个时钟周期。
//所以时钟周期还要除以GPU频率才能算出真正的时间
//https://blog.csdn.net/adi_1987/article/details/53465370
//为了方式编译器的优化，我们才里面申请一段共享内存，并且将共享内存的数据求和，并且将结果放到一个数组中
__global__ void gpu_function(clock_t* start, clock_t* stop, int* input, int* result_arr){
//每个线程块的第一个线程的记录线程
if(threadIdx.x == 0){
start[blockIdx.x] = clock();
}
//利用共享内存将在全局内存的输入数据取入
__shared__ int temp[192];
temp[threadIdx.x] = input[threadIdx.x];
__syncthreads();
//这里执行计算
int i = 0;
int result = 0;
for(i = 0; i < 192; i++){
result = result + temp[i];
}
//将数据放到全局内存中
if(threadIdx.x == 2){
result_arr[blockIdx.x] = result;
}
//这里执行同步
__syncthreads();
//这里记录结束的时间戳
if(threadIdx.x == 0){
stop[blockIdx.x] = clock();
}
}
//我们的显卡一共有15个SM，每个SM有192个sp，我们打算运行30个线程块，每个线程块开启192个线程。
int main(){
//这里激活内核函数
dim3 grid_dimention(30);
dim3 block_dimention(192);
clock_t* start = NULL;
clock_t* stop = NULL;
cudaMallocHost((void **)&start, 30*sizeof(clock_t));
cudaMallocHost((void **)&stop, 30*sizeof(clock_t));
//这里申请两个记录时间戳的数组来
clock_t* device_start = NULL;
clock_t* device_stop = NULL;
cudaMalloc((void **)&device_start, 30*sizeof(clock_t));
cudaMalloc((void **)&device_start, 30*sizeof(clock_t));
//输入数组指针
int* input;
//这里申请一个数组，这个数组里面有0-191个数字
cudaMallocHost((void **)&input, 192 * sizeof(int));
int i;
for(i = 0; i < 192; i++){
input[i] = i;
}
cudaEvent_t start_event;
cudaEvent_t stop_event;
cudaEventCreate(&start_event);
cudaEventCreate(&stop_event);
cudaEventRecord(start_event, 0);
//这里申请输出空间
int *device_result;
cudaMalloc((void **)&device_result, 30 * sizeof(int));
int *result;
cudaMallocHost((void **)&result, 30 * sizeof(int));
//将输入数组拷贝到显卡上
int* device_input;
cudaMalloc((void **)&device_input, 192 * sizeof(int));
cudaMemcpy(device_input, input, 192 * sizeof(int), cudaMemcpyHostToDevice);
gpu_function<<>>(device_start, device_start, device_input, device_result);
cudaThreadSynchronize();
//这里将数据拷出
cudaMemcpy(start, device_start, 30*sizeof(clock_t), cudaMemcpyDeviceToHost);
cudaMemcpy(stop, device_stop, 30*sizeof(clock_t), cudaMemcpyDeviceToHost);
cudaMemcpy(result, device_result, 30*sizeof(int), cudaMemcpyDeviceToHost);
cudaEventRecord(stop_event, 0);
cudaEventSynchronize(stop_event);
//记录时间
float time;
cudaEventElapsedTime(&time, start_event, stop_event);
printf("基于事件的计时方式：%f\n", time);
//这里获取GPU的主频，使用的方法就是cudaDeviceProp
cudaDeviceProp prop;
cudaGetDeviceProperties(&prop, 0);
int gpu_rate = prop.clockRate;
//查找最早的开始时间和最后结束的时间
clock_t min_clock = start[0];
clock_t max_clock = stop[0];
//选出第一个block的开始时间和最后一个block的结束时间
for(i = 1; i < 30; i++){
if(min_clock < start[i]){
min_clock = start[i];
}
if(max_clock > stop[i]){
max_clock = stop[i];
}
}
printf("%ld\n", (max_clock-min_clock)/gpu_rate);
cudaFreeHost(start);
cudaFreeHost(stop);
cudaFreeHost(result);
cudaFreeHost(input);
cudaEventDestroy(start_event);
cudaEventDestroy(stop_event);
cudaFree(device_start);
cudaFree(device_stop);
cudaFree(device_input);
cudaFree(device_result);
return 0;
}
这里我们采用了两套思路来计时，但是不是为什么在CUDA内部的clock的计时方式总是得出负数的结果，现在暂时还是没有发现问题的原因。而基于事件的方式得到的时间应该是正确的。在基于clock()的方法中，我们让每个块的第一个线程都记录这个块的开始和结束时间，通过第一个块的开始时间和最后一个块结束时间来得到GPU的运行时间。而基于事件的方法中，我们首先为开始和结束的时间打上锚点，cudaEventRecord()函数的参数分别是流和创建好的时间，而这个事件会在这个函数之前对应流的所有cuda操作完成后激活。而cudaEventSynchronize()，会卡死当前cpu线程，直到这个函数参数中的事件已经完成，这个函数的意义就是让stop_event事件完成之后在进行事件的计算，一般来讲，核函数的调用通常都是异步的，而数据的拷贝炒作都是同步的。
我们也可以在主机端来利用计时函数来进行GPU计算的计时。但是同时也要注意同步问题。在CUDA中我们即便没有调用sync相关的接口，与CUDA相关的函数通常都会被放到“0号流”中，这是一个缺省流，当用户不进行异步操作的时候，这些CUDA函数就会被发射到0号流中。0号流除了完成分配给自己的操作，还可以用来做全局的同步。当0号流和其他流同时存在的时候，0号流会在所有其他流还是之前开始，在其他流结束之后结束。所以在CUDA中，0号流、或者说缺省流在同步上的意义非比寻常。
3.2、存储器的访问优化
3.2.1、全局存储器的优化
全局存储器是就是显卡的显存，我们在全局存储器上有两点可以做访存优化，一个是负载均衡，一个是合并访存。在GPU中通常有多个内存控制器，每个内存控制器掌管显存的一部分，所以每个线程块全部集中访问全局内存的一个区域比均匀访问显存的各个区域会带来更大的性能劣化。所以我们在程序设计的时候要让不同的线程在显存在均匀访问。
此外我们可以根据合并访问来提升全局内存效能，在一个线程束的线程如果访问的位置在全局内存上相近，那么这个线程束需要的所有全局内存上的数据就可以仅通过一次访问得到。
大体的示意图如下：
图的上半部分就是一个个全局内存空间，而下半部分就是一个线程束中的一个个线程，因为这张图来自比较早的文献，所以一共是只有半个线程束的大小。但是随着硬件的发展，我们已经不再讲half-warp的概念了。无论是指令级同步的相邻线程数量，还是共享内存的bank conflict避免，还是全局内存合并读取，都是一个warp，也就是32个相邻线程。
实际上合并访问的触发条件在现在的设备中已经被大大放宽了，合并访问一要保证段对齐，也就是针对字长为8-128bit（对应的段大小就是，32个字节到128个字节）的访问粒度，线程束访问的最低地址要是这个段的首地址。那么剩下的线程针对这个段的访问就是合并的，可以不一一对应，也不用访问这个段的全部内容。即便地址没有对齐也没有关系，实际上也会按照128个字节为段来进行合并访问。所以说，我们有理由认为，只要同一个warp访问的范围没有超过128个字节，通常来讲都是可以进行合并访问的。
3.2.2、共享存储器的优化
共享存储器是全局内存的cache。我们需要防止一种叫做bank conflict的现象出来，以此来提升共享内存访问的性能。下面的示意图展现了bank conflict的由来。共享内存大概是这样的一个结构：
每一个绿色的方框代表一个32bit大小的存储单元，而一系列的存储单元会被串起来，共享32bit的位宽。被串起来的一组存储单元就被叫做一个bank。对于一个SM来说，一共有32个bank，这和线程束的大小的匹配的。并且在不同的bank中共享内存的地址是取模分布的，也就是说0-31号bank的地址分别是0-31，同时也是32-63；并且对于0号bank来说，它既存储了地址0的数据，也存储了地址32的数据。所以说在同一个warp的不同线程中，对于同一个bank的访问会导致两次串行的IO，从而降低共享内存的利用率。如果一个线程束中的不同线程访问的是不同bank的存储器，那么就可以带来完全并行的共享内存的访问。
当然还有一种情况是比较特殊的，那就是“广播字”。在共享内存中，单个的存储空间是可以进行广播的。也就是说如果不同线程对于同一个bank的访问是出现在同一个存储单元上（比如多个线程同时访问共享内存数组的32号地址，也就是0号bank），那么会可以就可以进行一次广播达到多次传输的效果，从而并不会产生bank conflict。
以上就是CUDA简明入门的所有内容。
扫一扫
1506
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
