CUDA入门

原创已于 2025-05-13 10:31:09 修改 · 754 阅读

10 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

于 2025-05-12 05:12:35 首次发布

CUDA入门

第一讲——CUDA介绍

一、CPU vs GPU

在这里插入图片描述

属性	CPU	GPU
核心数量	少（高性能）	多（低功耗）
延迟	低	高
吞吐量	低（适合顺序执行）	高（适合并行执行）
控制能力	强（分支预测、乱序执行）	弱
擅长任务类型	串行计算、控制密集型	并行计算、数据密集型

二、GPU加速方法与并行化工具

函数库（Libraries）

库名	主要功能
FFT	快速傅里叶变换（Fast Fourier Transform）
BLAS	基本线性代数子程序（Basic Linear Algebra Subprograms），如矩阵乘法
Thrust	类似于 C++ STL 的并行算法库，支持如 `sort`、`reduce`、`scan` 等并行操作
RAND	随机数生成库，适合蒙特卡洛模拟等场景
Statistics	提供统计运算功能，如均值、方差、直方图等

编译器指令

指令集	描述
OpenACC	一组高层指令（pragma）可用于在 C/C++/Fortran 中标注哪些代码应并行执行，并自动映射到 GPU

#pragma acc parallel loop
for (int i = 0; i < N; ++i) {
    a[i] = b[i] + c[i];
}

GPU编程语言

语言	描述
CUDA C/C++	NVIDIA 提供的编程模型，最常用的 GPU 编程语言
CUDA Fortran	Fortran 的 CUDA 扩展
HIP	AMD 推出的与 CUDA 类似的编程接口
OpenCL	一种跨平台的并行计算框架，支持多种硬件设备

三、CUDA编程模型基础

编程模型是底层计算机系统的抽象，用于表达算法和数据；程序语言和 API 用于上述抽象的实现。

CUDA 执行模型——异构主机+设备模型

角色	描述
Host	一般是 CPU，负责控制、分配任务
Device	指 GPU，负责大量数据并行计算

串行部分在主机的 C 程序中执行
并行部分由 GPU 上的 kernel（SIMD：Single Program Multiple Data）线程阵列执行

线程层次化结构

Grid(网格)
└── Block(线程块)  		  ← blockIdx
	└── Warp(线程束)	      ← 每 32 个线程组成一个 warp（隐式）
    	└── Thread(线程) 	   ← threadIdx

名称	描述	数量限制 / 特点
线程	最小执行单元，执行 kernel 中的代码	每个 Block最多1024 个线程
线程束	最小调度单元，包含最多32 个线程	同一warp内线程执行同一指令
线程块	一组线程的集合	block可以是一维、二维或三维，最多 $2^{10}$ 个线程
线程网格	包含所有Block，为每次kernel调用创建	Grid可以是一维、二维或三维，最多 $2^{16}$ 个线程块

内存层次化结构

内存类型	生命周期	访问范围	特点
全局内存	GPU 全局内存	所有线程可访问	容量大、速度慢
共享内存	每个 block 独享	同一block的线程共享	容量小、速度快（类似缓存）
寄存器	每线程独有	仅该线程访问	速度最快，但数量有限

第二讲——CUDA并行模型

一、CUDA执行模型

异构主机+设备模型

角色	描述
Host	一般是 CPU，负责控制、分配任务
Device	指 GPU，负责大量数据并行计算

串行部分在主机的 C 程序中执行
并行部分由 GPU 上的 kernel（SIMD）线程阵列执行

二、内存分配及数据传输API

在这里插入图片描述

内存分配

cudaMalloc和cudaFree

cudaError_t cudaMalloc(
    void **devPtr, 		// 二级指针（void**），用于接收分配的设备内存地址
    size_t size			// 需分配的字节数（size_t 类型）	
);
// cudaError_t表示操作状态（如cudaSuccess或错误码）
cudaError_t cudaFree(
    void* devPtr		// devPtr：指向待释放设备内存的指针
);

cudaMalloc 在 GPU 的全局内存中分配连续的线性空间，返回的指针 *devPtr 指向起始地址，分配的内存内容未初始化，需手动填充数据。主机（CPU）代码不能直接解引用 cudaMalloc 返回的指针进行读写，必须通过 cudaMemcpy 传输数据，必须调用 cudaFree 释放分配的内存，避免泄漏。

C/C++ 按值传递指针时，函数内部无法修改外部指针的值。通过传递 void**，cudaMalloc 可修改调用者的指针变量，使其指向设备内存地址。

cudaFree释放 GPU 全局内存中由 devPtr 指向的内存块，避免内存泄漏，调用后 devPtr 不会自动置空，需手动设为 nullptr 以防悬空指针。

cudaMallocHost和cudaFreeHost

cudaError_t cudaMallocHost(void **ptr, size_t size);
cudaError_t cudaFreeHost(void *ptr);

cudaMallocHost 在主机内存中分配连续的线性空间，必须调用 cudaFree 释放分配的内存，避免泄漏。

数据传输

cudaMemcpy

cudaError_t cudaMemcpy(
    void* dst,          	// 目标地址（主机或设备指针）
    const void* src,    	// 源地址（主机或设备指针）
    size_t count,       	// 复制的字节数
    cudaMemcpyKind kind 	// 传输方向标识
);

kind 枚举值	说明
`cudaMemcpyHostToHost`	主机内存 → 主机内存（极少使用）
`cudaMemcpyHostToDevice`	主机内存 → 设备显存（常见）
`cudaMemcpyDeviceToHost`	设备显存 → 主机内存（常见）
`cudaMemcpyDeviceToDevice`	设备显存 → 设备显存
`cudaMemcpyDefault`	自动判断方向（需统一内存支持）

cudaMemcpy 是 CUDA 编程中用于在主机（CPU）与设备（GPU）之间或设备内存内部复制数据的核心函数。默认同步执行，会阻塞 CPU 线程直到传输完成，异步版本需使用 cudaMemcpyAsync。

