CUDA编程模型

一、CUDA编程模型概述

并行计算的三层：

领域层，在编程和算法设计时考虑如何解析数据和函数；

逻辑层，在编程实现时确保线程和计算可以正确解决问题；

硬件层，通过理解线程如何映射到其核心从而提高性能。

1.CUDA编程结构

        1.分配GPU内存

        2.从CPU内存中拷贝数据到GPU内存

        3.调用CUDA的kernel函数完成运算

        4.将数据从GPU拷贝回CPU

        5.释放GPU空间 

2. 内存管理

cudaMalloc：向GPU设备分配一定字节的线性内存，并返回指针。

cudaMemcpy：负责主机和设备直之间的数据传输。此函数以同步方式执行，在cudaMemcpy函数返回以及传输操作完成之前主机应用程序阻塞。

主要需要注意cudaMemcpy函数

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

主机到设备：cudaMemcpy(d_A,h_A,nBytes,cudaMemcpyHostToDevice)

设备到主机：cudaMemcpy(h_A,d_A,nBytes,cudaMemcpyDeviceToHost)

主机到主机：cudaMemcpy(d_A,h_A,nBytes,cudaMemcpyHostToHost)

设备到设备：cudaMemcpy(h_A,d_A,nBytes,cudaMemcpyDeviceToDevice)

CUDA的内存主要分为两部分主机端和设备端，而且是具有层次，简单如下图所示，主要由全局内存和共享内存组成。

3. CUDA 执行模型

3.1 CUDA内核由一组线程执行

所有线程运行相同的代码（SPMD）， 每个线程都有一个ID，用于计算内存地址和做出控制决策。

 Block IDs 和Thread IDs。每个线程使用ID确定其作用于哪块数据。Block ID（1D，2D，3D），Thread ID（1D，2D，3D）。在核函数调用时需要 进行初始化i。

此处需要注意的是Block是一个组织线程的方式，并无明确硬件与其对应。

3.2 线程层次结构

Grid：一组由一个kernel启动所产生的线程block，简单说一个kernel会对应一个Grid，在全局内存中共享数据，在运行时动态调度。

Block： 协作线程阵列（CTA），线程之间的同步，在共享内存中共享数据，1D、2D、3D， 最多512个线程（看CUDA版本）。

kernel：在线程上运行的代码部分。

通常一个grid会被组成成二维，block会被组织为三维。 

一个block内的线程可以通过共享内存、原子操作、同步等进行协作，不同block的线程则无法同步。因为不同block的线程可能分布在不同的SM中 。

关于网格和块的大小设置见下一节。

3.3 线程的执行

        真正运行时一个warp的32个线程（类似于老式织布机的一排，织一次有32根线）运行在一个SM上，它们共享指令，每4个周期执行一条warp指令而且由SM进行动态调度。

 4. CUDA编程基础

4.1  启动CUDA核函数

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

其中第一个值为网格维度，也就是所启动的block数量，第二个值为每个块中的线程数量。通过指定这两个值配置内核中的线程数量和线程的布局。

核函数的执行是异步的，可以使用cudaDevicesSynchronize（void）来进行强制主机程序等待。 

 4.2 编写核函数

 __host__和__device__可以一起使用，这样函数在设备端、主机端均可调用。

 主机端的C向量加法：

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

核函数的向量加法：

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

二者相比，核函数中的循环体消失了，换用了内置的线程坐标进行索引。如果每个向量有32个元素可以考虑sumArrayOnGPU<<<1,32>>>(float *A, float *B, float *C)进行调用。

4.3 处理错误

验证核函数：

        （1）在核函数中使用printf（Fermi版本以上）

        （2）将执行参数设置为<<<1,1>>>，模拟了串行执行程序，可用于调试和验证结果

检查错误：

        #define CHECK(call)

如：CHECK(cudaDeviceSYnchronize())，检查核函数中哪一步就在哪一步外使用CHECK。

