“我还要学CUDA编程!”(二)——CUDA编程模型

FunnyWii
FunnyWii
发布于 2024-04-21 / 26 阅读
0
0

“我还要学CUDA编程!”(二)——CUDA编程模型

CUDA编程结构

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

  1. 把数据从CPU内存(HOST)拷贝至GPU内存(DEVICE)
  2. 调用该Kernel函数,对DEVICE中的数据进行操作
  3. 将数据从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使用 blockIdxthreadIdx来互相区分,这两个变量是CUDA的内置变量,为 unit3类型,即包含3个 uint的结构,可以使用 x y z三个字段进行指定。比如 blockIdx.xthreadIdx.y

GridBlock的维度依靠 dim3类型的 gridDimblockDim来指定,这种类型的数据也可以通过 x y z来获得。

CUDA网格结构.png

在一段 .cu程序中,会定义两组这种变量。一个是在Host中定义的 dim3类型的 GridBlock,该变量仅在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=10block.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);

threadIdx.png

由上图也能看出,数据在全局内存中是线性存储的,所以可以通过 blockIdxthreadIdx获取一个线程的唯一ID,同时建立线程和数据元素间的映射关系。一般情况下,可以通过1维 GridBlock来计算全局数据访问的唯一索引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)输出可读信息。

并行线程

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

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].



评论