CUDA 技巧与经验 关于block、thread

http://tech.it168.com/a2009/1116/809/000000809882.shtml

在使用CUDA进行GPGPU计算时，global + shared的黄金组合在速度上远远超过了texture，只有在以下两种情况下使用texture：

1.需要图像输出时
2.需要反复随机访问的，但内容不变的大块内容，如索引表，查找表
　　其他时候尽量使用global + shared

block、thread划分的简单原则：

用不同的block处理完全不相关的数据可以获得最好的性能。此时只需要在block内进行数据交换。例如可以用每个block处理一个独立的矩阵，或者大型矩阵中的一部分(如一行，或者一个子矩阵)。需要处理一列数据时尽量将其转化为按行存储后再处理，以避免bank conflict。
　　CUDA的thread是一种轻量级线程，比较适合处理小粒度问题，在block一级也仍然只适合处理中小粒度的问题。而使用MPI构造的系统，一般是在比较大的粒度上进行并行操作，以减少数据交换。因此即便是已经经过检验的传统并行算法，在CUDA上移植时也要注意是否有可能进行进一步的划分。
　　block与thread的划分原则与可以并行处理的最小子问题占用的存储器大小和子问题的数量都有关系。不同尺寸的同一问题，block和thread划分方法往往不同。例如子问题占用存储器资源有以下几种情况：
　　在同一block的shared memory可以储存多个并行子问题的数据，此时可以尽量用一个block处理相邻的并行子问题，或者若干个子问题合并成一个稍大的子问题。
　　每个block的shared memory只能储存一个并行子问题的数据，此时就是简单根据输入数据算输出
　　每个block的shared memory不足以处理一个并行子问题，此时可以使用texture和global作为存储器，或者将子问题进一步分拆为数个串行过程，用若干个kernel完成。
　　每个block的thread数量以192-256为宜，处理的矩阵或者图像尺寸最好长宽都是16的整数倍，如果可能尽量将其补成16的整数倍。
　　我们这里尽量只介绍使用CUDA进行计算的方法，而少涉及具体的算法。

核内的几种操作:

缩减问题(reduction)，推荐使用二叉树的方式来做。
　　例如：缩减加操作为可以用以下代码描述：

#define ADD     if (tid < 128) { temp[tid] += temp[tid + 128];} __syncthreads();\
    if (tid <  64) { temp[tid] += temp[tid +  64];} __syncthreads();\
    if (tid <  32) { temp[tid] += temp[tid +  32];} __syncthreads();\
    if (tid <  16) { temp[tid] += temp[tid +  16];} __syncthreads();\
    if (tid <   8) { temp[tid] += temp[tid +   8];} __syncthreads();\
    if (tid <   4) { temp[tid] += temp[tid +   4];} __syncthreads();\
    if (tid <   2) { temp[tid] += temp[tid +   2];} __syncthreads();\
    if (tid <   1) { temp[tid] += temp[tid +   1];}\

基于同样的树状结构进行缩减操作，还可以进行排序，累乘，寻找最大最小值，求方差等等
　　

邻域问题：

　　在图像处理中，经常可以遇到邻域问题，如4-邻域，8-邻域等。基本思想就是使用图像中的几个元素按照一定的算子计算输出矩阵中的一个元素
　　常见的8-邻域算子如:

-1 -1 -1 1 1 1
　　-1 8 -1 1 1 1
　　-1 -1 -1 高通滤波器 1 1 1 低通滤波器

等等，可以完成锐化、平滑、去噪，查找边缘等功能
　　cuda的imagedenoising例子使用的是texture，其实也可以使用shared memory来做。

基本思路是：
　　1.在原有图像周围扩展像素，向外扩展像素数由邻域决定，如4-邻域和8-邻域扩展1个像素，24-邻域扩展2个像素等
　　2.每个block处理输出图像中长宽为16x16的一小部分，每个thread处理输出图像中的1个像素。这样每个block的线程数是256，符合nv推荐的192-256的范围。
　　texture的版本可以直接参考imagedenoising
　　shared memory的版本的可以参考以下代码(仅作示例，不能运行，处理4-邻域问题，其中i_data为global中的原始图像,func代表使用算子计算的过程)

shared type s_data[18][18];
shared type o_data[16][16];
s_data[tid.x + 1] = i_data[bid.x*16+bid.y*640*18+tid.x-640];
s_data[18*17 + tid.x + 1] = i_data[bid.x*16+bid.y*640*18+tid.x];
s_data[18+tid.x*18] = i_data[bid.x*16+bid.y*640*18+tid.y*640-1];
s_data[35+tid.x*18] = i_data[bid.x*16+bid.y*640*18+tid.y*640+1];
s_data[(tid.x+1)*18+(tid.y)+1]=i_data[bid.x*16+bid.y*640*18+tid.y*640+tid.x]；//中间256个
o_data[tid.x*16+tid.y]=func(
s_data[tid.x*18+tid.y+1],//邻点1
s_data[(tid.x+1)*18+tid.y],//邻点4
s_data[(tid.x+1)*18+tid.y+1],//中心点
s_data[(tid.x+1)*18+tid.y+2],//邻点2
s_data[(tid.x+2)*18+tid.y+2],//邻点3)；

这种方法可以归纳为：以输出结果组织数据，将与输出相关的数据输入同一block
　

投影问题：

　　这里的投影问题特指通过投影进行维度压缩的问题，如医学成像中常见的radon变换(当然，医学影像处理中更常见的是radon逆变换，但逆变换用傅式变换处理更好，此处不再赘述)。

　　这一类问题有一些例子：给手照x光片，我们的手是三维的，而最后出来的x光片则是平面的。做核磁共振(CT)时，实际上传感器是一个线阵，一直绕着我们的身体旋转，在不同的角度将我们身体的一个二维的截面投影到一维的传感器上。
　　以CT为例，实质上是将输入矩阵进行旋转后投影到一行上。设图像与传感器夹角为theta，可以根据tid.x * sin(theta) + tid.y*cos(theta)求出点(tid.x, tid.y)在线阵上的投影。此时结果与输入的数据是一对多的关系，而且随着角度不同，很难确定是由哪些点决定最终输出。
　　此时可以采用这样的方法：按照输入组织数据，每个block处理一行或者一个子矩阵并各自输出结果，最终将各block结果进行累加。
　　用这样的方法解决CT造影问题，将输入图像按行输入各个block，各block将旋转后的行存入一个长度为floor(1.414*max(width,height))的行向量
　　然后将各行输出的结果累加，即可得到最终结果。