CUDA2.4-原理之性能优化及浮点运算

本部分来自于《大规模并行处理器编程实战》第六章、第七章。打算不再看这本书了，准备看《programming massively parallel processors 2nd》，即它的第二版，第一版是09年的，第二版是13年的，虽说第二版可是里面涉及的是cuda4.0 和5.0，然而现在2015年7月，cuda都7.0了，正所谓赶速度，完全赶不上啊。虽然说本书好，不过一个不小心，你费老大劲做的优化，发现其实新版本的cuda或者硬件完全不需要，果然有关cuda的最好的资料其实还是官方文档，因为这些完全赶不上速度啊。

一、性能优化

        CUDA kernel函数的执行速度很大程度上取决于每个设备的资源约束。而且不同的应用程序中，不同的约束可能决定并成为限制因素。在特定的CUDA设备上，可以通过一种资源代替另一种资源来提高英程序的性能。合理的策略有可能提高性能，也可能不起作用，所以需要测试，本章说的这些可以用来培养程序员对算法的直觉，如何来提高整体的性能。

1.1 线程执行问题

        相比较前面几个博文，没有过多的讨论每个块中线程的执行时间问题。在本书发售的时候，Nvidia公司是通过对块中的线程进行捆绑执行的，即执行的不是单个线程而是一个warp（包含32个线程），这样做降低硬件成本，而且一定程度上优化了存储器访问的服务。对于划分warp来说，如果是一维的，那么对于最后一个warp的划分如果不满32个，会将其他块中的线程拉过来补全成32个再执行。如果对于多维线程的块来说，划分warp前会把维度映射成一个线性顺序，y 和z 的坐标小的放前面，大的放后面。假如一个块有二维的线程，那么将所有threadIdx.y是1的线程放在threadIdx.y是0的线程后面（注意这里是y 不是x，即按照 Fortran语言的顺序，前面的变化快，后面的变化慢），以此类推（这里要注意这种访问形式，在后面的1.2很有用，相当于二维的矩阵是转置的，不是正常的那种形式，注意，不过如果将y想象成矩阵的行，x想像成矩阵的列，那就没问题了）。

图1

上图就是将二维的线程块映射成一维的顺序形式。（注意顺序）。而对于三维的块来说，首先将threadIdx.z为0的所有线程按照线性顺序排列，以此类推。将warp中所有线程执行玩一个指令，再去执行另一个指令，这种执行方式叫做SIMT（single Instruction multiple thread，单指令多线程）。如果在一个warp中所有的线程执行的是相同的指令，那么工作效率最高，如果说在if-then-else结构中，决策条件基于threadIdx的值，比如if（threadIdx>2）｛｝这种的，那么会导致线程分支，即按照两个控制流分支路径，这样线程0、1、2和3、4、5就不同；如果循环条件基于threadIdx的值，则循环也会引起线程分支。这种用法自然会导致很多重要的并行算法的产生。

        这里举个归约数组和的算法：即对于一个有着N个元素的数组来说，想要求的所有元素的和，那么归约采取的是两两结合，即假如10个元素（(1,2),(3,4),(5,6),(7,8),(9,10)），然后第二次（（1,2,3,4）,(5,6,7,8),(9,10)）。这种类型的，两两归约，因为前面计算一次之后就只有10/2=5个数值，所以第二次只需要计算5/2=2次（每次执行完一次，再重新分配，比如前面的第二次（9,10）因为缺少2个，就无法被4整除，所以第二次只计算了2次，然后（1,2,3,4,5,6,7,8）,(9,10)）接着最后得到最终答案；代码如下：

1.__shared__float partialSum[]

2.unsigned int t = threadIdx.x;
3.for(unsigned int stride = 1;
4.    stride < blockDim.x; stride *= 2)
5.{
6.__syncthreads();
7.if(t %( 2 * stride) == 0)
8.   partialSUm[t] += paritialSum[t + stride];
9.}

执行过程如下：

图2

注意：上述的归约其实导致块中有一般的线程从来不执行，是很浪费的，所以需要修改kernel函数，留待作业。

在上面的代码中，就是当threadIdx.x的值为偶数的时候才执行加法，所以会导致那些不执行第8行代码的线程需要通过另外一条路径。下面是修改的kernel函数，即不再采用两两相加的方法，这样的是具有较少线程分支的kernel函数：

1.__shared__float partialSum[]

2.unsigned int t = threadIdx.x;
3.for(unsigned int stride = blockDim.x>>1;
4.    stride > 0; stride >>= 1)
5.{
6.__syncthreads();
7.if(t < stride)
8.   partialSUm[t] += paritialSum[t + stride];
9.}

修改之后的代码的性能有所不同，上述代码中通过移位来代替除法操作，降低开销。执行过程如下：

图3

因为第一次迭代中，线程0-255都执行加法，而线程256-511不执行加法。由于warp中包含的32个线程对应的threadIdx.x值是连续的，因此都0-7个warp所有线程都执行加法，第8-15个warp则跳过加法，由于warp中所有线程都通过相同的路径，所以没有线程分支。可是还是因为有if的存在，kernel函数的分支并未完全消除。在执行第5次迭代的时候，第8行代码的线程个数低于32.也就是，最后5次迭代中分别只有16、8、4、2、1个线程执行加法运算，所以仍然存在分支。

1.2全局存储器的带宽
        制约CUDA kernel函数的一个重要因素就是全局存储器的访问数据。之前有讨论过如何减少访问的流量来达到加速的目的，这里接着讨论存储器合并技术。使得更加有效的将数据从全局存储器中移动到共享存储器和寄存器上。因为cuda系统采用的是DRAM的全局存储器，这种DRAM单元为了加快数据访问的速度，采用并行进程的方式，即当DRAM芯片中的传感器接收到请求的单元的索引的时候，会顺带把其附近的单元的电位一起传送过来，如果应用程序在访问单元改变前能够充分利用这种来自多个连续单元的数据，则会比真正的随机顺序的单元访问要块得多。所以需要kernel函数安排数据的访问顺序。（本书发布时，现在不知道cuda的设备还是不是采用DRAM）在G80/GT200中，考虑这样一个事实：同一个waro中的线程在任何给定的时间内都执行同一条指令，也就是说当同一个wari中所有线程执行同一条指令访问全局存储器中的连续单元时，这种访问模式是最好的。如线程0访问单