通于一而万事毕，无心得而鬼神服。

本文以数组加法为例，简述CUDA的编程模型

基本概念

确保了解《CUDA的线程结构》中的内容再进行下去，因为对于线程结构本章会省略一部分描述

编程模型

首先，编程模型是写程序时可被控制的部分，通过编程模型可以控制计算设备的行为；其次，编程模型对于编译器和库函数来说是一层通信抽象；

再次，对于GPU来说，编程模型可分为

• 核函数
• 内存管理
• 线程管理
• 流

最后，对于Cuda编程，可抽象出几个过程

• 领域层（业务主导要实现的目标）
• 逻辑层（基于逻辑组织线程代码）
• 硬件层（思考如何将线程映射到机器上，以及思考如何提升性能）

接下来，会从以下方面说明Cuda的编程结构

• 内存
• 线程

• 核函数

  • 启动
  • 编写
  • 验证
• 错误处理

内存管理

首先应该理解异构环境下内存的样貌：通过PCIe总线分割和通信的两种内存：主机内存和设备内存

一个异构环境，通常有多个CPU多个GPU，他们都通过PCIe总线相互通信，也是通过PCIe总线分隔开的。所以我们要区分一下两种设备的内存：

• 主机：CPU及其内存
• 设备：GPU及其内存

注意这两个内存从硬件到软件都是隔离的（CUDA6.0 以后支持统一寻址），我们目前先不研究统一寻址，我们现在还是用内存来回拷贝的方法来编写调试程序，以巩固大家对两个内存隔离这个事实的理解。

CUDA的内存管理是C风格的，列表如下：

另外，CUDA中设备内存也是分层次的，分为单网格中公用的全局内存和单线程块中公用的共享内存，大致如图：

线程管理

在《CUDA的线程结构》中，曾经用草图展示过cuda线程的组织形式，这里用更好的图展示一下：

另外，不同线程块的线程物理隔离无法相互影响，同一个线程块中的线程可以实现以下两种操作：

• 同步
• 共享内存

最后再重申一下《CUDA的线程结构》中提过的，如何标识一个线程：

• blockIdx（线程块在线程网格内的位置索引）
• threadIdx（线程在线程块内的位置索引）

这两个内置结构体基于 uint3 定义，包含三个无符号整数xyz的结构

还有blockIdx中三个字段的范围和threadIdx中三个字段的范围：

• gridDim
• blockDim

dim3类型(基于uint3定义的数据结构)的变量，也包含三个字段xyz

另外，线程块在3个轴上的个数和一个线程块的线程数是有限制的：

• 线程块在xyz三轴上分别最多：2^31-1,65535,65535
• 一个线程块线程总数最多1024，在三维上最多:1024,1024,64

核函数

之前hello world时已经用到了核函数，但是没有给一个明确定义，这里给出：

核函数就是在CUDA模型上诸多线程中运行的那段串行代码，这段代码在设备上运行，用NVCC编译，产生的机器码是GPU的机器码，所以我们写CUDA程序就是写核函数

• 第一步我们要确保核函数能正确的运行产生正切的结果
• 第二步优化CUDA程序的部分，无论是优化算法，还是调整内存结构，线程结构都是要调整核函数内的代码，来完成这些优化的。

接下来将以数组加法为例来实现和使用核函数，以及应用一些内存管理的函数

单机版数组加法

无需多言，直接上代码，明确我们要干什么：

