计算机视觉大型攻略 —— CUDA(1)编程基础

基于CUDA的并行编程在计算机视觉和机器学习邻域得到了广泛的应用。[1]Professional Cuda C Programming一书系统的介绍了CUDA的编程模型和各种优化的奇技淫巧，虽说文中GPU的系统架构有些过时，但是基础理论和架构仍然适合当前的主流GPU。准备写几篇blog总(fan)结(yi)一下书中第二章到第六章的部分，最后再举一个例子将这几章的内容贯穿起来。

参考文献

[1] PROFESSIONAL CUDA C Programming. John Cheng, Max Grossman, Ty McKercher.

CUDA编程基础

本文对应[1]的第二章，CUDA Programming Model。

上图取自[1]的Figure2-2，描述了一个典型CUDA程序的结构。从图中可知，一个CUDA程序由Host Code和Parallel code构成，他们分别执行在Host端(CPU)和Device端(GPU)。

• Host端：包括CPU和host内存。通常Host端的程序负责复杂的控制逻辑，以及定义如何将算法数据映射到Device端。Host code采用ANSI C语言。
• Device端: 包括GPU和显存。Device端的程序需要设计成轻量级的线程，处理并行任务，控制逻辑简单。Device端的函数通常也被称作核函数(kernel)。Device code采用CUDA C语言。

Host端与Device端是异步的，即host启动一个kernel后，立刻返回，因此CPU与GPU独立且并行的执行各自的程序。

一个CUDA程序通常包含三个过程，

1. 从CPU内存拷贝数据到显存。(Host memory->Device memory)
2. 调用Kernel函数处理显存数据。
3. 将结果从显存拷贝到内存。

Device端内存管理

Host端与Device端各自拥有独立的内存空间。Host端的内存采用标准C函数分配，拷贝和释放。而CUDA runtime提供了在Host端分配Device内存的函数。如下表，

可通过查阅CUDA编程手册进一步了解他们的函数原型。

其中cudaMemcpy实现了Host端和Device端内存之间的拷贝。

CUDA线程管理

NVIDIA采用Single Instruction, Multiple Thread(SIMT)的多核并行架构。主机端启动一个Kernel函数后，设备端启动大量的线程“同步”执行。（这里同步加上了引号，因为他们只是逻辑上同步，实质上由于计算资源的局限性，CUDA仍然会对线程进行调度）。

CUDA采用两级结构管理这些线程。

如上图所示，

• 一个Kernel生成的所有线程构成一个Grid。这些线程共享相同的全局内存空间(global memory space)。
• 一个Grid包含许多线程blocks。每个block包含许多线程，这些线程可以使用block内的同步机制，共享block的shared memory。

因此，对于Device上的一个线程，可由两个变量唯一确定：blockIdx, threadIdx。这两个变量是由编译器built-in的，在kernel函数中可以直接访问。每个线程都有自己独立的blockIdx, threadIdx, 使用这两个变量来确定自己的逻辑行为或者找到自己应当处理的数据。cuda为grid和block提供了最多三维的划分方式。分别使用blockIdx.x, blockIdx.y, blockIdx.z, threadIdx.x, threadIdx.y, threadIdx.z表示。上图是一个二维划分的例子，也可以理解为blockId.z和threadIdx.z均为1。

另一组kernel函数内嵌(build-in)的变量为blockDim, gridDim。分别表示block, grid每一维的线程和block数量。

__global__ kernel_name<<<grid, block>>>(argument list)

上面这个函数用于执行一个Kernel函数，grid代表了block的数量，block代表了一个block内的线程数量。__global__是关键字，表示该函数可由Host和Device调用，在Device上执行。CUDA C中还有一些其他类似的关键字，Table2-2给出了一个对比。

Grid与Block

Grid代表了一个Grid内的Block数量，Block代表了一个Block内的线程数量。对任何一个并行任务来讲，grid与block的值在设备允许的范围内是可以随意选取的。然而他们的取值一定程度上会影响到的算法执行速度，这取决于GPU性能与算法的设计。大多数时候需要通过实验来确定最终结果。 

