type
status
date
slug
summary
tags
category
icon
password
声明
本文修改自以下知乎专栏, 主要润色了语言风格, 改善了可读性.
前言
2006年, NVIDIA公司发布了CUDA, 即建立在NVIDIA的GPUs上的一个通用并行计算平台和编程模型, 基于CUDA编程可以利用GPUs的并行计算引擎来更加高效地解决比较复杂的计算难题. 近年来, GPU最成功的一个应用领域就是深度学习领域, 基于GPU的并行计算已经成为训练深度学习模型的标配.
GPU并不是一个独立运行的计算平台, 而是需要与CPU协同工作, 可以视作CPU的协处理器. 因此当我们在说GPU并行计算时, 其实指的是CPU+GPU的异构计算架构. 在该异构计算架构中, GPU与CPU通过PCIe总线连接在一起来协同工作, CPU所在位置称为为主机端(host), 而GPU所在位置称为设备端(device), 如下图所示.

可以看到GPU包括更多的运算核心, 特别适合(控制逻辑简单的)(数据大量并行的)计算密集型任务, 例如大型矩阵运算; 而CPU的运算核心较少, 但是可以实现复杂的逻辑运算, 因此更适合(控制逻辑复杂的)(数据串行的)控制密集型任务. 另外, CPU上的线程是重量级的, 上下文切换开销大, 但是GPU由于存在很多核心, 其线程是轻量级的, 上下文切换开销小. 因此, CPU+GPU的异构计算平台可以实现优势互补, CPU负责处理逻辑复杂的串行程序, 而GPU则重点处理数据密集型的并行计算程序, 从而发挥各自的优势.

CUDA是NVIDIA公司所开发的GPU编程模型, 它提供了简易的GPU编程接口, 基于CUDA编程可以构建基于GPU计算的应用程序. CUDA提供了对C/C++, Python, Fortran等编程语言的支持, 这里我们选择CUDA C/C++接口对CUDA编程进行讲解, 开发平台为Windows 10 + VS 2013.

CUDA编程模型基础
基本概念
在给出CUDA的编程实例之前, 这里先对CUDA编程模型中的一些概念及基础知识做个简单介绍.
CUDA编程模型是一个异构模型, 需要CPU和GPU协同工作. 在CUDA中, host和device这两个概念非常重要. 我们用host指代CPU及其内存, 用device指代GPU及其内存. CUDA程序中既包含host程序, 又包含device程序, 它们分别在CPU和GPU上运行. 同时, host与device之间可以进行通信和数据拷贝. 典型的CUDA程序的执行流程如下:
- 分配host内存, 并进行数据初始化
- 分配device内存, 并将数据从host拷贝到device
- 调用CUDA的核函数, 在device上完成指定的运算
- 将device的运算结果拷贝到host
- 释放device和host上分配的内存
上面流程中最重要的一个环节是调用CUDA的核函数来执行并行计算. 核(kernel)是CUDA中一个重要概念, 是指要在device上大量并行线程中并行执行的函数. 核函数声明时需要加上
__global__符号, 在调用时需要用<<<grid, block>>>来指定该核函数执行时所需要的线程数量. 在CUDA中, 每一个执行核函数的线程都具有一个唯一的线程号thread ID, 这个ID值可以通过(正在执行的?)核函数的内置变量threadIdx获得. 因为GPU是异构模型, 所以编码时需要区分期望在host或device上运行的代码. CUDA通过函数类型限定词来区别将在host或device上运行的函数. 以下是三个主要的函数类型限定词:
__global__: 期望由host调用(一些特定的GPU也可以由device调用), 在device上执行, 返回类型必须是void, 不支持可变参数, 不能成为类成员函数, 不支持static变量, 只能访问device端内存. 注意使用__global__定义的核函数是异步的, 这意味着host不等该核函数执行完就会执行下一步.
__device__: 只能由device调用, 在device上执行. 不可以和__global__同时用.
__host__(默认类型): 只能由host调用, 在host上执行. 一般省略不写, 不可以和__global__同时用, 但可和__device__同时使用, 此时函数会在device和host上各自编译一份.
CUDA的线程层次结构
要深入理解kernel, 必须要对kernel的线程层次结构有一个清晰的认识.
首先, GPU上运行的是很多并行化的轻量级线程. 在device上执行核函数时, 实际上是启动很多线程, 所有线程同时执行该核函数. 一个核函数所启动的所有线程称为一个网格(grid), 同一个网格上的所有线程共享相同的全局内存空间(global memory). 网格是线程结构的第一层次, 而网格内部又可以分为很多线程块(block), 一个线程块里面包含很多线程(thread), 这是第二个层次.
如下图所示, 这是一个
grid和block均为二维(2-dim)的线程组织. <<<grid, block>>>标识中, grid和block都是类型为dim3的变量. dim3可以看成是包含三个无符号整数成员 的结构体变量. 在定义时, 缺省值将自动初始化为1, 因此grid和block可以灵活地定义为1-dim, 2-dim或3-dim结构. 对于图中结构(取水平方向为x轴), 定义的grid和block如下所示. 核函数kernel_fun在调用时必须通过加上<<<grid, block>>>来指定要使用的线程数及结构. Debug小技巧: 调试时可以将核函数设置为单线程.
kernel_fun<<<1,1>>>(params...) 
所以, 一个线程(thread)需要两个内置的坐标变量(
blockIdx, threadIdx)来唯一标识. 它们都是类型为dim3的变量, 其中blockIdx指明线程所在的block在grid中的坐标, 而threaIdx指明线程本身在block中的坐标, 如图中的Thread(1, 1)满足:一个线程块(block)上的线程是放在同一个流式多处理器(SM)上执行的, 但是单个SM处理器的资源有限, 所以单个线程块中的线程数是有限制的(现代GPUs的线程块可支持的线程数可达1024个). 有时候, 我们需要知道一个线程在块中的全局ID, 此时就必须先要知道block的组织结构, 可以通过线程的内置变量
blockDim来获得, 它记录线程块各个维度的大小. 对于2-dim的块 , 线程的ID值为 ; 对于3-dim的块 , 线程 的ID值为 . 另外线程还有内置变量gridDim, 它记录网格各个维度的大小. kernel的这种线程组织结构天然适合向量和矩阵运算. 我们可以利用上图的2-dim结构实现两个矩阵的加法, 每个线程负责处理矩阵中一个位置的两个元素相加, 代码如下所示. 我们定义线程块大小为
(16, 16), 然后将N * N大小的矩阵均分为多个线程块来执行加法运算. CUDA的内存模型
这里简单介绍一下CUDA的内存模型. 如下图所示, 可以看到, 每个线程有自己私有的本地内存(Local Memory), 而每个线程块还包含共享内存(Shared Memory),可以被线程块中所有线程共享, 其生命周期与线程块一致. 此外, 同一网格内所有的线程都可以访问该网格的全局内存(Global Memory). 除此之外网格还拥有一些只读内存块, 例如常量内存(Constant Memory)和纹理内存(Texture Memory). 内存结构涉及到程序优化, 这里不作深入探讨.

