cuda开发进阶

1 线程组织

CUDA线程分成Grid和Block两个层次：

一个单独的kernel启动的所有线程组成一个grid，grid中所有线程共享global memory。一个grid由许多block组成，block由许多线程组成。grid和block都可以是一维、二维或者三维，上图是一个二维的grid和二维的block。相关CUDA内置变量如下：

blockIdx：block索引，blockIdx.x表示block的x坐标
threadIdx：线程索引，同理blockIdx
blockDim：block维度，上图中blockDim.x=5
gridDim：grid维度，同理blockDim

一般会把grid组织成2D,block为3D。grid和block都使用dim3作为声明：

dim3 block{3};
dim3 grid{(nElem+block.x-1)/block.x}

dim3仅为host端可见，其对应的device端类型为uint3。

2 kernel启动

CUDA kernel调用格式为：

kernel_name<<<grid,block>>>(arument list);

其中grid和block即为前文介绍的类型dim3的变量。通过两个变量可以配置一个kernel的线程的总和，以及线程的组织形式。例如：

kernel_name<<<4,8>>>(argument list);

该行代码表明grid为一维，有个block，block为一维，每个block有8个线程，故此共有4*8=32个线程。

所有CUDA kernel的启动都是异步的，当CUDA kernel被调用时，控制权会立即返回给CPU。

函数类型标识符：

__device__和__host__可以组合使用。

kernel的限制：

仅能获取device memeory
必须返回void类型
不支持可变数目参数
不支持静态变量
不支持函数指针
异步

3 线程索引

线程的组织形式对程序的性能影响至关重要。矩阵在memory中是row-major线性存储的：

在kernel中，线程的唯一索引非常有用，为例确定线程的索引，我们以2D为例：

线程和block索引
矩阵中元素的坐标
线性global memory的偏移

首先，可以将thread和block索引映射到矩阵坐标：

ix=threadIdx.x+blockIdx.x*blockDim.x
iy=threadIdx.y+blockIdx.y*blockDim.y

之后可以利用上述变量计算线性地址：

idx=iy*nx+ix

上图展示了block和thread索引，矩阵坐标以及线性地址之间的关系。相邻thread拥有连续的threadIdx.x，也就是索引为(0,0),(1,0),(2,0),(3,0)...的thread连续，而不是(0,0),(0,1),(0,2),(0,3)...连续，这跟线性代数中的矩阵不一样。

我们来验证下面的关系：

thread_id（2,1）block_id（1,0） coordinate（6,1） global index 14 ival 14

如下图：

4 device管理

掌握GPU信息查询信息很重要，这可以帮助设置kernel的执行配置。runtime api查询GPU信息：

cudaError_t cudaGetDeviceProperties(cudaDeviceProp *prop, int device);

如何决定最佳GPU？对于支持多GPU系统，需要从中选择一个来作为我们的device，选择出最佳计算性能GPU的一种方法就是由其拥有的处理器数量决定，可以使用下面的代码来选择最佳GPU。

int numDevices = 0;
cudaGetDeviceCount(&numDevices);
if (numDevices > 1) {
    int maxMultiprocessors = 0, maxDevice = 0;
    for (int device=0; device<numDevices; device++) {
        cudaDeviceProp props;
        cudaGetDeviceProperties(&props, device);
        if (maxMultiprocessors < props.multiProcessorCount) {
            maxMultiprocessors = props.multiProcessorCount;
            maxDevice = device;
        }
    }
    cudaSetDevice(maxDevice);
}

设置device：对于多GPU系统，使用nvidia-smi可以查看各GPU属性，每个GPU从0开始依次标注，使用环境变量CUDA_VISIBLE_DEVICES可以指定GPU而不用修改应用程序。可以设置环境变量CUDA_VISIBLE_DEVICES-2来屏蔽其他GPU，这样只有GPU2能被使用。当然也可以使用CUDA_VISIBLE_DEVICES-2,3来设置多个GPU，他们的device ID分别为0和1。

5 GPU架构

SM（Streaming Multiprocessor）是GPU架构中非常重要的部分，GPU硬件的并行性就是用SM决定的。以Fermi架构为例，其包含以下主要组成部分：

CUDA cores
Shared Memory/L1 Cache
Register File
Load/Store Units
Special Function Units
Warp Scheduler

GPU中，每个SM都设计成支持数以百计的线程并行执行，并且每个GPU都包含了很多的SM。当一个kernel启动后，thread会被分配到这些SM中执行。大量的thread可能会被分配到不同的SM，但是同一个block中的thread必须在同一个SM中并行执行。

CUDA采用Single Instruction Mutiple Thread（SIMT）的架构来管理和执行thread，这些thread以32个为单位组成一个单元，称为warps。warp中所有的线程并行的执行相同的指令。每个thread拥有自己的instruction address counter和状态寄存器，并且用该线程自己的数据执行指令。

SIMT和SIMD（Single Instruction Mutiple Data）类似，SIMT应该算是SIMD的升级版，更灵活，但效率略低，SIMT是NVIDIA提出的GPU新概念。二者都通过将同样的指令广播给多个执行单元来实现并行。一个主要的不同就是，SIMD要求所有的vector element在一个统一的同步组里同步执行，而SIMT运行线程们在一个warp中独立执行。SIMT有三个SIMD没有的主要特征：

每个thread拥有自己的instruction address counter
每个thread拥有自己的状态寄存器
每个thread可以有自己独立的执行路径

一个block只会由一个SM调度，block一旦被分配好SM，该block就会一直驻留在该SM中，直到执行结束。一个SM可以同时拥有多个block。下图显示了软硬件方面的术语：

需要注意的是，大部分thread只是逻辑上并行，并不是所有的thread可以在物理上同时进行。这会导致同一个block中的线程可能会有不同步调。并行thread之间的共享数据会导致竞态；多线程请求同一个数据会导致未定义行为。CUDA提供了API来同步一个block的thread以保证在进行下一步处理前，所有thread都到达某个时间点。不过，我们是没有什么原子操作来保证block之间同步的。

