概述

CUDA编程模型提供了一个计算机架构抽象作为应用程序和其可用硬件之间的桥梁。通信抽象是程序与编程模型实现之间的分界线，它通过专业的硬件原语和操作系统的编译器或库来实现。利用编程模型所编写的程序指定了程序的各组成部分是如何共享信息及相互协作的。
编程模型可以理解为，我们要用到的语法，内存结构，线程结构等这些我们写程序时我们自己控制的部分，这些部分控制了异构计算设备的工作模式，都是属于编程模型。

除了与其他并行编程模型共有的抽象外，CUDA编程模型还利用GPU架构的计算能力提供了以下几个特有功能

• 一种通过层次结构在GPU中组织线程的方法
• 一种通过层次结构在GPU中访问内存的方法

以程序员的角度可以从以下几个不同的层面来看待并行计算：

• 领域层
• 逻辑层
• 硬件层

在编程与算法设计的过程中，你最关心的应是在领域层如何解析数据和函数，以便在并行运行环境中能正确、高效地解决问题。
进入编程阶段，你的关注点应转向如何组织并发线程。在这个阶段，你需要从逻辑层面来思考，以确保你的线程和计算能正确地解决问题。
在硬件层，通过理解线程是如何映射到核心可以帮助提高其性能。

CUDA编程模型

在一个异构环境中包含多个CPU和GPU，每个GPU和CPU的内存都由一条PCI-Express总线分隔开。

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

CUDA编程模型主要是异步的，因此在GPU上进行的运算可以与主机-设备通信重叠。串行代码在主机CPU 上执行，而并行代码在GPU上执行。一个典型的CUDA程序包括由并行代码互补的串行代码。主机代码按照ANSI C标准进行编写，而设备代码使用CUDA C进行编写。

内存管理

内存管理在传统串行程序是非常常见的，寄存器空间，栈空间内的内存由机器自己管理，堆空间由用户控制分配和释放，CUDA程序同样，只是CUDA提供的API可以分配管理设备上的内存，当然也可以用CDUA管理主机上的内存，主机上的传统标准库也能完成主机内存管理。
下面表格有一些主机API和CUDA C的API的对比：

设备内存分配：cudaMalloc
cudaError_t cudaMalloc(void **devPtr, size_t size);

• 在 GPU 显存中分配一段连续的线性内存
• devPtr 用来返回设备内存地址
• size 表示要分配的字节数
• 返回值类型是 cudaError_t，用于表示是否成功

主机与设备数据传输：cudaMemcpy

cudaError_t cudaMemcpy(
    void *dst,
    const void *src,
    size_t count,
    enum cudaMemcpyKind kind
);
示例
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);

这个函数以同步方式执行，因为在cudaMemcpy函数返回以及传输操作完成之前主机应用程序是阻塞的。
这意味着下面的代码中，第二行必须等第一行的数据传输完成后才会执行：

cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
kernel<<<grid, block>>>(d_a);

CUDA 错误处理
CUDA 运行时 API 通常返回 cudaError_t 类型。
常见返回值：

可以用 cudaGetErrorString 将错误码转换为可读信息：

cudaError_t err = cudaMalloc((void**)&d_a, bytes);

if (err != cudaSuccess) {
    printf("CUDA error: %s\n", cudaGetErrorString(err));
}

内存层次结构
在GPU内存层次结构中，最主要的两种内存是全局内存和共享内存。全局类似于CPU 的系统内存，而共享内存类似于CPU的缓存。

示例

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

int main() {
    int n = 1024;
    size_t bytes = n * sizeof(float);

    float *h_a = (float*)malloc(bytes);
    float *d_a = nullptr;

    // 1. 初始化主机数据
    for (int i = 0; i < n; i++) {
        h_a[i] = i;
    }

    // 2. 分配设备内存
    cudaMalloc((void**)&d_a, bytes);

    // 3. 主机数据复制到设备
    cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);

    // 4. 这里通常会启动 kernel
    // kernel<<<grid, block>>>(d_a);

    // 5. 设备数据复制回主机
    cudaMemcpy(h_a, d_a, bytes, cudaMemcpyDeviceToHost);

    // 6. 释放内存
    cudaFree(d_a);
    free(h_a);

    return 0;
}

线程管理

Grid、Block、Thread 的关系

可以理解为：一个 Grid 里面有很多 Block，一个 Block 里面有很多 Thread。
线程块内外的协作关系
同一个线程块内的线程可以协作，主要方式有：

