易成 Institute of High Energy Physics

Presentation on theme: "易成 Institute of High Energy Physics"— Presentation transcript:

1 chengyi@nvidia.com 易成 Institute of High Energy Physics
CUDA C/C++ 编程介绍 易成 Institute of High Energy Physics

2 CUDA C/C++ 编程介绍：今天介绍的内容
基本概念 异构计算，主机和设备 __global__, __device__, __host__ 线程的两层组织结构，Grid 和 Block 内存组织架构 Shared memory/Global memory/Register Global memory/Constant memory/Texture memory Cache 内存管理 cudaMalloc(), cudaFree() cudaMemcpyAsync(), cudaMemcpy() __global__:主机端访问，设备端执行 __device__:设备端访问，设备端执行 __host__:主机端访问，主机端执行 学习资料链接：

3 CUDA? CUDA (Compute Unified Device Architecture) CUDA C/C++
使用 GPU 大量计算单元和高带宽实现通用计算的计算平台 提出了对硬件和软件要求 支持 CUDA 的GPU： CUDA: CUDA C/C++ 基于标准的 C/C++ CUDA C/C++ 是对它的一个扩展，实现异构计算 主要介绍 CUDA C/C++ 编程 CUDA是一个支持GPU计算核心进行计算的平台，是一个通用的高带宽计算平台。CUDA是一个支持GPU计算的非常重要的软件，现在的异构计算加速器有很多，包括像AMD的GPU，FPGA，以及国内的寒武纪的计算卡等等，但是他们没有类似于CUDA的计算平台。AMD是基于OpenCL进行编程开发，用这个平台进行开发非常繁琐，据我所知，苹果即将放弃使用OpenCL。 CUDA除了支持通用的Tesla系列GPU之外，还可以支持nvidia其他的GPU产品，包括quadro，GeForce，tegra等等。 CUDA C是一个对C/C++的扩展，可以基于CUDA编程，实现异构计算。

4 异构计算 线程块 线程 基本概念/术语 索引 内存组织 线程同步 异步执行

5 GPU 硬件架构 一个简单的例子： Hello, World! 内容框架 一个实际的例子： 矢量点乘 CUDA编程演示 总结

6 Heterogeneous Computing (异构计算)
术语: Host (主机) CPU 和它的内存 (host memory) Device (设备) GPU 和它的内存 (device memory) 介绍另个概念：主机和设备。 主机就是：服务器、笔记本、工作站等。有自己的CPU和内存； 设备：就是我们的GPU，GPU其实也是一个比较独立的计算部件，他有自己的计算核心和内存，但是没有独立的操作系统。 主机和设备共用同一套操作系统，在计算的时候，我们需要把数据从主机内存拷贝到设备，才能让GPU对数据进行计算。 Host Device

7 低延迟 VS 高吞吐量 CPU GPU 通过复杂的缓存体系结构，减小指令延迟 大量的晶体管处理逻辑任务 集成大量计算单元，获得高吞吐量
CPU的特点的是：有复杂的缓存体系，可以减小指令的延迟，但是占用了大量的晶体管，擅长处理复杂的逻辑任务。计算核心的主频高，但是数量少，现在Intel 铂金8180最多可到28核心，基准频率2.5GHz；AMD最多可以到32核心。 但是GPU是一个众核的架构，有大量的计算核心，相比较而言，缓存的比率比较小，因此，非常适合做数学计算。 在编程过程中，两种架构有什么区别呢？ 面对CPU架构，我们在编程的时候，我们是按照串行的思维模式进行编程的，去思考他的逻辑关系。即使我们使用MPI进行并行编程，也是把每个进程当做一个计算节点，互相通过消息传递进行通信，每个进程还是串行的执行。 面对GPU架构，我们编程的时候就需要考虑如果使用几百，甚至几千个核心进行并行计算。这也是我们今天要介绍的内容。 CPU计算的串行编程模式，就好比我们每个人，我们自己只需要考虑我们自己该做什么，怎么才能让自己的工作效率更高。 GPU计算的并行编程模式，就好比一家公司，需要考虑怎么给公司的每个员工分配任务，让公司的整体效率更高。 我们在考虑CPU平台上的编程模型时，实际上还是串行的思路，按照串行的执行顺序去编程。当然，现在CPU已经有很多核心了，但如果我们不额外添加并行语句，如MPI，或OpenMP，我们编写的C/C++代码还是串行的。但是在GPU平台，我们编程的时候就需要考虑多核并行。这是二者的区别。 那么，CPU架构与GPU架构又有什么不同呢？ 这里，我们给出了CPU与GPU的简易架构图，我们可以发现： CPU是缓存优化的处理器，缓存很大。而GPU是并行吞吐优化的处理器，计算核心很多 CPU将更多的晶体管用于缓存和流控，而GPU则更多的用于计算核心 CPU 通过复杂的缓存体系结构，减小指令延迟 大量的晶体管处理逻辑任务 GPU 集成大量计算单元，获得高吞吐量 更多的晶体管用于数学计算