同一个warp中的thread可以任意顺序执行，active warps被SM资源限制。当一个warp空闲时，SM可以调度驻留在该SM中一个可用warp。在并非的warp之间切换是没有什么消耗的，因为硬件资源早就被分配到所有的thread和block，所有该新调度的warp的状态已经存储在SM中了。

SM可以看作GPU的心脏，寄存器和共享内存是SM的稀缺资源。CUDA将这些资源分配给所有驻留在SM中的thread。因此，这些有限的资源就使每个SM中的active warps由非常严格的限制，也就限制了并行能力。所有，掌握部分硬件知识，有助于CUDA性能提升。

5.1 Fermi架构

Fermi是第一个完整的GPU计算架构：

512个accelerator cores 即所谓CUDA cores （包含ALU和FPU）
16个SM，每个SM包含32个CUDA core
6个384位GDDR5 DRAM，支持6GB global on-board memory
GigaThread engine将thread blocks分配给SM调度
768KB L2 cache
每个SM有16个load/store单元，允许每个clock cycle为16个thread(即所谓half-warp,不过现在不提这个概念)计算源地址和目的地址
Special function units（SFU）用来执行sin cosine等
每个SM两个warp scheduler，两个instruction dispatch unit，当一个block被分配到一个SM中后，所有该block中的thread会被分到不同的warp中
Fermi (compute capability 2.x) 每个SM同时可以处理48个warp共计1536个thread。

每个SM由以下几个部分组成：

执行单元（CUDA cores）
调度分配warp的单元
shared memory，register file，L1 cache

5.2 Kepler架构

Kepler相较于Fermi更快，效率更高，性能更好。

15个SM
6个64位memory controller
192个单精度CUDA cores，64个双精度单元，32个SFU，32个load/store单元（LD/ST）
增加register file到64K
每个Kepler的SM包含四个warp scheduler，8个instruction dispatches，使得每个SM可以同时issue和执行4个warp
Kepler K20x (compute capability 3.5) 每个SM可以同时调度64个warp攻击2048个thread

5.2.1 Dynamic Parallelism

Dynamic Parallelism是Kepler的新特性，运行GPU动态启动新的Grid。有了这个特性，任何kernel内都可以启动其它的kernel了。这样直接实现了Kernel的递归已经解决了kernel之间数据依赖问题。

5.2.2 Hyper-Q

Hyper-Q是Kepler的另一个新特征，增加了CPU和GPU之间硬件上的联系，使CPU可以在GPU上同时运行更多任务。这样就可以增加GPU的利用率减少CPU的闲置时间。Fermi依赖一个单独的硬件上的工作队列来从CPU传递任务给GPu，这样在某个任务阻塞时，会导致之后的任务无法得到处理，Hyper-Q解决了这个问题。相应的，Kepler为GPU和CPU提供了32个工作队列。

6 warp

逻辑上，所有thread是并行的，但是，从硬件的角度来说，实际上并不是所有的thread能够同一时刻执行，接下来我们来解释有关warp的一些本质。

6.1 warps 与 blocks

warp是SM的基本执行单元。一个warp包含32个并行的thread，这32个thread执行SIMT模式。也就是说，所有thread执行同一条指令，并且每个thread会使用各自的data执行该指令。

block可以是一维、二维或者三维的，但是从硬件角度看，所有thread都被组织成一维，每个thread都有唯一的ID。

每个block的warp数量可以由下面公式计算获得：

一个warp中的线程必然在同一个block中，如果block所含线程数目不是warp大小的整数倍，那么多出的哪些thread所在warp中，会剩余一些inactive的thread，也就是说，解释凑不够整数倍的thread，硬件也会为warp凑足，只不过哪些thread是inactive状态，需要注意的是，即使这部分thread是inactive，也会消耗SM资源。

6.2 warp divergence

控制流语句普遍存在各种编程中，GPU支持传统的，C-style，显示控制流结构，例如if..else,for,while等等。CPU有复杂的硬件设计可以很好的做分支预测，即预测应用程序走哪个path。如果预测正确，那么CPU只会有很小的消耗。和CPU对比来说，GPU没有那么复杂的分支预测（CPU和GPU差异引起的）。那么，因为所有同一个warp中的thread必须执行相同的指令，如果这些线程遇到控制流语句是，如果进入不同分支，那么同一时刻除了正在执行的分支之外，其余分支都被阻塞了，十分影响性能。这类问题就是warp divergence。

注意，warp divergence问题只会发生在同一个warp中。下图展示了warp divergence问题：

为了获得最好的性能，就需要避免同一个warp存在不同的执行路径。避免该问题的方法很多，比如这样一个情形，假设有两个分支，分支的决定条件是thread的唯一ID的奇偶性：

__global__ void mathKernel1(float *c) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;
    if (tid % 2 == 0) {
        a = 100.0f;
    } else {
        b = 200.0f;
    }
    c[tid] = a + b;
}

一种方法是，将条件改为以warp大小为步调，然后取奇偶，如下：

__global__ void mathKernel2(void) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;
    if ((tid / warpSize) % 2 == 0) {
        a = 100.0f;
    } else {
        b = 200.0f;
    }
    c[tid] = a + b;
}

代码：

int main(int argc, char **argv) {
// set up device
int dev = 0;
cudaDeviceProp deviceProp;
cudaGetDeviceProperties(&deviceProp, dev);
printf("%s using Device %d: %s\n", argv[0],dev, deviceProp.name);
// set up data size
int size = 64;
int blocksize = 64;
if(argc > 1) blocksize = atoi(argv[1]);
if(argc > 2) size = atoi(argv[2]);
printf("Data size %d ", size);
// set up execution configuration
dim3 block (blocksize,1);
dim3 grid ((size+block.x-1)/block.x,1);
printf("Execution Configure (block %d grid %d)\n",block.x, grid.x);
// allocate gpu memory
float *d_C;
size_t nBytes = size * sizeof(float);
cudaMalloc((float**)&d_C, nBytes);
// run a warmup kernel to remove overhead
size_t iStart,iElaps;
cudaDeviceSynchronize();
iStart = seconds();
warmingup<<<grid, block>>> (d_C);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
printf("warmup <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x, iElaps );
// run kernel 1
iStart = seconds();
mathKernel1<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
printf("mathKernel1 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps );
// run kernel 3
iStart = seconds();
mathKernel2<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel2 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps );
// run kernel 3
iStart = seconds ();
mathKernel3<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel3 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps);
// run kernel 4
iStart = seconds ();
mathKernel4<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel4 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps);
// free gpu memory and reset divece
cudaFree(d_C);
cudaDeviceReset();
return EXIT_SUCCESS;
}

