CUDA 学习笔记

Posted on Edited on

参考资料：

获得 GPU 加速的关键

总结

必要条件

  1. 数据传输比例较小
  2. 核函数的算术强度较高
  3. 核函数中定义的线程数目较多

提高性能的技巧

  1. 减少主机与设备之间的数据传输
  2. 提高核函数的算术强度
  3. 增大核函数的并行规模

SM 流多处理器

官方有时简称 Multiprocessor

SM 的构成

  1. 寄存器
  2. 共享内存
  3. 常量内存的缓存
  4. 纹理和表面内存的缓存
  5. L1 缓存
  6. （常见4个）线程束调度器
  7. 执行核心：INT32, FP32, FP64, 单精度浮点数超于函数的特殊函数单元(special function unit, SFUs), 混合精度的张量核心(tensor cores)

SM 的占有率

当并行规模较小，有些 SM 占有率为零，导致程序性能低下。当并行规模足够大时，也有可能得到非 100% 的占有率。

考虑指标（查询官方文档 CUDA_Occupancy_Calculator.xls 的图灵架构）：

  1. 一个 SM 最多拥有线程块个数$\ N_b = 16$。
  2. 一个 SM 最多拥有的线程个数$\ N_t = 1536$。
  3. 线程块大小最大为$\ 1024$。

当并行规模足够大（核函数配置的总线程数足够多），需要分几种情况讨论 SM 的理论占有率：

  1. 寄存器核共享内存使用量很少。SM 的占有率完全由线程块大小决定。首先，由于线程束大小是 32，线程块大小最好是 32 的倍数。其次，线程块大小不能小于$\ N_t/N_b$​，才可能利用最大线程数量。因此，96 的线程块大小就能获得 100% 的占有率。类似的，128 的线程块大小在开普勒架构下能获得满占有率。（我选择 128）
  2. 寄存器带来占有率瓶颈。一个 SM 最多有 65536（64K）个寄存器。如果令线程数最大化（1536），那么平均一个线程可用 42.7 个寄存器。当每个线程所用的寄存器个数大于 42 时，SM 的占有率将低于 50%。
  3. 有限的共享内存对占有率的约束。一个 SM 拥有 65536 字节共享内存。如果线程块大小为 128，且SM 是线程满载（1536），要令占有率为 100%，则网格大小为 12，一个线程块最多有 5461 字节的共享内存。

运行时API查询设备

参考deviceQuery.cpp

全局内存的合理使用

合并访问

全局内存具有最高的延迟，所以配置有缓存，每次 Cache Miss 默认读取 32 字节。由于缓存命中问题，需要提高内存访问的合并度，从而提高效率。简单起见，只考虑全局内存到 L2 缓存。可以定义合并度(degree of coalescing) 为一个线程束请求的字节数除以实际数据传输处理的字节数。合并度体现了显存带宽的利用率。

CUDA 运行时 API 函数分配的内存的首地址至少是256字节的整数倍。

只要确保每当线程束读取一次内存后，该内存对应的32字节范围都被当前线程束利用，则合并度就是100%。

对于一个线程束访问同一个全局内存地址的 4 字节，属于广播式的非合并访问（因为只用了 4 字节而非 32 字节，合并度为 12.5%）。如果内存只读，那么适合采用常量内存。

写优先于只读

有时候合并写入与合并读取无法兼顾，例如矩阵转置。此时应当优先保证写入是合并的。

这是因为从帕斯卡架构开始，若编译器判定一个全局内存变量在整个核函数的范围内都只读，则会自动用函数__ldg()读取全局内存，附带缓存效果，缓解非合并访问的影响。但对于全局内存的写入，没有类似的优化手段。

共享内存的合理使用

在核函数内，使用__shared__修饰的变量（数组）将使用共享内存。同一线程块使用同一个共享内存副本，而不同线程块的共享内存是互相独立的，不可见的。

利用共享内存，可以消除一些无法兼顾读写合并的全局内存访问，例如数组归约，先将要读取的数据从全局内存拷贝至共享内存（然后__syncthreads()同步一下），可以带来细微的性能提升。

对共享内存的访问越频繁，性能提升越明显。

动态共享内存

共享内存大小可以在核函数的执行配置中指定：

1
<<<grid_size, block_size, sizeof(real) * block_size>>>