4.4 示例程序：向量加法

#include <cuda_runtime.h>
#include <stdio.h>
 
//错误处理宏
#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);\
    }\
}
 
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 current %d\n", hostRef[i], gpuRef[i], i);
            break;
        }
    }
    if (match) 
        printf("Array match.\n\n");
}
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 sumArraysOnGPU(float *A, float *B, float *C)
{
    int i = threadIdx.x;
    C[i] = A[i]+B[i];
    
}
 
int main()
{
    //printf("%s Starting...\n", argv[0]);
    
    //设置设备，有多个GPU的情况下
    int dev = 0;
    cudaSetDevice(dev);
    
    //设置向量的大小，此处设为32
    int nElem = 32;
    printf("Vector size: %d\n", nElem);
    //分配主机内存
    size_t nBytes = nElem * sizeof(float);
    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);
    
    //对向量A、B进行初始化
    int i = 0;
    for (; i<nElem; i++){
        h_A[i] = i;
        h_B[i] = i;
    }
    memset(hostRef, 0, nBytes);
    memset(gpuRef, 0, nBytes);
    
    //分配GPU内存
    float *d_A, *d_B, *d_C;
    cudaMalloc((float**)&d_A, nBytes);
    cudaMalloc((float**)&d_B, nBytes);
    cudaMalloc((float**)&d_C, nBytes);
    
    //将数据从CPU内存传输至GPU内存
    cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);
    
    //设置kernel函数的块数和线程数,此处为一个块，32线程
    dim3 block(nElem);
    dim3 grid(nElem / block.x);
    
    //调用kernel函数进行执行
    printf("Execution configuration <<<%d, %d>>>\n", grid.x, block.x);
    sumArraysOnGPU<<<grid, block>>>(d_A, d_B, d_C);
    
    //执行结束后将结果从GPU内存传输回CPU内存
    cudaMemcpy(gpuRef, d_C, nBytes, cudaMemcpyDeviceToHost);
    
    //在CPU同样执行一次加法，二者对比观察是否准确
    sumArraysOnHost(h_A, h_B, hostRef, nElem);
    
    //检查二者结果是否相同
    checkResult(hostRef, gpuRef, nElem);
    
    //释放GPU内存
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    
    //释放CPU内存
    free(h_A);
    free(h_B);
    free(hostRef);
    free(gpuRef);
    
    return 0;
    
    
}

运行结果：

 4.5 运行时间计算

第一种方法是在CPU使用cpuSecond（）函数在核函数上下，同时加上cudaDeviceSynchronize（）来计算时间。

第二种方法是使用nvprof工具

//查看帮助
nvprof --help
 
//测试内核
nvprof ./01-hello

5. 组织并行线程

这一节介绍每一个线程是怎么确定唯一的索引，然后建立并行计算，并且不同的线程组织形式是怎样影响性能的：

• 二维网格二维线程块
• 一维网格一维线程块
• 二维网格一维线程块

5.1 使用块和线程建立矩阵索引

对于一个给定的线程，首先可以通过将线程和块索引映射到矩阵坐标上来获取线程块和线程索引的全局内存偏移量，然后将这些矩阵坐标映射到全局内存的存储单元中。

第一步，线程和块索引映射到矩阵坐标上

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

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

第二步，将矩阵坐标映射至全局内存的索引上

        idx = iy   *   nx   +   ix

 这样我们就得到了每个线程的唯一标号，并且在运行时kernel是可以访问这个标号的。前面讲过CUDA每一个线程执行相同的代码，也就是异构计算中说的多线程单指令，如果每个不同的线程执行同样的代码，又处理同一组数据，将会得到多个相同的结果，显然这是没意义的，为了让不同线程处理不同的数据，CUDA常用的做法是让不同的线程对应不同的数据，也就是用线程的全局标号对应不同组的数据。

全局内存中线性存储，6行8列：

我们要做管理的就是：