输出：

$ ./simpleDivergence using Device 0: Tesla M2070
Data size 64 Execution Configuration (block 64 grid 1)
Warmingup elapsed 0.000040 sec
mathKernel1 elapsed 0.000016 sec
mathKernel2 elapsed 0.000014 sec

我们也可以直接使用nvprof这个工具来度量性能：

$ nvprof --metrics branch_efficiency ./simpleDivergence

输出为：

Kernel: mathKernel1(void)
1 branch_efficiency Branch Efficiency 100.00% 100.00% 100.00%
Kernel: mathKernel2(void)
1 branch_efficiency Branch Efficiency 100.00% 100.00% 100.00%

Branch Efficiency的定义如下：

到这里应该奇怪为什么二者表现相同，实际上当我们代码很简单，可以被预测时，CUDA编译器会自动帮助我们优化代码。稍微提一下GPU分支预测，这里，一个被称为预测变量的东西会被设置成1或者0，所有分支都会得到执行，但是只有预测值为1时才会得到执行。当条件状态少于某一个阈值时，编译器会将分支指令替换为预测指令，因此，现在回到自动优化问题，一份较长的代码就会导致warp divergence了。

6.3 Resource Partitioning

一个warp的context包括以下三部分：

Program counter
Register
Shared memory

在同一个context中切换是没有消耗的，因为在整个warp的生命周期内，SM处理的每个warp执行的contex都是on-chip的。每个SM有一个32位register集合放在register file中，还有固定数量的shared memory，这些资源都被thread瓜分了。由于资源有限，所有，如果thread比较多，那么每个thread占用的资源就较少；thread少，占用资源就较多，这个需要根据自己的要求做出一个平衡。

资源限制了驻留在SM中block的数量，不同device，register和shared memory的数量也不同，比如前文介绍的Fermi和Kepler的差别。如果没有足够的资源，kernel的启动就会失败。

当一个block获得足够的资源时，就成为active block。block中的warp就称为active warp。active warp又可以分为下面三类：

Selected warp
Stalled warp
Eligible warp

SM中warp调度器每个cycle会挑选active warp送去执行，一个被选中的warp被称为selected warp，没被选中但已经做好准备被执行的称为Eligible warp，没准备好要执行的称为Stalled warp。warp适合执行需要满足下面两个条件：

32个CUDA core有空
所有当前指令的参数都准备就绪

例如：kepler任何时刻的active warp数目必须小于或等于64个。selected warp 数目必须小于或等于4个。如果一个warp阻塞了，调度器会挑选一个Eligible warp准备去执行。

CUDA编程中应该重视对计算资源的分配：这些资源限制了active warp的数量。因此，我们必须掌握硬件的一些限制，为了最大化GPU利用率，我们必须最大化active warp的数目。

6.4 Latency Hiding

指令从开始到结束消耗的clock cycle称为指令的latency。当每个cycle都有eligible warp被调度是，计算资源就会得到充分利用，基于此，我们可以将每个指令的latency隐藏于其它warp的指令过程中。

当涉及到指令latency时，指令可以被区分为下面两种：

Arithmetic instruction
Memory instruction

顾明思义，Arithmetic instruction latency是一个算数操作的始末间隔。另一个则是指load或store的始末间隔。二者的latency大约为：

10~20 cycles for arithmetic operations
400~800 cycles for global memory accesses

下图是一个简单的执行流程，当warp0阻塞时，执行其它的warp，当warp变为eligible时重新执行。

如何评估active warps的数量来hide latency？Little’s Law可以提供一个合理的估计：

对于Arithmetic operations来说，并行性可以表示为用来hide Arithmetic latency的操作的数目。下表显示了Fermi和Kepler相关数据，这里以(a+b*c)作为操作的例子。不同的算数指令，throughput(吞吐)也是不同的。

这里的throughput定义为每个SM每个cycle的操作数目。由于每个warp执行同一种指令，因此每个warp对应32个操作。所以，对于Feimi来说，每个SM需要640/32=20个warp来保持计算资源的充分利用。这也意味着，arithmetic operations的并行可以表达为操作的数目或warp的数目。二者的关系也对应了两种方式来增加并行性：

Instruction-level Parallelism(ILP)：同一个thread中更多的独立指令
Thread-level Parallelism(TLP)：更多并发的eligible threads

对于Memory operations，并行性可以表达为每个cylce的byte数目。

因为memory throughput总是以GB/Sec为单位，我们需要先作相应的转化。可以通过以下指令来查看device的memory frequency：

$ nvidia-smi -a -q -d CLOCK | fgrep -A 3 "Max Clocks" | fgrep "Memory"

以Fermi为例，其memory frequency可能是1.566GHz，Kepler的是1.6GHz。那么转换过程为

144GB/sec÷1.566GHz=92Bytes/Cycle

乘上这个92就得到上图的74，这里的数字是针对device的，而不是每个SM。

有了这些数据，我们就可以做一些计算了。以Fermi为例，加上每个thread的任务是将一个float(4bytes)类型的数据从global memory移至SM用来计算，你应该需要大约18500个thread，也就是597个warp来隐藏所有的memory latency。

74KB÷4 bytes/thread = 18500 threads

18500 threads ÷ 32 threads/warp = 579 warps

Fermi有16个SM，所以每个SM需要579÷16=36个warp来隐藏memory latency。

6.5 Occupancy