CUDA程序执行的物理层
还有一点很重要,你需要对GPU的硬件实现有一个基本的认识.
上面说到了kernel的线程组织层次, 从中可以看出调用一个核函数实际上会启动很多线程. 这些线程是逻辑上并行的, 但是在物理层却并不一定. 这其实和CPU的多线程有类似之处: 多线程如果没有多核支持, 在物理层也是无法实现并行的(只能分时调度). 好在GPU拥有很多CUDA核心, 只要能充分利用CUDA核心就可以充分发挥GPU的并行计算能力.
GPU硬件的核心组件叫做流式多处理器(Streaming Multiprocessor, SM). SM的核心组件包含CUDA核心, 共享内存和寄存器等. SM可以并发地执行数百个线程, 其并发能力就取决于它所拥有的资源数. 当一个核函数被执行时, 其网格中的各个线程块被分配到各个SM上, 一个线程块只能在一个SM上被调度, 而一个SM一般可以调度多个线程块(多对一). 一个核函数的各个线程块在实际运行时可能被分配给多个SM, 由此观之不难发现grid只是逻辑层, 而SM才是执行的物理层.
SM采用的是SIMT(Single-Instruction, Multiple-Thread, 单指令多线程)架构, 其基本执行单元是线程束(warps). 因此当线程块被分配到某个SM上时,它将进一步被划分为多个线程束, 因为这才是SM的基本执行单元. 一个线程束包含
32个线程, 这些线程同时执行相同的指令, 但是每个线程各自维护自己的指令地址计数器PC和寄存器状态, 也有自己独立的执行路径. 所以尽管线程束中的线程同时从同一程序地址执行,但是可能具有不同的行为, 比如遇到分支结构时, 一些线程可能进入其中一个分支执行, 但是期间另外一些线程有可能因为不进入该分支而不执行, 只能死等, 因为GPU规定线程束中所有线程在同一周期内必须执行相同的指令. 因此, 线程束的执行路径分化往往会导致性能下降. 一个SM上能够并行的线程束数量是有限的, 这是因为资源限制: SM要为每个线程块分配共享内存, 还要为各线程束中的每个线程分配独立的寄存器, 所以SM的性能会影响其所支持的并发运行的线程块和线程束数量. 总之, 网格和线程块只是逻辑划分, 实际上一个核函数的所有线程在物理层其实不一定是同时并发执行的. 所以核函数的grid和block的配置不同, 可能会导致性能出现差异, 这点需要特别注意. 另外, 由于SM的基本执行单元是包含
32个线程的线程束, 所以block的大小(即block中包含的线程数量)一般要设置为32的倍数. 
在进行CUDA编程前, 最好先检查一下自己的GPU的硬件配置, 这样才可以有的放矢. 可以通过下面的程序获得GPU的配置属性:
向量加法实例
知道了CUDA编程基础, 我们就来个简单的实战: 利用CUDA编程实现两个向量的加法. 在实现之前, 先简单介绍一下CUDA编程中内存管理的API. 首先是在device上分配内存的
cudaMalloc函数: 这个函数和C语言中的
malloc类似, 行为是在device上申请一定字节大小的显存, 其中devPtr是指向所分配内存的指针. 要释放分配的内存则使用cudaFree(void* devptr)函数, 和C语言中的free对应. 返回类型cudaError_t为一个整型枚举值, 编码CUDA API的执行结果, 可以通过将其传给以下API来获得人类可读的错误信息: 还有一个重要的函数
cudaMemcpy, 负责host和device之间的数据通信: 其中
src指向数据源, dst指向目标区域, count是复制的字节数, kind控制复制的方向, 类型有: cudaMemcpyHostToHost, cudaMemcpyHostToDevice, cudaMemcpyDeviceToHost及cudaMemcpyDeviceToDevice四种. 此外, 由于host端和device端的代码执行是异步的, 所以CUDA还设计了用于同步的API:
值得注意的是
cudaMemcpy函数带有隐式的同步语义, 即只有核函数执行完毕之后才能将结果拷贝回host端. 现在我们来实现一个向量加法的实例. 这里
grid和block都设计为1-dim. 首先定义kernel如下: 其中
stride是整个grid中包含的线程数. 有时候向量的分量数很多, 这时候可以在各个线程中执行多个元素(元素总数/线程总数)的加法, 相当于使用了多个grid来处理(实则也就是一个线程干多个线程的活…). 这是一种grid-stride loop方式. 不过下面的例子中, 一个线程只需要处理一个元素就能handle, 所以核函数中的for循环只会执行一轮就会退出. 下面我们给出计算向量加法的完整实现: 这里我们的向量大小为
1<<20, 而block大小为256,那么grid大小是4096. 故核函数的线程层级结构如下图所示: 
使用
nvprof工具可以分析核函数的运行情况, 结果如下所示, 可以看到调用核函数费时约1.5ms. 调整block的大小, 对比不同配置下的核函数运行情况, 这台机器上当
block为128时费时约1.6ms, 而block为512时费时约1.7ms, 当block为64时费时约2.3ms. 由此观之block并非越大越好, 而要根据实际情况适当设置. 在上面的实现中, 我们需要单独在host和device上进行各自的内存分配, 并且需要进行数据拷贝, 这很容易出错. 好在CUDA 6.0引入了统一内存(Unified Memory)来避免这种麻烦. 统一内存简单来说就是使用一个托管内存来同步管理host和device中的内存, 并且自动在host和device中进行数据传输. CUDA中使用
cudaMallocManaged函数分配托管内存: 利用统一内存,可以将上面的程序简化如下:
相比之前的代码, 使用统一内存的实现明显更简洁了🤠!
值得注意的是, device上核函数的执行与host端是异步的! 尽管托管内存会自动进行数据传输, 这里也要调用
cudaDeviceSynchronize()函数来保证device和host同步, 这样后面才可以正确访问核函数的计算结果. 矩阵乘法实例
最后我们再实现一个稍微复杂一些的例子: 两个矩阵的乘法. 设输入矩阵为 和 ,要得到 . 实现思路是每个线程计算 的一个元素值 . 对于矩阵运算, 应该选用2-dim的
grid和block. 首先定义矩阵的结构体:
然后我们来实现矩阵乘法的核函数. 这里我们定义了两个辅助的
__device__函数, 分别用于获取矩阵的元素值和为矩阵元素赋值. 具体代码如下: 最后我们采用统一内存编写矩阵相乘的测试实例:
这里矩阵大小为 , 设计的线程的
block大小为 , 那么可以算出grid大小为 . 最终测试结果如下: 当然, 这不是最高效的实现, 后面还可以再继续优化☺️.
小结
最后只有一句话: CUDA入门容易, 但是深入难! 希望你不是从入门到放弃...
参考资料
- John Cheng, Max Grossman, Ty McKercher. Professional CUDA C Programming, 2014.
- 作者:JAY
- 链接:https://notion-next-353pmh21m-jays-projects-ab02da23.vercel.app//article/2f63690b-3830-8045-9ac2-eb641b5436a0
- 声明:本文采用 CC BY-NC-SA 4.0 许可协议,转载请注明出处。