• 线程和块索引（来计算线程的全局索引）
• 矩阵中给定点的坐标（ix,iy）
• (ix,iy)对应的线性内存的位置

代码检查： 

#include <cuda_runtime.h>
#include <stdio.h>
 
#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);\
    }\
}
 
void printMatrix(int *C, const int nx, const int ny)
{
    int *ic = C;
    printf("\nMatrix: (%d.%d)\n", nx, ny);
 
    for (int iy = 0; iy < ny; iy++)
    {
        for (int ix = 0; ix < nx; ix++)
        {
            printf("%3d", ic[ix]);
 
        }
 
        ic += nx;
        printf("\n");
    }
 
    printf("\n");
    return;
}
 
__global__ void printThreadIndex(int *A, const int nx, const int ny)
{
    int ix = threadIdx.x + blockIdx.x * blockDim.x;
    int iy = threadIdx.y + blockIdx.y * blockDim.y;
    unsigned int idx = iy * nx + ix;
 
    printf("thread_id (%d,%d) block_id (%d,%d) coordinate (%d,%d) global index"
           " %2d ival %2d\n", threadIdx.x, threadIdx.y, blockIdx.x, blockIdx.y,
           ix, iy, idx, A[idx]);
}
 
int main(int argc, char **argv)
{
    printf("%s Starting...\n", argv[0]);
 
    // 获取可用的GPU设别信息
    int dev = 0;
    cudaDeviceProp deviceProp;
    CHECK(cudaGetDeviceProperties(&deviceProp, dev));
    printf("Using Device %d: %s\n", dev, deviceProp.name);
    CHECK(cudaSetDevice(dev));
 
    // 设置矩阵维度8X6
    int nx = 8;
    int ny = 6;
    int nxy = nx * ny;
    int nBytes = nxy * sizeof(float);
 
    // 分配主机内存
    int *h_A;
    h_A = (int *)malloc(nBytes);
 
    // 初始化数据
    for (int i = 0; i < nxy; i++)
    {
        h_A[i] = i;
    }
    printMatrix(h_A, nx, ny);
 
    // GPU分配内存
    int *d_MatA;
    CHECK(cudaMalloc((void **)&d_MatA, nBytes));
 
    // 从CPU内存传输数据至GPU内存
    CHECK(cudaMemcpy(d_MatA, h_A, nBytes, cudaMemcpyHostToDevice));
 
    // 设置执行配置。两行四列的block
    dim3 block(4, 2);
    dim3 grid((nx + block.x - 1) / block.x, (ny + block.y - 1) / block.y);
 
    // 执行kernel函数
    printThreadIndex<<<grid, block>>>(d_MatA, nx, ny);
    CHECK(cudaGetLastError());
 
    // 释放内存
    CHECK(cudaFree(d_MatA));
    free(h_A);
 
    // reset device
    CHECK(cudaDeviceReset());
 
    return (0);
}

运行结果：

以其中threadIdx = (0, 0), blockIdx=  (0,1)为例。

        ix = 0 + 0 * 2 = 0

        iy = 0 + 1 * 2 = 2

在矩阵中为（0，2）

将矩阵坐标映射至全局内存的索引上

        idx = 2 * 8 + 0 = 16

5.2 使用二维网格和二维块对矩阵求和

在5.1的代码上进行修改

    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);
    cudaDeviceSynchronize();
    iElaps = seconds() - iStart;

void sumMatrixOnHost(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;
    }
 
    return;
}

// 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];
}

5.3 使用一维网格和一维块对矩阵求和

配置如下：

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

__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];
        }
}

5.4使用二维网格一维块 

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

// 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];
}

 总结：用不同的线程组织形式会得到正确结果，但是效率有所区别。从此例子中可以看出：

• 改变执行配置（线程组织）能得到不同的性能
• 传统的核函数可能不能得到最好的效果
• 一个给定的核函数，通过调整网格和线程块大小可以得到更好的效果