https://mp.weixin.qq.com/s/KgK3ertk9XVTxWhynv2AgA
本系列是为了弥补教程和实际应用之间的空白,帮助大家理解 CUDA 编程并最终熟练使用 CUDA 编程。你不需要具备 OpenGL 或者 DirectX 的知识,也不需要有计算及图形学的背景。
1 CPU 和 GPU 的基础知识
2 CUDA 编程的重要概念
3 并行计算向量相加
4 实践
4.1 向量相加 CUDA 代码
4.2 实践向量相加
5 给大家的一点参考资料
提到处理器结构,有2个指标是经常要考虑的:延迟和吞吐量。所谓延迟,是指从发出指令到最终返回结果中间经历的时间间隔。而所谓吞吐量,就是单位之间内处理的指令的条数。
下图1是 CPU 的示意图。从图中可以看出 CPU 的几个特点:
图1:CPU 的示意图
所以综合以上三点,CPU 在设计时的导向就是减少指令的时延,我们称之为延迟导向设计,如下图3所示。
下图2是 GPU 的示意图,它与之前 CPU 的示意图相比有着非常大的不同。从图中可以看出 GPU 的几个特点 (注意紫色和黄色的区域分别是缓存单元和控制单元):
图2:GPU 的示意图
所以,GPU 在设计过程中以一个原则为核心:增加简单指令的吞吐。因此,我们称 GPU 为吞吐导向设计,,如下图3所示。
图3:CPU 是延迟导向设计,GPU 是吞吐导向设计
那么究竟在什么情况下使用 CPU,什么情况下使用 GPU 呢?
CPU 在连续计算部分,延迟优先,CPU 比 GPU ,单条复杂指令延迟快10倍以上。
GPU 在并行计算部分,吞吐优先,GPU 比 CPU ,单位时间内执行指令数量10倍以上。
适合 GPU 的问题:
CUDA (Compute Unified Device Architecture),由英伟达公司2007年开始推出,初衷是为 GPU 增加一个易用的编程接口,让开发者无需学习复杂的着色语言或者图形处理原语。
OpenCL (Open Computing Languge) 是2008年发布的异构平台并行编程的开放标准,也是一个编程框架。OpenCL 相比 CUDA,支持的平台更多,除了 GPU 还支持 CPU、DSP、FPGA 等设备。
下面我们将以 CUDA 为例,介绍 GPU 编程的基本思想和基本操作。
首先主机端 (host) 和设备端 (device),主机端一般指我们的 CPU,设备端一般指我们的 GPU。
一个 CUDA 程序,我们可以把它分成3个部分:
第1部分是: 从主机 (host) 端申请 device memory,把要拷贝的内容从 host memory 拷贝到申请的 device memory 里面。
第2部分是: 设备端的核函数对拷贝进来的东西进行计算,来得到和实现运算的结果,图4中的 Kernel 就是指在 GPU 上运行的函数。
第3部分是: 把结果从 device memory 拷贝到申请的 host memory 里面,并且释放设备端的显存和内存。
图4:一个 CUDA 程序可以分成3个部分
CUDA 编程中的内存模型
这里就引出了一个非常重要的概念就是 CUDA 编程中的内存模型。
从硬件的角度来讲:
CUDA 内存模型的最基本的单位就是 SP (线程处理器)。每个线程处理器 (SP) 都用自己的 registers (寄存器) 和 local memory (局部内存)。寄存器和局部内存只能被自己访问,不同的线程处理器之间呢是彼此独立的。
由多个线程处理器 (SP) 和一块共享内存所构成的就是 SM (多核处理器) (灰色部分)。多核处理器里边的多个线程处理器是互相并行的,是不互相影响的。每个多核处理器 (SM) 内都有自己的 shared memory (共享内存),shared memory 可以被线程块内所有线程访问。
再往上,由这个 SM (多核处理器) 和一块全局内存,就构成了 GPU。一个 GPU 的所有 SM 共有一块 global memory (全局内存),不同线程块的线程都可使用。
上面这段话可以表述为:每个 thread 都有自己的一份 register 和 local memory 的空间。同一个 block 中的每个 thread 则有共享的一份 share memory。此外,所有的 thread (包括不同 block 的 thread) 都共享一份 global memory。不同的 grid 则有各自的 global memory。
图5:CUDA 内存模型,硬件角度
从软件的角度来讲:
图6:CUDA 内存模型,软件角度
如下图6所示,所谓线程块内存模型在软件侧的一个最基本的执行单位,所以我们从这里开始梳理。线程块就是线程的组合体,它具有如下这些特点:
如下图7所示的线程块就是由256个线程组成的,它执行的任务就是一个最基本的向量相加的一个操作。在线程块内,这256个线程的计算是彼此互相独立的,并行的。下面的这个 [i],就是如何确定每个线程的索引 (在显存中的位置)。在计算完以后 (图中弯箭头的头部),会设置一个时钟,将这256个线程的计算结果进行同步。
图7:一个256个线程组成的线程块
以上就是一个256位向量的加的操作的并行处理方法,得到最终的向量加的结果。
所谓网格 (grid),其实就是线程块的组合体,如下图8所示。
CUDA 核函数由线程网格 (数组) 执行。每个线程都有一个索引,用于计算内存地址和做出控制决策。在计算完以后 (图中所有弯箭头的头部),会设置一个时钟,将这N个线程块的计算结果进行同步。
图8:网格就是线程块的组合体
线程块 id & 线程 id:定位独立线程的门牌号
核函数需要确定每个线程在显存中的位置,我们之前提到 CUDA 的核函数是要在设备端来进行计算和处理的,在执行核函数时需要访问到每个线程的 registers (寄存器) 和 local memory (局部内存)。在这个过程中需要确定每一个线程在显存上的位置。所以我们需要像图9那样使用线程块的 index 和线程的 index 来确定线程在显存上的位置。
图9:使用线程块的 index 和线程的 index 来确定线程在显存上的位置
如图9所示,图9中的线程块索引是2维的,每个网格都由2×2个线程块组成;线程索引是3维的,每个线程块都由2×4×2个线程组成,所以代码应该是:
图10:线程Id计算
图10中:M=N=2,P,Q,S=2,4,2。
每个线程x的那一维应该是线程块的索引×线程块的x维度大小+线程的索引。(设备端线程x的那一维的索引)。
每个线程y的那一维应该是线程块的索引×线程块的y维度大小+线程的索引。(设备端线程y的那一维的索引)。
线程束 (warp)
前面我们提到,如图11所示的每一行由1个控制单元加上若干计算单元所组成,这些所有的计算单元执行的控制指令是一个。这其实就是个非常典型的 "单指令多数据流机制"。
图11:一个线程束 (warp):采用单指令多数据流机制
单指令多数据流机制是说:执行的指令是一条,只不过不同的计算单元使用的数据是不一样的。而上面这一行,我们就称之为一个线程束 (warp)。
所以,SM 采用的 SIMT (Single-Instruction, Multiple-Thread,单指令多线程) 架构,warp (线程束) 是最基本的执行单元。一个 warp 包含32个并行 thread,这些 thread 以不同数据资源执行相同的指令。一个 warp 只包含一条指令,所以:warp 本质上是线程在 GPU 上运行的最小单元。
由于warp的大小为32,所以block所含的thread的大小一般要设置为32的倍数。
当一个 kernel 被执行时,grid 中的线程块被分配到 SM (多核处理器) 上,一个线程块的 thread 只能在一个SM 上调度,SM 一般可以调度多个线程块,大量的 thread 可能被分到不同的 SM 上。每个 thread 拥有它自己的程序计数器和状态寄存器,并且用该线程自己的数据执行指令,这就是所谓的 Single Instruction Multiple Thread (SIMT),如图12所示。
图12:Single Instruction Multiple Thread (SIMT)
下面我们就用一个实际的例子来看看 CUDA 编程具体是如何操作的。例子就是两个长度为N的张量相加,如下图13所示。
图13:两个张量相加
在 CPU 中完成相加的操作很简单:
// Compute vector sum C = A+B void vecAdd(float* A, float* B, float* C, int n) { for (i= 0, i< n, i++) C[i] = A[i] + B[i]; } int main() { // Memory allocation for A_h, B_h, and C_h // I/O to read A_hand B_h, N elements … vecAdd(A_h, B_h, C_h, N); }
要在 GPU 中完成这一操作,首先我们想一下它是否适合使用 GPU,我们当时总结了四个特点:
所以,向量相家的任务适合在 GPU 上编程。
再回顾下 GPU 运算步骤,如图4所示:
一个 CUDA 程序,我们可以把它分成3个部分:
第1部分是: 从主机 (host) 端申请 device memory,把要拷贝的内容从 host memory 拷贝到申请的 device memory 里面。
第2部分是: 设备端的核函数对拷贝进来的东西进行计算,来得到和实现运算的结果,图4中的 Kernel 就是指在 GPU 上运行的函数。
第3部分是: 把结果从 device memory 拷贝到申请的 host memory 里面,并且释放设备端的显存和内存。
如下:
#include <cuda.h> void vecAdd(float* A, float* B, float* C, int n) { int size = n* sizeof(float); float* A_d, B_d, C_d; … 1. // Allocate device memory for A, B, and C // copy A and B to device memory 2. // Kernel launch code –to have the device // to perform the actual vector addition 3. // copy C from the device memory // Free device vectors }
下面我们把这些内容细化到函数。
设备端代码:
主机端代码:
内存是插在主板上的内存插槽上的内存条,而显存是独立显卡上焊在显卡上的内存芯片。
申请显存的函数 cudaMalloc():
在主机端完成显存的申请,得到相应的指针。
图14:申请显存的函数 cudaMalloc()
释放显存的函数 cudaFree( ):
将指向显存的指针释放掉。
图15:释放显存的函数 cudaFree( )
内存和显存之间互相拷贝的函数 cudaMemcpy( ):
参数含义是:终点的指针,起点的指针,拷贝的大小,模式 (主机端到设备端,设备端到主机端,设备端之间的拷贝)
图16:内存和显存之间互相拷贝的函数 cudaMemcpy( )
以上三个函数是 CUDA 帮我们写好的,如果调用的话需要先:
# include cuda.h
下面就是具体的 C++ 代码实现:
申请内存的大小是 n *sizeof(float),定义3个指针 A_d,B_d,C_d。
cudaMalloc 函数需要传入 1. 指针的指针 (指向申请得到的显存的指针)。2. 申请显存的大小。 所以分别传入 &A_d 和 size。同理后面依次传入 &B_d 和 size,&C_d 和 size。
cudaMemcpy 函数需要传入 1. 终点的指针。2. 起点的指针。3. 拷贝的大小。4. 模式。 所以分别传入 A_d, A, size, cudaMemcpyHostToDevice。同理后面依次传入 B_d, B, size, cudaMemcpyHostToDevice 和 C, C_d, size, cudaMemcpyHostToDevice。
最后把设备端申请的显存都释放掉。cudaFree 函数需要传入设备端申请显存的指针,即 A_d,B_d,C_d。
void vecAdd(float* A, float* B, float* C, int n) { int size = n * sizeof(float); float* A_d, *B_d, *C_d; 1. // Transfer A and B to device memory cudaMalloc((void **) &A_d, size); cudaMemcpy(A_d, A, size, cudaMemcpyHostToDevice); cudaMalloc((void **) &B_d, size); cudaMemcpy(B_d, B, size, cudaMemcpyHostToDevice); // Allocate device memory for cudaMalloc((void **) &C_d, size); 2. // Kernel invocation code –to be shown later … 3. // Transfer C from device to host cudaMemcpy(C, C_d, size, cudaMemcpyDeviceToHost); // Free device memory for A, B, C cudaFree(A_d); cudaFree(B_d); cudaFree(C_d); }
下面我们进入最重要的部分,即:如何自己书写一个 kernel 函数。
核函数调用的注意事项
CUDA 编程的标识符号
不同的表示符号对应着不同的工作地点和被调用地点。核函数使用 __global__ 标识,必须返回 void。__device__ & __host__ 可以一起用。
图17:CUDA 编程的标识符号
下面,按照我们刚才的对核函数的介绍,我们展示了向量相加的代码。
代码讲解:
首先,看到 __global__ 标识,返回的是 void,就意味着 vecAddKernel 函数是一个在 host 端调用,在 device 端执行的核函数。它的三个参数就是我们之前申请好的指向三段显存的指针。
通过 int i= threadIdx.x+ blockDim.x* blockIdx.x; (线程的索引,线程块的索引,线程块维度的大小) 来计算好要访问的线程的索引的位置。
那么如何在主机端调用呢?我们使用尖括号**<<<网格 grid 维度,线程块 block 维度>>>**来包括:线程块数 ceil(n/256) 和一个线程块的线程数256。
图18:向量相加的代码
第1步主机端 __host__ 修饰:申请显存,内存。显存,内存的互相拷贝。内存,显存释放。比如图19中申请的网格是 ceil(n/256) 维的代表一个网格有 ceil(n/256) 个线程块;线程块是256维的,代表一个线程块有256个线程。
第2步设备端 __global__ 修饰:计算索引绝对位置,并行计算。
图19:主机端和设备端代码
详细地讲,核函数只能在主机端调用,调用时必须申明执行参数。调用形式如下:
Kernel<<<Dg,Db, Ns, S>>>(param list);
<<<>>> 运算符内是核函数的执行参数,告诉编译器运行时如何启动核函数,用于说明内核函数中的线程数量,以及线程是如何组织的。
<<<>>> 运算符对 kernel 函数完整的执行配置参数形式是 <<<Dg, Db, Ns, S>>>
最后我们简单介绍下 CUDA 编程如何执行编译的过程。因为我们之前在 CPU 上编程,使用 g++ 或 gcc 进行编译,再通过 link 生成可执行程序。那么在 GPU 端,编译器就是 NVCC (NVIDIA Cuda compiler driver)。
通常我们会把和 GPU 相关的头文件放在 .h 文件里,把设备端执行的程序 (__global__ 定义的函数) 放在 .cu 文件里,这些程序我们用 NVCC 来进行编译。主机端的程序放在 .h 和 .cpp 里面,这些程序我们可以继续用 g++ 或 gcc 来进行编译。
通常我们有这几种编译的方法:
图20:CUDA 编程如何执行编译的过程
CUDA 中 threadIdx,blockIdx,blockDim,gridDim 的使用
下面这张图21比较清晰的表示的几个概念的关系:
图21:几个变量的关系
cuda 通过<<< >>>符号来分配索引线程的方式,我知道的一共有15种索引方式。
这一节我们通过一个实例直观感受下 CUDA 并经计算究竟能使这些计算简单,并行度高的操作加速多少。
我们先看一下 CPU 执行向量相加的代码:
#include <iostream> #include <cstdlib> #include <sys/time.h> using namespace std; void vecAdd(float* A, float* B, float* C, int n) { for (int i = 0; i < n; i++) { C[i] = A[i] + B[i]; } } int main(int argc, char *argv[]) { int n = atoi(argv[1]); cout << n << endl; size_t size = n * sizeof(float); // host memery float *a = (float *)malloc(size); float *b = (float *)malloc(size); float *c = (float *)malloc(size); for (int i = 0; i < n; i++) { float af = rand() / double(RAND_MAX); float bf = rand() / double(RAND_MAX); a[i] = af; b[i] = bf; } struct timeval t1, t2; gettimeofday(&t1, NULL); vecAdd(a, b, c, n); gettimeofday(&t2, NULL); //for (int i = 0; i < 10; i++) // cout << vecA[i] << " " << vecB[i] << " " << vecC[i] << endl; double timeuse = (t2.tv_sec - t1.tv_sec) + (double)(t2.tv_usec - t1.tv_usec)/1000000.0; cout << timeuse << endl; free(a); free(b); free(c); return 0; }
注释:
float*a =(float*)malloc(size); 分配一段内存,使用指针 a 指向它。
for 循环产生一些随机数,并放在分配的内存里面。
vecAdd(float* A,float* B,float* C,int n) 要输入指向3段内存的指针名,也就是 a, b, c。
gettimeofday 函数来得到精确时间。它的精度可以达到微妙,是C标准库的函数。
最后的 free 函数把申请的3段内存释放掉。
编译:
g++ -O3 main_cpu.cpp -o VectorSumCPU
我们再看一下 CUDA 执行向量相加的代码:
#include <iostream> #include <cstdlib> #include <sys/time.h> #include <cuda_runtime.h> using namespace std; __global__ void vecAddKernel(float* A_d, float* B_d, float* C_d, int n) { int i = threadIdx.x + blockDim.x * blockIdx.x; if (i < n) C_d[i] = A_d[i] + B_d[i]; } int main(int argc, char *argv[]) { int n = atoi(argv[1]); cout << n << endl; size_t size = n * sizeof(float); // host memery float *a = (float *)malloc(size); float *b = (float *)malloc(size); float *c = (float *)malloc(size); for (int i = 0; i < n; i++) { float af = rand() / double(RAND_MAX); float bf = rand() / double(RAND_MAX); a[i] = af; b[i] = bf; } float *da = NULL; float *db = NULL; float *dc = NULL; cudaMalloc((void **)&da, size); cudaMalloc((void **)&db, size); cudaMalloc((void **)&dc, size); cudaMemcpy(da,a,size,cudaMemcpyHostToDevice); cudaMemcpy(db,b,size,cudaMemcpyHostToDevice); cudaMemcpy(dc,c,size,cudaMemcpyHostToDevice); struct timeval t1, t2; int threadPerBlock = 256; int blockPerGrid = (n + threadPerBlock - 1)/threadPerBlock; printf("threadPerBlock: %d \nblockPerGrid: %d \n",threadPerBlock,blockPerGrid); gettimeofday(&t1, NULL); vecAddKernel <<< blockPerGrid, threadPerBlock >>> (da, db, dc, n); gettimeofday(&t2, NULL); cudaMemcpy(c,dc,size,cudaMemcpyDeviceToHost); //for (int i = 0; i < 10; i++) // cout << vecA[i] << " " << vecB[i] << " " << vecC[i] << endl; double timeuse = (t2.tv_sec - t1.tv_sec) + (double)(t2.tv_usec - t1.tv_usec)/1000000.0; cout << timeuse << endl; cudaFree(da); cudaFree(db); cudaFree(dc); free(a); free(b); free(c); return 0; }
注释:
首先要用 __global__ 来修饰。
vecAdd(float* A,float* B,float* C,int n) 要输入指向3段显存的指针名,也就是 d_a, d_b, d_c。
float*da =NULL; 定义空指针。
cudaMalloc((void**)&da, size); 申请显存,da 指向申请的显存,注意 cudaMalloc 函数传入指针的指针 (指向申请得到的显存的指针)。
cudaMemcpy(da,a,size,cudaMemcpyHostToDevice) 把内存的东西拷贝到显存,也就是把 a, b, c 里面的东西拷贝到 d_a, d_b, d_c 中。
int threadPerBlock =256; int blockPerGrid =(n + threadPerBlock -1)/threadPerBlock; 计算线程块和网格的数量。
vecAddKernel <<< blockPerGrid, threadPerBlock >>> (da, db, dc, n); 调用核函数。
gettimeofday 函数来得到精确时间。它的精度可以达到微妙,是C标准库的函数。
最后的 free 函数把申请的3段内存释放掉。
编译:
/usr/local/cuda/bin/nvcc main_gpu.cu -o VectorSumGPU
编译之后得到可执行文件 VectorSumCPU 和 VectorSumGPU 之后,我们可以执行一下比较下运行时间 (注意要在 linux 下运行):
在 CPU 下,执行1000000000次加需要4.18秒。
./VectorSumCPU 1000000000 4.18261
在 GPU 下,执行1000000000次加只需要1.6e-05秒,哇。
(base) wjh19@iccv:~/mage/CUDA/db$ ./VectorSumGPU 1000000000 threadPerBlock: 256 blockPerGrid: 3906250 1.6e-05
GPU 对于计算简单,并行度高的计算果然可以大幅提速!!!
在 CPU 下,执行1000次加需要1e-06秒。
(base) wjh19@iccv:~/mage/CUDA/db$ ./VectorSumCPU 1000 1e-06
在 GPU 下,执行1000次加需要1.3e-05秒。
(base) wjh19@iccv:~/mage/CUDA/db$ ./VectorSumGPU 1000 threadPerBlock: 256 blockPerGrid: 4 1.3e-05
GPU 对于少量计算效率反倒不如 CPU。
参考
1.深蓝学院课程讲解:https://www.shenlanxueyuan.com/course/410
2. D. Kirk and W. Hwu, “Programming Massively Parallel Processors –A Hands-on Approach, Second Edition”
3. CUDA by example, Sanders and Kandrot
4. Nvidia CUDA C Programming Guide:https://docs.nvidia.com/cuda/cuda-c-programming-guide/
5. CS/EE217 GPU Architecture andProgramming
如果觉得有用,就请分享到朋友圈吧!
计算机视觉SLAM 分享SLAM/SFM定位建图技术(特征匹配、位姿估计、图优化等)、视觉定位,深度感知(双目视觉、TOF、结构光等)以及XR(VR、AR、MR等)领域干货,探讨以上领域在学术界与产业界研究与应用。 53篇原创内容 公众号