当一个warp阻塞了，SM会执行另一个eligible warp。理想情况是，每时每刻到保证cores被占用。Occupacy就是每个SM的active warp 占最大warp数目的比例：

我们可以使用cudaError_t cudaGetDeviceProperties(struct cudaDeviceProp *prop, int device)函数来获取warp最大数目。然后用maxThreadsPerMultiProcessor来获取具体数值。

grid和block的配置准则：

保证block中thread数目是32的倍数。
避免block太小：每个block最少128或256个thread
根据kernel需要的资源调整block。
保证block数目远大于SM的数目。
多做实验来挖掘出最好的配置。

Occupacy专注每个SM中可以并行的thread或者warp数目。不管怎样，Occupacy不是唯一的性能指标，当occupacy达到某个值时，再做优化就不可能有效果了。其实还有许多其它指标需要调节，将在后文继续探讨。

6.6 Synchronize

同步是并行编程的一个普遍的问题。在CUDA的世界里，有两种方式实现同步：

System-level：等待所有host和device的工作完成
Block-level：等待device中的block中所有thread执行到某个点

因为CUDA API和host代码是异步的，cudaDeviceSynchronize可以用来停住CPU等待CUDA中的操作完成。因为block中的thread执行顺序补丁，CUDA提供了一个function来同步block中的thread：__device__ void __syncthreads(void);。此函数被调用，block中的每个thread都会等待所有其它thread执行到某个点来实现同步。

7 Kernel性能调节

7.1 代码准备

kernel函数：

__global__ void sumMatrixOnGPU2D(float *A, float *B, float *C, int NX, int NY) {
    unsigned int ix = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int iy = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned int idx = iy * NX + ix;
    if (ix < NX && iy < NY) {
        C[idx] = A[idx] + B[idx];
    }
}

指定一个比较大的数据矩阵，包含16384个元素：

int nx = 1<<14;
int ny = 1<<14;

通过main函数配置block的维度

if (argc > 2) {
    dimx = atoi(argv[1]);
    dimy = atoi(argv[2]);
}
dim3 block(dimx, dimy);
dim3 grid((nx + block.x - 1) / block.x, (ny + block.y - 1) / block.y);

编译

 nvcc -O3 -arch=sm_20 sumMatrix.cu -o sumMatrix

7.2 使用nvprof检查warps的活动

在做各项数据比较的时候需要有个基准，这里使用四个block配置的时间消耗作为基准观察，分别为（32，32）（32,16）（16,32）和（16,16）。下面是几种配置的时间消耗输出结果：

$ ./sumMatrix 32 32
sumMatrixOnGPU2D <<< (512,512), (32,32) >>> elapsed 60 ms
$ ./sumMatrix 32 16
sumMatrixOnGPU2D <<< (512,1024), (32,16) >>> elapsed 38 ms
$ ./sumMatrix 16 32
sumMatrixOnGPU2D <<< (1024,512), (16,32) >>> elapsed 51 ms
$ ./sumMatrix 16 16
sumMatrixOnGPU2D <<< (1024,1024),(16,16) >>> elapsed 46 ms

比较这几个结果，不难发现，最慢的是第一个（32，32），最快的是第二个（32,16），这里可以猜测到的是，拥有更多的block并行性更好。这个猜测可以使用nvprof的ahieved_occupancy这个metirc参数来验证。该参数前文已经介绍过。实际上就是每个SM在每个cycle能够达到的最大active warp数目占总数的比例。下面是使用该参数后得到的结果：

$ nvprof --metrics achieved_occupancy ./sumMatrix 32 32
sumMatrixOnGPU2D <<<(512,512), (32,32)>>> Achieved Occupancy 0.501071
$ nvprof --metrics achieved_occupancy ./sumMatrix 32 16
sumMatrixOnGPU2D <<<(512,1024), (32,16)>>> Achieved Occupancy 0.736900
$ nvprof --metrics achieved_occupancy ./sumMatrix 16 32
sumMatrixOnGPU2D <<<(1024,512), (16,32)>>> Achieved Occupancy 0.766037
$ nvprof --metrics achieved_occupancy ./sumMatrix 16 16
sumMatrixOnGPU2D <<<(1024,1024),(16,16)>>> Achieved Occupancy 0.810691

从上面的输出可以得知两件事：

由于第二个配置比第一个配置有更多的block，device就会达到更多active warp。也就是第二个性能优于第一个的原因。
第四个的achieved Occupacy最高，但不是最快的。因此，较高的achieved Occupacy并不一定就意味着更好的性能，也就是说还有更多的因素影响GPU的性能。

7.3 用nvprof检查内存操作

对于C[idx] = A[idx] + B[idx]来说共有三个memory操作：两个memory load和一个memory store。要查看这些操作的效率可以使用nvprof的两个metric参数，如果想要查看memory的throughput，则可使用gld_throughput：

$ nvprof --metrics gld_throughput./sumMatrix 32 32
sumMatrixOnGPU2D <<<(512,512), (32,32)>>> Global Load Throughput 35.908GB/s
$ nvprof --metrics gld_throughput./sumMatrix 32 16
sumMatrixOnGPU2D <<<(512,1024), (32,16)>>> Global Load Throughput 56.478GB/s
$ nvprof --metrics gld_throughput./sumMatrix 16 32
sumMatrixOnGPU2D <<<(1024,512), (16,32)>>> Global Load Throughput 85.195GB/s
$ nvprof --metrics gld_throughput./sumMatrix 16 16
sumMatrixOnGPU2D <<<(1024,1024),(16,16)>>> Global Load Throughput 94.708GB/s

不难看出，第四个拥有最高的load throughput，但是却比第二个慢（第二个也就是第四个的一半），所以，较高的load throughput也不一定就有较高的性能。之后讲到memory transaction时会具体分析这种现象的原因。简单说，就是高load throughput有可能是一种假象，如果需要的数据在memory中存储格式未对齐、不连续，则会导致不必要的load操作，所以本文的efficiency会这么低。

然后，我们可以使用nvprof的gld_efficiency来度量load efficiency，该metric参数是指我们确切需要的global load throughput与实际得到的global load memory的比值。这个metric参数就可以让我们知道，APP的load操作利用device memory bandwidth的程度：