cudaMemcpyAsync

cudaError_t cudaMemcpyAsync(
    void* dst,                 // 目标地址（主机或设备指针）
    const void* src,           // 源地址（主机或设备指针）
    size_t count,              // 复制的字节数
    cudaMemcpyKind kind,       // 传输方向标识
    cudaStream_t stream = 0    // 关联的 CUDA 流（默认为0，即默认流）
);

cudaMemcpyAsync 是 CUDA 中用于异步内存复制的核心函数，允许在主机（CPU）与设备（GPU）之间或设备内存内部并行执行数据传输，同时不会阻塞主机线程。主机内存必须为页锁定内存（通过 cudaMallocHost 分配），否则性能会下降或无法异步执行。将 cudaMemcpyAsync与核函数执行分配到不同流中，实现并行，如计算与数据传输重叠。显式同步通过 cudaStreamSynchronize(stream) 等待流中所有操作完成。

向量加法主机代码

int main(){
    float *a, *b, *out;
    float *d_a, *d_b, *d_out; 
    a   = (float*)malloc(sizeof(float) * N);
    b   = (float*)malloc(sizeof(float) * N);
    out = (float*)malloc(sizeof(float) * N);

    for(int i = 0; i < N; i++){
        a[i] = 1.0f;
        b[i] = 2.0f;
    }

    cudaMalloc((void**)&d_a, sizeof(float) * N);
    cudaMalloc((void**)&d_b, sizeof(float) * N);
    cudaMalloc((void**)&d_out, sizeof(float) * N);

    cudaMemcpy(d_a, a, sizeof(float) * N, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, b, sizeof(float) * N, cudaMemcpyHostToDevice);

    int block_size = 256;
    int grid_size = (N + block_size - 1) / block_size;
    vector_add <<< grid_size, block_size >>> (d_out, d_a, d_b, N);

    cudaMemcpy(out, d_out, sizeof(float) * N, cudaMemcpyDeviceToHost);

    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_out);
    free(a); 
    free(b); 
    free(out);
}

三、数据并行化与线程

CUDA 线程的本质

每个 CUDA 线程可看作一个“虚拟化”的 Von-Neumann 处理器
每个线程独立运行同一份内核程序（SIMD），但使用不同的数据索引

并行线程阵列（线程网格）

CUDA kernel 的执行单位是 线程网格（Grid）
每个线程网格由多个 线程块（Block） 组成
所有线程运行相同代码，但通过索引处理不同数据

名称	CUDA 变量
线程	`threadIdx.x/y/z`
线程束	（隐式）
线程块	`blockIdx.x/y/z`、`blockDim.x/y/z`
网格	`gridDim.x/y/z`

__global__ void vector_add(float *out, float *a, float *b, int n) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if(tid < n) {
        out[tid] = a[tid] + b[tid];
    }
}

在这里插入图片描述

如向量加法代码和图解所示，当主机执行代码vector_add <<< grid_size, block_size >>> (d_out, d_a, d_b, N); 线程网格一共有grid_size个block，每个block有block_size个线程，那么共有grid_size * block_size个线程执行vector_add函数，对于每个线程可以通过tid = blockIdx.x * blockDim.x + threadIdx.x;得到其在整个网络的编号，当然我们只需要前n个线程就够用了，每个线程处理一个位置上的加法即可，这就是数据并行化与线程。

四、多维内核配置

前两节向量加法代码中，我们分配数据使用一维数组，获取线程tid是线性计算，第一讲中提到网格可以是一维、二维、三维，线程块也可以是一维、二维、三维，这就为处理多维数据提供了便利，比如可以让线程块为 $2^{3} * 2^{3} * 2^{3}$ 的三维线程块，同时网格也为 $2^{3} * 2^{3} * 2^{3}$ 的三维网格，那么一共有 $2^{6} * 2^{6} * 2^{6}$ 个线程，恰好就可以一个线程操作 $2^{6} * 2^{6} * 2^{6}$ 三维数组的一个值了。本节旨在通过二维图片处理展示多维线程网格和线程块的使用。

在这里插入图片描述

如图所示，我们需要配置大小为 $16 * 16$ 的线程块，并配置恰当的线程网格以覆盖 $62 * 76$ 的图片

// 图片大小n x m，宽为m，高为n
dim3 DimGrid((n - 1)/16 + 1, (m - 1)/ 16 + 1, 1);
dim3 DimBlock(16, 16, 1);
kernelFunction <<< DimGrid, DimBlock >>> (d_Pin, d_Pout, m, n);

灰度化处理核函数示例

在这里插入图片描述

#define CHANNELS 3

__global__ void colorConvert(unsigned char *grayImage, 
                             unsigned char *rgbImage, 
                             int width, int height) {
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    
    if (x < width && y < height) {
        int grayOffset = y * width + x;
        int rgbOffset = grayOffet * CHANNELS;
        
        unsigned char r = rgbImage[rgbOffset];
        unsigned char g = rgbImage[rgbOffset + 1];
        unsigned char b = rgbImage[rgbOffset + 2];
        
        grayImage[grayOffset] = 0.21f * r + 0.71f * g + 0.07f *  b;
    }
}

模糊图像处理核函数示例

在这里插入图片描述

__global__ void blurKernel(unsigned char *in, 
                           unsigned char *out, 
                           int width, int height) {
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    
    if (x < width && y < height) {
        int pixVal = 0;
        int pixels = 0;
        
        for (int blurX = -BLUR_SIZE; blurX < BLUR_SIZE + 1; ++blurX) {
            for (int blurY = -BLUR_SIZE; blurY < BLUR_SIZE + 1; ++blurY) {
                int curX = x + blurX;
                int curY = y + blurY;
                
                if (curX > -1 && curX < weight && curY > -1 && curY < height) {
                    pixVal += in[curY * width + curX];
                    pixels++;
                }
            }
        }
        
        out[y * width + x] = (unsigned char)(pixVal / pixels);
    }
}

