CUDA 笔记2

Reference:

https://face2ai.com/program-blog/#GPU%E7%BC%96%E7%A8%8B%EF%BC%88CUDA%EF%BC%89

4 共享内存

4.1 共享内存概述

GPU的共享内存在片上，延迟低带宽高。常用语：

  • 块内线程通信的通道
  • 用于全局内存数据的缓存
  • 用于转换数据来优化全局内存访问模式

回顾GPU的存储结构：

5-1

共享内存的生命周期是整个线程块执行的过程。每个线程对于共享内存的访问如果是独立进行的，则效率是比较高的，如果会产生竞态，就要解决冲突问题。如果线程束内32个线程访问同一地址，由一个线程访问完后进行广播。

一个SM上的所有的正在执行的线程块共同使用物理的共享内存，所以共享内存也成为了活跃线程块的限制，块使用的共享内存越小，线程块级别的并行度就越高。

4.1.1 共享内存分配

声明共享内存可以在核函数内，也可以在核函数外。也可以动态声明。

1
__shared__ float a[size_x][size_y];

size_x和size_y需要在编译期确定，如果需要动态声明，需要使用extern关键字：

1
2
extern __shared__ int tile[];
kernel<<<grid,block,isize*sizeof(int)>>>(...);

动态声明只支持1维数组。

4.1.2 共享内存存储与访问模式

共享内存是一维地址空间。其具有特殊的形式是共享内存分为32个内存模型，称为存储体，对应线程束的32个线程，如果线程访问共享内存时访问不同的存储体，就可以在无冲突的情况下一个事务完成访问。

线程束访问共享内存有三个模式：

  • 并行访问：多地址访问多存储体
  • 串行访问：多地址访问同一存储体
  • 广播访问：单一地址读取单一存储体

并行访问可能会有冲突，如果有小部分冲突，可以将冲突的部分分隔执行。如果完全冲突，就变成了串行访问。而广播访问相较于并行访问的带宽利用率差，延迟接近。

理想并行访问：

不规则但不冲突访问：

可能产生冲突的访问：如果是同一地址就通过广播解决，否则产生冲突

一个存储体的宽度与设备计算能力有关，3.x以上的为8字节。也就是说如果32个线程束访问32个地址上的小于8字节的数据，是可以并行进行的。

4字节的存储体：

5-5

8字节的存储体：

5-6

无冲突的两种情况：

5-7

5-8

冲突的情况：两个线程访问了一个存储体的不同地址

5-9

存储体冲突会严重影响共享内存访问效率。有时可以用填充来降低冲突：

假设只有4个存储体

5-11

上面的存储方式，如果声明为shared int a[5][4];就会有冲突，但是如果加入一行：

5-11-1

只有四个存储体，要访问的数据就会错开：

5-12

以下语句可以查询存储体宽度：

1
2
3
4
5
6
cudaError_t cudaDeviceGetSharedMemConfig(cudaSharedMemConfig * pConfig);
//返回值
cudaSharedMemBankSizeFourByte
cudaSharedMemBankSizeEightByte
//配置存储体带宽：
cudaError_t cudaDeviceSetShareMemConfig(cudaSharedMemConfig config);

在不同核函数之间修改配置需要进行一次同步。存储体大小不会影响共享内存的使用，但会对性能可能有较大的影响。

4.1.3 配置共享内存

共享内存的配置有两个函数，一个是针对设备的，另一个是针对核函数的：

1
2
cudaError_t cudaDeviceSetCacheConfig(cudaFuncCache cacheConfig);
cudaError_t cudaFuncSetCacheConfig(const void* func,enum cudaFuncCacheca cheConfig);

config有几种：

1
2
3
4
cudaFuncCachePreferNone: no preference(default)
cudaFuncCachePreferShared: prefer 48KB shared memory and 16 KB L1 cache
cudaFuncCachePreferL1: prefer 48KB L1 cache and 16 KB shared memory
cudaFuncCachePreferEqual: prefer 32KB L1 cache and 32 KB shared memory

如果共享内存使用多，更多的共享内存更好。如果寄存器使用更多，L1Cache越多越好。两种存储的行为是不一样的。L1通过缓存行进行数据的访问和删除。而共享内存是完全编程控制的。

4.1.3 同步

CUDA采用宽松内存模型，即内存访问不是一定按照程序中的顺序执行的。因此必须进行同步。显示的同步是__sycthreads()调用。保证同步后所有的全局内存，共享内存都是所有线程可见的。

如果只是需要内存同步，可以用内存栅栏来同步。内存栅栏有三个维度：块，网格，系统。

1
2
3
void __threadfence_block();
void __threadfence();		//网格
void __threadfence_system();

这一节仅仅是一个概述，因此就到此结束了。

4.2 共享内存布局

本节重点要解决设计核函数时的两个问题：

  • 跨内存存储体映射数据元素
  • 从线程索引到共享内存偏移的映射

对于一个二维的共享内存：

1
2
3
4
#define N 32
...
__shared__ int x[N][N];
...

习惯的访问方式可能是这样的：