CUDA 内置线程变量
CUDA 在核函数中提供了一些内置变量，用来标识当前线程的位置。

这些变量都可以通过 .x、.y、.z 访问三个维度。

threadIdx.x
blockIdx.x
blockDim.x
gridDim.x

一维线程索引计算

int i = blockIdx.x * blockDim.x + threadIdx.x;

每个线程通过不同的 i 处理数组中的不同元素。
为什么要计算全局线程 ID
GPU 中每个线程执行的是同一份代码，所以必须让每个线程知道：
我是谁？我应该处理哪一份数据？
这就是 blockIdx、threadIdx 和 blockDim 的作用。
设blockDim.x = 4;

dim3 类型
CUDA 使用 dim3 类型定义 Grid 和 Block 的维度。
dim3 block(256);
dim3 grid(4);
这表示：

• 每个 Block 有 256 个线程
• Grid 中有 4 个 Block
• 总线程数是 4 * 256 = 1024

dim3 可以是一维、二维或三维：

dim3 block1D(256);
dim3 block2D(16, 16);
dim3 block3D(8, 8, 4);

根据数据规模计算 Grid 大小
实际编程中，通常先确定 Block 大小，再根据数据总量计算 Grid 大小。

int n = 1000;
int blockSize = 256;
int gridSize = (n + blockSize - 1) / blockSize;

为什么这样写？
因为 n 不一定刚好能被 blockSize 整除。
这个公式可以向上取整，保证线程数量足够覆盖所有数据。

由于可能会产生多余线程，所以 kernel 中通常要加边界判断：

__global__ void kernel(float *data, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < n) {
        data[i] = data[i] * 2.0f;
    }
}

启动一个核函数

执行配置 <<<grid, block>>>

kernel<<<grid, block>>>(args);

我们可以使用dim3类型的grid维度和block维度配置内核，也可以使用int类型的变量，或者常量直接初始化：

kernel_name<<<4,8>>>(argument list);

上面这条指令的线程布局是：

不同线程布局示例
假设有 32 个数据元素需要处理，可以有多种线程组织方式。

方式一：4 个 Block，每个 Block 8 个 Thread
kernel<<<4, 8>>>();

方式二：1 个 Block，每个 Block 32 个 Thread
kernel<<<1, 32>>>();

方式三：32 个 Block，每个 Block 1 个 Thread
kernel<<<32, 1>>>();

CUDA 核函数启动是异步的
普通 C 函数调用通常是同步的：

foo();
bar();  // foo 执行完后才会执行 bar

但 CUDA 核函数启动是异步的：

kernel<<<grid, block>>>(args);

当 kernel 被启动后，控制权会很快返回给在这里插入代码片 CPU。
这并不代表 GPU 已经执行完 kernel。
如果希望 CPU 等待 GPU 执行完成，可以使用：

cudaDeviceSynchronize();

示例：

kernel<<<grid, block>>>(args);
cudaDeviceSynchronize();

cudaMemcpy 的隐式同步
虽然 kernel 启动是异步的，但某些 CUDA API 会产生隐式同步。
例如：

cudaMemcpy(h_data, d_data, bytes, cudaMemcpyDeviceToHost);

当从 GPU 拷贝数据回 CPU 时，主机端通常会等待之前的 kernel 执行完成，并等待数据传输完成。
也就是说：

kernel<<<grid, block>>>(d_data);
cudaMemcpy(h_data, d_data, bytes, cudaMemcpyDeviceToHost);

这里 cudaMemcpy 会等待 kernel 完成后再复制数据
完整示例

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

__global__ void addOne(float *data, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < n) {
        data[i] += 1.0f;
    }
}

int main() {
    int n = 1000;
    size_t bytes = n * sizeof(float);

    float *h_data = (float*)malloc(bytes);
    float *d_data = nullptr;

    for (int i = 0; i < n; i++) {
        h_data[i] = i;
    }

    cudaMalloc((void**)&d_data, bytes);
    cudaMemcpy(d_data, h_data, bytes, cudaMemcpyHostToDevice);

    int blockSize = 256;
    int gridSize = (n + blockSize - 1) / blockSize;

    addOne<<<gridSize, blockSize>>>(d_data, n);

    cudaMemcpy(h_data, d_data, bytes, cudaMemcpyDeviceToHost);

    cudaFree(d_data);
    free(h_data);

    return 0;
}