五、线程调度

线程调度基础

CUDA 使用硬件执行资源流式处理器SM来并行调度线程。每个线程块（Block）作为一个调度单位分配给 SM，CUDA 调度器以零开销方式调度线程束（warp），不需要显式上下文切换。线程块之间无依赖，可以以任意顺序被调度执行。CUDA 程序不依赖硬件中有多少 SM，因此具有很好的可扩展性。

SM资源限制

资源类型	例：Fermi上限
每个 SM 可分配线程块数	最多 8 个
每个 SM 可容纳线程数	最多 1536 线程

如果线程块有 256 个线程，则最多 1536 / 256 = 6 个线程块可并行执行
如果线程块有 512 个线程，则最多 1536 / 512 = 3 个线程块可并行执行

这里忽略了线程块所需共享内存、所需寄存器数量对容纳线程块的影响，实际需要考虑这些资源限制。

Warp调度

每个线程块被划分为多个线程束（Warp），每个Warp有32个线程。Warp是 SM 的基本调度单位。

示例：

一个线程块含256个线程，则包含256 / 32 = 8个 warp
若一个 SM 分配3个此线程块，则总共有8 × 3 = 24个 warp 被调度执行

上面说CUDA调度器以零开销方式调度线程束，不需要显式上下文切换，是因为若一个 SM 分配3个此线程块，SM同时驻留这3个线程块，每个线程块的上下文信息一直保留在硬件中，这些线程块被称为是 “驻留的”，Warp 则是 “活跃的”，切换Warp执行，不需要加载/保存寄存器或状态，只需改变指令发射目标，这个成本非常低，称为“零开销线程调度”。

例题

问题： 在执行矩阵乘法时，哪种线程块大小配置更好？
A. 8×8 B. 16×16 C. 32×32

每个线程块线程数：A: 64 B: 256 C: 1024
Fermi 每个 SM 最多 1536 个线程，最多 8 个线程块

配置	可同时分配的线程块数	SM 利用率（线程数）
A 8×8	最多 8 个（64×8）	512
B 16×16	最多 6 个（256×6）	1536（满载）
C 32×32	最多 1 个（1024×1）	1024

第三讲——内存与数据划分

一、内存访问效率

简单矩阵乘法示例

__global__ void MatrixMulKernel(float *M, float *N, float *P, int Width) {
    int Row = blockIdx.y * blockDim.y + threadIdx.y;
    int Col = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (Col < Width && Row < Width) {
        float Pvalue = 0.0;
        
        for (int k = 0; k < Width; ++k) {
            Pvalue += M[Row * Width + k] * N[k * Width + Col];
        }
        
        P[Row * Width + Col] = Pvalue;
    }
}

性能评估指标——CGMA

$\frac{FLOPs}{Global Memory Accesses}$ ，它衡量的是每一次全局内存访问带来多少次浮点运算。

高CGMA表示你“用一次内存访问换来了更多计算”，提高了资源利用率；内存带宽是GPU上的重要瓶颈，CGMA越大，越不容易成为瓶颈；高CGMA程序更有可能接近GPU的FLOPs峰值。因此CGMA越大越好。

简单矩阵乘法性能评估

for (int k = 0; k < Width; ++k) {
    Pvalue += M[Row * Width + k] * N[k * Width + Col];
}
P[Row * Width + Col] = Pvalue;

每次循环进行1次乘法和1次加法，一共循环 Width 次，因此每个线程执行 2 * Width 次浮点运算，每个线程读 M[Row * Width + k] 访问 Width 次，读 N[k * Width + Col] 访问 Width 次，写 P[Row * Width + Col] 写 1 次，合计 2 * Width + 1 次全局内存访问。
$\frac{FLOPs}{Global Memory Accesses} = \frac{2 * Width}{2 * Width + 1} \approx 1 FLOPs$

假设硬件参数	数值
浮点峰值性能	1500 GFLOPS
DRAM内存带宽	200 GB/s

矩阵乘法每次全局内存访问传输4 字节（float），那么每秒最多能够访问内存 $200GB/s4B=50G/s\frac{200GB/s}{4B}=50G/s$ 次，由于CGMA=1，浮点运算速度为 $50 G / s * C G M A = 50 G F L O P s$ ，理论上有1500 GFLOPs，实际只有50 GFLOPs，内存是瓶颈，执行速度限制为设备峰值浮点执行率的 3.3%。

CUDA 编程中，内存访问效率对程序性能的影响是决定性的，因为 GPU 擅长计算，但常常被数据访问延迟所“拖慢”。如果不优化内存访问，线程大部分时间都在“等数据”，计算资源就被浪费了，并且高并发线程下，若每个线程访问内存方式不合理，会造成“内存带宽争用”和“访存不合并”，大幅降低吞吐量。

内存访问效率

在这里插入图片描述

变量声明/类型	内存类型	作用域	生命周期	相对访问速度	说明
`int a;`	Register	Thread	Thread	最快	每个线程私有，位于 SM 的寄存器中，访问延迟最低
`__shared__ int a;`	Shared	Block	Block	很快	SM 内的片上存储，线程块内共享，用于数据复用和通信
`__device__ int a;`	Global	Grid	Grid	较慢	所有线程可读写，位于 DRAM，访问代价高
`__device__ __constant__ int a;`	Constant	Grid	Grid	快	所有线程共享，只读常量，小而快，适合广播常量

共享内存特点

属性	描述
本质	每个 SM 上的片上内存，也称为 Scratchpad Memory
速度	比全局内存快得多，延迟低（接近寄存器），吞吐高
作用域	限于线程块内（block-level）
生命周期	随线程块生命周期结束而消失
访问方式	使用普通指针变量（数组或变量名），通过 `__shared__` 声明
分配方式	静态（编译时大小固定）或动态（内核调用时指定大小）