$ nvprof --metrics gld_efficiency ./sumMatrix 32 32
sumMatrixOnGPU2D <<<(512,512), (32,32)>>> Global Memory Load Efficiency 100.00%
$ nvprof --metrics gld_efficiency ./sumMatrix 32 16
sumMatrixOnGPU2D <<<(512,1024), (32,16)>>> Global Memory Load Efficiency 100.00%
$ nvprof --metrics gld_efficiency ./sumMatrix 16 32
sumMatrixOnGPU2D <<<(1024,512), (16,32)>>> Global Memory Load Efficiency 49.96%
$ nvprof --metrics gld_efficiency ./sumMatrix 16 16
sumMatrixOnGPU2D <<<(1024,1024),(16,16)>>> Global Memory Load Efficiency 49.80%

从上述结果可知，最好两个的load efficiency只是前两个的一半。这也可以解释，为什么较高的throughput和较高的Occupacy却没有产生较好的性能。尽管最后两个的load操作数目要多不少（因为二者的throughput较高），但是他们的load effecitiveness却低不少（由于efficiency较低）

观察最后两个可以发现，他们的block的x象限配置是warp的一般，由前文推测可知，该象限应该保持未warp大小的整数倍。关于具体原因将在后续文中详细解释。

7.4 其它并行性

我们现在可以得出一个结论就是blockDim.x应该是warp大小的整数倍，这样做是很容易提升load efficiency。现在，我们还有其它疑惑，比如：

继续调整blockDim.x是否会继续增加load throughput？
还有其他方法能增大并行性吗？

现在，我们重新调整一个基准数据出来，这两个问题可以从这个基准分析个大概：

$ ./sumMatrix 64 2
sumMatrixOnGPU2D <<<(256,8192), (64,2) >>> elapsed 0.033567 sec
$ ./sumMatrix 64 4
sumMatrixOnGPU2D <<<(256,4096), (64,4) >>> elapsed 0.034908 sec
$ ./sumMatrix 64 8
sumMatrixOnGPU2D <<<(256,2048), (64,8) >>> elapsed 0.036651 sec
$ ./sumMatrix 128 2
sumMatrixOnGPU2D <<<(128,8192), (128,2)>>> elapsed 0.032688 sec
$ ./sumMatrix 128 4
sumMatrixOnGPU2D <<<(128,4096), (128,4)>>> elapsed 0.034786 sec
$ ./sumMatrix 128 8
sumMatrixOnGPU2D <<<(128,2048), (128,8)>>> elapsed 0.046157 sec
$ ./sumMatrix 256 2
sumMatrixOnGPU2D <<<(64,8192), (256,2)>>> elapsed 0.032793 sec
$ ./sumMatrix 256 4
sumMatrixOnGPU2D <<<(64,4096), (256,4)>>> elapsed 0.038092 sec
$ ./sumMatrix 256 8
sumMatrixOnGPU2D <<<(64,2048), (256,8)>>> elapsed 0.000173 sec
Error: sumMatrix.cu:163, code:9, reason: invalid configuration argument

从上面的数据，我们可以分析出来的是：

最后一个配置(256,8)不可行，block中总共的thread数目超过了1024，这是GPU的硬件限制。
最好的结果是第四个(128,2)。
第一个启动了最多block，但不是最快的。
因为第二个与第四个在一个block中拥有相同数目的thread，本应该猜测二者拥有相同的表现，但是实际却是第二个略逊色，所以blockDim.x的大小是关键。
剩下的相对第四个都拥有较少的block数目，所以并行规模也是影响性能的关键因素。

现在，我们推测，拥有block最少的应该会有一个最低的achieved Occupacy，而拥有最多的block应该会到达最高的achieved Occupacy：

$ nvprof --metrics achieved_occupancy ./sumMatrix 64 2
sumMatrixOnGPU2D <<<(256,8192), (64,2) >>> Achieved Occupancy 0.554556
$ nvprof --metrics achieved_occupancy ./sumMatrix 64 4
sumMatrixOnGPU2D <<<(256,4096), (64,4) >>> Achieved Occupancy 0.798622
$ nvprof --metrics achieved_occupancy ./sumMatrix 64 8
sumMatrixOnGPU2D <<<(256,2048), (64,8) >>> Achieved Occupancy 0.753532
$ nvprof --metrics achieved_occupancy ./sumMatrix 128 2
sumMatrixOnGPU2D <<<(128,8192), (128,2)>>> Achieved Occupancy 0.802598
$ nvprof --metrics achieved_occupancy ./sumMatrix 128 4
sumMatrixOnGPU2D <<<(128,4096), (128,4)>>> Achieved Occupancy 0.746367
$ nvprof --metrics achieved_occupancy ./sumMatrix 128 8
sumMatrixOnGPU2D <<<(128,2048), (128,8)>>> Achieved Occupancy 0.573449
$ nvprof --metrics achieved_occupancy ./sumMatrix 256 2
sumMatrixOnGPU2D <<<(64,8192), (256,2) >>> Achieved Occupancy 0.760901
$ nvprof --metrics achieved_occupancy ./sumMatrix 256 4
sumMatrixOnGPU2D <<<(64,4096), (256,4) >>> Achieved Occupancy 0.595197

实际上（64,2）的Achieved Occupacy竟然是最低，经过他有最多的block，它达到了硬件对block数量的限制。

第四个（128,2）和第七个（256,2）拥有差不多的Achieved Occupacy。我们对这两个再做一个试验，再次增大，将blockDim.y设置未了，也就减少了block的大小：

$ ./sumMatrix 128 1
sumMatrixOnGPU2D <<<(128,16384),(128,1)>>> elapsed 0.032602 sec
$ ./sumMatrix 256 1
sumMatrixOnGPU2D <<<(64,16384), (256,1)>>> elapsed 0.030959 sec

这次配置产生了最佳的性能，特别是（256,1）要比（128,1）要更好，再次检查achieved occupacy，load throughput和load efficiency：