编写核函数

我们会启动核函数了，但是核函数哪里来的？当然我们写的，核函数也是一个函数，但是声明核函数有一个比较模板化的方法：

__global__ void kernel_name(argument list);

CUDA 函数限定符

其中

__host__ __device__ float add(float a, float b) {
    return a + b;
}

核函数的限制
CUDA 核函数有一些限制：

为什么核函数里没有 for 循环
CPU 代码中：

for (int i = 0; i < N; i++) {
    C[i] = A[i] + B[i];
}

表示一个线程依次处理所有元素。
GPU 代码中：

int i = blockIdx.x * blockDim.x + threadIdx.x;
C[i] = A[i] + B[i];

表示每个线程只处理一个元素。
如果启动了 N 个线程，那么这 N 个线程就可以并行处理 N 个数组元素。

为什么要验证核函数

CUDA 核函数运行在 GPU 上，并且通常由大量线程并行执行。
因此，核函数写完后不能只看程序能否运行，还需要确认：

• 计算结果是否正确
• 线程索引是否正确
• 内存访问是否越界
• CPU 与 GPU 的计算结果是否一致

验证核函数的基本思路是：
用 CPU 端的串行结果作为参考结果，再与 GPU 端结果进行对比。

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

//GPU 端核函数：
__global__ void sumArraysOnGPU(float *A, float *B, float *C, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < N) {
        C[i] = A[i] + B[i];
    }
}

//然后比较 CPU 结果和 GPU 结果：
void checkResult(float *hostRef, float *gpuRef, int N) {
    double epsilon = 1.0E-8;

    for (int i = 0; i < N; i++) {
        if (abs(hostRef[i] - gpuRef[i]) > epsilon) {
            printf("Results do not match!\n");
            printf("host %f gpu %f at index %d\n", hostRef[i], gpuRef[i], i);
            return;
        }
    }

    printf("Check result success!\n");
}

浮点数计算可能存在微小误差，因此不要直接写：

if (hostRef[i] != gpuRef[i])

更推荐使用误差范围：

if (abs(hostRef[i] - gpuRef[i]) > epsilon)

处理错误

普通 C 程序里，函数调用通常是同步的，哪里错了比较容易定位。但 CUDA 不一样，kernel 启动后 CPU 会马上继续往下执行，GPU 端真正的错误可能过一会儿才暴露出来。比如 kernel 里访问越界，可能不会在 kernel 启动那一行立刻报错，而是在后面的 cudaMemcpy 或 cudaDeviceSynchronize 时才被发现。
定义一个错误处理宏封装所有的CUDA API调用，这简化了错误检查过程：

#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);                                                           \
    }                                                                      \
}

之后每次调用 CUDA API 时，不直接写：

cudaMalloc((void**)&d_A, bytes);

而是写成：

CHECK(cudaMalloc((void**)&d_A, bytes));

这样如果 cudaMalloc 失败，程序会输出错误发生的文件名、行号、错误码和可读的错误原因。

编译执行

nvcc test.cu -o test
./test

给和函数计时

CUDA 程序的性能分析不能只看整个程序运行时间，因为整个程序里包含很多部分：

