CUDA C编程（十九）重叠内核执行和数据传输

  重叠内核和数据传输表现出不同的行为，并且需要考虑一些与并发内核执行相比不同的因素。Fermi GPU和Kepler GPU有两个复制引擎队列：一个用于将数据传输到设备，另一个用于从设备中将数据提取出来。因此，最多可以重叠两个数据传输，并且只有当他们的方向不同并且被调度到不同的流时才能这样做，否则，所有的数据传输都将是串行的。在应用程序中，还需要检验数据传输和内核执行之间的关系，从而可以区分以下两种情况：1.如果一个内核使用数据A，那么对A进行数据传输必须安排在内核启动前，切必须位于相同的流中；2.如果一个内核完全不使用数据A，那么内核执行和数据传输可以位于不同的流中。在第二种情况下，实现内核和数据传输的并发执行是很容易的：将它们放置在不同的流中，这就已经向运行时表示了并发地执行它们是很安全地。然而，在第一种情况下，要实现数据传输和内核执行之间的重叠会更复杂，因为内核依赖数据作为输入。当内核和传输之间存在依赖关系时，可以使用向量加法示例来检验如何实现重叠数据传输和内核执行。

使 用 深 度 优 先 调 度 重 叠
  实现向量加法的CUDA程序，其基本结果包含3个主要步骤：1.将两个输入向量从主机复制到设备中；2.执行向量加法运算；3.将单一的输出向量从设备返回主机中。从这些步骤中也许不能明显看出计算和通信是如何被重叠的。为了在向量加法中实现重叠，需要将输入和输出数据集划分为子集，并将来自一个子集的通信与来自其他子集的计算进行重叠。具体对向量加法来说，需要将两个长度为N的向量加法问题划分为长度为N/M的向量相加的M个子问题。因为这里的每个子问题都是独立地，所以每一个都可以被安排在不同的CUDA流中，这样它们的计算就可以重叠的。

#include <stdio.h>
#include <cuda_runtime.h>

#define NSTREAM 4
#define BDIM 128
#define CHECK(call)                                                            \
{                                                                              \
    const cudaError_t error = call;                                            \
    if (error != cudaSuccess)                                                  \
    {                                                                          \
        fprintf(stderr, "Error: %s:%d, ", __FILE__, __LINE__);                 \
        fprintf(stderr, "code: %d, reason: %s\n", error,                       \
                cudaGetErrorString(error));                                    \
    }                                                                          \
}

void initialData(float *ip, int size)
{
    int i;

    for(i = 0; i < size; i++)
    {
        ip[i] = (float)(rand() & 0xFF) / 10.0f;
    }
}

void sumArraysOnHost(float *A, float *B, float *C, const int N)
{
    for (int idx = 0; idx < N; idx++)
        C[idx] = A[idx] + B[idx];
}

//加法内核
__global__ void sumArrays(float *A, float *B, float *C, const int N)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if (idx < N)
    {
        for (int i = 0; i < N; ++i)
        {
            C[idx] = A[idx] + B[idx];
        }
    }
}

void checkResult(float *hostRef, float *gpuRef, const int N)
{
    double epsilon = 1.0E-8;
    bool match = 1;

    for (int i = 0; i < N; i++)
    {
        if (abs(hostRef[i] - gpuRef[i]) > epsilon)
        {
            match = 0;
            printf("Arrays do not match!\n");
            printf("host %5.2f gpu %5.2f at %d\n", hostRef[i], gpuRef[i], i);
            break;
        }
    }

    if (match) printf("Arrays match.\n\n");
}