1
2
3
4
5
#define N 32
...
__shared__ int x[N][N];
...
int a=x[threadIdx.y][threadIdx.x];

这样的方式实际上也是最高效的，因为一个线程束中的线程是线程块中按行优先划分的，即线程束中的threadIdx.x是连续的。因此这样的访问方式，对应访问的也是共享内存的一行。并且每个线程访问一个存储体：(红色的是以上这种方式的访问模式)

5-12

因此数据要优先行主序访问，和线程束的线程对应。

4.2.1 Example1

用一个简单的程序读写共享内存，这个程序的核函数只做两个操作，将线程索引值写入二维共享内存，再按照行主序读取值存入全局内存：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include<cuda_runtime.h>
const int BLOCKDIM = 32*32;

__global__ void func(int *out){
    __shared__ int tile[32][32];
    int idx = threadIdx.y*blockDim.x + threadIdx.x;
    tile[threadIdx.y][threadIdx.x] = idx;
    __syncthreads();
    out[idx] = tile[threadIdx.y][threadIdx.x];
}

int main(){
    int *glob_h, *glob_d;
    int msize = sizeof(int)*BLOCKDIM*BLOCKDIM;
    glob_h = (int *)malloc(msize);
    cudaMalloc(&glob_d, msize);
    dim3 block(32,32);
    dim3 grid((BLOCKDIM-1)/block.x+1, (BLOCKDIM-1)/block.y+1);
    func<<<grid,block>>>(glob_d);
    cudaMemcpy(glob_h, glob_d, msize,cudaMemcpyDeviceToHost);
    cudaFree(glob_d);
    free(glob_h);
    return 0;
}

image-20241029210647380

如果按列主序访问，则：

image-20241029210818538

所有的访问都会触发bank conflict。此时的核函数运行时间是行主序的6.3倍。

还可以尝试动态分配共享内存，和静态分配只是方法上不一样。

1
2
3
4
5
6
7
8
__global__ void func(int *out){
    extern __shared__ int tile[];
    ...;
}
int main(){
    ...;
    func<<<grid,block,32*32*sizeof(int)>>>(glob_d);
}

上文已经提到过了，可以通过填充使需要引用的数据交错来消除访问冲突，只需要将第二维度添加一个填充，就可以观察到bank conflicts完全消失。动态分配的共享内存填充时需要注意索引，因为数据是一维存储的。

在上面的例子，数据是正方形且和block大小一致，如果是矩形，并且读和写的时候主序不一致。就要考虑索引问题，需要先将索引转换成线性，再重新计算行和列的坐标。这里跳过了这部分内容，待有需要时再仔细考虑。

4.3 减少全局内存访问

在CUDA执行模型一章节已经介绍过避免分支分化和循环展开等方式的核函数优化方式。本节从减少全局内存访问的角度优化核函数：

  • 重新安排数据访问模式避免线程束分化
  • 展开循环以保证有足够的操作使指令和内存带宽饱和

这一节使用了之前的完全展开的规约计算为例。完全展开的规约计算：

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
__global__ void reduceGmem(int * g_idata,int * g_odata,unsigned int n)
{
	//set thread ID
	unsigned int tid = threadIdx.x;
	unsigned int idx = blockDim.x*blockIdx.x+threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the
	int *idata = g_idata + blockIdx.x*blockDim.x;

	//in-place reduction in global memory
	if(blockDim.x>=1024 && tid <512)
		idata[tid]+=idata[tid+512];
	__syncthreads();
	if(blockDim.x>=512 && tid <256)
		idata[tid]+=idata[tid+256];
	__syncthreads();
	if(blockDim.x>=256 && tid <128)
		idata[tid]+=idata[tid+128];
	__syncthreads();
	if(blockDim.x>=128 && tid <64)
		idata[tid]+=idata[tid+64];
	__syncthreads();
	//write result for this block to global mem
	if(tid<32)
	{
		volatile int *vmem = idata;
		vmem[tid]+=vmem[tid+32];
		vmem[tid]+=vmem[tid+16];
		vmem[tid]+=vmem[tid+8];
		vmem[tid]+=vmem[tid+4];
		vmem[tid]+=vmem[tid+2];
		vmem[tid]+=vmem[tid+1];

	}

	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];

}