$ nvprof --metrics achieved_occupancy ./sumMatrix 256 1
$ nvprof --metrics gld_throughput ./sumMatrix 256 1
$ nvprof --metrics gld_efficiency ./sumMatrix 256 1

输出：

Achieved Occupancy 0.808622
Global Load Throughput 69.762GB/s
Global Memory Load Efficiency 100.00%

可以看出，最佳配置既不是拥有最高的achieved occupacy，也不是最高laod throughput的，所以不存在唯一metric来优化计算性能，而是需要从众多metric中寻求一个平衡。

总结：

在大多数情况下，并不存在唯一的metric可以精确的优化性能。
哪个metric或者event对性能影响大是由kernel具体代码决定的。
在众多的metric和event中寻求一个平台。
Grid/block heuristics（启发）为调节性能提供了不错的切入点。

8 分支分歧与循环展开

8.1 避免分支分歧

有时，控制依赖于thread索引。同一个warp中，一个条件分支可能告知很差的性能。通过重新组织数据获取模式可以减少或避免warp divergence。

8.2 并行规约问题

我们现在计算一个数组N个元素的和。这个过程用CPU编程很容易实现：

int sum = 0;
for (int i = 0; i < N; i++)
    sum += array[i];

那么，如果array的元素非常多呢？应用并行计算可以大大提升这个过程的效率。鉴于甲方的交换律等性质。这个求和过程可以以元素的任意顺序来进行：

将输入数组切割成很多小的块。
用thread来计算每个块的和。
对这些块的结果再求和得最终结果。

数组的切割主旨是，用thread求数组中按一定规律配对的两个元素和，然后将所有结果组合成一个新的数组，然后再次求配对两元素和，多次迭代，直到数组中只有一个结果。

比较直观的两种实现方式：

Neighbored pair：每次迭代都是相邻两个元素求和。
Interleaved pair：按一定跨度配对两个元素。

下图展示了两种方法的求解过程，对于有N个元素的数组，这个过程需要N-1此求和,Log(N)步。Interleaved pair的跨度是半个数组长度。

下面是用递归实现的Interleaved pair代码：

int recursiveReduce(int *data, int const size) {
    // terminate check
    if (size == 1) return data[0];
    // renew the stride
    int const stride = size / 2;
    // in-place reduction
    for (int i = 0; i < stride; i++) {
        data[i] += data[i + stride];
    }
    // call recursively
    return recursiveReduce(data, stride);
}

上述讲的这类问题术语叫规约问题。Parallel reduction（并行规约）是指迭代减少操作，是并行算法中非常关键的一种操作。

8.3 并行规约中的发散

这部分以neighbored pair 为参考研究：

在这kernel里面，有两个global memory array，一个用来存放数组所有数据，另一个用来存放部分和。所有block独立的执行求和操作。__syncthreads用来保证每次迭代，所有的求和操作都做完，然后进入下一步迭代。

__global__ void reduceNeighbored(int *g_idata, int *g_odata, unsigned int n) {
    // set thread ID
    unsigned int tid = threadIdx.x;
    // convert global data pointer to the local pointer of this block
    int *idata = g_idata + blockIdx.x * blockDim.x;
    // boundary check
    if (idx >= n) return;
        // in-place reduction in global memory
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        if ((tid % (2 * stride)) == 0) {
            idata[tid] += idata[tid + stride];
        }
        // synchronize within block
        __syncthreads();
    }
    // write result for this block to global mem
    if (tid == 0) g_odata[blockIdx.x] = idata[0];
}

因为没有办法让所有的block同步，所以最后将所有block的结果送回host来进行串行计算，如下图所示：

main代码：

int main(int argc, char **argv) {
// set up device
int dev = 0;
cudaDeviceProp deviceProp;
cudaGetDeviceProperties(&deviceProp, dev);
printf("%s starting reduction at ", argv[0]);
printf("device %d: %s ", dev, deviceProp.name);
cudaSetDevice(dev);
bool bResult = false;
// initialization
int size = 1<<24; // total number of elements to reduce
printf(" with array size %d ", size);
// execution configuration
int blocksize = 512; // initial block size
if(argc > 1) {
blocksize = atoi(argv[1]); // block size from command line argument
}
dim3 block (blocksize,1);
dim3 grid ((size+block.x-1)/block.x,1);
printf("grid %d block %d\n",grid.x, block.x);
// allocate host memory
size_t bytes = size * sizeof(int);
int *h_idata = (int *) malloc(bytes);
int *h_odata = (int *) malloc(grid.x*sizeof(int));
int *tmp = (int *) malloc(bytes);
// initialize the array
for (int i = 0; i < size; i++) {
// mask off high 2 bytes to force max number to 255
h_idata[i] = (int)(rand() & 0xFF);
}
memcpy (tmp, h_idata, bytes);
size_t iStart,iElaps;
int gpu_sum = 0;
// allocate device memory
int *d_idata = NULL;
int *d_odata = NULL;
cudaMalloc((void **) &d_idata, bytes);
cudaMalloc((void **) &d_odata, grid.x*sizeof(int));
// cpu reduction
iStart = seconds ();
int cpu_sum = recursiveReduce(tmp, size);
iElaps = seconds () - iStart;
printf("cpu reduce elapsed %d ms cpu_sum: %d\n",iElaps,cpu_sum);
// kernel 1: reduceNeighbored
cudaMemcpy(d_idata, h_idata, bytes, cudaMemcpyHostToDevice);
cudaDeviceSynchronize();
iStart = seconds ();
warmup<<<grid, block>>>(d_idata, d_odata, size);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
cudaMemcpy(h_odata, d_odata, grid.x*sizeof(int), cudaMemcpyDeviceToHost);
gpu_sum = 0;
for (int i=0; i<grid.x; i++) gpu_sum += h_odata[i];
printf("gpu Warmup elapsed %d ms gpu_sum: %d <<<grid %d block %d>>>\n",
iElaps,gpu_sum,grid.x,block.x);
// kernel 1: reduceNeighbored
cudaMemcpy(d_idata, h_idata, bytes, cudaMemcpyHostToDevice);
cudaDeviceSynchronize();
iStart = seconds ();
reduceNeighbored<<<grid, block>>>(d_idata, d_odata, size);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
cudaMemcpy(h_odata, d_odata, grid.x*sizeof(int), cudaMemcpyDeviceToHost);
gpu_sum = 0;
for (int i=0; i<grid.x; i++) gpu_sum += h_odata[i];
printf("gpu Neighbored elapsed %d ms gpu_sum: %d <<<grid %d block %d>>>\n",
iElaps,gpu_sum,grid.x,block.x);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
cudaMemcpy(h_odata, d_odata, grid.x/8*sizeof(int), cudaMemcpyDeviceToHost);
gpu_sum = 0;
for (int i = 0; i < grid.x / 8; i++) gpu_sum += h_odata[i];
printf("gpu Cmptnroll elapsed %d ms gpu_sum: %d <<<grid %d block %d>>>\n",
iElaps,gpu_sum,grid.x/8,block.x);
/// free host memory
free(h_idata);
free(h_odata);
// free device memory
cudaFree(d_idata);
cudaFree(d_odata);
// reset device
cudaDeviceReset();
// check the results
bResult = (gpu_sum == cpu_sum);
if(!bResult) printf("Test failed!\n");
return EXIT_SUCCESS;
}