“Scratchpad” 是计算机体系结构中的术语，指临时高速缓存区域，由程序员显式管理，与自动缓存（如 L1/L2）不同，CUDA 的共享内存类似一个 用户可控的高速缓存，可以明确写代码指示何时读/写，提高数据复用性，减少全局内存访问。

提高内存访问效率的方法总结

方法	描述
提高计算密度（CGMA）	提高每次内存访问所执行的计算量，减少访问频率，提高吞吐率
使用共享内存（shared）	减少全局内存访问，局部数据复用（如 block tile 矩阵乘法）
内存访问对齐（coalescing）	使同一 warp 内线程访问连续内存，减少访问次数
使用常量内存	小量只读共享数据，用常量内存加速广播

二、分块并行算法

分块并行概念

在这里插入图片描述

“分块并行算法”是并行编程中一种常用的优化策略，特别适用于 CUDA 等 GPU 编程，通过将数据划分为多个小块（Tile/Block），以更好地利用 共享内存 和 线程块结构，从而提升性能。

分块并行同步

在 CUDA 分块并行算法中，同步操作主要是为了确保所有线程在共享内存读写时保持一致性与数据安全。

操作场景	是否需要同步？	原因
所有线程向共享内存写入后	是	防止读取未完成的数据
所有线程从共享内存读取后	是	防止其他线程提前覆盖数据
每个线程只读写自己的共享内存区域	否	无数据竞争，不需要同步

分块并行算法中的同步策略

分块识别：根据线程 ID 计算当前分块索引
加载共享内存：从全局内存读取 tile → 写入共享内存
同步：__syncthreads() 确保 tile 完全加载完毕
计算：所有线程并行计算当前 tile 中的乘积累加
同步：再次 __syncthreads() 准备进入下一轮 tile 加载
重复：处理所有分块直到计算完成

分块矩阵乘法示例

在这里插入图片描述

#define TILE_WIDTH 16  // 每个线程块处理 16x16 子矩阵

__global__ void MatrixMulKernel(float *M, float *N, float *P, int Width) {
	// 假设Width % TILE_WIDTH == 0
    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    int bx = blockIdx.x;  int by = blockIdx.y;
    int tx = threadIdx.x; int ty = threadIdx.y;

    int Row = by * TILE_WIDTH + ty;
    int Col = bx * TILE_WIDTH + tx;

    float Pvalue = 0.0;

    // 每轮迭代处理一个 tile（子矩阵）
    for (int k = 0; k < Width / TILE_WIDTH; ++k) {
        // 加载 M、N 的 tile 到共享内存
        Mds[ty][tx] = M[Row * Width + (k * TILE_WIDTH + tx)];
        Nds[ty][tx] = N[(k * TILE_WIDTH + ty) * Width + Col];

        __syncthreads();  // 等待所有线程加载完毕

        // 对当前 tile 执行部分乘法累加
        for (int l = 0; l < TILE_WIDTH; ++l)
            Pvalue += Mds[ty][l] * Nds[l][tx];

        __syncthreads();  // 确保计算完再加载下一个 tile
    }

    // 将结果写回全局内存
    P[Row * Width + Col] = Pvalue;
}

在这里插入图片描述

并行本质

上面的代码本质就是上图的公式，这里可以看到 $Cij=∑k=1tAikBkjC_{ij} = \sum^{t}_{k=1}A_{ik}B_{kj}$ ，其中 $C_{ij}$ ，体现在当前是 $(i, j)$ 线程块，就是代码中的(by, bx)，这隐含在当前线程的背景中，接下来的关于k的for循环对应公式中的累加部分，代码中的Mds Nds对应了 $A_{ik}$ $ B_{kj} $，接下来的循环就是对$ A_{ik}$ $ B_{kj} $进行矩阵乘法，最后$ C_{ij}$的(ty,tx)的值Pvalue进行k次点乘，且累加，实现了矩阵块的计算。

同步问题

关于分块矩阵乘法核函数的同步问题，可以看到代码中使用了两个__syncthreads()进行同步，第一个同步函数是为等待此线程所在线程块的所有线程将它们对应在矩阵中分块的数据写入共享内存，因为下面进行值计算时，一个线程不仅要用到自己加载的数据也要用到其它线程加载的数据，因此等待防止读取未完成的数据；第二个同步函数，确保计算完再加载下一个块，因此这是一个循环，到k+1又要加载一边数据，如果不同步，下一个加载进去了，这边还在计算，就会导致错误，因此等待防止其他线程提前覆盖数据。

分块(线程块)大小注意事项

每个线程块应该包含很多线程

TILE_WIDTH为16的线程块包含256个线程
TILE_WIDTH为32的线程块包含1024个线程

当线程块宽度为16时，在每个阶段，每个线程块从全局内存中进行 $2 * 256 = 512$ 次浮点载入，用于 $256 * (2 * 16) = 8, 192$ 次乘法/加法运算， $\frac{8192}{512}=16$ 。

当线程块宽度为32时，在每个阶段，每个线程块从全局内存中进行 $2 * 1024 = 2048$ 次浮点载入，用于 $1024 * (2 * 32) = 65, 536$ 次乘法/加法运算， $\frac{65536}{2048}=32$ 。

假设每个SM都有16KB的共享内存，最多驻留8个线程块，最多1536个线程：

当TILE_WIDTH = 16时，每个线程块使用的共享内存大小为 $2 * 256 * 4 B = 2 K B$ ，对于16KB的共享内存，可以满足最多8个线程块执行，对于最多1536个线程，最多满足6个线程块执行，综上可以驻留6个线程块。

当TILE_WIDTH = 32时，每个线程块使用的共享内存大小为 $2 * 1024 * 4 B = 8 K B$ ，对于16KB的共享内存，可以满足最多2个线程块执行，对于最多1536个线程，最多满足1个线程块执行，综上可以驻留1个线程块。