8 Heterogeneous Computing
#include <iostream> #include <algorithm> using namespace std; #define N #define RADIUS 3 #define BLOCK_SIZE 16 __shared__ int temp[BLOCK_SIZE + 2 * RADIUS]; __global__ void stencil_1d(int *in, int *out) { int gindex = threadIdx.x + blockIdx.x * blockDim.x; int lindex = threadIdx.x + RADIUS; temp[lindex] = in[gindex]; // Read input elements into shared memory if (threadIdx.x < RADIUS) { temp[lindex - RADIUS] = in[gindex - RADIUS]; temp[lindex + BLOCK_SIZE] = in[gindex + BLOCK_SIZE]; } // Synchronize (ensure all the data is available) __syncthreads(); int result = 0; // Apply the stencil for (int offset = -RADIUS ; offset <= RADIUS ; offset++) result += temp[lindex + offset]; out[gindex] = result; // Store the result void fill_ints(int *x, int n) { fill_n(x, n, 1); int main(void) { int *in, *out; // host copies of a, b, c int size = (N + 2*RADIUS) * sizeof(int); int *d_in, *d_out; // device copies of a, b, c // Alloc space for host copies and setup values in = (int *)malloc(size); fill_ints(in, N + 2*RADIUS); out = (int *)malloc(size); fill_ints(out, N + 2*RADIUS); // Alloc space for device copies cudaMalloc((void **)&d_in, size); cudaMalloc((void **)&d_out, size); // Copy to device cudaMemcpy(d_in, in, size, cudaMemcpyHostToDevice); cudaMemcpy(d_out, out, size, cudaMemcpyHostToDevice); stencil_1d<<<N/BLOCK_SIZE,BLOCK_SIZE>>>(d_in + RADIUS, d_out + RADIUS); // Launch stencil_1d() kernel on GPU // Copy result back to host cudaMemcpy(out, d_out, size, cudaMemcpyDeviceToHost); // Cleanup free(in); free(out); cudaFree(d_in); cudaFree(d_out); return 0; parallel fn 一段异构计算的程序代码，代码中会同时有串行代码和并行代码，串行代码在 CPU 上执行；并行代码在 GPU 上并行执行。一般来说，我们会把代码中最消耗时间的部分并行化，在GPU上去执行。 serial code parallel code serial code

9 CUDA 编程的“三部曲” 申请GPU 内存，将数据从 CPU 内存拷贝到 GPU 内存 PCI Bus 三部曲的第一步：
拷贝数据的通道是PCI总线，因为现在Intel CPU还不支持nvlink，所以只能使用PCI-E总线进行数据拷贝。 申请GPU 内存，将数据从 CPU 内存拷贝到 GPU 内存

10 CUDA 编程的“三部曲” 申请GPU 内存，将数据从 CPU 内存拷贝到 GPU 内存
PCI Bus 指令和数据： 指令：加载GPU的计算程序，也就是kernel函数； 数据：从GPU内存读取数据进行计算。 申请GPU 内存，将数据从 CPU 内存拷贝到 GPU 内存 加载 GPU 程序 (kernel)，完成 GPU 上的并行计算

11 CUDA 编程的“三部曲” 申请 GPU 内存，将数据从 CPU 内存拷贝到 GPU 内存
PCI Bus 申请 GPU 内存，将数据从 CPU 内存拷贝到 GPU 内存 加载 GPU 程序 (kernel)，完成 GPU 上的并行计算。默认为异步执行模式 将结果从 GPU 内存拷回 CPU 内存，释放GPU 内存

12 GPU 内存组织 (Pascal) SM-0 SM-1 SM-N Registers Registers Registers L2
SMEM Tex/L1 C SMEM Tex/L1 C SMEM Tex/L1 CUDA编程与优化经常会与不同内存打交道，因此大家还需要对内存架构非常熟悉 整个GPU的内存体系包括了三个层次，与CPU类似。 片下的全局内存，全局内存存储空间很大，但是延迟很高，因此我们在片上和片下之间加入了L2的缓存。当GPU线程要在全局内存中读取数据的时候，首先是查看L2是否已经缓存了这些数据，如果有了，那么直接从L2中读取，速度更快。我们也称刚才的方式为缓存命中。 片上的寄存器、L1缓存、共享内存、常量和纹理的缓存。 还有多个SM共享的L2缓存 不同的内存有不同的特性，因此对不同的内存有不同的使用方法和优化方案，我们在编程优化中会讲到。 L2 Global Memory