同时要修改共享内存变量的声明方式

1
2
3
4
5
// before
__shared__ real a[128];

// after
extern __shared__ real a[];

动态和静态的声明方式在性能上几乎没有差别。

避免共享内存的bank冲突

共享内存在物理上被分为 32 个（恰好是线程束大小）个 bank，多个线程同时访问同一个 bank 会导致冲突；多个线程同时访问不同的 bank 能取得更高的性能。

除了开普勒架构（8 字节）外，每 4 个字节被划分到一个 bank。例如 0000 属于 bank0，0004 属于 bank1 …… 0128 属于 bank0。可以看出，每个 bank 的相邻层的地址相差 128。

一个线程束试图同时访问同一个 bank 中的 n 层数据将导致 n 次 内存事务(memory transaction)。亦称为 n 路 bank 冲突。n 很大的 bank 冲突是要尽量避免的。

在矩阵转置中，写入共享内存时易产生 bank 冲突，可以通过修改数组大小解决。

记得双精度浮点数占 8 字节。

原子函数

例如在数组归约（加法）中，要令GPU完成求和，需要使用原子操作。

特性

原子函数在全局内存和共享内存上体现原子性。

原子函数不依赖内存栅栏，不需要线程同步或顺序约束（参考官方文档）。

原子函数只能用于设备函数。

原子性粒度

系统原子性：在任意 CPU 和 GPU 的任意线程上保持原子性。函数名有后缀_system，如atomicAdd_system。

设备原子性：在当前 GPU 的任意线程上保持原子性。函数名无附加后缀。

线程块原子性：在当前 GPU 的同一线程块中的任意线程上保持原子性。函数名有后缀_block，如atomicAdd_block。

原子函数速查

注意：所有的函数都返回旧值

支持的类型可参考书第 92 页。

1
2
3
4
5
6
7
8
9
10
11
12
13
// T in {int, unsigned int, unsigned long long int}
using UI = unsigned int;
T atomicAdd(T *address, T val); // 加法，支持浮点
T atomicSub(T *address, T val); // 减法，支持浮点
T atomicExch(T *address, T val); // new = val，支持float
T atomicMin(T *address, T val); // 最小值
T atomicMax(T *address, T val); // 最大值
UI atomicInc(UI *address, UI val); // new = (old >= val) ? 0 : (old+1)
UI atomicDec(UI *address, UI val); // new = ((old == 0) || (old > val)) ? val : (old-1)
T atomicAnd(T *address, T val); // 按位与
T atomicOr(T *address, T val); // 按位或
T atomicXor(T *address, T val); // 按位异或
T atomicCAS(T *address, T compare, T val); // new = old == compare ? val : old。支持unsigned short int

线程束与协作组

SIMT

单指令-多线程（single instruction multiple thread, SIMT）：不同线程共享同一个 PC。SIMT 的通病是分支发散（branch divergence），所有的分支会产生串行的时间开销。

在伏特架构之前，一个线程束共享一个 PC。不同线程束之间没有分支发散问题。

从伏特架构开始，引入了独立线程调度（independent thread scheduling）机制。每个线程有自己的 PC，使得需要用户自己控制线程束内的同步。另一个代价是每个线程需要用两个寄存器来做 PC。如果旧代码出现线程束内不安全的问题，可以指定虚拟架构为低于伏特架构的计算能力。

线程束内的线程同步

__syncwarp()比线程块同步函数__syncthreads()更加廉价。

1
2
// 函数原型
void __syncwarp(unsigned mask = 0xffffffff);

其中掩码表示参与同步的线程，默认 32 个线程全部参与。

在分治时，当问题规模缩小到一个线程束内，可以不使用__syncthreads而用__syncwarp。

更多线程束内的基本函数

注意：若当前线程不参与，则函数的返回值是无定义的。

线程束表决函数（warp vote functions）

线程束匹配函数（warp match functions）待续

线程束洗牌函数（warp shuffle functions）

线程束矩阵函数（warp matrix functions）待续

线程束表决

1
2
3
4
5
6
7
8
unsigned __ballot_sync(unsigned mask, int predicate);
// 若当前线程参与，则同步表决下一次参与意愿（当 predicate 非零，则令返回值对应位置 1， 否则置 0。）。相当于从旧掩码表决出一个新掩码。