初始化数组，使其包含16M元素：

int size = 1<<24;

kernel配置为1D grid和1D block：

dim3 block (blocksize, 1);
dim3 block ((siize + block.x – 1) / block.x, 1);

编译：

$ nvcc -O3 -arch=sm_20 reduceInteger.cu -o reduceInteger

运行：

$ ./reduceInteger starting reduction at device 0: Tesla M2070
with array size 16777216 grid 32768 block 512
cpu reduce elapsed 29 ms cpu_sum: 2139353471
gpu Neighbored elapsed 11 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
Improving Divergence in Parallel Reduction

考虑上节if判断条件：

if ((tid % (2 * stride)) == 0)

因为这表达式只对偶数ID的线程为true，所以其导致很高的divergent warps。第一次迭代只有偶数ID的线程执行了指令，但是所有线程都要被调度；第二次迭代，只有四分之一的thread是active的，但是所有的thread仍然要被调度。

我们可以重新组织每个线程对应的数组索引来强制id相邻的thread来处理求和操作。如下图所示（主要图中的thread id与上一个图的区别）：

新的代码：

__global__ void reduceNeighboredLess (int *g_idata, int *g_odata, unsigned int n) {
    // set thread ID
    unsigned int tid = threadIdx.x;
    unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
    // convert global data pointer to the local pointer of this block
    int *idata = g_idata + blockIdx.x*blockDim.x;
    // boundary check
    if(idx >= n) return;
    // in-place reduction in global memory
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        // convert tid into local array index
        int index = 2 * stride * tid;
        if (index < blockDim.x) {
            idata[index] += idata[index + stride];
        }    
        // synchronize within threadblock
        __syncthreads();
    }
    // write result for this block to global mem
    if (tid == 0) g_odata[blockIdx.x] = idata[0];
}

注意这行代码：

int index = 2 * stride * tid;

因为步调乘以了2，下面的语句使用block的前半部分thread来执行求和：

if (index < blockDim.x)

对于一个有512个thread的block来说，前八个warp执行第一轮reduction，剩下八个warp什么也不干；第二轮，前4个warp执行，剩下12个warp什么也不干。因此，就彻底不存在divegence了（重申：divergence只发生于同一个warp）。最后的五轮还是会导致divergence，因为这个时候需要执行thread已经凑不够一个warp了。

// kernel 2: reduceNeighbored with less divergence
cudaMemcpy(d_idata, h_idata, bytes, cudaMemcpyHostToDevice);
cudaDeviceSynchronize();
iStart = seconds();
reduceNeighboredLess<<<grid, block>>>(d_idata, d_odata, size);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
cudaMemcpy(h_odata, d_odata, grid.x*sizeof(int), cudaMemcpyDeviceToHost);
gpu_sum = 0;
for (int i=0; i<grid.x; i++) gpu_sum += h_odata[i];
printf("gpu Neighbored2 elapsed %d ms gpu_sum: %d <<<grid %d block %d>>>\n",iElaps,gpu_sum,grid.x,block.x);

运行结果：

$ ./reduceInteger Starting reduction at device 0: Tesla M2070
vector size 16777216 grid 32768 block 512
cpu reduce elapsed 0.029138 sec cpu_sum: 2139353471
gpu Neighbored elapsed 0.011722 sec gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu NeighboredL elapsed 0.009321 sec gpu_sum: 2139353471 <<<grid 32768 block 512>>>

新的实现比原来的快了1.26。我们也可以使用nvprof的inst_per_warp参数来查看每个warp上执行的指令数目的平均值。

$ nvprof --metrics inst_per_warp ./reduceInteger

输出，原来的是新kernel的两倍还多，因为原来的有许多不必要的操作也执行了：

Neighbored Instructions per warp 295.562500
NeighboredLess Instructions per warp 115.312500

在查看throughput：

$ nvprof --metrics gld_throughput ./reduceInteger

输出，新的kernel拥有更大的throughput，因为虽然I/O操作数目相同，但是其消耗时间短：

Neighbored Global Load Throughput 67.663GB/s
NeighboredL Global Load Throughput 80.144GB/s
Reducing with Interleaved Pairs

Interleaved pair模式的初始步调是block大小的一半，每个thread处理像半个block的两个数据之和。和之前的图示相比，工作的thread数目没有变化。但是每个thread的load/store global memory的位置是不同的。

Interleaved pair的kernel实现：

/// Interleaved Pair Implementation with less divergence
__global__ void reduceInterleaved (int *g_idata, int *g_odata, unsigned int n) {
// set thread ID
unsigned int tid = threadIdx.x;
unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
// convert global data pointer to the local pointer of this block
int *idata = g_idata + blockIdx.x * blockDim.x;
// boundary check
if(idx >= n) return;
// in-place reduction in global memory
for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
if (tid < stride) {
idata[tid] += idata[tid + stride];
}
__syncthreads();
}
// write result for this block to global mem
if (tid == 0) g_odata[blockIdx.x] = idata[0];
}

