CUDA编程结构

在GPU上执行的函数称为CUDA核函数（Kernel Function)，核函数会被GPU上多个线程执行。典型的CUDA程序遵循如下模式：

把数据从CPU内存（HOST）拷贝至GPU内存（DEVICE）
调用该Kernel函数，对DEVICE中的数据进行操作
将数据从DEVICE传送至HOST

CUDA内存管理

GPU内存分配：cudaMalloc

cudaError_t cudaMalloc (void** devPtr,  size_t size)

函数返回一个 devPtr，指向所分配的内存。

HOST<->DEVICE数据传输：cudaMemcpy

cudaError_t cudaMemcpy ( void* dst, const void* src, size_t count, cudaMemcpyKind kind )

函数将 src指向的数据，传输至 dst，复制方向由 kind制定，包含如下4种：

cudaMemcpyHostToHost = 0
    Host -> Host 
cudaMemcpyHostToDevice = 1
    Host -> Device 
cudaMemcpyDeviceToHost = 2
    Device -> Host 
cudaMemcpyDeviceToDevice = 3
    Device -> Device

需要注意的是，这个函数是以同步的方式执行，即 cudaMemcpy完成操作并返回值前，程序是阻塞的。

此外，除了Kernel启动外的CUDA函数调用均会返回一个错误枚举类型 cudaError_t，若GPU内存分配成功返回 cudaSuccess，否则返回 cudaErrorMemoryAllocation，可以通过 cudaGetErrorString(cudaError_t error)得到可读的错误消息。

简单的数据传输示例

nvcc中封装了多种编译工具

#include <stdlib.h>
#include <string.h>
#include <time.h>

// A、B、C都是数组，其本质是指向数组首元素的指针。
// 对于一个数组A，A[i] 等价于 *(A + i) 。
__global__ void sumArrayGPU(float *A, float *B, float *C, const int N){
    for(int i = 0;i < N; i++){
        C[i] = A[i] + B[i];
    }
}

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

void initData(float *ip, int size){
    time_t t;
    srand((unsigned int) time(&t));
    for(int i = 0;i < size;i++){
        ip[i] = (float) (rand() &0xFF) / 10.0f;
    }
}

int main(int argc, char** argv){
    const int N = 1024;
    auto memSize = N * sizeof(float);
    float *h_A, *h_B, *h_C;
    h_A = (float *)malloc(memSize);
    h_B = (float *)malloc(memSize);
    h_C = (float *)malloc(memSize);
    initData(h_A, N);
    initData(h_B, N);
    sumArrayCPU(h_A, h_B, h_C, N);

    float *d_A, *d_B, *d_C;
    cudaMalloc((void**) &d_A, memSize);
    cudaMalloc((void**) &d_B, memSize);
    cudaMalloc((void**) &d_C, memSize);
    cudaMemcpy(d_A, h_A, memSize, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, memSize, cudaMemcpyHostToDevice);
    sumArrayGPU<<<16,16>>>(d_A, d_B, d_C, N);
    cudaMemcpy(h_C, d_C, memSize, cudaMemcpyDeviceToHost);

    free(h_A);
    free(h_B);
    free(h_C);
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    return 0;
}

线程管理

由一个Kernel启动所产生的所有线程统称为一个 Grid，此 Grid中所有线程共享全局内存。

一个 Grid中包含多个 Block，一个 Block又包含多个 Thread。不同 Block内的 Thread不能协作。

不同 Thread使用 blockIdx和 threadIdx来互相区分，这两个变量是CUDA的内置变量，为 unit3类型，即包含3个 uint的结构，可以使用 x y z三个字段进行指定。比如 blockIdx.x 、threadIdx.y

Grid和 Block的维度依靠 dim3类型的 gridDim和 blockDim来指定，这种类型的数据也可以通过 x y z来获得。

CUDA网格结构.png

在一段 .cu程序中，会定义两组这种变量。一个是在Host中定义的 dim3类型的 Grid和 Block，该变量仅在Host中可见。另一个是在Kernel函数中定义的 uint3类型，仅在Device端可见。下面代码中会看到这两种变量定义位置的区别。

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

__global__ void checkDim(){
    //内置的预初始化Thread和Grid变量，仅在Device访问
    printf("threadIdx:(%d,%d,%d)  blockIdx:(%d,%d,%d)  blockDim:(%d,%d,%d)" 
    "gridDim:(%d,%d,%d) \n" ,threadIdx.x, threadIdx.y, threadIdx.z, blockIdx.x, blockIdx.y, blockIdx.z,
    blockDim.x ,blockDim.y, blockDim.z, gridDim.x, gridDim.y, gridDim.z);
}