/**
 * 单机版:两数组相加
 *
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
void add(const double *x,const double *y,double *z,const int N) {
    for (int i=0;i<N;i++) {
        z[i] = x[i]+y[i];
    }
}

// 函数2: 检查结果
void check(const double *z,const int N) {
    bool has_error = false;
    for (int i=0;i<N;i++) {
        if (fabs(z[i]-c) > EPSILON) {
            has_error = true;
            break;
        }
    }
    printf("%s\n",has_error?"Has errors":"No errors");
}

int main() {

    //创建数组xyz
    double* x = new double[N];
    double* y = new double[N];
    double* z = new double[N];

    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = a;
        y[i] = b;
        z[i] = 0;
    }

    //计算数组
    add(x,y,z,N);
    //检查结果
    check(z,N);

    //删除数组指针
    delete[] x;
    delete[] y;
    delete[] z;

    return 0 ;
}

CUDA程序和内存管理

先看代码，比起hello world，这是一个完整的CUDA程序，它包含了下图所示的结构：

/**
 * cuda版:两数组相加
 *
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
__global__ void add(const double *x,const double *y,double *z) {
    const int n  = blockDim.x*blockIdx.x + threadIdx.x;
    z[n] = x[n] + y[n];
}

// 函数2: 检查结果
void check(const double *z,const int N) {
    bool has_error = false;
    for (int i=0;i<N;i++) {
        if (fabs(z[i]-c) > EPSILON) {
            has_error = true;
            break;
        }
    }
    printf("%s\n",has_error?"Has errors":"No errors");
}

int main() {

    //主机: 创建数组xyz
    double* x = new double[N];
    double* y = new double[N];
    double* z = new double[N];

    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = a;
        y[i] = b;
        z[i] = 0;
    }

    //主机: 创建设备数组xd,yd

    double* xd;
    double* yd;

    cudaMalloc((void **)&xd,M);
    cudaMalloc((void **)&yd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
    cudaMemcpy(yd,y,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    //const dim3 block_size(128,1,1);
    const int block_size=128;
    const dim3 grid_size(N/block_size,1,1);

    //设备: 计算数组
    add<<<grid_size,block_size>>>(xd,yd,yd);

    //主机: 复制设备结果到主机
    cudaMemcpy(z,yd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(z,N);


    //删除数组指针
    cudaFree(xd);
    cudaFree(yd);

    delete[] x;
    delete[] y;
    delete[] z;

    return 0 ;
}

先把核函数的部分放在一边，代码里有以下几点可以注意：

• 隐式设备初始化：了解过opengl或d3d程序的可以明显感觉，CUDA是没有显式初始化设备(即GPU)的函数。在第一次调用一个和设备管理及版本查询功能无关的运行时 API函数时，设备将自动初始化。

• 内存管理一节已经大致了解了内存管理的相关函数，这里就是应用的地方，需要注意的是cudaMalloc和cudaMemcpy

  • cudaMalloc的第一个参数是指针的指针
  • cudaMemcpy的最后一个参数规定了数据传输的方式（从名称就可以看出）

CUDA程序和核函数

然后来看核函数，可分为三个方面：

• 启动核函数
• 编写核函数
• 验证核函数

首先，CPU是一个控制者，运行核函数，要从CPU发起，第一，核函数是异步的，第二，如《CUDA的线程结构》所说，三个尖括号<<<grid,block>>>内是对设备代码执行的线程结构的配置：

• 网格的布局（线程块）
• 线程块的布局（线程）

其次，同步函数有两种

• cudaDeviceSynchronize显式同步
• 涉及到设备端不执行完，主机没办法进行，比如内存拷贝函数，隐式同步

然后看到核函数的写法，cuda中不仅有核函数，还有一些其他类型，通过不同的标识表明，如核函数是__global__：

有些函数可以同时定义为 device 和 host ，这种函数可以同时被设备和主机端的代码调用，主机端代码调用函数很正常，设备端调用函数与C语言一致，但是要声明成设备端代码，告诉nvcc编译成设备机器码，同时声明主机端设备端函数，那么就要告诉编译器，生成两份不同设备的机器码。

其他方面核函数基本和普通函数一致，但有一些限制：

1. 除了统一寻址，只能访问设备内存
2. 必须有void返回类型
3. 不支持可变数量的参数
4. 不支持静态变量
5. 显示异步行为
6. 不可作为成员，但可以包装
7. 在动态并行的前提下，可以递归调用

另外设备函数可用__noinline__ __forceinline__来作为内联函数和非内联函数

CUDA小技巧，当我们进行调试的时候可以把核函数配置成单线程的：

下面是一些例子：

• 函数的使用
• 斐波那契数列，本程序意义不大，旨在展示__device__ __host__的共用
• 斐波那契数列，本程序意义不大，旨在展示递归的使用

/**
 * 函数的使用
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
//TODO: 注意到三种写法效率是一样的
__device__ double addfunc1(double x,double y) {
    return x+y;
}
__device__ double addfunc2(double x,double y,double *z) {
    *z = x+y;
}
__device__ double addfunc3(double x,double y,double &z) {
    z= x+y;
}
__global__ void add(const double *x,const double *y,double *z) {
    const int n  = blockDim.x*blockIdx.x + threadIdx.x;
    if (n<N) {
        z[n] = x[n] + y[n];
        //z[n] = addfunc1(x[n],y[n]);
        //addfunc2(x[n],z[n],&z[n]);
        //addfunc3(x[n],y[n],z[n]);
    }
}

/**
 * 斐波那契数列，本程序意义不大，旨在展示__device__ __host__的共用
 *
*/
#include<iostream>


const int N = 10;
const int M = sizeof(int)*N;

// 函数1: 斐波那契数列
__device__ __host__ void fibN(int* A, const int N) {
    if (N<0) return;
    if (A[N]==-1) {
        if (N>1) {
            fibN(A,N-1);
            fibN(A,N-2);
            A[N] = A[N-1] + A[N-2];
        }
        if (N<=1) {
            A[N]=N;
        }
    }
}

