CUDA心得体会

访问global虽然仍然需要最多数百个时钟周期，但是总的来说比cpu的内存还是快多了...如果使用得当的话。CUDA是符合SIMD-PRAM模型的，也就是n台功能相同的处理机（block看成处理机,每个thread看成处理机的一部分更好些），一个容量无限大的共享处理器M。其中M就是global memory，虽然在CUDA中也是有限的，但是比起可怜的shared还是多很多了。

global作为共享存储器，最重要的功能就是用来进行block之间的通信。这里能不能叫做通信其实有问题，因为block之间实际上没有通信...CUDA上的通信实际上就是刷新，运行完一个kernel，把所有的需要与其他block交换的数据都写进global，然后再launch一个kernel读入上次生成的数据...如此往复循环。
如果kernel的输入与输出是同一个矩阵，就是一个原地（in-place）操作，下个kernel还是读取这个矩阵。
如果输入矩阵与输出矩阵不同，而且需要反复操作这两个矩阵，就构成了一个乒乓操作，例如第一次kernel使用ping为输入，pong为输出，下次就用pong输入，ping输出。这种方法在迭代等数值计算中经常使用。

texure也是全局存储器，速度比global还要更快，但是我不推荐大家在GPGPU中使用texture。这是因为texture是只读的，array中的数据不随对应的矩阵改变而改变，需要反复绑定，而绑定占用的时间不一定比使用texture带来的提高小。大多数情况下，使用global + shared比使用texture要快。当然，在需要大块的只读数据如查找表，索引表时可以使用texture。

关于bankconflict
简单的说，矩阵中的数据是按照bank存储的，第i个数据存储在第i%16个bank中。一个block要访问shared memory，只要能够保证以其中相邻的16个线程一组访问thread，每个线程与bank是一一对应就不会产生bank conflict。否则会产生bankconflict，访存时间成倍增加，增加的倍数由一个bank最多被多少个thread同时访问决定。有一种极端情况，就是所有的16个thread同时访问同一bank时反而只需要一个访问周期，此时产生了一次广播。

关于coalsced reading
global存储器需要按照索引连续访问
以下的访问方式可以获得不错的速度：
s_data[ tid] = d_data[ tid];
s_data[ tid] = d_data[ tid + n * 16]; //索引开始地址为半warp整数倍
s_data[ tid] = ( tid == 0)？0.0f : d_data[ tid]；//允许几个不读
以下的方式则很慢：
s_data[ tid] = d_data[ 128 - tid];
s_data[ tid] = d_data[ 15 + tid];
s_data[ tid] = tid%2? d_data[ tid+1]: d_data[ tid -1];

下面有一些小技巧可以避免bank conflict 或者提高global存储器的访问速度
1. 尽量按行操作，需要按列操作时可以先对矩阵进行转置
2. 划分子问题时，使每个block处理的问题宽度恰好为16的整数倍，使得访存可以按照 s_data[tid]=i_data[tid]的形式进行
3. 使用对齐的数据格式，尽量使用nvidia定义的格式如float3,int2等，这些格式本身已经对齐。
4. 当要处理的矩阵宽度不是16的整数倍时，将其补为16的整数倍，或者用malloctopitch而不是malloc。
5. 利用广播，例如s_odata[tid] = tid%16 < 8 ? s_idata[tid] : s_idata[15];会产生8路的块访问冲突
而用：s_odata[tid]=s_idata[15];s_odata[tid]= tid%16 < 8 ? s_idata[tid] : s_data[tid]; 则不会产生块访问冲突
6. 需要使用结构体时，使用矩阵的结构体，不要使用结构体的矩阵。一些结构体类型可以通过指针转换提高访显存速度
7. 需要大量随机读时，用texture

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的整数倍。

以下是我用thread数为256的block自适应处理不同尺寸问题的方法，对一些问题有效。
对尺寸为2-16（2,4,8,16）的写一类核，一个block处理128-16个问题。
对尺寸为32-128（32，64，（96），128）的核，一个block处理8-2个问题，但与小尺寸问题使用不同的优化方法（子问题大于warp size），
对尺寸为256的核，使用一个核处理
对尺寸256以上的核，使用多个核或者一个核内进行循环处理。
对尺寸不是上述2的n次方或者32整数倍的问题，拓展后按照上述方法处理。

我们这里尽量只介绍使用CUDA进行计算的方法，而少涉及具体的算法。
核内的几种操作:

缩减问题（reduction），推荐使用二叉树的方式来做。
例如：缩减加操作为可以用以下代码描述：
#define ADD     if (tid < 128) { temp[tid] += temp[tid + 128];} __syncthreads();\
    if (tid     if (tid     if (tid     if (tid     if (tid     if (tid     if (tid

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

 

邻域问题：
在图像处理中，经常可以遇到邻域问题，如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))的行向量
然后将各行输出的结果累加，即可得到最终结果。

 

源自： 张舒 Blog

来自 “ ITPUB博客 ” ，链接：http://blog.itpub.net/22785983/viewspace-619729/，如需转载，请注明出处，否则将追究法律责任。

转载于:http://blog.itpub.net/22785983/viewspace-619729/