int main(int argc, char** argv){
    const int nElem = 6;
    //手动定义的dim3类型的Grid和Block Size，仅在Host访问
    dim3 block(3);
    dim3 grid((nElem + block.x -1) / block.x);

    printf("grid.x %d grid.y %d grid.z %d \n" ,grid.x, grid.y , grid.z);
    printf("block.x %d block.y %d block.z %d \n" ,block.x, block.y , block.z);

    checkDim<<<grid,block>>>();
    cudaDeviceSynchronize();
    return 0;
}

在 dim3 grid((nElem + block.x -1) / block.x);处， +block.x -1的原因很多人可能不懂，这是为了避免最后一个块不能完全包含一个 Block。比如当 nElem=10，block.x=3，若不这样处理，10/3 = 3，3个 Block，每个 Block包含3个 Thread，这样一共处理9个任务，相比10，少了1个；通过 +block.x -1，(10+3-1)/3 = 4，多出1个 Block来处理第10个任务。简单来说，就是向上取整。

调用Kernel函数

Kernel_Function<<<grid,block>>>(argument list);

<<<grid,block>>>中便是kernel函数的执行配置。

假设有4个 block，每个 block中又有8个 thread，那么核函数为 Kernel_Function<<<4,8>>>(argument list);

由上图也能看出，数据在全局内存中是线性存储的，所以可以通过 blockIdx和 threadIdx获取一个线程的唯一ID，同时建立线程和数据元素间的映射关系。一般情况下，可以通过1维 Grid和 Block来计算全局数据访问的唯一索引：int i = blockIdx.x * blockDim.x + threadIdx.x

Kernel函数的调用与Host的线程是异步的，需要调用下面函数来强制Host等待所有kernel函数执行结束：
cudaError_t cudaDeviceSynchronize(void);

此外，cudaMemcpy()隐式地包含同步功能，当调用 cudaMemcpy在Host和Device进行拷贝时，Host线程将等待 cudaMemcpy操作执行完毕才会继续执行程序。

Kernel函数限定符

先提一嘴，Kernel函数必须是 void返回类型。

有三种限定符：

__global__ Host调用，Device执行
__device__ Device调用，Device执行
__host__ Host调用，Host执行（可省略）

当 __device__ 和 __host__ 一起使用时，函数可以同时在Host和Device上编译。

Kernel函数和C函数的区别

__global__ void sumArrayGPU(float *A, float *B, float *C, const int N){
	int i = threadIdx.x;
        C[i] = A[i] + B[i];
}

从上面这段代码可以看出，for循环消失了。因为内置的 thread坐标替代了原本的数组索引。

错误处理宏

#define CHECK(call)                                                    
{                                                                      
   const cudaError_t error = call;                                     
   if (error != cudaSuccess)                                           
   {                                                                   
      printf("Error: %s:%d, ", __FILE__, __LINE__);                    
      printf("code:%d, reason: %s\n", error, cudaGetErrorString(error)); 
      exit(1);                                                         
   }                                                                   
}

CHECK(cudaMemcpy(d_A, h_A, memSize, cudaMemcpyHostToDevice));如果产生错误会报告一个错误代码，并通过 cudaGetErrorString(error)输出可读信息。

并行线程

由于矩阵用行优先的方式在全局内存中进行存储，对一个给定的线程，可以通过 Thread和 Block索引映射到矩阵坐标上：

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

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

然后将矩阵坐标映射到全局内存的索引上：

idx = iy * nx +ix

矩阵求和 - 2D Grid & 2D Block

使用2D Grid 和 2D Block 计算矩阵加法的代码如下：

#include "common.h"

void SumMatCPU(float *A, float *B, float *C, const int nx, const int ny){
    float *ia = A;
    float *ib = B;
    float *ic = C;

    for(int iy = 0; iy < ny; iy++){
        for(int ix = 0; ix < nx; ix++){
            ic[ix] = ia[ix] + ib[ix];
        }
        ia += nx;
        ib += nx;
        ic += nx;
    }
}

