三、 一些并行化技术应用实例

　　(1)CUDA平台下的box filter

　　Box filter是图像处理中非常常用的滤波处理方式，而且像素间是无关的。在并行算法设计中，主要有下面两种并行思想：细粒度和粗粒度。Box filter采用细粒度的并行方式，是和串行处理方式一样，将一帧图像细分为不同的处理过程。很明显，每一个点的处理需要n2-1步.这里可以有一个提高，分别沿着X和Y轴，完成两次一维的平均计算从而实现二维平均，这个并行过程可以在GPU上实现，这样就只需要n-1步了。

　　而粗粒度的并行实现方式，让每一个线程处理一整行和一整列，这样的方式比起前面的细粒度方式，能获得更好的性能。这是因为，粗粒度的处理方式，在两个GPU遍里面，每个线程能处理一个单一的点。然而，沿着X方向遍历的时候，对于全局存储器的读取方式是不满足合并访问规则的的，这是因为每个线程读取的是不同行的数据。采取纹理存储器能很好的规避这个问题，而对于Y方向的遍历就是满足合并访问规则的了。

　　(2)CUDA平台下的直方图生成

　　传统的图像直方图算法是通过逐点统计实现的，基本程序算法如下：

　　　　由于要确定的是整个图像的灰度分布情况，所以即使在并行处理环境中，不进行迭代运算也是无法生成图像的直方图的。要利用并行处理环境进行算法优化可以考虑将输入数据分成多个块，让某一个线程去处理其中一块形成一个局部直方图，最后归并成一个全局直方图。但是对于GPU平台这个算法需要进行调整，这是因为读取图像数据是顺序的而访问结果直方图数组是数据依赖的，考虑到GPU平台的存储器结构最优的是采用共享存储器而非全局存储器。GPU硬件上共享存储器的大小为16KB，共享存储区没有原子操作的支持。还有就是一个线程块应该包含128～256个线程以获得高效的运行。那么如果采用一个线程直接对应于一个子直方图的话，一个明显的限制就是16KB平均到最大一个包含192条线程的块，每个线程只有85bytes，因此这种方法最大情况共享存储区能适合每线程子直方图可达64单字节针脚数。字节计数器也引入了255字节的限制于单个执行线程所处理的数据大小，在执行线程间分割数据时需要考虑。也就是单个线程不能生成针脚数大于64的子直方图。

　　由于硬件提供了执行方式为SIMT的warp块(块内32 个并行线程为一组来创建、管理、调度和执行线程)，利用这一重要性质就可以设计共享存储区的原子操作。使用的方法是线程间不再有独立的子直方图而是每个warp(32个线程)共享同一存储区域生成子直方图，用这种方式存在warp内部共享存储区冲突。但硬件平台并未提供相应原子操作，所以这里需要利用一点程序设计技巧来解决冲突。下面就是解决冲突的方法：

　　　　这里图像数据值是从全局存储区读入的，取值范围是[0，255]，每个warp线程根据输入的图像数据使数组s_WarpHist[]相应位置增加1，数组s_WarpHist[]是与块内当前warp对应的块子直方图s_Hist[]的一部分。为了解决warp内的线程冲突，通过最末一个执行写入的线程将直方图标识，标识为直方图频数变量的高5位，因为warp的大小为32，所以只需要5位。这样的话，对于每一个线程读取当前直方图针脚值高5位就被替代为当前线程的标识符，然后将增加的频数写回子直方图的共享存储区。如果线程间要写入的不是同一针脚位置，就没有额外的操作，但当两个或多个线程在写同一位置时硬件会执行共享存储区写结合，使线程中被标识的频数被接受而拒绝所有其它未定的线程修改。在写尝试后每个线程读取来自同一个共享内存的位置。而那个可以写回自己修改的频数的线程，则退出循环并空闲等待warp中其余的线程。Warp继续它的执行，直到所有线程退出循环。由于每个warp使用自己的子直方图并且warp线程总是同步的，所以不依赖于warp的调度顺序，这个循环执行的迭代次数不会超过32次，而且只有在warp内所有线程读到的灰度值相同时才会发生。