idata是全局内存，并且通过同步和volatile保证内存写入的有序性，但是这也导致了每次写入都必须写回到全局内存，而不能使用缓存行。将idata的访问替换为对共享内存的访问：

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
__global__ void reduceSmem(int * g_idata,int * g_odata,unsigned int n)
{
	//set thread ID
    __shared__ int smem[DIM];
	unsigned int tid = threadIdx.x;
	//unsigned int idx = blockDim.x*blockIdx.x+threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the
	int *idata = g_idata + blockIdx.x*blockDim.x;

    smem[tid]=idata[tid];
	__syncthreads();
	//in-place reduction in global memory
	if(blockDim.x>=1024 && tid <512)
		smem[tid]+=smem[tid+512];
	__syncthreads();
	if(blockDim.x>=512 && tid <256)
		smem[tid]+=smem[tid+256];
	__syncthreads();
	if(blockDim.x>=256 && tid <128)
		smem[tid]+=smem[tid+128];
	__syncthreads();
	if(blockDim.x>=128 && tid <64)
		smem[tid]+=smem[tid+64];
	__syncthreads();
	//write result for this block to global mem
	if(tid<32)
	{
		volatile int *vsmem = smem;
		vsmem[tid]+=vsmem[tid+32];
		vsmem[tid]+=vsmem[tid+16];
		vsmem[tid]+=vsmem[tid+8];
		vsmem[tid]+=vsmem[tid+4];
		vsmem[tid]+=vsmem[tid+2];
		vsmem[tid]+=vsmem[tid+1];

	}

	if (tid == 0)
		g_odata[blockIdx.x] = smem[0];

}

除了将idata替换为共享内存，没有其他变化。执行后性能会有显著提升。由于时间有限，本节没有动手实验，直接搬了参考博客的代码。再进一步结合循环展开，充分利用带宽，性能还会有更大的提升。

4.4 合并全局内存访问

在之前的矩阵转置的例子，读取和写入总有一个是不连续的，存在非合并的访问。本节将学习如何使用共享内存合并访问。

为了避免交叉访问，可以用二维共享内存缓存矩阵数据，然后从共享内存中读取列存储到全局内存，因为共享内存按列读取的延迟更低，不过会遇到上述的冲突问题。这样读取是按照行进行的，写入也是按行进行的，原来的对全局内存的不连续读转换为了对共享内存的按列读取。

re-5

由于矩阵块可能是非方形，所以需要对索引重处理，具体而言，最开始从全局内存读取时线程是在行上连续的，因此可以用threadIdx.x作为读取的列号，连续读取行；而从共享内存读取列时，由于blockDim.y和blockDim.x可以不相等，不能用threadIdx.x作为行号了，要重新计算一个行号和列号，所以先将原来的线程线性ID算出来，再根据块的维度重新分配读取的数据索引。

image-20241031154324228

对于原矩阵中的一个数据点，其ix为threadIdx.x+blockDim.x*blockIdx.x，iy=threadIdx.y+blockDim.y*blockIdx.y。转置之后，其行号列号为：

1
2
3
4
5
6
7
unsigned int bidx,irow,icol;
bidx=threadIdx.y*blockDim.x+threadIdx.x;
irow=bidx/blockDim.y;
icol=bidx%blockDim.y;

ix=blockIdx.y*blockDim.y+icol;
iy=blockIdx.x*blockDim.x+irow;

bidx是线程在块内的线性ID，所以借助共享内存转置的核函数为：

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
__global__ void transformSmem(float * in,float* out,int nx,int ny)
{
	__shared__ float tile[BDIMY][BDIMX];
	unsigned int ix,iy,transform_in_idx,transform_out_idx;
// 1
	ix=threadIdx.x+blockDim.x*blockIdx.x;
    iy=threadIdx.y+blockDim.y*blockIdx.y;
	transform_in_idx=iy*nx+ix;
// 2
	unsigned int bidx,irow,icol;
	bidx=threadIdx.y*blockDim.x+threadIdx.x;
	irow=bidx/blockDim.y;
	icol=bidx%blockDim.y;
// 3
	ix=blockIdx.y*blockDim.y+icol;
	iy=blockIdx.x*blockDim.x+irow;
// 4
	transform_out_idx=iy*ny+ix;
	if(ix<nx&& iy<ny)
	{
		tile[threadIdx.y][threadIdx.x]=in[transform_in_idx];
		__syncthreads();
		out[transform_out_idx]=tile[icol][irow];
	}
}

不过实际运行结果，核函数的执行时间(3.43ms)比连续读不连续写快(3.8ms)，比连续写不连续读还是更慢(2.65ms)，由于块设置为了32*32，产生了冲突问题。

因此尝试做填充，只用了一个元素的填充就完全消除了冲突，核函数执行时间为2.14。

32*32不连续读，连续写入：

32*32借助共享内存连续读连续写入

由于之前尝试过，不连续读，连续写入在块比较小的时候，读取的缓存行内数据可以很快重用，所以可以减小块的大小以充分利用缓存，因此采用了8x32的块，保持写入的连续性，使整体块小一些，核函数执行时间只有1.56ms，现在对于使用共享内存的方法，也调整块的大小。在8x32得到了最佳性能1.57ms。其实在那一节已经通过nsys看过带宽已经达到理论峰值，无法提升，所以本章节使用共享内存并不会超过当时已经得出的带宽。但是8x32的共享内存版本存在bank conflicts，但没找到合适的填充完全消除。

4.5 常量内存

常量内存是核函数只读的，只有主机可以写常量内存。核函数通过将DRAM的常量内存缓存到片上的常量缓存来读取。对于常量内存，最佳的访问方式是线程束访问同一个位置，读取成本和线程束中线程读取常量内存的地址个数呈线性。

常量内存的声明前缀如下：

1
__constant