int __all_sync(unsigned mask, int predicate);
// 若当前线程参与，则同步投票，一票否决，返回投票是否通过。所有参与线程都同意才通过。

int __any_sync(unsigned mask, int predicate);
// 若当前线程参与，则同步投票，一票通过，返回投票是否通过。

线程束洗牌

定义“束内指标” int lane_id = threadIdx.x % w。

w 是逻辑线程束大小，只能取 2、4、8、16、32 中的一个。

1
2
3
4
5
6
7
8
9
10
11
T __shfl_sync(unsigned mask, T v, int srcLane, int w=warpSize);
// 参与线程返回标号为 srcLane 的变量 v 的值。这是一种广播式数据交换。 

T __shfl_up_sync(unsigned mask, T v, unsigned d, int w=warpSize);
// 标号为 t 的参与线程返回标号为 t - d 的线程中的 v 的值；若 t - d < 0 则返回原来的 v。相当于数据向上平移。

T __shfl_down_sync(unsigned mask, T v, unsigned d, int w=warpSize);
// 标号为 t 的参与线程返回标号为 t + d 的线程中的 v 的值；若 t + d >= w 则返回原来的 v。相当于数据向下平移。

T __shfl_xor_sync(unsigned mask, T v, int laneMask, int w=warpSize);
// 标号为 t 的参与线程返回标号为 t ^ laneMask 的线程中的 v 的值；相当于对应的两个线程交换数据。

使用线程束洗牌优化数组归约：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
constexpr unsigned int FULL_MASK = 0xffffffff;

void __global__ reduce_shfl(const real *d_x, real *d_y, const int N) {
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int n = bid * blockDim.x + tid;
    extern __shared__ real s_y[];
    s_y[tid] = (n < N ? d_x[n] : 0.0);
    __syncthreads();

    for (int offset = blockDim.x >> 1; offset >= 32; offset >>= 1) {
        if (tid < offset) {
            s_y[tid] += s_y[tid + offset];
        }
        __syncthreads();
    }

    real y = s_y[tid]; // 寄存器优化

    for (int offset = 16; offset > 0; offset >>= 1) {
        y +=  __shfl_down_sync(FULL_MASK, y, offset);
    }

    if (tid == 0) {
        atomicAdd(d_y, y);
    }
}

协作组

协作组（cooperative groups）可以看作线程块和线程束同步机制的推广。范围涵盖线程块内部，线程块之间（网格级）及设备之间的同步与协作。

引入头文件和命名空间

1
2
#include <cooperative_groups.h>
namespace cg = cooperative_groups; // 仅供参考

线程块级别的协作组

基本类型thread_group。有如下成员：

  1. void sync() 同步组内所有线程
  2. unsigned size() 返回组内总的线程数目
  3. unsigned thread_rank() 返回当前线程的组内标号
  4. bool is_valid() 检查定义的组是否违反 CUDA 的任何限制

导出类型thread_block，额外函数：

  1. dim3 group_index() 返回当前线程的线程块指标，相当于 blockIdx
  2. dim3 thread_index() 相当于 threadIdx

于是可以抽象出当前线程块：

1
2
thread_block g = this_thread_block();
g.sync(); // 等价于 __syncthreads()

可以将 thread_block 进行多次分割（但一组线程的数量只能是 2 的幂次）：

1
2
3
4
thread_group g32 = tiled_partition(this_thread_block(), 32); // 相当于线程束
thread_group g4 = tiled_partition(g32, 4); // 再次分割
// ---------- 模板化版本 ----------
thread_block_tile<32> g32 = tiled_partition<32>(this_thread_block());

这样的“线程块片”可以模仿线程束的一些行为，比如

1
2
3
4
5
6
7
8
9
unsigned __ballot_sync(int predicate);
int __all_sync(int predicate);
int __any_sync(int predicate);
T __shfl_sync(T v, int srcLane);
T __shfl_up_sync(T v, unsigned d);
T __shfl_down_sync(T v, unsigned d);
T __shfl_xor_sync(T v, int laneMask);

// 相对线程束的不同点：1. 不允许掩码，所有线程必须参与 2. 洗牌函数不再需要宽度参数，宽度由线程块片大小确定

再次优化数组归约