　　子直方图(即warp直方图)计算算法中的第一步count = s_WarpHist[data] & 0x07FFFFFFU;是将当前该针脚的频数放入count变量中，并将高位清零。

　　第二步增加相应针脚频数，并写入线程标识符，这里特别说明的是threadTag = threadIdx.x << (32 - WARP_LOG_SIZE);也就是将线程号经过移位操作放到高5位作为标识符。

　　第三步则是将增加后的频数写回到数组s_WarpHist[]。然后判断当前写入的频数是不是最后修改频数的线程所写入的，是就接受，不是则通过循环重新执行上述操作。

　　在完成了子直方图计算后就是归并每一个子直方图，这个过程应该分两步首先归并各个warp的子直方图到块直方图，然后将块直方图归并为全局直方图。但因为对于全局存储区存在原子操作，所以块内或块间的线程可同时更新全局直方图而不会写冲突，所以不需要每个块输出自己的子直方图。

　　(3)CUDA平台下的Fan-in及优化分析

　　所谓Fan-in算法求和就是利用多处理器实现树型求和，图3.1是一个树型求和的实例：

图3.1 Fan-in算法处理实例

　　可以看到该算法是利用多处理器的并行计算能力，将计算时间缩短，也可以很清楚的看到并行处理程序的特色“浅而宽”。

　　在本算法中Fan-in算法求和方法主要应用在直方图的求和计算中，全局直方图的针脚数为256。使用Fan-in算法求和步骤如下：

　　第一步将0和128，1和129，…，127和255求和，将结果保存在0，1，…，127中;

　　第二步将0和64，1和65，…，63和127求和，将结果保存在0，1，…，63中;

　　继续重复以上计算方式，直到只剩一个元素为止，最后结果保存在元素0中，这里要特别说明的是结合CUDA程序特性线程在取数据的执行方式是线程0→线程1→……→线程n，有局部性原理可以知道，最优的数据寻址方式是相邻寻址，所以这里采用的是0和128这样间隔求和而非0和1这样的相邻求和叠加模型，即图3.3

　　这种Fan-in算法的另外一种实现形式是纵横加工算法[30]，其核心思想与Fan-in算法一致，但在实现上略有不同。算法如下：

　　假设，N个数a1，a2，…，an，其中N=P2

　　第一步：a1+a2，…，ap+1+ap+2，…,aN-P+1+aN-P+2

　　第二步：(a1+a2)+a3,…，(ap+1+ap+2)+aP+3,…,(aN-P+1+aN-P+2)+aN-P+3

　　第P-1步：(a1+a2+…+aP-1)+ap,…,(aP+1+aP+2+…+a2P-1)+a2P,…,(aN-P+1+aN-P+2+…+aN-1)+aN

　　第P步到第P-1+log2P步按Fan-in算法处理。

　　由此我们对以上两种算法做简要的理论分析，Fan-in算法使用的处理线程个数为N/2，加速比为Sf和效率Ef分别为

　　而使用纵横加工算法使用的线程个数为 ，加速比Sz和效率Ez分别为

　　从上面的式子可以看到，传统Fan-in算法的加速比高于纵横加工算法，但后者为常数效率的并行算法，也就是效率不随问题规模和处理器个数的改变而改变，也就是说当问题规模较小，线程切换延时也比较小的时候采用Fan-in算法，程序运行时间更快，当超过CUDA的kernel函数承担线程极限，函数切换延时较大时则应采用纵横加工算法获得的执行效率更高。本文根据所要解决的实际问题采用的是传统Fan-in算法的基本改进方式。

图3.2 CUDA程序中直方图求和的Fan-in算法应用