在实际应用中，不能仅看 CGMA，也必须结合 SM 的资源限制（共享内存、线程数等）分析。TILE_WIDTH = 32 拥有更高 CGMA、更少内存带宽压力，但并发性差；TILE_WIDTH = 16 CGMA 适中，但支持更多线程块并发，更可能接近峰值计算性能。综上较小且更多的线程块可能是有利的。

任意矩阵维度的分块矩阵乘法

到目前为止，提出的分块矩阵乘法内核只能处理尺寸是分块宽度倍数的方阵，然而，实际应用程序需要处理任意大小的矩阵。可以将行和列填充成分块大小的倍数，但会产生大量的空间和数据传输时间开销。我们将采取一种不一样的方法。

使用 numTiles = ceil(M_WIDTH / TILE_WIDTH) 控制每阶段分块，使用 边界检查（if (...)）避免访问越界内存，输出时再次检查 Row 和 Col 是否在有效矩阵内，防止写越界。

#define TILE_WIDTH 16

__global__ void MatrixMulKernel(float* M, float* N, float* P, 
                                int M_HEIGHT, int M_WIDTH, int N_WIDTH) {
    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    int Row = blockIdx.y * TILE_WIDTH + threadIdx.y;
    int Col = blockIdx.x * TILE_WIDTH + threadIdx.x;

    float Pvalue = 0.0f;

    int numTiles = (M_WIDTH + TILE_WIDTH - 1) / TILE_WIDTH;

    for (int t = 0; t < numTiles; ++t) {
        if (Row < M_HEIGHT && (t * TILE_WIDTH + threadIdx.x) < M_WIDTH)
            Mds[threadIdx.y][threadIdx.x] = 
            	M[Row * M_WIDTH + (t * TILE_WIDTH + threadIdx.x)];
        else
            Mds[threadIdx.y][threadIdx.x] = 0.0f;

        if ((t * TILE_WIDTH + threadIdx.y) < M_WIDTH && Col < N_WIDTH)
            Nds[threadIdx.y][threadIdx.x] = 
            	N[(t * TILE_WIDTH + threadIdx.y) * N_WIDTH + Col];
        else
            Nds[threadIdx.y][threadIdx.x] = 0.0f;

        __syncthreads();

        for (int k = 0; k < TILE_WIDTH; ++k) {
            Pvalue += Mds[threadIdx.y][k] * Nds[k][threadIdx.x];
        }

        __syncthreads();
    }

    if (Row < M_HEIGHT && Col < N_WIDTH) {
        P[Row * N_WIDTH + Col] = Pvalue;
    }
}

第四讲——性能问题

一、Wraps and SIMD

在这里插入图片描述

Warps

每个线程块拆分为包含32线程的线程束，具体实现不是 CUDA 编程模型的一部分，不归程序员管。线程束是 SM中的调度单元，线程束中的线程以单指令多数据 (SIMD) 方式执行。

在这里插入图片描述

如图所示是一个线程块，在多维线程块中，线程块按照x、y、z维的顺序线性化为一维，线程块在线性化后被划分线程束。线程束内或线程束之间的任何顺序都不具参考性，如果线程之间存在任何依赖关系，则必须执行 __syncthreads() 才能获得正确的结果。

SIMD

在这里插入图片描述

流式多处理器是SIMD处理器，用于指令获取、解码和控制的控制单元在多个处理单元之间共享，前面提到线程束是SM的基本调度单位，可以理解为SM内有1个控制单元和32个处理单元。当线程束中的线程通过判断分支执行不同的代码流路径时，就会发生控制分支现象，一些线程采用then路径，另一些采用else路径。采用不同路径的线程的执行在当前GPU中被串行化，当分支或循环条件是线程索引的函数时，可能会出现控制分支现象如下图所示，CUDA 使用分支掩码机制，吞吐量约为 1/N（N 为分支数量）。

在这里插入图片描述

示例

__global__ void vector_add(float *out, float *a, float *b, int n) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if(tid < n) {
        out[tid] = a[tid] + b[tid];
    }
}

假设n为1,000，每个线程块的大小是 256 个线程

那么共有4个线程块，不妨标记为0、1、2、3，每个线程块包含8个线程束

显然线程块0、1、2的所有线程均在n的范围内，那么它们的24个线程束不会发生控制分支现象，因为所有线程的条件一致

对于线程块3，线程束0-6均在有效范围内，不会发生控制分支现象，线程束7中线程992-999在有效范围内，线程1000-1023在有效范围外，会发生控制分支现象

那么最终32个线程束只有一个会发生控制分支现象，串行化（分支掩码顺序执行各个分支）对控制分支的影响很小，对性能影响较小

二、控制分支现象对性能影响

边界条件检查对于并行代码的完整功能和稳健性至关重要，平铺矩阵乘法内核有许多边界条件检查，然而，这些检查可能会导致性能显着下降，例如以下分块加载代码中：

if (Row < M_HEIGHT && (t * TILE_WIDTH + threadIdx.x) < M_WIDTH)
    Mds[threadIdx.y][threadIdx.x] = 
        M[Row * M_WIDTH + (t * TILE_WIDTH + threadIdx.x)];
else
    Mds[threadIdx.y][threadIdx.x] = 0.0f;

在这里插入图片描述

载入矩阵M分块时的两类线程块。Type 1直到最后一个阶段之前，其分块都在有效范围内的线程块；Type 2其分块部分都在有效范围之外的线程块。

假设有16x16的分块和线程块，每个线程块有8个线程束
假设有100x100的方阵

每个线程会历经7个阶段 $(100 / / 16 + 1)$ ，有 49 个线程块。

有42个Type 1的线程块，即前6行的线程块，一共包含336个线程束，共有7次加载M矩阵分块的过程，因此共有2,352个线程束-阶段，但是只有最后一次加载存在控制分支现象，因此只有226个线程束-阶段存在控制分支现象