13 GPU 内存组织 Global Memory Shared Memory Register
位于片外；延迟~100 时钟周期；所有线程可见 Shared Memory 位于片上；延迟~10 时钟周期；以线程块为单位分配；线程块内所有线程共享 Register 位于片上；延迟~1 时钟周期；以线程为单位分配；线程独享 Global/Constant/Texture Memory 都位于片外 经过不同的 L1 cache L1/constant/Texture Cache Conatant/Texture: Read-Only Memory

14 GPU P100 SM (Streaming Multiprocessor )
SM/GPU 56 FP32 units 64 FP64 units 32 TensorCore - Register/SM 256 KB Shared Memory/SM 64 KB L1 Cache 24 KB CUDA编程与优化经常会与不同内存打交道，因此大家还需要对内存架构非常熟悉 整个GPU的内存体系包括了三个层次，与CPU类似。 片下的全局内存，全局内存存储空间很大，但是延迟很高，因此我们在片上和片下之间加入了L2的缓存。当GPU线程要在全局内存中读取数据的时候，首先是查看L2是否已经缓存了这些数据，如果有了，那么直接从L2中读取，速度更快。我们也称刚才的方式为缓存命中。 片上的寄存器、L1缓存、共享内存、常量和纹理的缓存。 还有多个SM共享的L2缓存 不同的内存有不同的特性，因此对不同的内存有不同的使用方法和优化方案，我们在编程优化中会讲到。

15 Unified Shared Memory/
GPU V100 SM (Streaming Multiprocessor ) GV100 SM/GPU 80 FP32 units 64 FP64 units 32 TensorCore 8 Register/SM 256 KB Unified Shared Memory/ L1 Cache 128 KB CUDA编程与优化经常会与不同内存打交道，因此大家还需要对内存架构非常熟悉 整个GPU的内存体系包括了三个层次，与CPU类似。 片下的全局内存，全局内存存储空间很大，但是延迟很高，因此我们在片上和片下之间加入了L2的缓存。当GPU线程要在全局内存中读取数据的时候，首先是查看L2是否已经缓存了这些数据，如果有了，那么直接从L2中读取，速度更快。我们也称刚才的方式为缓存命中。 片上的寄存器、L1缓存、共享内存、常量和纹理的缓存。 还有多个SM共享的L2缓存 不同的内存有不同的特性，因此对不同的内存有不同的使用方法和优化方案，我们在编程优化中会讲到。

16 GPU CUDA 编程模型 Software GPU Thread blocks are executed on SM ...
Threads are executed by cuda core Thread CUDA Core 这里我们要讲的是：编程模型和GPU的对应关系。 首先概括一下这几个概念。其中SM（Streaming Multiprocessor）和SP（streaming Processor）是硬件层次的，其中一个SM可以包含多个SP。thread是一个线程，多个thread组成一个线程块block，多个block又组成一个线程网格grid。 现在就说一下一个kenerl函数是怎么执行的。一个kernel程式会有一个grid，grid底下又有数个block，每个block是一个thread群组。在同一个block中thread可以通过共享内存（shared memory）来通信，同步。而不同block之间的thread是无法通信的。 CUDA的设备在实际执行过程中，会以block为单位。把一个个block分配给SM进行运算；而block中的thread又会以warp（线程束）为单位，对thread进行分组计算。目前CUDA的warp大小都是32，也就是说32个thread会被组成一个warp来一起执行。同一个warp中的thread执行的指令是相同的，只是处理的数据不同。 基本上warp 分组的动作是由SM 自动进行的，会以连续的方式来做分组。比如说如果有一个block 里有128 个thread 的话，就会被分成四组warp，第0-31 个thread 会是warp 1、32-63 是warp 2、64-95是warp 3、 是warp 4。而如果block 里面的thread 数量不是32 的倍数，那他会把剩下的thread独立成一个warp；比如说thread 数目是66 的话，就会有三个warp：0-31、32-63、64-65 。由于最后一个warp 里只剩下两个thread，所以其实在计算时，就相当于浪费了30 个thread 的计算能力；这点是在设定block 中thread 数量一定要注意的事！ 一个SM 一次只会执行一个block 里的一个warp，但是SM 不见得会一次就把这个warp 的所有指令都执行完；当遇到正在执行的warp 需要等待的时候（例如存取global memory 就会要等好一段时间），就切换到别的warp来继续做运算，借此避免为了等待而浪费时间。所以理论上效率最好的状况，就是在SM 中有够多的warp 可以切换，让在执行的时候，不会有「所有warp 都要等待」的情形发生；因为当所有的warp 都要等待时，就会变成SM 无事可做的状况了。 实际上，warp 也是CUDA 中，每一个SM 执行的最小单位；如果GPU 有16 组SM 的话，也就代表他真正在执行的thread 数目会是32*16 个。不过由于CUDA 是要透过warp 的切换来隐藏thread 的延迟、等待，来达到大量并行化的目的，所以会用所谓的active thread 这个名词来代表一个SM 里同时可以处理的thread 数目。而在block 的方面，一个SM 可以同时处理多个thread block，当其中有block 的所有thread 都处理完后，他就会再去找其他还没处理的block 来处理。假设有16 个SM、64 个block、每个SM 可以同时处理三个block 的话，那一开始执行时，device 就会同时处理48 个block；而剩下的16 个block 则会等SM 有处理完block 后，再进到SM 中处理，直到所有block 都处理结束。 Thread blocks are executed on SM Thread Block SM ... A kernel is launched as a grid of thread blocks Grid Device

