CUDA程序设计

主要内容 GPGPU及CUDA介绍 CUDA编程模型 多线程及存储器硬件

GPGPU及CUDA介绍

多核时代 多个适当复杂度、相对低功耗内核并行工作 核心时钟频率基本不变 配置并行硬件资源提高处理能力 Quad-core Opteron IBM Cell Broadband Engine nVidia GT200

GPU与CPU硬件架构的对比 CPU：更多资源用于缓存及流控制 GPU：更多资源用于数据计算 适合具备可预测、针对数组的计算模式 CPU Cache ALU Control DRAM DRAM CPU GPU

应用范围 CPU: control processor 不规则数据结构 不可预测存取模式 递归算法 分支密集型算法 单线程程序 GPU: data processor 规则数据结构 可预测存取模式 油气勘探、金融分析、医疗成像、有限元、基因分析、地理信息系统、…

GPGPU (General Purpose Computing on GPU)

GPGPU 核心思想 缺点 用图形语言描述通用计算问题 把数据映射到vertex或者fragment处理器 硬件资源使用不充分 存储器访问方式严重受限 难以调试和查错 高度图形处理和编程技巧

CUDA (Compute Unified Device Architecture) CUDA有效结合CPU+GPU编程 串行部分在CPU上运行 并行部分在GPU上运行 CPU Serial Code Grid 0 . . . GPU Parallel Kernel KernelA<<< nBlk, nTid >>>(args); CPU Serial Code Grid 1 . . . GPU Parallel Kernel KernelB<<< nBlk, nTid >>>(args);

CUDA极大提高了现有应用的效果 MRI Reconstruction Cartesian Scan Data Spiral Scan Data Gridding1 (b) (b) (c) (a) FFT Iterative Reconstruction Spiral scan data + Gridding + FFT Reconstruction requires little computation Based on Fig 1 of Lustig et al, Fast Spiral Fourier Transform for Iterative MR Image Reconstruction, IEEE Int’l Symp. on Biomedical Imaging, 2004

Advanced MRI Reconstruction Cartesian Scan Data Spiral Scan Data Gridding (b) (c) (a) (b) FFT Iterative Reconstruction Spiral scan data + Iterative recon Reconstruction requires a lot of computation

Advanced MRI Reconstruction Compute Q Q只和扫描参数有关 FHd是数据相关的 使用线性求解器计算ρ More than 99.5% of time Acquire Data Compute FHd The least-squares reconstruction algorithm poses the reconstruction problem in the form FHF times p = FHd. 1. P is the desired image. 2. FHF is a special data structure that can be derived very quickly from another data structure, Q. Q depends only on the scanner configuration and can therefore be computed before the scan data is even acquired. So, while computing Q is very expensive, this cost can be amortized over many reconstructions, and does not factor into the reconstruction time for any image. 3. FHd is a matrix that depends on both the scanner configuration and the data. Computing FHd is fairly expensive. Finally, given FHF and FHd, a linear solver is used to find the image. This last step is quite fast. In fact, 99.5% of the reconstruction time for a single image is devoted to computing FHd, and computing Q is even more expensive, so those are the computations that we chose to accelerate on the G80. Find ρ Haldar, et al, “Anatomically-constrained reconstruction from noisy data,” MR in Medicine.

