《CUDA编程：基础与实践》读书笔记(4)：CUDA流

1. CUDA流

一个CUDA流指的是由主机发出的在一个设备中执行的CUDA操作序列。除主机端发出的流之外，还有设备端发出的流，但本文不考虑后者。一个CUDA流中的各个操作按照主机发布的次序执行；但来自两个不同CUDA流的操作不一定按照某个次序执行，有可能是并发或者交错地执行。

任何CUDA操作都存在于某个CUDA流中，如果没有明确指定CUDA流，那么所有CUDA操作都是在默认流中执行的。非默认CUDA流由cudaStream_t类型的变量表示，它由如下CUDA运行时API产生与销毁：

cudaError_t cudaStreamCreate(cudaStream_t* pStream);
cudaError_t cudaStreamDestroy(cudaStream_t stream);

为了检查CUDA流中的所有操作是否都在设备中执行完毕，可以使用如下函数：

//阻塞主机直到stream中的所有操作都执行完毕
cudaError_t cudaStreamSynchronize(cudaStream_t stream);
//不阻塞主机，只检查stream中的所有操作是否都执行完毕，若是则返回cudaSuccess，否则返回cudaErrorNotReady
cudaError_t cudaStreamQuery(cudaStream_t stream);

为了产生多个相互独立的CUDA流、实现不同CUDA流之间的并发，主机在向某个CUDA流中发布命令后必须马上获取程序控制权，不等待该CUDA流中的命令在设备中执行完毕。下文将介绍主机如何在向某个CUDA流发布命令后马上取得控制权。此外，也可以在主机端使用多个线程控制多个CUDA流。

2. 核函数与主机的重叠执行

下面是默认CUDA流中数组相加的例子：

cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice);
cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice);
add<<<grid_size, block_size>>>(d_x, d_y, d_z);
cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost);

从设备的角度看，以上4个CUDA操作将在默认CUDA流中按顺序依次执行。从主机的角度看，数据传输是同步的（或者说是阻塞的），比如说主机在执行前两个cudaMemcpy语句时，会等待该命令执行完毕再继续往下走，所以在进行数据传输时，主机是闲置的，不能进行其它操作。不同的是，核函数的启动是异步的（或者说是非阻塞的），意思是主机发出调用核函数的命令后，不会等待命令执行完毕，而会立刻取得程序控制权，然后紧接着发出最后一个cudaMemcpy命令，但是该命令不会立即被执行，因为这是默认流中的CUDA操作，必须等待前一个CUDA操作（即核函数的调用）执行完毕才会开始执行。

根据上述分析可知，主机在发出核函数调用命令后会立刻继续执行接下来的命令。如果下一条命令是主机的某个计算任务，那么就可以实现核函数与主机计算任务的并行计算。

3. 核函数与核函数的重叠执行

因为同一个CUDA流中的CUDA操作在设备中是顺序执行的，所以要实现多个核函数之间的并行就必须使用多个CUDA流。在使用的多个CUDA流中可以有一个默认流，但此时各个流之间并不完全独立，本文不讨论这种情况，只讨论使用多个非默认流的情况。在非默认流中调用核函数时，执行配置必须包含一个流对象，一个名为my_kernel(...)的核函数只能用如下三种调用方式之一：

//N_grid是网格大小，最一般的情形是一个dim3类型的结构体，简单情况下可以是一个整数
//N_block是线程块大小，最一般的情形是一个dim3类型的结构体，简单情况下可以是一个整数
//N_shared是核函数中使用的动态共享内存的字节数，如果没有则设为0
//stream是cudaStream_t类型的CUDA流对象
my_kernel<<<N_grid, N_block>>>(...);
my_kernel<<<N_grid, N_block, N_shared>>>(...);
my_kernel<<<N_grid, N_block, N_shared, stream>>>(...);

下面的例子简单展示了如何使用非默认CUDA流重叠执行多个核函数：

#include "cuda_runtime.h"

void __global__ my_kernel()
{
    // do some calculations
}