int main(int argc, char **argv)
{
    printf("> %s Starting...\n", argv[0]);

    int dev = 0;
    cudaDeviceProp deviceProp;
    CHECK(cudaGetDeviceProperties(&deviceProp, dev));
    printf("> Using Device %d: %s\n", dev, deviceProp.name);
    CHECK(cudaSetDevice(dev));

    // check if device support hyper-q
    if (deviceProp.major < 3 || (deviceProp.major == 3 && deviceProp.minor < 5))
    {
        if (deviceProp.concurrentKernels == 0)
        {
            printf("> GPU does not support concurrent kernel execution (SM 3.5 "
                    "or higher required)\n");
            printf("> CUDA kernel runs will be serialized\n");
        }
        else
        {
            printf("> GPU does not support HyperQ\n");
            printf("> CUDA kernel runs will have limited concurrency\n");
        }
    }

    printf("> Compute Capability %d.%d hardware with %d multi-processors\n",
           deviceProp.major, deviceProp.minor, deviceProp.multiProcessorCount);

    // set up max connectioin
    char * iname = "CUDA_DEVICE_MAX_CONNECTIONS";
    setenv (iname, "1", 1);
    char *ivalue =  getenv (iname);
    printf ("> %s = %s\n", iname, ivalue);
    printf ("> with streams = %d\n", NSTREAM);

    // set up data size of vectors
    int nElem = 1 << 18;
    printf("> vector size = %d\n", nElem);
    size_t nBytes = nElem * sizeof(float);

    // malloc pinned host memory for async memcpy
    // 在前面向量加法的程序中，数据传输是通过同步复制函数来实现的
    // 要重叠数据传输和内核执行，必须使用异步复制函数
    // 因为异步复制函数需要固定的主机内存，所以需要使用cudaHostAlloc函数，在固定主机内存中修改主机数组的分配
    float *h_A, *h_B, *hostRef, *gpuRef;
    CHECK(cudaHostAlloc((void**)&h_A, nBytes, cudaHostAllocDefault));
    CHECK(cudaHostAlloc((void**)&h_B, nBytes, cudaHostAllocDefault));
    CHECK(cudaHostAlloc((void**)&gpuRef, nBytes, cudaHostAllocDefault));
    CHECK(cudaHostAlloc((void**)&hostRef, nBytes, cudaHostAllocDefault));

    // initialize data at host side
    initialData(h_A, nElem);
    initialData(h_B, nElem);
    memset(hostRef, 0, nBytes);
    memset(gpuRef,  0, nBytes);

    // add vector at host side for result checks
    sumArraysOnHost(h_A, h_B, hostRef, nElem);

    // malloc device global memory
    float *d_A, *d_B, *d_C;
    CHECK(cudaMalloc((float**)&d_A, nBytes));
    CHECK(cudaMalloc((float**)&d_B, nBytes));
    CHECK(cudaMalloc((float**)&d_C, nBytes));

    cudaEvent_t start, stop;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop));

    // invoke kernel at host side
    dim3 block (BDIM);
    dim3 grid  ((nElem + block.x - 1) / block.x);
    printf("> grid (%d, %d) block (%d, %d)\n", grid.x, grid.y, block.x,
            block.y);

    // sequential operation
    CHECK(cudaEventRecord(start, 0));
    CHECK(cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice));
    CHECK(cudaEventRecord(stop, 0));
    CHECK(cudaEventSynchronize(stop));
    float memcpy_h2d_time;
    CHECK(cudaEventElapsedTime(&memcpy_h2d_time, start, stop));

    CHECK(cudaEventRecord(start, 0));
    sumArrays<<<grid, block>>>(d_A, d_B, d_C, nElem);
    CHECK(cudaEventRecord(stop, 0));
    CHECK(cudaEventSynchronize(stop));
    float kernel_time;
    CHECK(cudaEventElapsedTime(&kernel_time, start, stop));

    CHECK(cudaEventRecord(start, 0));
    CHECK(cudaMemcpy(gpuRef, d_C, nBytes, cudaMemcpyDeviceToHost));
    CHECK(cudaEventRecord(stop, 0));
    CHECK(cudaEventSynchronize(stop));
    float memcpy_d2h_time;
    CHECK(cudaEventElapsedTime(&memcpy_d2h_time, start, stop));
    float itotal = kernel_time + memcpy_h2d_time + memcpy_d2h_time;

    printf("\n");
    printf("Measured timings (throughput):\n");
    printf(" Memcpy host to device\t: %f ms (%f GB/s)\n",
           memcpy_h2d_time, (nBytes * 1e-6) / memcpy_h2d_time);
    printf(" Memcpy device to host\t: %f ms (%f GB/s)\n",
           memcpy_d2h_time, (nBytes * 1e-6) / memcpy_d2h_time);
    printf(" Kernel\t\t\t: %f ms (%f GB/s)\n",
           kernel_time, (nBytes * 2e-6) / kernel_time);
    printf(" Total\t\t\t: %f ms (%f GB/s)\n",
           itotal, (nBytes * 2e-6) / itotal);

    // grid parallel operation
    //接下来，需要在NSTREAM流中平均分配该问题的任务。每一个流要处理的元素数量使用以下代码定义：
    int iElem = nElem / NSTREAM;
    size_t iBytes = iElem * sizeof(float);
    grid.x = (iElem + block.x - 1) / block.x;

    cudaStream_t stream[NSTREAM];

    for (int i = 0; i < NSTREAM; ++i)
    {
        CHECK(cudaStreamCreate(&stream[i]));
    }

    CHECK(cudaEventRecord(start, 0));

    // initiate all work on the device asynchronously in depth-first order
    // 接下来，可以使用一个循环来为几个流同时调度iElem个元素的通信和计算
    for (int i = 0; i < NSTREAM; ++i)
    {
        int ioffset = i * iElem;
        CHECK(cudaMemcpyAsync(&d_A[ioffset], &h_A[ioffset], iBytes,
                              cudaMemcpyHostToDevice, stream[i]));
        CHECK(cudaMemcpyAsync(&d_B[ioffset], &h_B[ioffset], iBytes,
                              cudaMemcpyHostToDevice, stream[i]));
        sumArrays<<<grid, block, 0, stream[i]>>>(&d_A[ioffset], &d_B[ioffset],
                &d_C[ioffset], iElem);
        CHECK(cudaMemcpyAsync(&gpuRef[ioffset], &d_C[ioffset], iBytes,
                              cudaMemcpyDeviceToHost, stream[i]));
    }

    CHECK(cudaEventRecord(stop, 0));
    CHECK(cudaEventSynchronize(stop));
    float execution_time;
    CHECK(cudaEventElapsedTime(&execution_time, start, stop));

    printf("\n");
    printf("Actual results from overlapped data transfers:\n");
    printf(" overlap with %d streams : %f ms (%f GB/s)\n", NSTREAM,
           execution_time, (nBytes * 2e-6) / execution_time );
    printf(" speedup                : %f \n",
           ((itotal - execution_time) * 100.0f) / itotal);

    // check kernel error
    CHECK(cudaGetLastError());

    // check device results
    checkResult(hostRef, gpuRef, nElem);

    // free device global memory
    CHECK(cudaFree(d_A));
    CHECK(cudaFree(d_B));
    CHECK(cudaFree(d_C));

    // free host memory
    CHECK(cudaFreeHost(h_A));
    CHECK(cudaFreeHost(h_B));
    CHECK(cudaFreeHost(hostRef));
    CHECK(cudaFreeHost(gpuRef));

    // destroy events
    CHECK(cudaEventDestroy(start));
    CHECK(cudaEventDestroy(stop));

    // destroy streams
    for (int i = 0; i < NSTREAM; ++i)
    {
        CHECK(cudaStreamDestroy(stream[i]));
    }

    CHECK(cudaDeviceReset());
    return(0);
}

  下图显示了Tesla K40设备典型的时间轴，使用了8个硬件工作队列和4个CUDA流来重叠内核执行和数据传输。相对于阻塞的默认流执行，该流执行实现了近40%的性能提升。下图显示了三种重叠：1.不同流中内核的互相重叠；2.内核与其他流中的数据传输重叠；3.在不同流以及不同方向上的数据传输互相重叠。下图还展示了两种阻塞行为：1.内核被同一流中先前的数据传输所阻塞；2.从主机到设备的数据传输被同一方向上先前的数据传输所阻塞。

  虽然从主机到设备的数据传输是在4个不同的流中执行的，但时间轴显示它们是按顺序执行的，因为实际上它们是通过相同的复制引擎队列来执行的。接下来，可以尝试将硬件工作队列的数量减少至一个，然后重新运行，测试一下其性能。下图显示了在Tesla K40设备上产生的时间轴。与八个工作1队列的相比，没有显著差异。因为每个流只执行单一的一个内核，所以减少工作队列的数目并没有增加虚假依赖关系，同样，现存的虚假依赖关系(由主机到设备的复制队列引起的)也没有减少。
  减少K40中工作队列的数量，可以创造一个类似于Fermi GPU的环境：一个工作队列和两个复制队列。如果在Fermi GPU运行相同的测试，会发现虚假的依赖关系是确实存在的。这是由Kepler的工作调度机制导致的，在网格管理单元(Grid Management Unit,GMU)中实现。GMU负责对发送到GPU中的工作进行管理和排序。通过对GMU的分析有助于减少虚假的依赖关系。

  Kepler引入了一个新的网格管理和调度控制系统，即网格管理单元(GMU)。GMU可以暂停新网格的调度，使网格排队等待且暂停网格直到它们准备好执行，这样就使运行时变得非常灵活强大，动态并行就是一个很好的例子。在Fermi设备上，网格直接从流队列被传到CUDA工作分配器(CUDA Work Distributor,CWD)中。在Kepler设备上，网格被发送到GMU上，GMU在对GPU上执行的网格进行管理和优先级排序。GMU创建了多个硬件工作队列，从而减少或消除了虚假的依赖关系。通过GMU，流可以作为单独的工作流水线。即GMU被限制只能创建一个单一的硬件工作队列，根据以上测试结果证实，通过GMU进行的网格依赖性分析也可以帮助消除虚假的依赖关系。