初始化常量也就是将常量读取到片上通过以下函数完成：kind的默认参数是传输到设备。

1
cudaError_t cudaMemcpyToSymbol(const void *symbol, const void * src,  size_t count, size_t offset, cudaMemcpyKind kind)

本节以一维stencil计算为例，使用常量内存实现一维stencil。

计算的值为函数的导数，使用差分近似。

计算时c_i是固定的，因此可以存储到常量内存，便于线程束读取。而每个输入数据都会被使用8次，因此使用共享内存缓存输入数据。

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
__global__ void stencilKernel(float * in,float * out,int N)
{
    //开辟包含两端padding的共享内存
    __shared__ float smem[32+2*SIZE+2];
    //线程全局id
    int idx=threadIdx.x+blockDim.x*blockIdx.x;
    //线程在共享内存中的数据id 并将数据复制到共享内存
    int sidx=threadIdx.x+SIZE;
    if(idx>=N) return;
    smem[sidx]=in[idx];
    //确认是否可以填充
    if (threadIdx.x<SIZE)
    {
        if(idx>SIZE)		    		//有左端填充
            smem[sidx-SIZE]=in[idx-SIZE];
        if(idx<N-32)	//有右端填充
            smem[sidx+32]=in[idx+32];
    }
    __syncthreads();
    //无法计算的位置
    if (idx<SIZE||idx>=N-SIZE){
        out[idx] = 0.0;
        return;
    }

    float temp=0;
    //stencil
    #pragma unroll
    for(int i=1;i<=SIZE;i++)
    {
        temp+=coef[i-1]*(smem[sidx+i]-smem[sidx-i]);
    }
    out[idx]=temp;
}

11.78us。如果不使用共享内存，而是使用L1缓存来缓存数据为11.3us，这种情况下两种内存是差不多的。

除了常量内存，还可以使用只读缓存。使用只读缓存有两种方法，需要传入主机内存的指针。

1
2
3
4
5
6
7
8
9
10
__global__ void kernel(float* output, float* input) {
...
output[idx] += __ldg(&input[idx]);
...
}

void kernel(float* output, const float* __restrict__ input) {
...
output[idx] += input[idx];
}

常量缓存的优势是统一读取，而只读缓存适合分散读取。这两部分对于核函数都是只读的。需要视情况使用。使用方式比较简单，略过实现了。

4.6 线程束shuffle指令

线程束shuffle指令是CUDA提供的特殊机制，允许线程束内两个线程互相访问对方的寄存器。由于寄存器在硬件上是临近的，直接交换数据能够提供比任何其他形式内存更高效的数据访问。

shuffle指令是基于线程束提出的。每个线程都可以通过threadIdx.x确认线程束的ID和线程束内索引。

1
2
unsigned int LaneID=threadIdx.x%32;
unsigned int warpID=threadIdx.x/32;

4.6.1 线程束洗牌指令的不同形式

线程束洗牌指令有两组，分别用于整型和浮点。

线程束内交换整型变量的函数如下：

1
int __shfl(int var,int srcLane,int width=warpSize);

var就是需要得到的变量名。而srclane和width则用于确认获取数据的线程位置。默认的width是线程束大小，当width为默认值，srclane就是束内线程。

另一个函数调用时获取当前束内编号-delta编号线程的var值：

1
int __shfl_up(int var,unsigned int delta,int with=warpSize);

还有一个对称的：

1
int __shfl_down(int var,unsigned int delta,int with=warpSize);

最后一个是一个比较特别的指令：

1
int __shfl_xor(int var,int laneMask,int with=warpSize);

lanemask与当前线程索引得到的就是目标线程的编号。

浮点数的函数版本是上述函数的重载，只需要将var传入浮点类型。

4.6.2 使用shuffle指令进行数据交换

对于前几个常规的shuffle指令，使用都比较简单，例如：

1
2
3
4
5
__global__ void test_shfl_broadcast(int *in,int*out,int const srcLans){
    int value=in[threadIdx.x];
    value=__shfl(value,srcLans);
    out[threadIdx.x]=value;
}

__shfl_down和__shfl_up的使用也是一样的。另外，这两个调用超出范围的部分会保持原值，而shfl调用是会对width取模的：

1
2
3
4
5
__global__ void test_shfl_wrap(int *in,int*out,int const offset){
    int value=in[threadIdx.x];
    value=__shfl(value,threadIdx.x+offset);
    out[threadIdx.x]=value;
}

重点需要关注__shfl_xor，该操作足够灵活以组合出任何变换。先假设mask是1，这种情况下，最低位会翻转，因此就是完成FIGURE-5-23的蝶式变换。要实现指定的交换，需要仔细考虑mask的选取。因此接下来使用规约为例，使用__shfl_xor实现规约。

4.6.3 使用线程束洗牌指令的并行规约

在共享内存一节已经实现过用共享内存来传递线程规约的结果，本节使用shuffle指令来传递线程规约的结果，从而进一步提升效率，避免任何内存写入。