CPU 初始化数据
主机内存分配
GPU 显存分配
Host -> Device 数据传输
kernel 执行
Device -> Host 数据传输
结果验证
内存释放`

而给核函数计时关注的是其中最核心的一段：
kernel 在 GPU 上真正执行了多久

第一种方式是用 CPU 计时器

double cpuSecond() {
    struct timeval tp;
    gettimeofday(&tp, NULL);  //获取当前系统时间
    return ((double)tp.tv_sec + (double)tp.tv_usec * 1.e-6);
}


double iStart = cpuSecond();

sumArraysOnGPU<<<grid, block>>>(d_A, d_B, d_C, nElem);

cudaDeviceSynchronize();

double iElaps = cpuSecond() - iStart;
printf("sumArraysOnGPU elapsed %f sec\n", iElaps);

这里最重要的是：
cudaDeviceSynchronize();
因为 CUDA kernel 启动是异步的。也就是说，CPU 执行完：
sumArraysOnGPU<<<grid, block>>>(...);
之后会立刻往下走，并不会自动等待 GPU 把 kernel 执行完。
如果不加cudaDeviceSynchronize()，你测到的可能只是kernel 启动耗时，而不是kernel 执行耗时。

示例

//sumArraysOnGPU-timer.cu
#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <chrono>

#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);                                                           \
    }                                                                      \
}

double cpuSecond() {
    using namespace std::chrono;
    return duration<double>(high_resolution_clock::now().time_since_epoch()).count();
}

void initialData(float *ip, int size) {
    for (int i = 0; i < size; i++) {
        ip[i] = (float)(rand() & 0xFF) / 10.0f;
    }
}

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

void checkResult(float *hostRef, float *gpuRef, const int N) {
    double epsilon = 1.0E-8;

    for (int i = 0; i < N; i++) {
        if (fabs(hostRef[i] - gpuRef[i]) > epsilon) {
            printf("Results do not match!\n");
            printf("host %f gpu %f at index %d\n", hostRef[i], gpuRef[i], i);
            return;
        }
    }

    printf("Check result success!\n");
}

__global__ void sumArraysOnGPU(float *A, float *B, float *C, const int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < N) {
        C[i] = A[i] + B[i];
    }
}

int main() {
    int nElem = 1 << 24;   // 16M 个元素
    size_t nBytes = nElem * sizeof(float);

    printf("Vector size: %d\n", nElem);

    float *h_A = (float *)malloc(nBytes);
    float *h_B = (float *)malloc(nBytes);
    float *hostRef = (float *)malloc(nBytes);
    float *gpuRef = (float *)malloc(nBytes);

    initialData(h_A, nElem);
    initialData(h_B, nElem);

    memset(hostRef, 0, nBytes);
    memset(gpuRef, 0, nBytes);

    double iStart = cpuSecond();
    sumArraysOnHost(h_A, h_B, hostRef, nElem);
    double iElaps = cpuSecond() - iStart;

    printf("sumArraysOnHost elapsed %f sec\n", iElaps);

    float *d_A = NULL;
    float *d_B = NULL;
    float *d_C = NULL;

    CHECK(cudaMalloc((void **)&d_A, nBytes));
    CHECK(cudaMalloc((void **)&d_B, nBytes));
    CHECK(cudaMalloc((void **)&d_C, nBytes));

    CHECK(cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice));

    int blockSize = 1024;
    int gridSize = (nElem + blockSize - 1) / blockSize;

    dim3 block(blockSize);
    dim3 grid(gridSize);

    iStart = cpuSecond();

    sumArraysOnGPU<<<grid, block>>>(d_A, d_B, d_C, nElem);

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

    iElaps = cpuSecond() - iStart;

    printf("sumArraysOnGPU <<<%d, %d>>> elapsed %f sec\n",
           grid.x, block.x, iElaps);

    CHECK(cudaMemcpy(gpuRef, d_C, nBytes, cudaMemcpyDeviceToHost));

    checkResult(hostRef, gpuRef, nElem);

    CHECK(cudaFree(d_A));
    CHECK(cudaFree(d_B));
    CHECK(cudaFree(d_C));

    free(h_A);
    free(h_B);
    free(hostRef);
    free(gpuRef);

    CHECK(cudaDeviceReset());

    return 0;
}

nvcc ./sumArraysOnGPU-timer.cu -o sumArraysOnGPU-timer
./sumArraysOnGPU-timer
输出：
Vector size: 16777216
sumArraysOnHost elapsed 0.032108 sec
sumArraysOnGPU <<<16384, 1024>>> elapsed 0.000966 sec
Check result success!

第二种方式是用 nvprof 计时
普通运行程序：
./sumArraysOnGPU-timer
使用 nvprof：
nvprof ./sumArraysOnGPU-timer
nvprof 可以显示更详细的信息，比如：

kernel 执行时间
Host -> Device 拷贝时间
Device -> Host 拷贝时间
CUDA API 调用时间

组织并行线程

前面向量加法是一维数据，所以用一维 grid 和一维 block 就够了。现在进入矩阵加法，数据变成二维结构，于是就要思考：二维矩阵应该用二维线程组织，还是仍然用一维线程组织？

这一节通过矩阵加法比较三种方式：

1. 二维网格 + 二维线程块
2. 一维网格 + 一维线程块
3. 二维网格 + 一维线程块、

先看矩阵在内存中的存储方式。虽然矩阵是二维的，比如：
A[row][col]
但是在 GPU 全局内存中，它仍然是线性存储的，通常采用行优先方式：
index = row * nx + col
其中：

nx = 矩阵列数
ny = 矩阵行数

如果用变量表示：

int ix = 列索引;
int iy = 行索引;
int idx = iy * nx + ix;

这里 idx 就是二维矩阵元素在线性内存中的位置。

第一种方式：二维网格 + 二维线程块。
这是最直观的方式，因为矩阵是二维的，所以线程也按二维组织。
线程索引映射如下：

unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
unsigned int idx = iy * nx + ix;

核函数大概是：

__global__ void sumMatrixOnGPU2D(float *A, float *B, float *C, 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) {
        C[idx] = A[idx] + B[idx];
    }
}

执行配置类似：

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

sumMatrixOnGPU2D<<<grid, block>>>(d_A, d_B, d_C, nx, ny);

第二种方式：一维网格 + 一维线程块

这种方式只在 x 方向组织线程，每个线程负责矩阵中的一列，在线程内部用循环处理这一列的所有行。
线程索引：

unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;

核函数：

__global__ void sumMatrixOnGPU1D(float *A, float *B, float *C, 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;
            C[idx] = A[idx] + B[idx];
        }
    }
}

执行配置：

dim3 block(256);
dim3 grid((nx + block.x - 1) / block.x);

sumMatrixOnGPU1D<<<grid, block>>>(d_A, d_B, d_C, nx, ny);

这种方式的特点是：线程数量少了，每个线程做的事情变多了。
它不再是“一个线程处理一个元素”，而是“一个线程处理一列中的多个元素”。
这种方式有时候性能并不差，但它的并行度明显受限，因为总线程数量只和 nx 有关，而不是和 nx * ny 有关。

第三种方式：二维网格 + 一维线程块。

它使用二维 grid，但每个 block 仍然是一维的。可以理解为：

• blockIdx.x 负责列方向上的块划分
• blockIdx.y 负责行方向
• threadIdx.x 负责当前行中的列偏移

unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int iy = blockIdx.y;
unsigned int idx = iy * nx + ix;

核函数：

__global__ void sumMatrixOnGPUMix(float *A, float *B, float *C, 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) {
        C[idx] = A[idx] + B[idx];
    }
}

执行配置：

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

sumMatrixOnGPUMix<<<grid, block>>>(d_A, d_B, d_C, nx, ny);

这种方式仍然是一个线程处理一个矩阵元素，但是 block 只在 x 方向展开，grid 在 x 和 y 两个方向展开。
它可以看作是二维网格、二维块的一种特殊情况：block.y = 1。

设备管理

程序内部查询 GPU
最常用的是这几个 Runtime API：

//获取当前系统中 CUDA 可见设备的数量。
cudaError_t cudaGetDeviceCount(int *count); 

//获取指定 GPU 的设备属性。
cudaError_t cudaGetDeviceProperties(cudaDeviceProp *prop, int device);
    prop.name                  // GPU 名称
    prop.major                 // 计算能力主版本
    prop.minor                 // 计算能力次版本
    prop.totalGlobalMem        // 全局显存大小
    prop.multiProcessorCount   // SM 数量
    prop.maxThreadsPerBlock    // 每个线程块最大线程数
    prop.maxGridSize           // Grid 每个维度的最大尺寸
    prop.maxThreadsDim         // Block 每个维度的最大尺寸
    prop.sharedMemPerBlock     // 每个 Block 可用共享内存
    prop.warpSize              // warp 大小，通常是 32

//设置当前 CPU 线程后续 CUDA 操作使用哪一张 GPU。
cudaError_t cudaSetDevice(int device);

//获取当前 CPU 线程正在使用的 GPU 设备编号
cudaError_t cudaGetDevice(int *device);

如何选择 GPU
选择标准包括：

计算能力 compute capability
SM 数量
显存容量 totalGlobalMem
是否支持统一内存 managedMemory
是否支持并发 kernel concurrentKernels
是否支持异步拷贝 asyncEngineCount
当前 GPU 是否空闲
多 GPU 之间的拓扑关系

示例

int bestDevice = 0;
int bestScore = -1;

for (int dev = 0; dev < deviceCount; dev++) {
    cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, dev);

    int score = prop.multiProcessorCount * 100 + prop.major * 10 + prop.minor;

    if (score > bestScore) {
        bestScore = score;
        bestDevice = dev;
    }
}

cudaSetDevice(bestDevice);

命令行查询 GPU

nvidia-smi

列出所有 GPU：
nvidia-smi -L

查看某张卡的详细信息：
nvidia-smi -i 0 -q