应该用share memory还是global memory

system · 2013 年7 月 7 日 13:46

斑竹你好，在写图像相关的代码中遇到个问题。代码如下：

__global__ void CorrelationGPU(float *Img1 ,float *Img2,int Height,int Width,int *SubMatrix,int M,int m_nRow,int m_nCol)
{
int bky=blockIdx.y;        
int bkx=blockIdx.x;
int tdy=threadIdx.y; 
int tdx=threadIdx.x;
int i,j;
int y=SubMatrix[bky*m_nRow+bkx+m_nRow*m_nCol]; //原图block的中心点
int x=SubMatrix[bky*m_nRow+bkx];        
__shared__        float x1,x2,x3;
__shared__ float shared6[4847];        //大小是(M/2*2+1)*(M/2*2+1)*3,存储两个模板的值
M=M/2;
i=tdy-M;
j=tdx-M; //以块的中心为坐标
int M2=(2*M+1)*(2*M+1); //让块取成奇数*奇数,方便用

//////////////////////////////////////////////////////////////////////////
//求以(x,y)为中心的,(2M+1)*(2M+1)为模板的的相关系数        
shared6[0*M2+tdy*blockDim.x+tdx]=Img1[(y+i)*Width+x+j]; //f
shared6[1*M2+tdy*blockDim.x+tdx]=Img2[(y+i)*Width+x+j]; //g

shared6[2*M2+tdy*blockDim.x+tdx]=shared6[0*M2+tdy*blockDim.x+tdx]; //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=shared6[1*M2+tdy*blockDim.x+tdx];
__syncthreads();

for (i=M2/2;i>0;i=i>>1){ //奇数做规约,这里用i节省空间
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
x1=shared6[2*M2+0]; 
x2=shared6[3*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1]; //做的奇数规约x1~fm,x2~gm
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x1=x1/(M2);
x2=x2/(M2);

//用规约求相关系数cor
shared6[2*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[0*M2+tdy*blockDim.x+tdx]-x1); //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=(shared6[1*M2+tdy*blockDim.x+tdx]-x2)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2);
shared6[4*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2); 
__syncthreads();

for (i=M2/2;i>0;i=i>>1){//规约,得到fsum,gsum,fgsum
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
shared6[4*M2+tdy*blockDim.x+tdx]+=shared6[4*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
x1=shared6[2*M2+0]; //x1~fsum/cor,x2~gsum,x3~fgsum
x2=shared6[3*M2+0];
x3=shared6[4*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1];
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x3+=shared6[4*M2+blockDim.y*blockDim.x-1];
x1=x3/sqrt(x1*x2);
}

代码说明如下：
1：其中Img1,Img2是两幅图，高和宽分别为Height、Width，把两个图分割成m_nRowm_nCol份，每份的大小是MM，中心点存储在数组SubMatrix中，这里M=31。
2：在核函数中申请了一个共享数组，用来存储Img1、Img2的数据，以便用来计算相关系数；
3：因为线程是块的左上角为原点的，为了方便这里用i=ty-M;j=tx-M,把坐标原点放到了块的中心位置；
4：做相关的公式如下：
[attach]3248[/attach]

具体描述在原来这个帖子中有：http://cudazone.nvidia.cn/forum/forum.php?mod=viewthread&tid=7001
5：因此先用规约求均值：

shared6[0*M2+tdy*blockDim.x+tdx]=Img1[(y+i)*Width+x+j]; //f
shared6[1*M2+tdy*blockDim.x+tdx]=Img2[(y+i)*Width+x+j]; //g

shared6[2*M2+tdy*blockDim.x+tdx]=shared6[0*M2+tdy*blockDim.x+tdx]; //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=shared6[1*M2+tdy*blockDim.x+tdx];
__syncthreads();

这里第一部分是传入数据，第二部分是把数据拷贝到共享数据的另一块，因为在后面求相关系数的时候还要用，但是在递归求和的时候要改变数组里面的数据，所以这样是为了保护原有数据；
6：这部分就是用规约求均值：

for (i=M2/2;i>0;i=i>>1){ //奇数做规约,这里用i节省空间
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){ //这里是把和值给所有的线程,如果给share,后面再除的时候要出问题(出现多次除法)
x1=shared6[2*M2+0]; 
x2=shared6[3*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1]; //做的奇数规约x1~fm,x2~gm
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x1=x1/(M2);
x2=x2/(M2);

因为数组为奇数，所以在后面补加了一项。在这里得到的x1和x2就是两个块的均值。
7：但是在程序这里出了问题：

shared6[2*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[0*M2+tdy*blockDim.x+tdx]-x1); //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=(shared6[1*M2+tdy*blockDim.x+tdx]-x2)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2);
shared6[4*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2); 
__syncthreads();

for (i=M2/2;i>0;i=i>>1){//规约,得到fsum,gsum,fgsum
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
shared6[4*M2+tdy*blockDim.x+tdx]+=shared6[4*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
x1=shared6[2*M2+0]; //x1~fsum/cor,x2~gsum,x3~fgsum
x2=shared6[3*M2+0];
x3=shared6[4*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1];
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x3+=shared6[4*M2+blockDim.y*blockDim.x-1];
x1=x3/sqrt(x1*x2);

在第6步中x1,x2的值还是正确的，这里是58.064，但是在往下面走，到第7步中的__syncthreads值就变成了0.11869，这里明显出现错误。但是如果把__shared__ float x1,x2,x3定义成float x1,x2,x3就不会出现上面的错误。
8：我想是不是在这里因为x1,x2,x3是共享的，所以x1=x1/(M2);x2=x2/(M2) 被执行了多次，但是只是我猜测而已，还请版主大人指教。
9：我这里要用共享空间的原因是往后面做我是要用上面的数据得到一个方程的解，比如AX=B，其中A、B都是共享空间上的，A是66的，B是61的，通过求解出来的X值判断一个块内循环是否继续。即所有块内的线程就以这个共享的空间x为判断标准。（这个可能没见多明白，但是感觉和问题关系不大。。谢谢）

system · 2013 年7 月 7 日 15:19

LZ您好：

if (tdyblockDim.x+tdx==0){ //这里是把和值给所有的线程,如果给share,后面再除的时候要出问题(出现多次除法)
x1=shared6[2M2+0];
x2=shared6[3M2+0];
}
x1+=shared6[2M2+blockDim.yblockDim.x-1]; //做的奇数规约x1~fm,x2~gm
x2+=shared6[3M2+blockDim.y*blockDim.x-1];
x1=x1/(M2);
x2=x2/(M2);

您第六步这里已经出了问题。
x1和x2是shared memory中的变量，整个block都只有一份的，而您的kernel却是每个线程都跑一份。因此您后面对奇数规约的扫尾工作和除法求平均值的操作，其行为与您预期并不相同。

（以及，根据手册的说法，如果一个warp都向shared memory的同一位置普通写入（不是atomic操作），那么只会有一个线程写入成功，并且不保证是哪个线程。）

所以不管是累加还是相除，这里行为都不正确的。

您可以向前面“x1=shared6[2M2+0];
x2=shared6[3M2+0];”的时候一样，让0号线程干这个活，就能保证只干一次了。
以及，这就是您将后面4行写入前面if()的{}里面即可。

至于您直接float x1，x2也可以正确，那是因为各个线程都干了相同的活而已。
是正确的，也可以回写入shared memory，但是多了很多无用功不是么？

大致如此，请您思考一下shared memory的操作，并做相应修改。
祝您好运~

system · 2013 年7 月 8 日 01:00

嗯。。谢谢斑竹的解答。修改过后果然对了。。。我对斑竹解释的理解如下：
1：出现了多次操作的情况，在这里如果我块里面的线程数小于一个warp（32）,最初这种写法就应该是对的。因为“如果一个warp都向shared memory的同一位置普通写入（不是atomic操作），那么只会有一个线程写入成功，并且不保证是哪个线程”但是这里不论是那个线程写入都是一样的，因为他们是操作的是share memory。
2：但是在一个SM上同时只能运行一个warp，所以如果线程数大于了warp数（32），就会出现多个warp串行操作，就是会进行多次的累加和累除操作？那这里进行的次数是(块内线程数/32)次么？谢谢斑竹

system · 2013 年7 月 8 日 02:50

LZ您好：

不客气的，我和横扫斑竹将竭诚为您服务。

您的理解大体上是对的，我将稍微补充说明并总结如下：

1：如果您是一个warp向shared memory同一个地方写入相同的数据，那么只会有一个线程成功，并且按照手册说明，并不保证是哪个线程成功。但是因为这些线程干的都是同样的活，那么只要有一个成功即可，无需确切说明是哪个成功。

2：如果您是多个warp像shared memory同一个地方写入相同的数据，那么每个warp会有一个线程成功，究竟是warp内哪个线程成功并无法保证。但由于所有的线程干的都是同样的活，所以结果也是正确的，一共会写入warp数目次相同的结果。

3：上述1：，2：两种情况，都可以通过限定某个特定线程来干活（比如指定thread 0干活），也都可以不指定，在逻辑上都是正确的。虑及效率，建议单warp的情况直接写入即可，多warp的情况使用判断，将写入线程规定为特定线程或者特定warp均可。

因为CUDA的最小执行单元就是warp，一个warp里面指定一个线程写入和一个warp直接写入，硬件处理为一个线程写入，在写入上是一样快的，但是指定特定线程写入的话，会增加if的判断指令开销。

对于多个warp写入，建议if判断+一次写入的形式。

以及，为了逻辑清晰，您也可以始终加上if指定特定的线程干活，这个开销并不大的。

4：如果算法逻辑是如同您代码中那样，类似的累加等只需要进行一次，那么建议您if指定特定线程干活。
如果此时是单warp或者指定为单warp的话，因为warp内各个线程干的活一样，以及最后会随机有一个线程写入结果（warp内各线程结果相同），结果也是正确的。但为了不给您的编程和以后的理解代码造成困扰，建议直接按算法逻辑书写。
如果此时是多个warp的话，那么必须指定特定的线程（推荐）干活或者特定的warp（也可以）干活。而如您之前的写法，会导致该操作重复执行warp数目的次数。

另外需要说明的是，此时并非您说的那样“但是在一个SM上同时只能运行一个warp”，实际上新的GPU核心上一个SMX的SP数目远多于warp宽度32。以及，一个SM/SMX上的shared memory同时只能为一个warp服务。（请勿较真这里的“同时”是几个周期）
但是，无论一个SM上同时有几个warp，按照您之前kernel的写法，每个warp都会去操作一次shared memory。这是您算法实现规定的，而和具体的硬件细节无关。

5：如果您的多个线程（无论是单一warp内还是多个warp内）的结果都需要累加到同一个shared memroy位置，此时需要使用原子操作，能确保这些线程强制串行地完成累加或者其他积累性的操作。

相关内容大致总结如上，供您参考。

祝您编码顺利~

system · 2013 年7 月 9 日 01:13

谢谢斑竹耐心详细的解答，现在基本明白了。。谢谢

system · 2013 年7 月 9 日 02:26

不客气的，欢迎您常来论坛~

祝您编码顺利~