上面的例子只是简单交换了相邻的两个数据，shuffle的用法还有多种。最简单的来说，设mask是2的幂次方，那么交换的就是(lane + mask)%lane，因为mask位会翻转，其他位不变。对于规约，需要从16开始交换，即0和16线程交换值相加，依次类推，32个线程束内数据的规约如下：

1
2
3
4
5
6
7
8
9
__inline__ __device__ int warpReduce(int localSum)
{
    localSum += __shfl_xor(localSum, 16);
    localSum += __shfl_xor(localSum, 8);
    localSum += __shfl_xor(localSum, 4);
    localSum += __shfl_xor(localSum, 2);
    localSum += __shfl_xor(localSum, 1);
    return localSum;
}

与之前不同的是，每个线程都得到了规约的结果。线程束规约之后，再对线程束的结果进行规约就可以得到整个线程块的结果。(由于线程块的大小取了1024，下面的下标计算是比较简陋的)。以下是包含循环展开的使用shuffle指令的reduce，N=1 << 24所以边界没有仔细处理。

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
__global__ void reductShuffle(int *in, int *out, int size)
{
	//set thread ID
    __shared__ int smem[1024];
    unsigned int idx = blockDim.x*blockIdx.x*blockPerthread+threadIdx.x;
    if(idx+(blockPerthread-1) * blockDim.x >= size) return;
    int tmp0;
    tmp0 = in[idx];
    tmp0 += in[idx + blockDim.x];
    tmp0 += in[idx + blockDim.x * 2];
    tmp0 += in[idx + blockDim.x * 3];
    tmp0 += in[idx + blockDim.x * 4];
    tmp0 += in[idx + blockDim.x * 5];
    tmp0 += in[idx + blockDim.x * 6];
    tmp0 += in[idx + blockDim.x * 7];

	int mySum=tmp0;
	int laneIdx=threadIdx.x%warpSize;
	int warpIdx=threadIdx.x/warpSize;

	mySum=warpReduce(mySum);
	if(laneIdx==0)
		smem[warpIdx]=mySum;
	__syncthreads();

	mySum=(threadIdx.x<1024/warpSize)?smem[laneIdx]:0;
	if(warpIdx==0)
		mySum=warpReduce(mySum);

  if(threadIdx.x==0)
		out[blockIdx.x]=mySum;
}

最终的核函数时间如下：

Version T
基本规约 5.61ms
交错归约 2.77ms
交错归约+循环展开 510us
shuffle+循环展开 370us

从性能指标上来看，最大的提升是memory throughput，从交错循环+循环展开的141.53Gbyte/s提升到了181.53，因为块内规约的过程中没有任何访存操作，只有规约前累积8个块的数据，而这部分读取是连续的，因此整体访存效率很高，比之前的转置的带宽以及stream测试的带宽还要更高。

5 流和并发

本章主要介绍流和事件的概念，理解网格级并发，核函数执行与CPU执行/数据传输的OVERLAP等内容。

5.1 流和事件概述

5.1.1 流

使用CUDA求解问题的操作步骤包含主机端设备主存分配，数据传输，核函数启动，复制数据回主机。CUDA流能封装这些操作，并保持操作顺序。一个流中的操作是严格有序的，但不同流之间没有限制，多个流同时启动多个内核，就能实现网格级别并行，并且同时CPU还可以执行其他指令。

流在CUDA的API调用可以实现流水线和双缓冲。CUDA的API除了之前的同步接口，也有异步的，保证主机可以执行其他指令，执行多个流。不过流的并行还是受限于设备，只能提高设备的利用率，设备被占用时，其他流也只能等待。

之前的CUDA操作其实也是在流中运行的，只是隐式执行。流可以据此分为：

  • 隐式：空流
  • 显式：非空流

如果没有特别声明一个流，所有操作都是在默认空流中完成的。空流是无法管理的。所以控制流时必须使用非空流，基于流的异步内核支持以下类型的并发：

  • 主机和设备计算
  • 主机计算和主机设备数据传输
  • 主机设备传输和设备计算
  • 多个设备计算

异步执行自然有异步的函数，例如数据传输：

1
cudaError_t cudaMemcpyAsync(void* dst, const void* src, size_t count,cudaMemcpyKind kind, cudaStream_t stream = 0);

最后一个参数就是流，流的相关接口如下，destroy会等流执行完毕后才回收资源。

1
2
3
cudaStream_t a;
cudaError_t cudaStreamCreate(cudaStream_t* pStream);
cudaError_t cudaStreamDestroy(cudaStream_t stream);

执行异步数据传输时，主机内存必须是固定的。在全局内存一章已经说明了，主机内存可能随时移动，为了保证传输时内存是非分页的，需要使用固定内存分配：

1
2
cudaError_t cudaMallocHost(void **ptr, size_t size);
cudaError_t cudaHostAlloc(void **pHost, size_t size, unsigned int flags);

非空流执行内核时需要加入流：

1
kernel_name<<<grid, block, sharedMemSize, stream>>>(argument list);

在流的执行过程中，可以查询流的执行状态：