Code CPU GPU for (p = 0; p < numP; p++) { for (d = 0; d < numD; d++) { exp = 2*PI*(kx[d] * x[p] + ky[d] * y[p] + kz[d] * z[p]); cArg = cos(exp); sArg = sin(exp); rFhD[p] += rRho[d]*cArg – iRho[d]*sArg; iFhD[p] += iRho[d]*cArg + rRho[d]*sArg; } __global__ void cmpFhD(float* gx, gy, gz, grFhD, giFhD) { int p = blockIdx.x * THREADS_PB + threadIdx.x;   // register allocate image-space inputs & outputs x = gx[p]; y = gy[p]; z = gz[p]; rFhD = grFhD[p]; iFhD = giFhD[p]; for (int d = 0; d < SCAN_PTS_PER_TILE; d++) { // s (scan data) is held in constant memory float exp = 2 * PI * (s[d].kx * x + s[d].ky * y + s[d].kz * z); cArg = cos(exp); sArg = sin(exp); rFhD += s[d].rRho*cArg – s[d].iRho*sArg; iFhD += s[d].iRho*cArg + s[d].rRho*sArg; } grFhD[p] = rFhD; giFhD[p] = iFhD;

性能提升情况 S.S. Stone, et al, “Accelerating Advanced MRI Reconstruction using GPUs,” ACM Computing Frontier Conference 2008, Italy, May 2008.

计算结果对比

CUDA成功案例 广泛应用于生命科学、机械、石油、金融、数学、天文和通信等行业

医疗成像 Computed Tomography (CT) GE MRI (磁共振成像) GRAPPA 自动校准 加速网格化 快速重建 Digisens SnapCT Stone, UIUC Batenburg, Sijbers et al

量子化学 双电子积分 RI-MP2 correlation energy in Q-Chem 3.1 K Yasuda, Nagoya U, Japan Leslie Vogt, Harvard

分子动力学 现有的分子动力学软件 NAMD / VMD (alpha release) OpenMM: 分子建模 GROMACS (alpha release) HOOMD OpenMM: 分子建模 https://simtk.org/home/openmm

金融 Monte Calo模拟 投资组合优化 期权及衍生品定价 对冲基金 风险分析 CUDA中的 随机数发生器 SciFinance的Monte Calo定价模型 SciComp Co.

生物信息学和生命科学 序列对比 蛋白质对接 生物系统的随机仿真(SSA) 人体视觉皮层的自组织计算模型 分析基因表达的DNA微阵列工具 Schatz et al, U Maryland

流体动力学 3D Lattice-Boltzman解算器 基于Lattice-Boltzman的PDE解算器 Navier-Stokes解算器 等离子体湍流建模 Thibault and Senocak Tolke and Krafczy

电磁学和电磁力学 GPMAD: 离子束动力学模拟 FDTD法进行的光散射模拟 Acceleware的解算器 FDTD加速 Accelerware

天气, 大气, 海洋科学与空间建模 天气研究与预测模型 (WRF) 25% ~ 30%的性能提升 海啸模拟

加密编码

模式匹配

CUDA编程模型

CUDA设备与线程 计算设备（device） 计算密集部分使用大量线程并行的kernel GPU与CPU线程的区别 作为CPU（host）的协处理器 有独立的存储设备（device memory） 同时启动大量线程 计算密集部分使用大量线程并行的kernel GPU与CPU线程的区别 GPU的线程非常轻量，线程切换~1 cycle，而CPU需要~1000 cycle GPU上的线程数>1000时才能有效利用GPU的计算能力

Streaming Processor（SP） A fully pipelined, single-issue, inorder microprocessor 2 ALUs and a FPU Register file 32-bit scalar processing No instruction fetch and scheduling No cache

Streaming Multiprocessor（SM） An array of SPs 8 streaming processor 2 Special Function Units (SFU) A 16KB read/write shared memory Not a cache But a software-managed data store Multithreading issuing unit Instruction and constant cache

CUDA 程序基本结构 串行部分在CPU上运行（host） 并行部分在GPU上运行（device） CPU Serial Code (host) Grid 0 . . . GPU Parallel Kernel (device) KernelA<<< nBlk, nTid >>>(args); CPU Serial Code (host) Grid 1 . . . GPU Parallel Kernel(device) KernelB<<< nBlk, nTid >>>(args);

C扩展 Declspecs Keywords Intrinsics Runtime API Function launch global, device, shared, local, constant Keywords threadIdx, blockIdx Intrinsics __syncthreads Runtime API Memory, symbol, execution management Function launch __device__ float filter[N]; __global__ void convolve (float *image) { __shared__ float region[M]; ... region[threadIdx] = image[i]; __syncthreads() image[j] = result; } // Allocate GPU memory void *myimage = cudaMalloc(bytes) // 100 blocks, 10 threads per block convolve<<<100, 10>>> (myimage);

CUDA程序的编译 使用nvcc编译工具 nvcc <filename>.cu [-o excutable] 调试选项：-g（debug）、-deviceemu（CPU模拟GPU）

并行线程组织 并行性的维度 一维 y = a + b 二维 P = M  N 三维 CT or MRI

并行线程组织结构 Thread: 并行的基本单位 Thread block: 互相合作的线程组 Grid: 一组thread block Cooperative Thread Array (CTA) 允许彼此同步 通过快速共享内存交换数据 以1维、2维或3维组织 最多包含512个线程 Grid: 一组thread block 共享全局内存 Kernel: 在GPU上执行的核心程序 One kernel ↔ one grid

线程层次

Block and Thread IDs Blocks 和 Threads 具有IDs threadIdx, blockIdx Block ID: 1D or 2D Thread ID: 1D, 2D or 3D 由此决定相应处理数据

CUDA线程组织 … CUDA kernel函数由一系列线程组成 线程可划分为不同的Block 单指令多数据流（SPMD） 通过IDs确定处理的数据 线程可划分为不同的Block 在同一个block中，可以通过share memory、atomic operation和barrier synchronization进行协同 Thread Block 0 Thread Block 1 Thread Block N - 1 … float x = input[threadID]; float y = func(x); output[threadID] = y; threadID 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 … float x = input[threadID]; float y = func(x); output[threadID] = y; … float x = input[threadID]; float y = func(x); output[threadID] = y; …

一个简单的例子——Increment Array Elements //CPU program void inc_cpu(float *a, float b, int N) { for (intidx = 0; idx<N; idx++) a[idx] = a[idx] + b; } void main() … inc_cpu(a, b, N); //CUDA program __global__ void inc_gpu(float *a, float b, int N) { intidx =blockIdx.x* blockDim.x+ threadIdx.x; if (idx < N) a[idx] = a[idx] + b; } void main() … dim3 dimBlock (blocksize); dim3 dimGrid( ceil( N / (float)blocksize) ); inc_gpu<<<dimGrid, dimBlock>>>(a, b, N);

CUDA线程的同步 void __syncthreads(); Barrier synchronization 同步thread block之内的所有线程 避免访问共享内存时发生RAW/WAR/WAW 冒险 __shared__ float scratch[256]; scratch[threadID] = begin[threadID]; __syncthreads(); int left = scratch[threadID -1]; 在此等待，直至所有线程到达才开始执行下面的代码

存储器模型与内存分配 R/W per-thread registers R/W per-thread local memory 1-cycle latency R/W per-thread local memory Slow – register spilling to global memory R/W per-block shared memory But bank conflicts may drag down R/W per-grid global memory ~500-cycle latency But coalescing accessing could hide latency Read only per-grid constant and texture memories But cached

GPU Global Memory分配 cudaMalloc() cudaFree() 分配显存中的global memory 两个参数 对象数组指针 数组尺寸 cudaFree() 释放显存中的global memory 一个参数

GPU Global Memory分配 代码实例 分配64  64单精度浮点数组 数组指针Md 建议用“d”表示GPU显存数据结构 int BLOCK_SIZE = 64; float* Md; int size = BLOCK_SIZE * BLOCK_SIZE * sizeof(float); cudaMalloc((void**)&Md, size); cudaFree(Md);

Host - Device数据交换 cudaMemcpy() 在存储器直接传输数据 四个参数 目的对象数组指针 源对象数组指针 数组尺寸 传输方向 Host到Host Host到Device Device到Host Device到Device

Host - Device数据交换 代码实例 M.elements: CPU主存 Md: GPU显存 符号常数: cudaMemcpyHostToDevice和cudaMemcpyDeviceToHost cudaMemcpy(Md, M.elements, size, cudaMemcpyHostToDevice); cudaMemcpy(M.elements, Md, size, cudaMemcpyDeviceToHost);

CUDA变量与函数 CUDA引入的变量修饰词 __device__ __constant__ __shared__ 无修饰（Local变量） 储存于GPU上的global memory空间 和应用程序具有相同的生命期(lifetime) 可被grid中所有线程存取, CPU代码通过runtime函数存取 __constant__ 储存于GPU上的constant memory空间 __shared__ 储存于GPU上thread block内的共享存储器 和thread block具有相同的生命期(lifetime) 只能被thread block内的线程存取 无修饰（Local变量） 储存于SM内的寄存器和local memory 和具有相同的生命期(lifetime) Thread私有

Built-in dim3 Type 定义grid和thread block的组织 dim3 dimGrid(2, 2); dim3 dimBlock(4, 2, 2); kernelFunction<<< dimGrid, dimBlock>>>(…);

CUDA函数定义 __global__定义kernel函数 __device__和__host__ 可以组合使用 必须返回void Executed on the: Only callable from the: __device__ float DeviceFunc() device __global__ void KernelFunc() host __host__ float HostFunc() __global__定义kernel函数 必须返回void __device__和__host__ 可以组合使用 则被定义的函数在CPU和GPU上都被编译

CUDA函数定义 __device__ 函数不能用&运算符取地址 限制 不支持递归调用 不支持静态变量(static variable) 不支持可变长度参数函数调用 type va_list(stdarg.h) double average(int count, ...)

Kernel函数调用 调用时必须给出线程配置方式 __global__ void KernelFunc(...); dim3 DimGrid(100, 50); // 5000 thread blocks dim3 DimBlock(4, 8, 8); // 256 threads per block size_t SharedMemBytes = 64; // 64 bytes of shared memory KernelFunc<<< DimGrid, DimBlock, SharedMemBytes >>>(...);

CUDA数学函数 pow, sqrt, cbrt, hypot, exp, exp2, expm1, log, log2, log10, log1p, sin, cos, tan, asin, acos, atan, atan2, sinh, cosh,tanh, asinh, acosh, atanh, ceil, floor, trunc, round, etc. 只支持标量运算 许多函数有一个快速、较不精确的对应版本 以“__”为前缀，如__sin() 编译开关-use_fast_math强制生成该版本的目标码 每个多处理器包含两个超越函数计算单元

CUDA程序设计实例——方阵相乘 P = M * N （长宽均为WIDTH） 计算策略 每个线程计算矩阵P中的一个元素 N M P

第一步：CPU实现 // Matrix multiplication on the (CPU) host in double precision void MatrixMulOnHost(float* M, float* N, float* P, int Width) { for (int i = 0; i < Width; ++i) for (int j = 0; j < Width; ++j) { double sum = 0; for (int k = 0; k < Width; ++k) { double a = M[i * width + k]; double b = N[k * width + j]; sum += a * b; } P[i * Width + j] = sum; N k j WIDTH M P i WIDTH k WIDTH WIDTH

第二步：将矩阵数据传给显存 void MatrixMulOnDevice(float* M, float* N, float* P, int Width) { int size = Width * Width * sizeof(float); float* Md, Nd, Pd; … 1. // Allocate and Load M, N to device memory cudaMalloc(&Md, size); cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice); cudaMalloc(&Nd, size); cudaMemcpy(Nd, N, size, cudaMemcpyHostToDevice); // Allocate P on the device cudaMalloc(&Pd, size);

第三步：将计算结果传回内存 2. // Kernel invocation code – to be shown later … 3. // Read P from the device cudaMemcpy(P, Pd, size, cudaMemcpyDeviceToHost); // Free device matrices cudaFree(Md); cudaFree(Nd); cudaFree (Pd); }

第四步：kernel函数 // Matrix multiplication kernel – per thread code __global__ void MatrixMulKernel(float* Md, float* Nd, float* Pd, int Width) { // 2D Thread ID int tx = threadIdx.x; int ty = threadIdx.y; // Pvalue is used to store the element of the matrix // that is computed by the thread float Pvalue = 0;

第四步：kernel函数（续） k tx ty ty tx k for (int k = 0; k < Width; ++k) { float Melement = Md[ty * Width + k]; float Nelement = Nd[k * Width + tx]; Pvalue += Melement * Nelement; } Pd[ty * Width + tx] = Pvalue; Nd k tx WIDTH Md Pd ty ty WIDTH tx k WIDTH WIDTH

第五步：调用kernel函数 2. // Kernel invocation code // Setup the execution configuration dim3 dimBlock(Width, Width); dim3 dimGrid(1, 1); // Launch the device computation threads! MatrixMulKernel<<<dimGrid, dimBlock>>>(Md, Nd, Pd);

局限性 Nd 每个线程都需要读： 矩阵规模受限于每个block允许的thread数目 Md Pd Md矩阵的一行 Nd矩阵的一列 Grid 1 Block 1 48 Thread (2, 2) WIDTH Md Pd Nd 每个线程都需要读： Md矩阵的一行 Nd矩阵的一列 计算与访存比约为1:1 矩阵规模受限于每个block允许的thread数目

参考资料——CUDA SDK SDK中包含许多CUDA范例

多线程及存储器硬件

Streaming Multiprocessor执行Thread Blocks 线程以block为单位分配到SM 视资源需求, 一个SM分配至多8个block SM in G80可以接受768个线程 256 (threads/block) * 3 blocks 128 (threads/block) * 6 blocks, etc 线程并发运行 SM分配并维护线程ID SM管理并调度线程

Thread Block Size Considerations 对于矩阵乘法, 哪个thread block尺寸最好: 8X8, 16X16或者32X32? 8X8: 64 threads/block. 每个SM至多接受768 threads, 即12 blocks。但是, SM至多接受8 blocks, 所以实际上仅有512 threads 16X16: 256 threads/block.每个SM至多接受768 threads, 即3 blocks → 只要其它计算资源许可，可以满负荷工作 32X32: 1024 threads/block. SM无法处理

线程调度和执行 Thread block内部线程组织为32-thread warps Warp是SM调度的基本单位 An implementation decision - not part of CUDA Warp是SM调度的基本单位 Warp就是一条32路SIMD指令 Half-warp是warp的前一半或后一半 访问存储器的基本单位

Warp调度和执行 下条指令中全部操作数就位的warps拥有执行资格 Warp中全部线程执行同一指令 不同指令路径的线程只能顺序执行 每次执行warp中一条可能的路径 N条指令路径→1/N throughput 应尽量避免在同一warp内出现分支

SM存储器资源 Register and local memory: per-thread Shared memory: per-block 线程私有 编译器自行分配 e.g. float a; Shared memory: per-block Block内所有线程共享 使数据尽量靠近处理器 动态分配到blocks e.g. __shared__ float region[M]; Constant cache Texture cache

寄存器阵列 G80中每个SM配置8192个寄存器 实例：假设每个block有1616个thread 当前设计选择，不属于CUDA 寄存器动态分配到划归SM的blocks中 一旦分配到某一block, 不能被其它blocks访问 同一block内部的线程只能使用分配给该线程的寄存器 实例：假设每个block有1616个thread 每个thread使用10个寄存器，那么每个block需要使用2560个寄存器，因此每个SM能容纳3个block，也就是768个thread 假如每个thread多使用1个寄存器，那么每个block需要使用2816个寄存器，SM就只能容纳2个block，造成并行度下降1/3 应该综合考虑并行度与访存开销的影响。假如在上面的情况下多使用1个寄存器能够使访存次数减少一半，那么实际性能反而有所提高

利用Shared Memory提高性能 存储器模型回顾 性能优化思路 R/W per-block shared memory 1-cycle latency But bank conflicts may drag down R/W per-grid global memory ~500-cycle latency But coalescing accessing could hide latency 性能优化思路 Shared Memory比Global Memory快几百倍 线程之间通过Shared Memory合作 使用一个或少量线程装载和计算thread block内全部线程共享的数据

实例：矩阵乘法性能优化 每个元素都需要被多个线程重复使用 将元素存入shared memory供线程共享 分块计算 ty tx N M P WIDTH M P ty WIDTH tx WIDTH WIDTH

矩阵乘法的分块计算 每个block计算一块小方矩阵Pd 每个thread读入Pd的一个元素 假设M和N的大小是小矩阵大小的整数倍 bx tx Md Nd Pd Pdsub TILE_WIDTH WIDTH bx tx 1 TILE_WIDTH-1 2 by ty TILE_WIDTHE 矩阵乘法的分块计算 每个block计算一块小方矩阵Pd 每个thread读入Pd的一个元素 假设M和N的大小是小矩阵大小的整数倍

线程结构配置 // Setup the execution configuration dim3 dimBlock(TILE_WIDTH, TILE_WIDTH); dim3 dimGrid(Width / TILE_WIDTH, Width / TILE_WIDTH);

kernel函数代码——整体结构 // Block index int bx = blockIdx.x; int by = blockIdx.y; // Thread index int tx = threadIdx.x; int ty = threadIdx.y; // Pvalue stores the element of the block sub-matrix // that is computed by the thread – automatic variable! float Pvalue = 0; // Loop over all the sub-matrices of M and N // required to compute the block sub-matrix for (int m = 0; m < Width/TILE_WIDTH; ++m) { code for calculating the sub-sum; };

地址计算 在第m次循环中每个线程需要读取的元素地址为 M矩阵： Md Nd Pd Pdsub TILE_WIDTH WIDTH bx tx 1 TILE_WIDTH-1 2 by ty TILE_WIDTHE m 地址计算 在第m次循环中每个线程需要读取的元素地址为 M矩阵： (by*TILE_WIDTH+ty)*width+m*TILE_WIDTH+tx N矩阵： (m*TILE_WIDTH+ty)*width+bx*TILE_WIDTH+tx P矩阵： (by*TILE_WIDTH+ty)*width+bx*TILE_WIDTH+tx

将数据读入Shared Memory并计算sub-sum __shared__ float Mds[TILE_WIDTH][TILE_WIDTH]; __shared__ float Nds[TILE_WIDTH][TILE_WIDTH]; Mds[ty][tx] = Md[(by*TILE_WIDTH+ty)*width+m*TILE_WIDTH+tx]; Nds[ty][tx] = Nd[(m*TILE_WIDTH+ty)*width+bx*TILE_WIDTH+tx]; __syncthreads(); for (int k = 0; k < TILE_WIDTH; ++k) Pvalue += Mds[ty][k] * Nds[k][tx];

将结果写回P矩阵 Pd[(by*TILE_WIDTH+ty)*width+bx*TILE_WIDTH+tx] = Pvalue; 大约能达到峰值性能的15%

Thank you！