__global__ void fibNN(int* A, const int N) {
    fibN(A, N);
}


// 函数2: 检查结果
//0、1、1、2、3、5、8、13、21、34
void check(const int *A,const int N) {
    for (int i=0;i<N;i++) {
        std::cout<<A[i]<<" ";
    }
    std::cout<<std::endl;
}

int main() {

    //主机: 创建数组xyz
    int* x = new int[N];


    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = -1;
    }

    //主机: 创建设备数组xd,yd

    int* xd;
    cudaMalloc((void **)&xd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    const dim3 block_size(128,1,1);
    const dim3 grid_size(128,1,1);

    //设备: 计算数组
    fibNN<<<grid_size,block_size>>>(xd,N-1);


    //主机: 复制设备结果到主机
    cudaMemcpy(x,xd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(x,N);


    for (int i=0;i<N;i++) {
        x[i] = -1;
    }
    fibN(x,N-1);
    check(x,N);

    //删除数组指针
    cudaFree(xd);
    delete[] x;


    return 0 ;
}

/**
 * cuda版:两数组相加
 *
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
__global__ void add(const double *x,const double *y,double *z) {
    const int n  = blockDim.x*blockIdx.x + threadIdx.x;
    if (n<N) {
        z[n] = x[n] + y[n];
    }
}

// 函数2: 检查结果
void check(const double *z,const int N) {
    bool has_error = false;
    for (int i=0;i<N;i++) {
        if (fabs(z[i]-c) > EPSILON) {
            has_error = true;
            break;
        }
    }
    printf("%s\n",has_error?"Has errors":"No errors");
}

int main() {

    //主机: 创建数组xyz
    double* x = new double[N];
    double* y = new double[N];
    double* z = new double[N];

    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = a;
        y[i] = b;
        z[i] = 0;
    }

    //主机: 创建设备数组xd,yd

    double* xd;
    double* yd;

    cudaMalloc((void **)&xd,M);
    cudaMalloc((void **)&yd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
    cudaMemcpy(yd,y,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    //const dim3 block_size(128,1,1);
    const int block_size=128;
    const dim3 grid_size(N/block_size,1,1);

    std::cout << "Size: " << N/block_size <<"=="<<block_size <<std::endl;

    //设备: 计算数组
    //TODO: 注意到这里使用了会超过N的执行参数
    add<<<800000,block_size>>>(xd,yd,yd);

    //主机: 复制设备结果到主机
    cudaMemcpy(z,yd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(z,N);


    //删除数组指针
    cudaFree(xd);
    cudaFree(yd);

    delete[] x;
    delete[] y;
    delete[] z;

    return 0 ;
}

/**
 * 斐波那契数列，本程序意义不大，旨在展示递归的使用
 *
*/
#include<iostream>


const int N = 10;
const int M = sizeof(int)*N;

// 函数1: 斐波那契数列
__global__ void fibN(int* A, const int N) {
    if (N<0) return;
    if (A[N]==-1) {
        if (N>1) {
            fibN<<<1,1>>>(A,N-1);
            fibN<<<1,1>>>(A,N-2);
            A[N] = A[N-1] + A[N-2];
        }
        if (N<=1) {
            A[N]=N;
        }
    }
}



// 函数2: 检查结果
//0、1、1、2、3、5、8、13、21、34
void check(const int *A,const int N) {
    for (int i=0;i<N;i++) {
        std::cout<<A[i]<<" ";
    }
    std::cout<<std::endl;
}

int main() {

    //主机: 创建数组xyz
    int* x = new int[N];


    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = -1;
    }

    //主机: 创建设备数组xd,yd

    int* xd;
    cudaMalloc((void **)&xd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    const dim3 block_size(128,1,1);
    const dim3 grid_size(128,1,1);

    //设备: 计算数组
    fibN<<<grid_size,block_size>>>(xd,N-1);


    //主机: 复制设备结果到主机
    cudaMemcpy(x,xd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(x,N);


    //删除数组指针
    cudaFree(xd);
    delete[] x;


    return 0 ;
}

//0 1 1 0 -2 -2 -2 -4 -2 -6 结果错误 但是可见可以这么用

验证和错误检测

所有编程都需要对错误进行处理，早期的编码错误，编译器会帮我们搞定，内存错误也能观察出来，但是有些逻辑错误很难发现，甚至到了上线运行时才会被发现，而且有些厉害的bug复现会很难，不总出现，但是很致命，而且CUDA基本都是异步执行的，当错误出现的时候，不一定是哪一条指令触发的，这一点非常头疼；这时候我们就需要对错误进行防御性处理了

错误检测宏如下：

#pragma once
#include <cstdio>

#define CHECK(call)                                   \
do                                                    \
{                                                     \
    const cudaError_t error_code = call;              \
    if (error_code != cudaSuccess)                    \
    {                                                 \
        printf("CUDA Error:\n");                      \
        printf("    File:       %s\n", __FILE__);     \
        printf("    Line:       %d\n", __LINE__);     \
        printf("    Error code: %d\n", error_code);   \
        printf("    Error text: %s\n",                \
            cudaGetErrorString(error_code));          \
        exit(1);                                      \
    }                                                 \
} while (0)

• cudaGetErrorString可把error_code转为相应的提示信息

用这个宏函数包装大部分CUDA运行时API函数。有一个例外是cudaEventQuery()函数，因为它很有可能返回cudaErrorNotReady,但又不代表程序出错了

#include "error.cuh"
#include <cmath>
#include <cstdio>

const double EPSILON = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;
void __global__ add(const double *x, const double *y, double *z, const int N);
void check(const double *z, const int N);

int main(void)
{
    const int N = 100000000;
    const int M = sizeof(double) * N;
    double *h_x = (double*) malloc(M);
    double *h_y = (double*) malloc(M);
    double *h_z = (double*) malloc(M);

    for (int n = 0; n < N; ++n)
    {
        h_x[n] = a;
        h_y[n] = b;
    }

    double *d_x, *d_y, *d_z;
    CHECK(cudaMalloc((void **)&d_x, M));
    CHECK(cudaMalloc((void **)&d_y, M));
    CHECK(cudaMalloc((void **)&d_z, M));
    //TODO:此处错误复制方向搞反了
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyDeviceToHost));
    CHECK(cudaMemcpy(d_y, h_y, M, cudaMemcpyDeviceToHost));

    const int block_size = 128;
    const int grid_size = (N + block_size - 1) / block_size;
    add<<<grid_size, block_size>>>(d_x, d_y, d_z, N);

    CHECK(cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost));
    check(h_z, N);

    free(h_x);
    free(h_y);
    free(h_z);
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_z));
    return 0;
}