1
2
cudaError_t cudaStreamSynchronize(cudaStream_t stream);
cudaError_t cudaStreamQuery(cudaStream_t stream);

Sychronize会阻塞主机到流完成，而Query则会立刻返回，如果流执行结束会返回cudaSuccess，否则返回cudaErrorNotReady。

下面是常见的多流执行的模式：

1
2
3
4
5
6
7
8
9
for (int i = 0; i < nStreams; i++) {
    int offset = i * bytesPerStream;
    cudaMemcpyAsync(&d_a[offset], &a[offset], bytePerStream, streams[i]);
    kernel<<grid, block, 0, streams[i]>>(&d_a[offset]);
    cudaMemcpyAsync(&a[offset], &d_a[offset], bytesPerStream, streams[i]);
}
for (int i = 0; i < nStreams; i++) {
    cudaStreamSynchronize(streams[i]);
}

5.1.2 流调度

在Fermi架构上所有流都是在单一硬件上串行执行的。当要执行某个网格的时候CUDA会检测任务依赖关系，如果其依赖于其他结果，那么要等结果出来后才能继续执行。单一流水线可能会导致虚假依赖关系：

这里只有红圈内部的操作是真正并行的，因为A执行时因为B依赖于A队列被阻塞，B执行时也是，只有到C时才发现可以并行执行P，然而事实上P和A B之间也没有依赖，可以在更早的时间并行进行。

解决上述虚假依赖的方式是建立多个工作队列。Hyper-Q技术就是用32个硬件队列同时执行多个流，最大化并发：

5.1.3 流的优先级

3.5以上的设备可以给流设定优先级，数值上更小的核函数先执行。

1
2
cudaError_t cudaStreamCreateWithPriority(cudaStream_t* pStream, unsigned int flags,int priority);   //创建有优先级的流
cudaError_t cudaDeviceGetStreamPriorityRange(int *leastPriority, int *greatestPriority);	//查询优先级

5.1.4 事件

事件用于同步流的执行，检查和记录流的进度。

1
2
3
cudaEvent_t event;
cudaError_t cudaEventCreate(cudaEvent_t* event);
cudaError_t cudaEventDestroy(cudaEvent_t event);

主要用途之一是记录事件之间的时间间隔。事件通过下面指令添加到CUDA流：

1
2
cudaError_t cudaEventRecord(cudaEvent_t event, cudaStream_t stream = 0);
cudaError_t cudaEventElapsedTime(float* ms, cudaEvent_t start, cudaEvent_t stop);

可以等待或查询事件完成状态：

1
2
cudaError_t cudaEventSynchronize(cudaEvent_t event);
cudaError_t cudaEventQuery(cudaEvent_t event);

5.1.5 流同步

既然流可以并行执行，也必然可以进行同步。以处理通信，避免竞争等。非空流对于所有CPU操作都是非阻塞的，之前已经提到了流的两种分类：

  • 异步流（非空流）
  • 同步流（空流/默认流）

显式声明的流都是异步流，异步流不会阻塞主机。但非空流并不都是非阻塞的。也可以分为阻塞流和非阻塞流。果一个非空流被声明为非阻塞的，那么不存在任何阻塞情况，如果声明为阻塞流，则会被空流阻塞。接下来详细说明这两种流。

阻塞流会被空流阻塞，例如下面的三个流：

1
2
3
kernel_1<<<1, 1, 0, stream_1>>>();
kernel_2<<<1, 1>>>();
kernel_3<<<1, 1, 0, stream_2>>>();

空流会等待stream1，stream2会等待空流，但是对于主机来说三个核都是异步的。默认创建的流是阻塞的，非阻塞流用以下函数接口创建：

1
2
3
4
cudaError_t cudaStreamCreateWithFlags(cudaStream_t* pStream, unsigned int flags);
//两种flag
cudaStreamDefault;// 默认阻塞流
cudaStreamNonBlocking: //非阻塞流，对空流的阻塞行为失效。

之前已经有cudaMemcpy和cudaDeviceSynchronize等进行同步。CUDA还可以通过事件实现跨流同步：

1
2
3
cudaError_t cudaEventSynchronize(cudaEvent_t event);
cudaError_t cudaEventQuery(cudaEvent_t event);
cudaError_t cudaStreamWaitEvent(cudaStream_t stream, cudaEvent_t event);

这条命令的含义是，指定的流要等待指定的事件，事件完成后流才能继续，这个事件可以在这个流中，也可以不在，当在不同的流的时候，就实现了跨流同步。

事件还可以进行一些配置：

1
2
3
4
5
6
cudaError_t cudaEventCreateWithFlags(cudaEvent_t* event, unsigned int flags);
//possible flag:
cudaEventDefault
cudaEventBlockingSync
cudaEventDisableTiming
cudaEventInterprocess