在计算机视觉和深度学习算法中，几乎所有的并行任务都是基于数据的并行，即有大量的数据需要采用相同的方式处理。我们对线程的划分，多数情况下等效于我们对数据的划分。比如当处理图像数据时，图像数据通常是一个二维矩阵，我们在划分Grid和Block时，也可采用上图(Figure 2-5)所示的二维划分方法。

下图为[1]中的一个二维划分的示意，

我们以一个例子来解释一下。假设我们要逐元素处理一幅nx*ny大小的单通道图像（比如计算像素的局部特征census或者做卷积等）。现在计划为每一个像素分配一个线程。

第一步， 需要将图像从Host端拷贝到Device端(即global memory，这里还没有涉及global memory的含义, 后续文章会写)。在CUDA编程中，多维数据通常会先在Host端转成一维数据，再将其拷贝到Device端，也就是说Device端看到的是一维内存空间。这样做可以降低程序的复杂度（我们只需要调用一次cudaMalloc函数，减少该函数的调用可以略微改善系统性能，当然主要是改善了程序猿的心情，从而提高了编程效率)。

通常情况下，图像处理中，Host端将多维数组转成一维数组并不会带来额外到的开销，以Opencv的数据结构Mat为例，在连续存储(绝大多数，只要不是自己找麻烦)的情况下，Mat本身就是按照Row-major方式在内存中以一维的形式存放数据，可以通过成员函数.data拿到数据指针。

下图是一个8*6的矩阵的一维Row-Major存放示意。

第二步，在Host端定义线程的划分，

dim3 block(32, 32);
dim3 grid((nx+block.x-1)/block.x, (ny+block.y-1)/block.y);

block定位成了32*32，共1024个线程。在多数的Nvidia架构下, 1024是一个block所支持的线程数的最大值。block的值确定后，grid可以由nx, ny和block.x, block.y计算出来。

第三步，在Device端(Kernel函数), 我们需要为每一个线程找到其对应的像素坐标。 通过内嵌变量blockIdx, blockDim, threadIdx，我们可以获取每个线程对应的图像像素的二维坐标（ix, iy）。

ix=threadIdx.x+blockIdx.x*blockDim.x
iy=threadIdx.y+blockIdx.y*blockDim.y

以ix为例，blockIdx.x表示x方向上block的索引，blockDim.x为每个block在x方向上的线程数量, threadIdx.x为该线程在当前block内x方向上的索引。通过上面的公式可以计算出在当前Grid上，x方向上的线程索引ix，y方向上的线程索引iy, 以此来唯一映射图像的像素。

由于图像在Device端(global memory)中采用一维row-major数组的存放方式，因此，Device内存的索引为

idx=iy*nx+ix;

异步与同步

由于Device与Host是异步的，CUDA还提供了一个同步函数。

cudaError_t cudaDeviceSynchronize(void);

该函数可以让Host端程序阻塞，直到Device端将所有线程执行完。另外有些CUDA runtime API也会隐式的同步Host端与Device端，如cudaMemcpy，该函数只有在前面所有Kernel均执行完采开始执行拷贝，直到拷贝结束后，程序才返回。

CUDA Kernel函数的特点和局限

• 只能访问Device Memory
• 返回类型：void
• 不支持可变数量的参数
• 不支持静态变量
• 不支持函数指针
• 与Host程序异步执行

一个例子: 两个矩阵求和

[1]第二章给出了代码示例，他们均可以从www.wrox.com/go/procudac下载。以sumMatrixOnGPU-2D-grid-2D-block.cu为例，将本节所有知识点再串一下。

该程序将两个二维矩阵按位相加。