17 第一个例子：Hello World!

18 Hello World! (C 代码) 标准的纯 C 代码，在 CPU 上执行 printf("Hello World!\n");
int main(void) { printf("Hello World!\n"); return 0; } 标准的纯 C 代码，在 CPU 上执行

19 Hello World! (CUDA C 代码) 两个新的语法 __global__ void mykernel(void) {
printf("Hello World!\n"); } int main(void) { mykernel<<<1,1>>>(); return 0; 两个新的语法 知识点： 新语法__global__ 新语法<<<,>>> CUDA中的同步和异步的概念。必须在return 0之前加一个同步的语句。 #include <stdio.h> #include <cuda_runtime.h> __global__ void helloworld(void) { printf("hello world !!!\n"); } int main() helloworld<<<1,2>>>(); cudaDeviceReset(); #显示的释放GPU资源，或者在此处添加同步函数 return 0; cudaDeviceReset();这句话如果没有，则不能正常的运行，因为这句话包含了隐式同步，GPU和CPU执行程序是异步的，核函数调用后成立刻会到主机线程继续，而不管GPU端核函数是否执行完毕，所以上面的程序就是GPU刚开始执行，CPU已经退出程序了，所以我们要等GPU执行完了，再退出主机线程。

20 __global__ CUDA C/C++ 关键字 __global__ 修饰的函数 其它的函数修饰符
__global__ void mykernel(void) { printf("Hello World!\n"); } CUDA C/C++ 关键字 __global__ 修饰的函数 从主机端调用 (CC3.x, 也可以从设备端调用) 在设备端执行 Pascal GPU (CC6.X, P100: CC6.0, P40: CC6.1) 其它的函数修饰符 __device__: 设备端调用，设备端执行 __host__ : 主机端调用，主机端执行 __device__ 和 __host__ 可同时使用 (同时编译主机端和设备端两个版本) __global__是CUDA C的一个关键字，表示这个函数可以被主机端调用，在GPU上执行。 还有其他类似的修饰符：

21 <<<…>>>
mykernel<<<1,1>>>(); <<<X1, X2>>> 标记从主机端调用设备端函数以及相应的并行配置 也叫做 “kernel launch” (kernel 启动) X1 和 X2 分别为 kernel 的 grid (栅格) 和 block (线程块) 的设置 完成了在主机端调用、设备端执行一个 __global__ 函数过程 三个尖括号，表示kernel函数的并行配置。也就是说，这个函数需要用多少核心来，并行执行计算。 X1代表参数grid，X2代表参数block。 调用，启动内核函数

22 Hello World! (CUDA C 代码) Output: $ nvcc hello.cu $ a.out Hello World!
__global__ void mykernel(void) { printf("Hello World!\n"); } int main(void) { mykernel<<<1,1>>>(); return 0; mykernel 执行的任务很简单，“三部曲”只有“一部曲” Output: $ nvcc hello.cu $ a.out Hello World! $

23 一个实际的例子：矢量点乘 下面我们在看一个例子： 归约求和 sum a b c

24 一个实际的例子：矢量点乘 (只有一个分量)
GPU 上实现两个数相乘 __global__ void mul(int *a, int *b, int *c) { *c = (*a) * (*b); } __global__ 修饰的 mul 函数 mul() 在主机端调用 mul() 在设备端执行 我们来看一个例子： 两个矢量的点积计算

25 一个实际的例子：矢量点乘 (只有一个分量)
a, b, c 都是指针 __global__ void mul(int *a, int *b, int *c) { *c = (*a) * (*b); } mul() 在设备端执行, 因此 a, b 和 c 必须指向设备内存 因此，事先需要在设备端分配内存

26 内存管理 主机内存和设备内存位于不同的物理空间 内存管理的 CUDA API 设备 指针只能指向设备 内存 主机 指针只能指向主机 内存
可以在设备与主机之间传递 不能在主机端引用 (不能取地址里的内容) 主机 指针只能指向主机 内存 可以在主机与设备之间传递 不能在设备端引用 (不能取地址里的内容) 内存管理的 CUDA API cudaMalloc(), cudaFree(), cudaMemcpy() 与主机端的 malloc(), free(), memcpy()类似 知识点： 1、内存指针 2、内存空间的分配、释放、拷贝

27 矢量点乘 (只有一个分量) : kernel 回到 mul() kernel 设备端内存的申请在 main() 完成…
__global__ void mul(int *a, int *b, int *c) { *c = (*a) * (*b); } 设备端内存的申请在 main() 完成…

28 矢量点乘 (只有一个分量): main int main(void) { int a, b, c; // host copies of a, b, c int *d_a, *d_b, *d_c; // device copies of a, b, c int size = sizeof(int); // in bytes // Allocate space for device copies of a, b, c cudaMalloc((void **)&d_a, size); cudaMalloc((void **)&d_b, size); cudaMalloc((void **)&d_c, size); // Copy inputs to device cudaMemcpy(d_a, &a, size, cudaMemcpyHostToDevice); cudaMemcpy(d_b, &b, size, cudaMemcpyHostToDevice); Step 1: 申请 GPU 内存，CPU 到 GPU 数据拷贝 void **可以赋值给任何类型的变量 但是需要进行强制转换，不要(void **)也不会出错。 cudaMalloc((void **)&d_a, size);中其实省略（void **）也是可以运行的，因为定义的就是*d_a，&d_a就是 **

29 矢量点乘 (只有一个分量): main Step 2: 启动 kernel，完成设备端代码的计算
// Launch mul() kernel on GPU mul<<<1,1>>>(d_a, d_b, d_c); // Copy result back to host cudaMemcpy(&c, d_c, size, cudaMemcpyDeviceToHost); // Cleanup cudaFree(d_a); cudaFree(d_b); cudaFree(d_c); return 0; } Step 2: 启动 kernel，完成设备端代码的计算 Step 3: GPU 到 CPU 的数据拷贝，释放 GPU 内存

30 前面只用了一个block，一个线程计算，并没有充分发挥GPU的性能。下面我们来看一下怎么多核心的并行计算。
并行执行

31 矢量点乘 (多分量, N=512) : kernel GPU 的计算适合大量并发的任务 在 GPU 上执行 N 次
mul<<< 1, 1 >>>(); mul<<< N, 1 >>>(); 在 GPU 上执行 N 次 如何确保每次执行的是不同的数据呢？ 注意CUDA中一个GRID包含多个BLOCK。每个BLOCK包含多个THREAD。另外还有一个概念叫做wrap。WRAP是线程束的意思，也就是GPU实际执行运算的时候是以wrap单位的。比如wrap=32或者wrap=16.假设你的wrap是16，每次运行一定是16个线程一起运行。即使你使用了一个block里面的1个thread.GPU也会凑足16个thread，只是这些thread处于不活跃状态，此时就浪费了15个thread线程的资源。

32 矢量点乘 (多分量, N=512) : block, grid
<<<N, 1>>>: 并发执行 kernel mul N 次 术语: 在此，每次并发调用的 kernel 称为 一个 block 所有 block 的集合称为 grid <<<X1, X2>>> 中的 X1 设置 grid 的大小 每个 block 都有一个唯一标记的 ID blockIdx.x 可以用 blockIdx.x 去数组中索引对应的元素，以确保对应的数据参与计算 __global__ void mul(int *a, int *b, int *c) { c[blockIdx.x] = a[blockIdx.x] * b[blockIdx.x]; }

33 矢量点乘 (多分量, N=512) : kernel 设备端，所有的 block 并行执行:
__global__ void mul(int *a, int *b, int *c) { c[blockIdx.x] = a[blockIdx.x] * b[blockIdx.x]; } 设备端，所有的 block 并行执行: Block 0 Block 1 Block 2 Block 3 c[0] = a[0] * b[0]; c[1] = a[1] * b[1]; c[2] = a[2] * b[2]; c[3] = a[3] * b[3];

34 矢量点乘 (多分量, N=512) : 不能完成归约 在设备端完成了相应分量的乘积。
block 之间无法共享数据，不能在设备端完成归约求和操作 __global__ void mul(int *a, int *b, int *c) { c[blockIdx.x] = a[blockIdx.x] * b[blockIdx.x]; } 最终的归约求和在主机端完成。完整的 main 函数…

35 矢量点乘 (多分量, N=512) : main #define N 512 int main(void) {
int *a, *b, *c, sum=0; // host copies of a, b, c int *d_a, *d_b, *d_c; // device copies of a, b, c int size = N * sizeof(int); // Alloc space for device copies of a, b, c cudaMalloc((void **)&d_a, size); cudaMalloc((void **)&d_b, size); cudaMalloc((void **)&d_c, size); // Alloc space for host copies of a, b, c and setup input values a = (int *)malloc(size); random_ints(a, N); b = (int *)malloc(size); random_ints(b, N); c = (int *)malloc(size); void **可以赋值给任何类型的变量 但是需要进行强制转换，不要(void **)也不会出错。

36 矢量点乘 (多分量, N=512) : main // Copy inputs to device
cudaMemcpy(d_a, a, size, cudaMemcpyHostToDevice); cudaMemcpy(d_b, b, size, cudaMemcpyHostToDevice); // Launch mul() kernel on GPU with N blocks mul<<<N,1>>>(d_a, d_b, d_c); // Copy result back to host cudaMemcpy(c, d_c, size, cudaMemcpyDeviceToHost); for(int i=0; i<N; i++) sum += c[i]; // Reduce on Host // Cleanup free(a); free(b); free(c); cudaFree(d_a); cudaFree(d_b); cudaFree(d_c); return 0;}

37 矢量点乘 (多分量, N=512) : 设备端归约？ 能不能在设备端完成归约呢？<<<X1, X2>>>
X1: 设置 kernel 的 block 数目；block 之间不能通信 X2: 设置 kernel 的 block 内 thread (线程) 数； block 内线程间可以通过 shared memory (共享内存) 实现数据共享 mul<<<1, 1 >>>(); mul<<<1, N >>>(); block 内不同的线程用线程 ID 区分 threadIdx.x __shared__：声明共享内存 注意CUDA中一个GRID包含多个BLOCK。每个BLOCK包含多个THREAD。另外还有一个概念叫做wrap。WRAP是线程束的意思，也就是GPU实际执行运算的时候是以wrap单位的。比如wrap=32或者wrap=16.假设你的wrap是16，每次运行一定是16个线程一起运行。即使你使用了一个block里面的1个thread.GPU也会凑足16个thread，只是这些thread处于不活跃状态，此时就浪费了15个thread线程的资源。

38 设备端归约：kernel 设备端，所有的 thread 并行执行:
__global__ void mul(int *a, int *b, int *sum) { __shared__ int c[N]; c[threadIdx.x] = a[threadIdx.x] * b[threadIdx.x]; // Reduce code } 设备端，所有的 thread 并行执行: thread 0 thread 1 thread 2 thread 3 c[0] = a[0] * b[0]; c[1] = a[1] * b[1]; c[2] = a[2] * b[2]; c[3] = a[3] * b[3];

39 设备端归约实现 shared data thread ID shared data thread ID shared data

40 设备端归约代码 // Reduce code for (int i=N/2; i>0; i=i/2) {
if (threadIdx.x < i) { c[threadIdx.x] += c[threadIdx.x+i]; } if (threadIdx.x == 0) *sum = c[threadIdx.x]

41 设备端归约实现：问题？ shared data thread ID shared data thread ID shared data

42 设备端归约: 在哪同步呢? __syncthreads(); // Reduce code
for (int i=N/2; i>0; i=i/2) { if (threadIdx.x < i) { c[threadIdx.x] += c[threadIdx.x+i]; } if (threadIdx.x == 0) *sum = c[threadIdx.x]

43 设备端归约: kernel 完整代码 __global__ void mul(int *a, int *b, int *sum) { __shared__ int c[N]; c[threadIdx.x] = a[threadIdx.x] * b[threadIdx.x]; __syncthreads(); // Reduce code for (int i=N/2; i>0; i=i/2) { if (threadIdx.x < i) c[threadIdx.x] += c[threadIdx.x+i]; } if (threadIdx.x == 0) // shared memory to global memory *sum = c[threadIdx.x];

44 矢量点乘 (多分量, N=512) : main #define N 512 int main(void) {
int *a, *b, sum; // host copies of a, b int *d_a, *d_b, *d_sum; // device copies of a, b, sum int size = N * sizeof(int); // Alloc space for device copies of a, b, c cudaMalloc((void **)&d_a, size); cudaMalloc((void **)&d_b, size); cudaMalloc((void **)&d_sum, sizeof(int)); // Alloc space for host copies of a, b, c and setup input values a = (int *)malloc(size); random_ints(a, N); b = (int *)malloc(size); random_ints(b, N); 最大shared memory为49152 bytes，cuda C中int类型占4个字节，因此最大N为12288

45 矢量点乘 (多分量, N=512) : main // Copy inputs to device
cudaMemcpy(d_a, a, size, cudaMemcpyHostToDevice); cudaMemcpy(d_b, b, size, cudaMemcpyHostToDevice); // Launch mul() kernel on GPU with N threads mul<<<1, N>>>(d_a, d_b, d_sum); // Copy result back to host cudaMemcpy(&sum, d_sum, sizeof(int), cudaMemcpyDeviceToHost); // Cleanup free(a); free(b); cudaFree(d_a); cudaFree(d_b); cudaFree(d_sum); return 0; } mul<<<1, N>>>(d_a, d_b, d_sum);因为吃代码中，在设备端归约，最后if (threadIdx.x == 0) *sum = c[threadIdx.x],由于每个block都有一个threadIdx.x == 0，所以只能有一block参与计算，即gridDim.x=1；

46 小结 (1 of 2) 主机和设备 __global__ __device__ __host__ 主机 Host CPU & 内存
设备 Device GPU & 内存 __global__ __device__ __host__ __global__: 主机端调用、设备端执行 __device__: 设备端调用、设备端执行 __host__: 主机端调用、主机端执行 __host__ __device__ : 同时在设备端和主机端生成两个版本

47 小结 (2 of 2) 编写 CUDA 程序的三步曲 线程的两层组织结构 在设备端申请内存，并将数据从主机内存拷贝到设备内存
调用 kernel，完成设备端的并发计算 将计算结果从设备内存拷贝到主机内存，并释放内存 线程的两层组织结构 grid: block 的集合，由 blockIdx 标识。不同 block 之间不能通信 block: thread 的集合，由 threadIdx 标识 同一个 block 内的线程可通过 shared memory 通信 由 __sycnthreads() 同步

48 两种实现有什么不同呢？ 两者的性能地有巨大的差异 为什么呢？ 性能：多个线程块，每个只有一个线程 << 一个线程块，包含多个线程
Warp： 同一个线程块内的连续 32 个线程 Warp 是 GPU 上执行和调度单元 同一个 warp 内的所有线程执行相同的指令

49 Block 大小的设计 E.g. blockDim = 160 E.g. blockDim = 161 对于第一种情况 对于第二种情况
5 个 warp E.g. blockDim = 161 6 个 warp 对于第一种情况 N/32 个 warp, 没有浪费 对于第二种情况 N 个 warp, 31 / 32个thread 浪费 blockDim = 160表示每个block有160个thread

50 同时使用多个线程和线程块

51 矢量点乘：同时使用多个线程和线程块 已经介绍了 考虑采用多个 blocks 和 多个 threads 处理 N 很非常大的情形
多个 blocks, 每个 block 只有一个线程：无法归约 一个 block, 每个block 包含多个线程 ：block 中线程数目有限制，不能处理 N 非常大的情形 考虑采用多个 blocks 和 多个 threads 处理 N 很非常大的情形 首先考虑此时数据的索引，也就是线程与待处理数据之间的对应关系

52 多个 blocks 和 threads 时数据的索引
blockIdx.x 和 threadIdx.x 考虑有4个 blocks、每个有 8 threads线程与待处理数据的索引关系 threadIdx.x 1 7 2 3 4 5 6 Identical to finding offset in 1-dimensional storage of a 2-dimensional matrix: int index = x + width * y; blockDim.x=8,每个block有8个线程； mul<<< BLOCKS_NUM,THREADS_PER_BLOCK>>>，即mul<<<4,8>>> blockIdx.x = 0 blockIdx.x = 1 blockIdx.x = 2 blockIdx.x = 3 index = threadIdx.x + blockIdx.x * 8; 1 7 2 3 4 5 6

53 数组索引: 实例 哪个线程块中的线程会计算这个红色的数据？
1 31 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 blockDim.x M = 8 threadIdx.x = 5 由于blockIdx.x从0开始编号 1 7 2 3 4 5 6 blockIdx.x = 2 int index = threadIdx.x + blockIdx.x * 8; = * 8; = 21;

54 矢量点乘：同时使用多线程块和多线程 除 blockIdx.x 和 threadIdx.x 外，还要采用 blockDim.x 确定线程与待处理数据之间的对应关系 int index = threadIdx.x + blockIdx.x * blockDim.x; 在线程块内，使用 shared memory 实现线程间通信，完成线程块内部分和的计算 线程块之间，无法通信。因此，部分和的归约(即总和)，只能在主机端完成

55 同时使用多线程块和多线程：完整 kernel
__global__ void mul(int *a, int *b, int *sub_sum) { __shared__ int c[THREADS_PER_BLOCK]; int index = threadIdx.x + blockIdx.x * blockDim.x; c[threadIdx.x] = a[index] * b[index]; __syncthreads(); // Reduce code for (int i=THREADS_PER_BLOCK/2; i>0; i=i/2) { if (threadIdx.x < i) c[threadIdx.x] += c[threadIdx.x+i]; } if (threadIdx.x == 0) // shared memory to global memory sub_sum[blockIdx.x] = c[threadIdx.x] 假设THREADS_PER_BLOCK值为2的指数

56 同时使用多线程块和多线程：完整 main #define N (2048*2048)
#define THREADS_PER_BLOCK 512 #define BLOCKS_NUM (N / THREADS_PER_BLOCK ) int main(void) { int *a, *b, *sub_sum, sum=0; // host copies of a, b, sub_sum int *d_a, *d_b, *d_sub_sum; // device copies of a, b, sub_sum int size = N * sizeof(int); // Alloc space for device copies of a, b, sub_sum cudaMalloc((void **)&d_a, size); cudaMalloc((void **)&d_b, size); cudaMalloc((void **)&d_sub_sum, sizeof(int)* BLOCKS_NUM ); // Alloc for host copies of a, b, sub_sum and setup input values a = (int *)malloc(size); random_ints(a, N); b = (int *)malloc(size); random_ints(b, N); sub_sum = (int *)malloc(sizeof(int)* BLOCKS_NUM );

57 同时使用多线程块和多线程：完整 main // Copy inputs to device
cudaMemcpy(d_a, a, size, cudaMemcpyHostToDevice); cudaMemcpy(d_b, b, size, cudaMemcpyHostToDevice); // Launch mul() kernel on GPU mul<<< BLOCKS_NUM,THREADS_PER_BLOCK>>>(d_a, d_b, d_sub_sum); // Copy result back to host cudaMemcpy(sub_sum, d_sub_sum, sizeof(int)* BLOCKS_NUM , cudaMemcpyDeviceToHost); // Reduce on Host for(int i=0; i< BLOCKS_NUM; i++) sum += sub_sum[i]; // Cleanup free(a); free(b); free(sub_sum); cudaFree(d_a); cudaFree(d_b); cudaFree(d_sub_sum)；return 0;}

58 小结 多个线程和多个线程块协同 内建变量(built-in) 利用内建变量 blockDim.x，实现线程与待处理数据的索引关系
每个线程块完成部分和的归约 部分和归约成总和在主机端完成 内建变量(built-in) threadIdx : 线程 ID blockIdx : 线程块 ID blockDim : block的维度 gridDim : grid 的维度 这四个内建变量都是 dim3 型变量 (threadIdx.x; threadIdx.y; threadIdx.z)

59 问题?

60 举例 某高中总共有学生3000人，包含高一，高二，高三三个年级；每个年级10个班，每个班100人，按照5列20行排列座位。
按照CUDA的原理，可对应写为如下表达： dim3 grid(10,3) , block(5,20) gridDim.x=10,gridDim.y=3 blockDim.x=5, blockDim.y=20

61 一维和多维thread 1 7 2 3 4 5 6 如果定义：dim3 block(8,2)； 1 7 2 3 4 5 6 1 7 2 3
threadIdx.x blockDim.x= 8 blockDim.y= 1 1 7 2 3 4 5 6 blockIdx.x = 0 blockIdx.x = 1 blockIdx.x = 2 blockIdx.x = 3 如果定义：dim3 block(8,2)； 默认情况为一维thread，一维block 如果需要使用二维thread，需要定义dim3 block(2,3), 这样定义的话，blockDim.x=2,blockDim.y=3，即threadIdx.x为0~1，threadIdx.y为0~2， threadIdx.x threadIdx.y=0 1 7 2 3 4 5 6 blockDim.y=2 1 7 2 3 4 5 6 threadIdx.y=1 blockIdx.x = 0 blockIdx.x = 1 blockIdx.x = 2 blockIdx.x = 3 blockDim.x= 8

62 多维block 如果定义： int gridDimx=4; int gridDimy=2;
dim3 grid(gridDimx,gridDimy); 1 7 2 3 4 5 6 1 7 2 3 4 5 6 blockIdx.x = 0 blockIdx.x = 1 blockIdx.x = 2 blockIdx.x = 3 gridDim.y=2 blockIdx.y = 0 1 7 2 3 4 5 6 1 7 2 3 4 5 6 blockIdx.x = 0 blockIdx.x = 1 blockIdx.x = 2 blockIdx.x = 3 blockIdx.y = 1 gridDim.x=4

63 矢量点乘：同时使用多线程块和多线程