其中cudaEventBlockingSync指定使用cudaEventSynchronize同步会造成阻塞调用线程。cudaEventSynchronize默认使用cpu周期不断重复查询事件状态，而当指定了事件是cudaEventBlockingSync的时候，会将查询放在另一个线程中，而原始线程继续执行，直到事件满足条件，才会通知原始线程，这样可以减少CPU的浪费，但是由于通讯的时间，会造成一定的延迟。cudaEventDisableTiming表示事件不用于计时，可以减少系统不必要的开支也能提升cudaStreamWaitEvent和cudaEventQuery的效率。cudaEventInterprocess表明可能被用于进程之间的事件。

5.2 并发内核执行

本节介绍并发内核执行的一些基本问题：

  • 深度优先或广度优先调度
  • 调整硬件工作队列
  • 避免虚假依赖
  • 检查默认流的阻塞
  • 流之间添加依赖关系
  • 检查资源使用对并发的影响

5.2.1 流的并发与工作队列

本节直接采用了reference博客的代码运行测试。首先是并发流的创建和运行：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
cudaStream_t *stream=(cudaStream_t*)malloc(n_stream*sizeof(cudaStream_t));
for(int i=0;i<n_stream;i++)
{
    cudaStreamCreate(&stream[i]);
}
dim3 block(1);
dim3 grid(1);
cudaEvent_t start,stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start);
for(int i=0;i<n_stream;i++)
{
    kernel_1<<<grid,block,0,stream[i]>>>();
    kernel_2<<<grid,block,0,stream[i]>>>();
    kernel_3<<<grid,block,0,stream[i]>>>();
    kernel_4<<<grid,block,0,stream[i]>>>();
}
cudaEventRecord(stop);
CHECK(cudaEventSynchronize(stop));
float elapsed_time;
cudaEventElapsedTime(&elapsed_time,start,stop);
printf("elapsed time:%f ms\n",elapsed_time);

每个流都按顺序执行4个kernel。最后的同步是必要的。原代码计算量有点不够，增加任意一些计算量就可以在nvvp看到了。此外上面创建流的操作，可以用openmp多线程并行完成。

对于之前的虚假依赖问题，由于现在已经和很少有Fermi架构的GPU，不需要考虑这个问题了。对应的解决方式是用BFS的方式启动kernel，即所有流先启动第一个kernel，再继续启动第二个。

上面的nvvp结果，并发的流只有8个，这是因为Kepler支持的最大Hyper-Q 工作队列数为32，默认只开启8个，因为工作队列有资源消耗。修改工作队列配置的方式是修改环境变量：

1
2
3
4
#For Bash or Bourne Shell:
export CUDA_DEVICE_MAX_CONNECTIONS=32
#For C-Shell:
setenv CUDA_DEVICE_MAX_CONNECTIONS 32

5.2.2 并发资源的限制与流的阻塞和依赖

并发的最大化是受硬件资源限制的，上面的例子的核函数资源很小，因此多个工作队列上的多个流完全可以并行，但是如果增加核函数使用的资源，流的并发性就会降低。例如将block和网格都增加到32。执行时间会显著变长，如果进一步增加，nvvp也可以看到流的并发性降低。

默认情况下，空流会被阻塞，非空流也会被空流阻塞，例如下图的核函数：

1
2
3
4
5
6
7
for(int i=0;i<n_stream;i++)
{
    kernel_1<<<grid,block,0,stream[i]>>>();
    kernel_2<<<grid,block,0,stream[i]>>>();
    kernel_3<<<grid,block>>>();
    kernel_4<<<grid,block,0,stream[i]>>>();
}

上面的阻塞以及虚假依赖是需要避免的，但有时是需要创建流之间的依赖关系的，这时就要使用事件。事件不需要计时，所以创建时，可以声明为：

1
2
3
4
5
cudaEvent_t * event=(cudaEvent_t *)malloc(n_stream*sizeof(cudaEvent_t));
for(int i=0;i<n_stream;i++)
{
    cudaEventCreateWithFlag(&event[i],cudaEventDisableTiming);
}

通过事件添加依赖：

1
2
3
4
5
6
7
8
9
for(int i=0;i<n_stream;i++)
{
    kernel_1<<<grid,block,0,stream[i]>>>();
    kernel_2<<<grid,block,0,stream[i]>>>();
    kernel_3<<<grid,block,0,stream[i]>>>();
    kernel_4<<<grid,block,0,stream[i]>>>();
    cudaEventRecord(event[i],stream[i]);
    cudaStreamWaitEvent(stream[n_stream-1],event[i],0);
}

这样最后一个流就会等待前面所有的流执行完毕后才可以完成。

5.3 重叠内核执行和数据传输

Fermi架构和Kepler架构都有两个数据传输队列，从设备到主机和从主机到设备是分离的，因此可以重叠完成。只要kernel的计算和特定数据传输无关，二者就可以在不同的流中重叠执行。

本节使用向量加法为例，N_REPEAT用于控制计算强度，便于nvvp查看性能数据：

1
2
3
4
5
6
7
8
9
10
__global__ void sumArraysGPU(float*a,float*b,float*res,int N)
{
    int idx=blockIdx.x*blockDim.x+threadIdx.x;
    if(idx < N)
    //for delay
    {
        for(int j=0;j<N_REPEAT;j++)
            res[idx]=a[idx]+b[idx];
    }
}