有7个Type 2的线程块，即最后一行的线程块，一共包含56个线程束，共有7次加载M矩阵分块的过程，因此有392个线程束阶段，每个Type 2直到最后一个阶段之前，其线程块的前2个线程束处于有效范围内，后6个线程束处于有效范围外，且最后一个阶段，只有前两个线程束存在即有线程在有效范围内又有线程在有效范围外，因此，最后只有 $2 * 7 = 14$ 个线程束阶段存在控制分支现象

因此总共有2744个线程束-阶段，有350个线程束-阶段会发生控制分支现象

预计性能影响为$\frac{336+14}{2352+392} * 0.5=\frac{350}{2744} * 0.5 \approx 0.06375 = 6.375% $

总结

预计的性能影响取决于数据，对于更大的矩阵，影响会明显更小。一般来说，控制分支对大型输入数据集的边界条件检查的影响应该是不重要的，应该毫不犹豫地使用边界检查来确保完整的功能，内核充满控制流结构这一事实并不意味着控制分支会大量发生。稍后我们将介绍一些自然会导致控制分支的算法模式。

三、并行规约/归约

并行归约是一种常见的并行计算模式，广泛用于将一组输入数据归约为一个结果，例如求和、最大值、最小值、乘积等。它在图像处理、科学计算、深度学习等领域非常常见。

阶段1: [a0+a1, a2+a3, a4+a5, a6+a7]
阶段2: [(a0+a1)+(a2+a3), (a4+a5)+(a6+a7)]
阶段3: (((a0+a1)+(a2+a3)) + ((a4+a5)+(a6+a7))) → 1个值

并行归约常将数据集划分为更小的块，使用归约树将每个块的结果汇总为最终答案，时间复杂度为 $O (l o g N)$ 。

归约树每一层的操作数为：
$\frac{1}{2}N + \frac{1}{4}N + \frac{1}{8}N + \cdots + \frac{1}{N}N = (1 - \frac{1}{N})N = N - 1$
这和串行归约（N-1 次操作）一样，因此它是 工作高效（work-efficient） 的。
$\text{Average Parallelism} = \frac{\text{Total Work}}{\text{Critical Path Length}} = \frac{N - 1}{\log_2 N}$
举例当 $N = 1{,}000{,}000$ ，平均并行度为 $999,99920≈50,000\frac{999{,}999}{20} \approx 50{,}000$ ，这是理论上可以支持的平均并行度，实际上第一层最多需要 $N / 2$ 个线程同时进行操作，对于 $N = 1{,}000{,}000$ ，峰值为 500,000 个线程。但在后续阶段，线程数减半、减半……许多线程资源将闲置，这不是资源高效（resource-efficient） 的使用方式。

使用共享内存进行in-place归约

原始向量在设备全局内存中，共享内存用于保存部分求和向量，一开始，部分求和向量只是原始向量，每一步都使部分求和向量更接近总和，最终总和将在部分求和向量的元素0中，由于读取和写入部分的求和值，减少了全局内存流量。

在这里插入图片描述

__global__ void reduce_sum(float *input, float *globalSum, int N) {
    __shared__ float partialSum[2 * BLOCK_SIZE];

    unsigned int tid = threadIdx.x;
    unsigned int start = 2 * blockIdx.x * blockDim.x;
    
    // 边界检查避免越界
    partialSum[tid] = 
        (start + tid < N) ? input[start + tid] : 0.0f;
    partialSum[blockDim.x + tid] = 
        (start + blockDim.x + tid < N) ? input[start + blockDim.x + tid] : 0.0f;

    // In-place 归约
    for (unsigned int stride = 1; stride <= blockDim.x; stride <<= 1) {
        __syncthreads();	// 确保在进行下一步之前已经生成了每个版本的部分求和的所有元素
        if (tid % stride == 0) {
            partialSum[2 * tid] += partialSum[2 * tid + stride];
        }
    }

    // 每个 block 的线程 0 使用原子加将结果写入全局求和变量
    if (tid == 0) {
        atomicAdd(globalSum, partialSum[0]);
    }
}

每次迭代的stride变大，满足 tid % stride == 0 的线程数会越来越少，那么同一个warp（线程束）内，只有部分线程满足 if 条件执行加法，一部分线程不满足条件，即使没执行实际计算，也发生了分支发散。

比stride=32时，每个 warp 中最多只有一个线程在做事，其它线程不满足条件发散。当stride = 64, 128, 256, 512, 1024，一个些warp中的所有线程都跳过了 if，所以虽然不发散，但白白浪费资源。

可以改变索引的使用以改善发散行为，始终将部分求和结果压缩到partialSum[ ]数组的前面位置，保持活动线程连续。

在这里插入图片描述

__global__ void reduce_sum(float *input, float *globalSum, int N) {
    __shared__ float partialSum[2 * BLOCK_SIZE];

    unsigned int tid = threadIdx.x;
    unsigned int start = 2 * blockIdx.x * blockDim.x;
    
    // 边界检查避免越界
    partialSum[tid] = 
        (start + tid < N) ? input[start + tid] : 0.0f;
    partialSum[blockDim.x + tid] = 
        (start + blockDim.x + tid < N) ? input[start + blockDim.x + tid] : 0.0f;

    // In-place 归约
    for (unsigned int stride = blockDim.x; stride > 0; stride >>= 1) {
        __syncthreads();	// 确保在进行下一步之前已经生成了每个版本的部分求和的所有元素
        if (tid % stride == 0) {
            partialSum[2 * tid] += partialSum[2 * tid + stride];
        }
    }

    // 每个 block 的线程 0 使用原子加将结果写入全局求和变量
    if (tid == 0) {
        atomicAdd(globalSum, partialSum[0]);
    }
}

对于一个包含 1024 个线程的线程块，前 6 步中没有发散现象，每一步有1024、512、256、128、64、32 个连续线程处于活动状态，每个线程束中的所有线程都处于活动状态或全部处于非活动状态，后 5 步中仍将有发散现象。

