在进行基于体数据的PDE计算时,总是会涉及到邻接单元(neighbir cell)的访问,想要提高计算效率就必须尽量共享邻域数据,减少访问全局内存的次数.不同于二维的情况,尤其是需要多次迭代计算时,三维纹理的效率很多时候差强人意,且需要
在下一步迭代开始前进行大量数据的复制操作.而如果用二维纹理,cache命中率让人不敢恭维,且同样需要内存复制操作.虽然线型层次化纹理不需要复制操作,但最好的情况下也只能获得二维纹理同样的缓寸奖励. 这里介绍基于共享缓寸提高邻域数据共享率以减少全局内存访问次数的共享内存和寄存器循环切片法（SRC:Slice of Shared Memory-Refister Cycle Method）. 由于我没有有效的画图功能工具，所以虽然不情愿而且可能还会影响你对我在这里所说的话的理解，但是也只能如此，如有疑问可以联系我--
QQ:295553381
Email:cyrosly@163.
不多废话。首先，让我们在脑中设想一个立方体作为对体数据的抽象，然后将该立方体在三个维方向上均匀的分割，然后我们 将所有分割的单元（cell)影射到平铺的thread grid上。注意：这里cell = sub volume.
grid policy: 每个线程块（thread-block)处理一个cell,所以这里整个grid的线程的数量并不等于整个体数据的元素数量。
从每个cell 的第一层开始，首先将第一遍的沿Z轴正反方向的邻接点载入寄存器。然后在后面的循环中循环交换前中后
处理单元的次序，并更新共享内存。 共享内存的大小=(blockDimY+2)*(blockDimX+2)。
不同的机器最高效的块尺寸的选择不同。在我的机器上16x8可以获得最好的性能。
虽然如何分配共享缓寸可能有N种策略。当然如果你愿意。不过这样的划分可以最小化使用量（当然可以用register代
替，但是那样可并发的thread-block数量就会因为平均可用寄存器的数量增加而减少反而影响性能（实验过）。
另外也不要疑惑，在这里（后面的程序）的共享内存访问不会有bank conflicts(所以不用在X维度使用多余的辅助空
间来避免冲突）。 为了便于计算线程对数据的索引（注意，由于每个thread-block-plane对应的是一个体数据，而不是简单的平面映射
所以计算索引时要格外小心，你 也可能会因为在大脑中组织空间的逻辑结构时的”纠缠不清“而崩溃o(^!^)o。所以不要抓狂，
你可以简单的将它设置成对应的体数据的X-Y层大小，而在程序中沿Z轴顺序后推，事实上很多时候循环体的层数对性能影响
不大，例如256x256x256与512x512x64相同大小的体数据但决定不同的grid划分策略时，效率几乎没有差别,当然也只是有些时候.后面的测试你将会看到差距.
然后的计算中每个循环遍中都有一slice共享内存和寄存器中的的数据被替换和交换：计算完后当前slice的后置slice
成为下一遍的当前slice，而前置slice将连同后置slice的邻接单云的的数据一并更新。并成为后置slice。而当前slice
将成为前置slice.
最后还需要注意边界条件的处理，需要一个单独的内核来计算边界单元的值。
算法描述：
<1>
global/tex
load------------>register:front neghbors
global/tex
load------------>register:middle
global/tex
load------------>register:back neghbors
global/tex
<2>
__LOOP__[ 0 : layers-1 ]
<2|0>
smem
store<--------registers
/registers
compute
\shared mem
sync
<2|1>
update shared mem slice
sync
swap order of slice
<2|2>
store the result to global memory which computed in step <2|0>
kernel code:
__constant__ uint gridSizeU;
__constant__ uint gridSizeV;
__constant__ uint slice;
__constant__ uint subvolume;
__constant__ uint sublayers;
extern "C"
__global__
void kernel_volflo_pressure( float* target,
const float* source,
const float* div,
float coeff )
{
__shared__ float cell[ DIMY+2 ][ DIMX+2 ];
const uint tidx=IM( DIMX, blockIdx.x )+threadIdx.x;
const uint iU=tidx&( gridSizeU-1 )+1;
const uint iV=IM( DIMY, blockIdx.y )+threadIdx.y+1;
const uint subdomainID=tidx/gridSizeU;
uint idx=slice+IM( subdomainID, subvolume )+IM( gridSizeU, iV )+iU;
float cell010=source[ idx ];
float cell100=source[ idx-slice ];
float cell001=source[ idx+slice ];
const uint slotU=threadIdx.x+1;
const uint slotV=threadIdx.y+1;
const uint cc0=threadIdx.x==0;
const uint cc1=threadIdx.y==0;
for( uint layer=0; layer
cell[ slotV ][ slotU ]=cell010;
if( cc0 ){
cell[ slotV ][ 0 ]=source[ idx-1 ];
cell[ slotV ][ DIMX+1 ]=source[ idx+DIMX ];
}
if( cc1 ){
cell[ 0 ][ slotU ]=source[ idx-gridSizeU ];
cell[ DIMY+1 ][ slotU ]=source[ idx+IM( DIMY, gridSizeU ) ];
} __syncthreads();
target[ idx ]=0.16666667f*( cell[ slotV ][ slotU-1 ]+
cell[ slotV ][ slotU+1 ]+
cell[ slotV-1 ][ slotU ]+
cell[ slotV+1 ][ slotU ]+
cell100+cell001-coeff*div[ idx ] );
cell100=cell[ slotV ][ slotU ];
cell010=cell001; idx+=slice;
cell001=source[ idx ];
}
}
性能比较：在512x512x64的网格上（我的机器上一次可以分配的显存的最大尺寸）
配置：
GPU:8800 GTS 320MB
CPU:DualCore Intel Core 2 Duo E6650,2333 MHz(7x333)
效率对比：
CPU:只计算一次大约14秒，
GPU:迭代计算500次大概3秒
Profile片段:
GTS 8800 320MB:
Reference:
<
--Tobias Brandvik and Graham Pullan
--Whittle Laboratory , Department of Engineering,
--University of Cambridge , Cambridge , CB3 0DY,UK
<
--Robert Bridson1:University of British Columbia
--Matthias M¨ uller-Fischer: AGEIA Inc.
本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/Cyrosly/archive/2009/04/13/4069279.aspx
fig0.jpg来自 “ ITPUB博客 ” ，链接：https://blog.itpub.net/23056049/viewspace-627794/ ，如需转载，请注明出处，否则将追究法律责任。