int main(int argc, char **argv)
{
    printf("%s Starting...\n", argv[0]);
  
    // set up device
    int dev = 0;
    cudaDeviceProp deviceProp;
    CHECK(cudaGetDeviceProperties(&deviceProp, dev));
    printf("Using Device %d: %s\n", dev, deviceProp.name);
    CHECK(cudaSetDevice(dev));

cudaGetDeviceProperties函数返回当前GPU的信息。

当机器上有多个GPU的时候，cudaSetDevice函数指定使用哪个GPU。

    // set up data size of matrix
    int nx = 1 << 14;
    int ny = 1 << 14;

    int nxy = nx * ny;
    int nBytes = nxy * sizeof(float);
    printf("Matrix size: nx %d ny %d\n", nx, ny);

    // malloc host memory
    float *h_A, *h_B, *hostRef, *gpuRef;
    h_A = (float *)malloc(nBytes);
    h_B = (float *)malloc(nBytes);
    hostRef = (float *)malloc(nBytes);
    gpuRef = (float *)malloc(nBytes);

    // initialize data at host side
    double iStart = seconds();
    initialData(h_A, nxy);
    initialData(h_B, nxy);
    double iElaps = seconds() - iStart;
    printf("Matrix initialization elapsed %f sec\n", iElaps);

初始化两个大小为nx*ny的矩阵h_A, h_B。一般情况下，前缀“h_”代表指向Host端内存的指针，“d_”代表指向Device端内存的指针。指向Device端内存的指针不能在Host代码中访问，只能在Device端代码(Kernel函数)中访问。

    // malloc device global memory
    float *d_MatA, *d_MatB, *d_MatC;
    CHECK(cudaMalloc((void **)&d_MatA, nBytes));
    CHECK(cudaMalloc((void **)&d_MatB, nBytes));
    CHECK(cudaMalloc((void **)&d_MatC, nBytes));

    // transfer data from host to device
    CHECK(cudaMemcpy(d_MatA, h_A, nBytes, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(d_MatB, h_B, nBytes, cudaMemcpyHostToDevice));

分配Device内存。并且将Host端内存h_A, h_B拷贝到Device端d_MatA，d_MatB。CHECK()函数检查函数返回值。

    // invoke kernel at host side
    int dimx = 32;
    int dimy = 32;
    dim3 block(dimx, dimy);
    dim3 grid((nx + block.x - 1) / block.x, (ny + block.y - 1) / block.y);

    iStart = seconds();
    sumMatrixOnGPU2D<<<grid, block>>>(d_MatA, d_MatB, d_MatC, nx, ny);
    CHECK(cudaDeviceSynchronize());

    iElaps = seconds() - iStart;
    printf("sumMatrixOnGPU2D <<<(%d,%d), (%d,%d)>>> elapsed %f sec\n", grid.x,
           grid.y,
           block.x, block.y, iElaps);
    // check kernel error
    CHECK(cudaGetLastError());

为矩阵每一个元素分配一个线程，总共nx*ny个线程。分配block和grid的大小。block划分为32*32，并根据线程数计算出必要的block数。并调用Kernel函数sumMatrixOnGPU2D。

Host端调用Kernel函数后立刻返回，在cudaDeviceSynchronize()阻塞，直到Device端所有线程执行完。

由于Kernel函数调用后立刻返回，无法通过函数返回值检查是否出错。需要调用cudaGetLastError()检查Kernel函数执行中是否存在错误。

    // copy kernel result back to host side
    CHECK(cudaMemcpy(gpuRef, d_MatC, nBytes, cudaMemcpyDeviceToHost));

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

    // free device global memory
    CHECK(cudaFree(d_MatA));
    CHECK(cudaFree(d_MatB));
    CHECK(cudaFree(d_MatC));

程序结果在d_MatC中。将其从Device端内存拷贝到Host端。最后使用cudaFree释放Device端内存。

Kernel函数：

// grid 2D block 2D
__global__ void sumMatrixOnGPU2D(float *MatA, float *MatB, float *MatC, int nx,
                                 int ny)
{
    unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
    unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
    unsigned int idx = iy * nx + ix;

    if (ix < nx && iy < ny)
        MatC[idx] = MatA[idx] + MatB[idx];
}

首先计算每个线程对应的二维坐标ix, iy，然后拿到矩阵Row-Major的一维索引。求和，将结果放到MatC中。

总的来说，掌握了这一章基本上就可以去读写CUDA程序，优化项目了。然而，为了最大化的利用好CUDA的计算能力，仍需要对其内在机制（内存管理，线程管理）进行更深入的理解。