int main(void)
{
    const int NUM_STREAMS = 16;
    const int block_size = 128;
    const int grid_size = 8;
    cudaStream_t streams[NUM_STREAMS];

    for (int n = 0; n < NUM_STREAMS; ++n)
    {
        cudaStreamCreate(&(streams[n]));
    }

    for (int n = 0; n < NUM_STREAMS; ++n)
    {
        my_kernel<<<grid_size, block_size, 0, streams[n]>>>();
    }

    for (int n = 0; n < NUM_STREAMS; ++n)
    {
        cudaStreamDestroy(streams[n]);
    }

    return 0;
}

利用CUDA流并发执行多个核函数可以提升GPU硬件的利用率，减少闲置的SM，从而整体上获得性能提升。但当所有CUDA流中对应核函数的线程数总和超过一定阈值后，再增加CUDA流的数量就不会带来更高的加速比了，反而可能使程序的性能下降。制约加速比的因素是GPU的计算资源。

4. 核函数与数据传输的重叠执行

要实现核函数与数据传输的并发，必须让这两个操作处于不同的非默认流，而且数据传输必须使用cudaMemcpy的异步版本，即cudaMemcpyAsync函数。如果使用同步的数据传输函数，主机向一个流发出输出传输命令后就必须等待数据传输完毕，这样核函数与数据传输的重叠也就无法实现。异步传输函数的原型是：

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

在使用异步数据传输函数时，需要将主机内存定义为不可分页内存，这样在程序运行期间操作系统就不会改变主机内存的物理地址。如果给cudaMemcpyAsync函数传入的主机内存是可分页内存，那么函数就会退化到cudaMemcpy，从而导致同步传输，无法达到核函数与数据传输重叠执行的效果。不可分页主机内存的分配与释放可以用如下函数：

cudaError_t cudaMallocHost(void **ptr, size_t size);
cudaError_t cudaFreeHost(void *ptr);

下面给出一个使用CUDA流重叠执行核函数和数据传输的例子：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

const int N = 1 << 22;
const int M = sizeof(float) * N;
const int NUM_STREAMS = 64;
cudaStream_t streams[NUM_STREAMS];

void __global__ add(const float* x, const float* y, float* z, int N)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        z[n] = x[n] + y[n];
    }
}

int main(void)
{
    float *h_x, *h_y, *h_z;
    cudaMallocHost(&h_x, M);
    cudaMallocHost(&h_y, M);
    cudaMallocHost(&h_z, M);
    for (int n = 0; n < N; ++n)
    {
        h_x[n] = 1.23f;
        h_y[n] = 2.34f;
    }

    float *d_x, *d_y, *d_z;
    cudaMalloc(&d_x, M);
    cudaMalloc(&d_y, M);
    cudaMalloc(&d_z, M);

    for (int i = 0; i < NUM_STREAMS; i++)
    {
        cudaStreamCreate(&(streams[i]));
    }

    int N1 = N / NUM_STREAMS;
    int M1 = M / NUM_STREAMS;
    for (int i = 0; i < NUM_STREAMS; i++)
    {
        int off = i * N1;
        cudaMemcpyAsync(d_x + off, h_x + off, M1, cudaMemcpyHostToDevice, streams[i]);
        cudaMemcpyAsync(d_y + off, h_y + off, M1, cudaMemcpyHostToDevice, streams[i]);
        int block_size = 128;
        int grid_size = (N1 - 1) / block_size + 1;
        add<<<grid_size, block_size, 0, streams[i]>>>(d_x + off, d_y + off, d_z + off, N1);
        cudaMemcpyAsync(h_z + off, d_z + off, M1, cudaMemcpyDeviceToHost, streams[i]);
    }


    for (int i = 0; i < NUM_STREAMS; i++)
    {
        cudaStreamDestroy(streams[i]);
    }

    cudaFreeHost(h_x);
    cudaFreeHost(h_y);
    cudaFreeHost(h_z);
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);

    return 0;
}