这个过程，数据的传输和计算是不能重叠的，但是向量的各个位是独立的，因此可以把向量分块，分为N-SEGMENT个流去执行，就能重叠计算和数据传输：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
cudaStream_t stream[N_SEGMENT];
for(int i=0;i<N_SEGMENT;i++)
{
    CHECK(cudaStreamCreate(&stream[i]));
}
cudaEvent_t start,stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start,0);
for(int i=0;i<N_SEGMENT;i++)
{
    int ioffset=i*iElem; CHECK(cudaMemcpyAsync(&a_d[ioffset],&a_h[ioffset],nByte/N_SEGMENT,cudaMemcpyHostToDevice,stream[i])); CHECK(cudaMemcpyAsync(&b_d[ioffset],&b_h[ioffset],nByte/N_SEGMENT,cudaMemcpyHostToDevice,stream[i]));
    sumArraysGPU<<<grid,block,0,stream[i]>>>(&a_d[ioffset],&b_d[ioffset],&res_d[ioffset],iElem);
    CHECK(cudaMemcpyAsync(&res_from_gpu_h[ioffset],&res_d[ioffset],nByte/N_SEGMENT,cudaMemcpyDeviceToHost,stream[i]));
}
//timer
CHECK(cudaEventRecord(stop, 0));
CHECK(cudaEventSynchronize(stop));

数据传输需要调用异步接口，还需要注意主机内存要用cudaHostAlloc声明为固定内存。可以观察到核函数和数据在主机和设备之间的传输重叠。这里因为计算强度不够，所以无法完全cover数据传输的时间。

5.4 GPU和CPU并行

在GPU执行流的同时CPU也可以完成一部分工作。可以通过事件来确认流是否执行结束：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
cudaEvent_t start,stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start,0);
for(int i=0;i<N_SEGMENT;i++)
{
    int ioffset=i*iElem;
    CHECK(cudaMemcpyAsync(&a_d[ioffset],&a_h[ioffset],nByte/N_SEGMENT,cudaMemcpyHostToDevice,stream[i]));
    CHECK(cudaMemcpyAsync(&b_d[ioffset],&b_h[ioffset],nByte/N_SEGMENT,cudaMemcpyHostToDevice,stream[i]));
    sumArraysGPU<<<grid,block,0,stream[i]>>>(&a_d[ioffset],&b_d[ioffset],&res_d[ioffset],iElem);
    CHECK(cudaMemcpyAsync(&res_from_gpu_h[ioffset],&res_d[ioffset],nByte/N_SEGMENT,cudaMemcpyDeviceToHost,stream[i]));
}
//timer
CHECK(cudaEventRecord(stop, 0));
int counter=0;
while (cudaEventQuery(stop)==cudaErrorNotReady)
{
    counter++;
}
printf("cpu counter:%d\n",counter);

5.5 流回调

流回调可插入流中，在前面的任务完成后就会调用这个函数。回调函数不可以调用CUDA API，不可以执行同步。流函数有固定的参数：

1
2
3
void CUDART_CB my_callback(cudaStream_t stream, cudaError_t status, void *data) {
    printf("callback from stream %d\n", *((int *)data));
}

使用以下接口插入流：

1
cudaError_t cudaStreamAddCallback(cudaStream_t stream,cudaStreamCallback_t callback, void *userData, unsigned int flags);

简单的回调插入：

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
void CUDART_CB my_callback(cudaStream_t stream,cudaError_t status,void * data)
{
    printf("call back from stream:%d\n",*((int *)data));
}
int main(int argc,char **argv)
{
    //asynchronous calculation
    int iElem=nElem/N_SEGMENT;
    cudaStream_t stream[N_SEGMENT];
    for(int i=0;i<N_SEGMENT;i++)
    {
        CHECK(cudaStreamCreate(&stream[i]));
    }
    cudaEvent_t start,stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start,0);
    for(int i=0;i<N_SEGMENT;i++)
    {
        int ioffset=i*iElem;
        CHECK(cudaMemcpyAsync(&a_d[ioffset],&a_h[ioffset],nByte/N_SEGMENT,cudaMemcpyHostToDevice,stream[i]));
        CHECK(cudaMemcpyAsync(&b_d[ioffset],&b_h[ioffset],nByte/N_SEGMENT,cudaMemcpyHostToDevice,stream[i]));
        sumArraysGPU<<<grid,block,0,stream[i]>>>(&a_d[ioffset],&b_d[ioffset],&res_d[ioffset],iElem);
        CHECK(cudaMemcpyAsync(&res_from_gpu_h[ioffset],&res_d[ioffset],nByte/N_SEGMENT,cudaMemcpyDeviceToHost,stream[i]));
        CHECK(cudaStreamAddCallback(stream[i],my_callback,(void *)(stream+i),0));
    }
    //timer
    CHECK(cudaEventRecord(stop, 0));
    int counter=0;
    while (cudaEventQuery(stop)==cudaErrorNotReady){
        counter++;
    }
}