这次我们采用的技巧有：

  1. 使用静态全局内存
  2. 调用两次核函数，舍弃原子加法函数（这会提高精确度）
  3. 提高线程利用率（在规约之前进行小部分求和）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#include <real.hpp>
#include <iostream>
#include <algorithm>
#include <widgets.hpp>

constexpr unsigned int FULL_MASK = 0xffffffff;

void __global__ reduce_cp(const real *d_x, real *d_y, const int N) {
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    extern __shared__ real s_y[];

    real y = 0.0;
    const int stride = blockDim.x * gridDim.x; // 以网格大小为跨度
    for (int n = bid * blockDim.x + tid; n < N; n += stride) {
        y += d_x[n]; // 确保一个网格能覆盖所有数据
    }
    s_y[tid] = y;
    __syncthreads();

    // 线程块内，跨线程束折半归约
    for (int offset = blockDim.x >> 1; offset >= 32; offset >>= 1) {
        if (tid < offset) {
            s_y[tid] += s_y[tid + offset];
        }
        __syncthreads();
    }

    y = s_y[tid];

    for (int offset = 16; offset > 0; offset >>= 1) {
        y += __shfl_down_sync(FULL_MASK, y, offset);
    }

    if (tid == 0) {
        d_y[bid] = y; // 返回线程块结果
    }
}

constexpr int N = 1e8;
constexpr int GRID_SIZE = 10240;
constexpr int BLOCK_SIZE = 128;
__device__ real d_input[N];
__device__ real d_output[GRID_SIZE];

real reduce(const real *d_x) {
    real *d_y;
    CHECK(cudaGetSymbolAddress((void**)&d_y, d_output));
    constexpr int shared_size = sizeof(real) * BLOCK_SIZE;
    reduce_cp<<<GRID_SIZE, BLOCK_SIZE, shared_size>>>(d_x, d_y, N);
    reduce_cp<<<1, 1024, sizeof(real) * 1024>>>(d_y, d_y, GRID_SIZE);

    real h_y[1] = {0};
    CHECK(cudaMemcpy(h_y, d_y, sizeof(real), cudaMemcpyDeviceToHost));
    // CHECK(cudaMemcpyFromSymbol(h_y, d_output, sizeof(real)));

    return h_y[0];
}

int main() {
    using namespace std;
    static real input[N];
    std::fill_n(input, N, 1.23f);
    CHECK(cudaMemcpyToSymbol(d_input, input, N * sizeof(real)));
    real *d_x;
    CHECK(cudaGetSymbolAddress((void**)&d_x, d_input));
    printf("%f\n", reduce(d_x));
    return 0;
}

// 运行结果：
// 123000064.000000
// 对比 CPU 裸归约方法的加速比：67

CUDA 流

核函数外部的并行

主要有以下情形：

  1. 核函数计算与数据传输之间的并行
  2. 主机计算与数据传输之间的并行
  3. 不同的数据传输之间的并行
  4. 核函数计算与主机计算之间的并行
  5. 不同核函数之间的并行

若两个任务的运行时间相近，则尽量令他们并行。反之，则并行和串行的性能差距不大。

使用 CUDA 流的主要目的是尽量取得核函数外部的并行能力。

使用 CUDA 流

CUDA 流的使用相对简单，只要记住一些规则就好。

  1. CUDA 程序有一个默认流。若 API 不需要指定流句柄，可推测它用于默认流。

  2. 创建、销毁 CUDA 流

    1
    2
    3
    4
    5
    6
    7
    8
    cudaStream_t stream[5];
    for (auto &s : stream) {
        cudaStreamCreate(&s); // 创建流，获得句柄
    }
    // ...
    for (auto &s : stream) {
        cudaStreamDestroy(s); // 销毁流，句柄不抹除
    }

  3. 检测 CUDA 流是否空闲 

    1
    2
    cudaError_t cudaStreamSynchronize(cudaStream_t stream); // 阻塞主机，直到 CUDA 流完成所有操作
    cudaError_t cudaStreamQuery(cudaStream_t stream); // 返回 cudaSuccess（已完成）/ cudaErrorNotReady（未完成）

  4. cudaMemcpy 会阻塞主机。但是核函数的调用总是异步的。应该尽量早地调用核函数，随后再进行主机的计算任务。

  5. 在核函数配置中指派 CUDA 流：

    1
    my_kernel<<<N_grid, N_block, N_shared, stream_id>>>(...); // 四个模板参数必须齐全

  6. 异步的数据传输需要：

    1. 使用 cudaMemcpyAsync 函数。
    2. 涉及的主机内存必须是不可分页的，否则 API 将会退化为同步版本。

    异步传输的过程将由 GPU 的 DMA（direct memory access）接管。

  7. 可以用 CUDA 运行时 API 来申请不可分页的内存（non-pageable memory），又称固定内存（pinned memory）。不可分页意味着操作系统无权修改虚拟地址所对应的物理地址。

    1
    2
    3
    cudaError_t cudaMallocHost(void **ptr, size_t size);
    cudaError_t cudaHostAlloc(void **ptr, size_t size, size_t flags);
    cudaError_t cudaFreeHost(void *ptr);

  8. 良好的 CUDA 流指派会从流水线重叠中获得并行加速。但是 CUDA 流是有启动开销，且硬件环境有限，过多的 CUDA 流会拉低性能。