__global__ void SumMatGPU2D(float *matA, float *matB, float * matC, int nx, int ny){
    unsigned int ix = threadIdx.x + blockDim.x * blockIdx.x;
    unsigned int iy = threadIdx.y + blockDim.y * blockIdx.y;
    unsigned int idx = ix + iy * nx;

    if(ix < nx && iy < ny){
        matC[idx] = matA[idx] + matB[idx];
    }
}

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

    int nx = 1<<14; // 1<<14 16384; 1<<16 65535
    int ny = 1<<14;
    int nxy = nx * ny;
    int memSize = nxy * sizeof(float);

    printf("Matrix size : %d, %d\n", nx ,ny);

    //Malloc host mem
    float *h_A, *h_B, *hostRef , * gpuRef;
    h_A = (float *)malloc(memSize);
    h_B = (float *)malloc(memSize);
    hostRef = (float *)malloc(memSize);
    gpuRef = (float *)malloc(memSize);

    double start = seconds();
    initData(h_A, nxy);
    initData(h_B, nxy);
    double epls = seconds() - start;
    printf("Init Data eplased time: %f \n", epls);

    memset(hostRef,0,memSize);
    memset(gpuRef,0,memSize);
    start = seconds();
    SumMatCPU(h_A, h_B, hostRef, nx, ny);
    epls = seconds() - start;
    printf("Sum Mat on CPU cost : %f \n", epls);

    float *d_MatA, *d_MatB, *d_MatC;
    CHECK(cudaMalloc((void**)&d_MatA,memSize));
    CHECK(cudaMalloc((void**)&d_MatB,memSize));
    CHECK(cudaMalloc((void**)&d_MatC,memSize));

    cudaMemcpy(d_MatA, h_A, memSize, cudaMemcpyHostToDevice);
    cudaMemcpy(d_MatB, h_B, memSize, cudaMemcpyHostToDevice);

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

    start = seconds();
    SumMatGPU2D<<<grid, block>>>(d_MatA, d_MatB, d_MatC, nx, ny);
    cudaDeviceSynchronize();
    epls = seconds() - start;
    printf("SumMatGPU2D cost time : %f \n", epls);
    CHECK(cudaGetLastError());

    CHECK(cudaMemcpy(gpuRef, d_MatC, memSize, cudaMemcpyDeviceToHost));

    checkResult(hostRef, gpuRef, nxy);
    cudaFree(d_MatA);
    cudaFree(d_MatB);
    cudaFree(d_MatC);
    free(h_A);
    free(h_B);
    free(hostRef);
    free(gpuRef);

    cudaDeviceReset();
    return 0;
}

矩阵求和 - 1D Grid & 1D Block

对于一个Size为(nx,ny)的矩阵，使用1维的Grid和1维的Block意味着每个Thread要处理 ny个数据，相比上面2D Grid和2D Block的情况，全局索引映射会有很大不同。由于只有 threadIdx.x有用，因此还需要一个 for循环来处理每个线程中 ny元素：

__global__ void sumMatrixOnGPU1D(float *MatA, float *MatB, float *MatC, int nx, int ny){
    unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
    if (ix < nx ){
        for (int iy = 0; iy < ny; iy++){
            int idx = iy * nx + ix;
            MatC[idx] = MatA[idx] + MatB[idx];
        }
    }
}

矩阵求和 - 2D Grid & 1D Block

在这种情况下，可知 threadIdx.y = 0即数据都在第0行，且 blockDim.y = 1即1维Block，因此Kernel函数可写为：

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

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

参考文章

[1] NVIDIA 技术博客. (2022). CUDA 编程手册系列第二章: CUDA 编程模型概述. [online] Available at: https://developer.nvidia.com/zh-cn/blog/cuda-model-intro-cn/ [Accessed 14 Apr. 2024].
[2]docs.nvidia.com. (n.d.). CUDA Runtime API :: CUDA Toolkit Documentation. [online] Available at: https://docs.nvidia.com/cuda/cuda-runtime-api [Accessed 14 Apr. 2024].

‌
‌

Menu

Share

“我还要学CUDA编程！”（二）——CUDA编程模型

CUDA编程结构

CUDA内存管理

简单的数据传输示例

线程管理

调用Kernel函数

Kernel函数限定符

Kernel函数和C函数的区别

错误处理宏

并行线程

矩阵求和 - 2D Grid & 2D Block

矩阵求和 - 1D Grid & 1D Block

矩阵求和 - 2D Grid & 1D Block

参考文章

Comment

Ubuntu显卡驱动以及CUDA、cuDNN、TensorRT的安装

高达系列作品时间线

我还在学CUDA编程（四）——内存管理及访问模式

单目相机像素坐标系到世界坐标系的转换

Ubuntu下的常用软件安装和环境配置（持续更新）

单目相机的深度图，以及绝对深度和相对深度

点云数据 (Point Cloud Data)

继续学CUDA编程（三）——CUDA执行模型（上）

多目标跟踪算法的评价指标

多目标跟踪中的目标匹配算法