四、内存并行性

DRAM

在这里插入图片描述

DRAM一个地址只能从内存单元核心阵列（二维阵列）中取出一个bit，平常内存是将多个内存单元核心阵列组合在一起，例如8个，那么一个地址就能取出一个Byte，当然也可以是核心阵列的位线组成一个Byte，需要看实际实现。实际的 DRAM 模块是将多个这样的阵列组成一个“总线宽度”（通常是 8-bit、16-bit、32-bit 或 64-bit）。比如 DDR3/DDR4 DIMM 通常是 64-bit = 8 Byte 的总线宽度。

DRAM迸发

DRAM 的核心结构是一个二维的内存单元矩阵，每次访问前，需要激活一个“行”，把这整行的数据加载到行缓冲区，这个操作由于开行延迟比较慢，但一旦行缓冲加载完，就可以连续读多列数据，DRAM的迸发就是通过设置的迸发长度，顺序读取该行中连续的几个列，迸发读取不代表从 DRAM 核心一次就“并行”取出多个bit，而是先取到行缓冲中，然后依次或流水线式地输出多个bit，再结合上面提到的总线宽度为64，那么就可以依次或流水线式地输出多个8Byte，如下图所示将一个红方块看作一个8Byte。

在这里插入图片描述

又因为即使迸发，核心阵列访问延迟也不会减少，因此，实际DRAM会有多个内存通道，可以理解为多根总线，不同总线流水线迸发，就可以做到减少核心阵列访问延迟。

在这里插入图片描述

示例

NVIDIA GTX280 GPU全局内存带宽峰值 = 141.7GB/s

全局内存接口对于典型的64位接口带宽17.6 8GB/s，我们需要更多的带宽（141.7 GB/s），因此需要 8 个内存通道。

每个地址空间被划分为多个迸发部分，每当访问一个位置时，同一部分中的所有其他位置也会传送给处理器，示例16 字节的地址空间，4 字节的迸发部分，如下图所示。事实上，我们至少有 4GB 地址空间，迸发部分大小为 128 字节或更多。

在这里插入图片描述

合并访存

在这里插入图片描述

当一个线程束的所有线程都执行一条加载指令时，如果所有访问的位置都落入同一个迸发部分，那么只会发出一个DRAM请求，访问将被全部合并。

不合并访存

在这里插入图片描述

当访问的位置跨越迸发段边界时，合并失败，访问和传输的某些字节未被线程使用。

合并访存方法

如果数组访问中的索引采用以下形式，则线程束中的访问是对连续位置的访问。

在这里插入图片描述

那么当矩阵乘法读入矩阵M一行时，访问合并，而读入N的一列时，访问不合并，分块矩阵乘法分块也可以如此分析。

第五讲——并行计算算法

一、并行直方图

并行直方图计算模式是一种重要且有效的计算方法，就每个线程的输出行为而言，与我们目前介绍的所有模式都有很大不同：输出可以被所有参与的线程修改。

接下来以文本直方图为例分析，其它直方图可以类比。定义bins为划分字母表的若干个部分，每个部分包含四个字母：a-d、e-h、i-l、m-p、q-t、u-x、y-z。对于输入字符串中的每个字符，相应的bin计数增加。在短语 “Programming Massively Parallel Processors”中，输出直方图如下所示：

在这里插入图片描述

数据划分

分段分区

在这里插入图片描述

分段分区导致内存访问效率低下，相邻线程不访问相邻内存位置，访问未被合并，DRAM 带宽利用率低。

在这里插入图片描述

交错分区

在这里插入图片描述

交错分区所有线程处理元素的连续部分，他们都移动到下一部分并重复处理，内存访问被合并，DRAM 带宽利用率高。

在这里插入图片描述

数据竞态条件

当多个线程并发访问相同的内存位置，且至少有一个线程是写操作，又没有正确同步机制时，就会发生数据争用。这和操作系统里面进程的数据竞态相同，这里不展开。

原子操作是一种不可中断的读-改-写操作，确保多个线程对同一个内存地址的并发访问是安全的。

CUDA 提供了一系列原子操作的内建函数，这里介绍原子加法atomicAdd(int* address, int val)从全局或共享内存中地址指向的位置读取 32 位变量 old，计算 (old + val)，并将结果存储回同一地址的内存中。该函数返回变量 old。还有其它原子操作可以查看官方文档。

基本的直方图内核

__global__ void histo_kernel(unsigned char *buffer, long size, unsigned int *histo) {
    unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
    // 交错划分
    int stride = blockDim.x * gridDim.x;
    
    while(i < size) {
        int alphabet_position = buffer[i] - "a";
        if (alphabet_position >= 0 && alphabet_position <26)
            // 原子加法
            atomicAdd(&(histo[alphabet_position/4]), 1);
        i += stride;
    }
}

原子操作的内存层次结构与性能影响

虽然原子操作可以保证数据一致性，但它们不是无代价的，通常有较高的延迟和较低的吞吐率，尤其是在高并发场景，如果多个线程访问同一个地址（同一个 bin），原子操作会产生冲突，硬件必须进行序列化，导致性能下降。

全局内存/DRAM上的原子操作

最通用，但性能最差。
延迟：几百个时钟周期（300-800+ cycles）。
如果多个线程竞争访问同一个地址，会严重拖慢吞吐量。

L2 缓存层的原子操作

L2 缓存可以帮助降低全局原子访问延迟。
CUDA 11 之后，大多数全局内存原子操作可以直接发生在 L2，绕过 DRAM。
比纯全局内存快，但仍较共享内存慢。

共享内存上的原子操作

最快的原子操作，适用于block 内线程之间的并发更新。
延迟低，吞吐量高（尤其在现代 GPU 架构如 Ampere 上）。
缺点是只能用于单个 block 内部，不能跨 block。