指令速查

编译器选项

1
2
--ptxas-options=-v # 报道每个核函数的寄存器使用数量
--maxrregcount= # 限制所有核函数的寄存器使用量

修饰符

1
2
3
4
5
6
7
8
__global__ # 核函数修饰符
__device__ # 设备函数、变量修饰符
__host__ # 主机函数修饰符，一般只与__device__同时出现
__noinline__ # 建议函数非内联
__forceinline__ # 建议函数内联
__launch_bounds__() # 修饰核函数，限制寄存器使用量
__shared__ # 修饰全局变量，使它成为共享内存
__managed__ # 修饰全局变量，使它成为统一内存，必须与__device__同时出现

常用代码

CMake模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# 2021/08/02
cmake_minimum_required(VERSION 3.20)

project(hello LANGUAGES CUDA VERSION 0.1)

option(USE_DOUBLE "Use real as double, otherwise float" OFF)

find_package(CUDAToolkit REQUIRED)

# set(CMAKE_CUDA_COMPILER_ID NVIDIA)

# set(ENV{PATH} "C:/Program\ Files\ (x86)/Microsoft\ Visual\ Studio/2019/Community/VC/Tools/MSVC/14.29.30037/bin/Hostx64/x64:$ENV{PATH}")

# report registers used in each core function
add_compile_options(--ptxas-options=-v)

include_directories(./include)

add_executable(hello main.cu)

# compute capability, see https://developer.nvidia.com/zh-cn/cuda-gpus#compute
set_property(TARGET hello PROPERTY CUDA_ARCHITECTURES 86)
# set_property(TARGET hello PROPERTY CUDA_ARCHITECTURES 75)

# install(TARGETS hello DESTINATION .)

报错、计时工具

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
// widgets.hpp
#pragma once
#include <cstdio>
#include <chrono>

#define CHECK(callee)                                            \
    do {                                                         \
        const cudaError_t err = callee;                          \
        if (err == cudaSuccess) break;                           \
        printf("CUDA error at %s(%d)\n", __FILE__, __LINE__);    \
        printf("    Function:   %s\n", __FUNCTION__);            \
        printf("    Error code: %d\n", err);                     \
        printf("    Error hint: %s\n", cudaGetErrorString(err)); \
        exit(1);                                                 \
    } while (0)

template<typename F>
void cudaTiming(const F &func) {
    cudaEvent_t start, stop;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop));
    CHECK(cudaEventRecord(start));
    cudaEventQuery(start);

    func();

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    float elapsed_time;
    CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));
    printf("CUDA Time = %g ms.\n", elapsed_time);
    CHECK(cudaEventDestroy(start));
    CHECK(cudaEventDestroy(stop));
}

template<typename F>
void hostTiming(const F &func) {
    auto start = std::chrono::steady_clock::now();

    func();

    auto end = std::chrono::steady_clock::now();
    std::chrono::duration<double, std::milli> duration = end - start;
    printf("Host Time = %g ms.\n", duration.count());
}

精度控制模板

1
2
3
4
5
6
7
8
9
10
// real.hpp
#pragma once

#ifdef USE_DOUBLE
using real = double;
constexpr real EPS = 1e-15;
#else
using real = float;
constexpr real EPS = 1e-6f;
#endif