void __global__ add(const double *x, const double *y, double *z, const int N)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        z[n] = x[n] + y[n];
    }
}

void check(const double *z, const int N)
{
    bool has_error = false;
    for (int n = 0; n < N; ++n)
    {
        if (fabs(z[n] - c) > EPSILON)
        {
            has_error = true;
        }
    }
    printf("%s\n", has_error ? "Has errors" : "No errors");
}

用上述方法不能捕捉调用核函数的相关错误，因为核函数不返回任何值(回顾一下，核函数必须用void修饰)。有一个方法可以捕捉调用核函数可能发生的错误

#include "error.cuh"
#include <cmath>
#include <cstdio>

const double EPSILON = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;
void __global__ add(const double *x, const double *y, double *z, const int N);
void check(const double *z, const int N);

int main(void)
{
    const int N = 100000000;
    const int M = sizeof(double) * N;
    double *h_x = (double*) malloc(M);
    double *h_y = (double*) malloc(M);
    double *h_z = (double*) malloc(M);

    for (int n = 0; n < N; ++n)
    {
        h_x[n] = a;
        h_y[n] = b;
    }

    double *d_x, *d_y, *d_z;
    CHECK(cudaMalloc((void **)&d_x, M));
    CHECK(cudaMalloc((void **)&d_y, M));
    CHECK(cudaMalloc((void **)&d_z, M));
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice));


    //TODO:错误在线程块的大小太大
    const int block_size = 1280;
    const int grid_size = (N + block_size - 1) / block_size;
    add<<<grid_size, block_size>>>(d_x, d_y, d_z, N);
    CHECK(cudaGetLastError());
    CHECK(cudaDeviceSynchronize());

    CHECK(cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost));
    check(h_z, N);

    dim3 a;
    free(h_x);
    free(h_y);
    free(h_z);
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_z));
    return 0;
}

void __global__ add(const double *x, const double *y, double *z, const int N)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        z[n] = x[n] + y[n];
    }
}

void check(const double *z, const int N)
{
    bool has_error = false;
    for (int n = 0; n < N; ++n)
    {
        if (fabs(z[n] - c) > EPSILON)
        {
            has_error = true;
        }
    }
    printf("%s\n", has_error ? "Has errors" : "No errors");
}

CHECK(cudaGetLastError());
CHECK(cudaDeviceSynchronize());

其中，第一条语句的作用是捕捉第二个语句之前的最后一个错误，第二条语句的作用是同步主机与设备。之所以要同步主机与设备，是因为核函数的调用是异步的，即主机发出调用核函数的命令后会立即执行后面的语句，不会等待核函数执行完毕。