使 用 广 度 优 先 调 度 重 叠
  先前的例子表明，当采用广度有限的方式调度内核时，Fermi GPU可以实现最好的效果。现在，将在重叠数据传输和计算内核中，检验广度优先排序产生的效果，下面的代码演示了使用广度优先的方法来调度流间的计算和通信：

 // initiate all asynchronous transfers to the device
    for (int i = 0; i < NSTREAM; ++i)
    {
        int ioffset = i * iElem;
        CHECK(cudaMemcpyAsync(&d_A[ioffset], &h_A[ioffset], iBytes,
                              cudaMemcpyHostToDevice, stream[i]));
        CHECK(cudaMemcpyAsync(&d_B[ioffset], &h_B[ioffset], iBytes,
                              cudaMemcpyHostToDevice, stream[i]));
    }

    // launch a kernel in each stream
    for (int i = 0; i < NSTREAM; ++i)
    {
        int ioffset = i * iElem;
        sumArrays<<<grid, block, 0, stream[i]>>>(&d_A[ioffset], &d_B[ioffset],
                &d_C[ioffset], iElem);
    }

    // enqueue asynchronous transfers from the device
    for (int i = 0; i < NSTREAM; ++i)
    {
        int ioffset = i * iElem;
        CHECK(cudaMemcpyAsync(&gpuRef[ioffset], &d_C[ioffset], iBytes,
                              cudaMemcpyDeviceToHost, stream[i]));
    }

  下图显示了在K40设备上只使用一个硬件工作队列时的时间轴。与深度优先的方法相比它没有明显的差异，因为Kepler的双向调度机制有助于消除虚假的依赖关系。但如果在Fermi设备上运行相同的测试，在整体性能方面会发现，使用广度优先方法不如使用深度优先方法。由主机到设备复制队列上的争用导致的虚假依赖关系，在主机到设备间的传输完成时，将阻止所有的内核启动。因此，对于Kepler设备而言，在大多数情况下无需关注其工作调度顺序。而在Fermi设备上，要注意这些问题，并且队不同的调度方案做出评估，使工作负载找到最佳的任务调度顺序。