注意下面的语句，步调被初始化为block大小的一半：

for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {

下面的语句使得第一次迭代时，block的前半部分thread执行相加操作，第二次是前四分之一，以此类推：

if (tid < stride)

下面是加入main的代码：

cudaMemcpy(d_idata, h_idata, bytes, cudaMemcpyHostToDevice);
cudaDeviceSynchronize();
iStart = seconds();
reduceInterleaved <<< grid, block >>> (d_idata, d_odata, size);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
cudaMemcpy(h_odata, d_odata, grid.x*sizeof(int), cudaMemcpyDeviceToHost);
gpu_sum = 0;
for (int i = 0; i < grid.x; i++) gpu_sum += h_odata[i];
printf("gpu Interleaved elapsed %f sec gpu_sum: %d <<<grid %d block %d>>>\n",iElaps,gpu_sum,grid.x,block.x);

输出：

$ ./reduce starting reduction at device 0: Tesla M2070
with array size 16777216 grid 32768 block 512
cpu reduce elapsed 0.029138 sec cpu_sum: 2139353471
gpu Warmup elapsed 0.011745 sec gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu Neighbored elapsed 0.011722 sec gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu NeighboredL elapsed 0.009321 sec gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu Interleaved elapsed 0.006967 sec gpu_sum: 2139353471 <<<grid 32768 block 512>>>

这次相对第一个kernel又快了1.69，比第二个也快了1.34.这个效果主要由global memory的load/store模式导致的。

8.4 循环展开

loop unrolling 是用来优化循环，减少分支的方法，该方法简单说就是把本应该在多次loop中完成的操作，尽量压缩到一次loop。循环体展开的程度称为loop unrolling factor（循环展开因子），loop unrolling对于顺序数组操作性能影响很大，考虑如下代码：

for (int i = 0; i < 100; i++) {
    a[i] = b[i] + c[i];
}

如下重复一次循环体操作，迭代数目将减少一半：

for (int i = 0; i < 100; i += 2) {
    a[i] = b[i] + c[i];
    a[i+1] = b[i+1] + c[i+1];
}

从高级语言层面是无法看出性能提升的原因的，需要从low-level instruction层面去分析，第二阶段循环次数减少了一半，而循环体两句语句的读写操作的执行在CPU上市可以同时执行、互相独立的，所以相对第一段，第二段性能要更好。

Unrolling在CUDA编程中意义更重。我们的目标依然是通过减少指令执行消耗，增加更多的独立指令来提高性能。这样就会增加更多的并行操作从而产生更高的智力和内地带宽（bandwith）。也就提供了更多的eligible warps来帮助hide instruction/memory latency。

8.5 规约展开

在前文的reduce Interleaved中，每个block处理一部分数据，我们给这数据起名data block。下面的代码是reduceInterleaved的修正版本，每个block都是以两个data block作为源数据进行操作（前文中，每个block处理一个data block）。这是一种cyclic partitioning：每个thread作用于多个data block，并且从每个data block中取出一个元素处理。

__global__ void reduceUnrolling2 (int *g_idata, int *g_odata, unsigned int n) {
    // set thread ID
    unsigned int tid = threadIdx.x;
    unsigned int idx = blockIdx.x * blockDim.x * 2 + threadIdx.x;

    // convert global data pointer to the local pointer of this block
    int *idata = g_idata + blockIdx.x * blockDim.x * 2;

    // unrolling 2 data blocks
    if (idx + blockDim.x < n) g_idata[idx] += g_idata[idx + blockDim.x];
    __syncthreads();

    // in-place reduction in global memory
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            idata[tid] += idata[tid + stride];
        }
        // synchronize within threadblock
        __syncthreads();
    }

    // write result for this block to global mem
    if (tid == 0) g_odata[blockIdx.x] = idata[0];
}

注意下面的语句，每个thread从相邻的data block中取数据，这一步实际上就是将两个data block规约成一个。

if (idx + blockDim.x < n) g_idata[idx] += g_idata[idx+blockDim.x];

global array index也要相应的调整。因为，相对之前的版本，同样的数据，我们只需要原来一半的thread就能解决问题。要注意的是，这样做也会降低warp或block的并行性（因为thread少了）：

main增加代码如下：

cudaMemcpy(d_idata, h_idata, bytes, cudaMemcpyHostToDevice);
cudaDeviceSynchronize();
iStart = seconds();
reduceUnrolling2 <<< grid.x/2, block >>> (d_idata, d_odata, size);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
cudaMemcpy(h_odata, d_odata, grid.x/2*sizeof(int), cudaMemcpyDeviceToHost);
gpu_sum = 0;
for (int i = 0; i < grid.x / 2; i++) gpu_sum += h_odata[i];
printf("gpu Unrolling2 elapsed %f sec gpu_sum: %d <<<grid %d block %d>>>\n",iElaps,gpu_sum,grid.x/2,block.x);

posted @ 2026-01-11 09:43 钟齐峰阅读(10) 评论(0) 收藏举报

刷新页面返回顶部

cuda开发进阶

cuda开发进阶

1 线程组织

2 kernel启动

3 线程索引

4 device管理

5 GPU架构

5.1 Fermi架构

5.2 Kepler架构

5.2.1 Dynamic Parallelism

5.2.2 Hyper-Q

6 warp

6.1 warps 与 blocks

6.2 warp divergence

6.3 Resource Partitioning

6.4 Latency Hiding

6.5 Occupancy

6.6 Synchronize

7 Kernel性能调节

7.1 代码准备

7.2 使用nvprof检查warps的活动

7.3 用nvprof检查内存操作

7.4 其它并行性

8 分支分歧与循环展开

8.1 避免分支分歧

8.2 并行规约问题

8.3 并行规约中的发散

8.4 循环展开

8.5 规约展开

公告