原子操作位置	延迟（cycles）	并发吞吐量（线程间争用时）
Shared Memory	10-50	非常高（低冲突时接近无损）
L2 Cache	100-300	中等（取决于访问模式）
Global Memory	300-800+	最低（冲突时吞吐极低）

私有化直方图内核

在这里插入图片描述

私有化是一种用于并行化应用程序的强大且常用的技术

私有化的成本

创建和初始化私有副本的开销
将私有副本的内容累积到最终副本的开销

私有化的收益

访问私有副本和最终副本时的数据争用和串行化要少得多
整体性能往往可以提升10倍以上

私有化的局限

私有直方图大小需要很小，适合共享内存

#define NUM_BINS 7

__global__ void private_histo_kernel (unsigned char *buffer, 
                                       long size, 
                                       unsigned int *histo) {
    __shared__ int private_histo[NUM_BINS];

    if (threadIdx.x < NUM_BINS) private_histo[tid] = 0;
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;
    
    __syncthreads();
    while(i < size) {
        int alphabet_position = buffer[i] - "a";
        if (alphabet_position >= 0 && alphabet_position <26)
            // 共享内存原子加法
            atomicAdd(&(private_histo[alphabet_position/4]), 1);
        i += stride;
    }
    __syncthreads();

    // Step 3: 合并私有直方图到全局直方图
    if (threadIdx.x < NUM_BINS) {
        atomicAdd(&histo[tid], private_histo[tid]);
    }
}

二、并行扫描算法

扫描/前缀和算法

定义

给定一个长度为 $n$ 的数组 $[x0,x1,…,xn−1][x_0, x_1, \ldots, x_{n-1}]$ ，和一个二元结合运算（如加法、乘法、最大值等），扫描操作返回一个新数组：
$[y_0, y_1, \ldots, y_{n-1}] = [x_0, x_0 \oplus x_1, \ldots, x_0 \oplus x_1 \oplus \ldots \oplus x_{n-1}]$
例如，对于加法操作：

输入： $3, 1, 7, 0, 4, 1, 6, 3$
输出： $3, 4, 11, 11, 15, 16, 22, 25$

串行代码

y[0] = x[0];
for (int i = 1; i < n; i++) {
    y[i] = y[i - 1] + x[i];
}

时间复杂度：O(N)，需要 N-1 次加法

朴素并行算法

每个线程 $T_i$ 负责计算 $y_i$ ，但需要遍历 $x_0$ 到 $x_i$ ，时间复杂度为 O(N)（总共需要 $N(N+1)2\frac{N(N+1)}{2}$ 次加法），并没有比串行更快。

重要性

扫描操作是许多高级并行算法的核心步骤。比如：

并行排序（如 radix sort）
图算法（如连通分量分析、BFS）
稀疏矩阵运算（如压缩行存储 CSR 的构建）
流处理和分段处理
动态分配任务/资源（如内存空间、数组分区等）

并行扫描是并行计算中的“积木”，不是最终目标，但几乎所有大厦都要用它来搭建。

低效并行扫描内核

__global void inefficient_scan_kernel(float *X, float *Y, int InputSize) {
    __shared__ float XY[SECTION_SIZE];  // SECTION_SIZE = blockDim.x
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    XY[threadIdx.x] = (i < InputSize) ? X[i] : 0.0f;
    __syncthreads();

    for (unsigned int stride = 1; stride <= threadIdx.x; stride <<= 1) {
        float in1 = XY[threadIdx.x - stride];
        __syncthreads();
        XY[threadIdx.x] += in1;
        __syncthreads();
    }

    if (i < InputSize) Y[i] = XY[threadIdx.x];
}

在这里插入图片描述

这个并行算法的本质，是将当前位置及之前的数据进行归约，因此这个核函数只能进行一个线程块内的扫描。

此扫描算法执行 $l o g (n)$ 次并行迭代，迭代过程分别进行了 $(n - 1), (n - 2), (n - 4), \dots, (n - n / 2)$ 次加法，共计加法次数 $n * l o g (n) - (n - 1)$ 是 $O (n * l o g (n))$ 复杂度。此扫描算法运行效率不高，顺序扫描算法才执行n次加法。虽然理论加法次数少，但由于同步多、依赖重、并行度差，在实际 GPU 上常常慢于预期。

当执行资源因工作效率低而饱和时，并行算法可能比顺序算法慢。

高效并行扫描内核

__global void efficient_scan_kernel(float *X, float *Y, int InputSize) {
    __shared__ float XY[2 * BLOCK_SIZE];
    unsigned int tid = threadIdx.x;
    int globalOffset = blockIdx.x * 2 * BLOCK_SIZE;
    
    if (globalOffset + tid < InputSize) {
    	XY[tid] = input[globalOffset + tid];
    } else {
        XY[tid] = 0.0f;
    }
    if (globalOffset + blockDim.x + tid < InputSize) {
    	XY[BLOCK_SIZE + tid] = input[globalOffset + BLOCK_SIZE + tid];
    } else {
        XY[BLOCK_SIZE + tid] = 0.0f;
    }
    __syncthreads();
    
    for (unsigned int stride = 1; stride <= BLOCK_SIZE; stride <<= 1) {
        int index = (tid + 1) * stride * 2 - 1;
        if (index < 2 * BLOCK_SIZE) 
            XY[index] += XY[index - stride];
        __syncthreads();
    }
    
    for (unsigned int stride = BLOCK_SIZE / 2; stride > 0; stride >>= 1) 
        int index = (tid + 1) * stride * 2 - 1;
        if (index + stride < 2 * BLOCK_SIZE) 
            XY[index + stride] += XY[index];
        __syncthreads();
    }

    if (globalOffset + tid < InputSize) 
        Y[globalOffset + tid] = XY[tid];
	if (globalOffset + BLOCK_SIZE + tid < InputSize) 
        Y[globalOffset + BLOCK_SIZE + tid] = XY[BLOCK_SIZE + tid];
}