万能的斑竹大人，帮我看看这个程序那里有问题。我是在原来那个图像相关的例子上修改的，原来帖子链接如下：http://cudazone.nvidia.cn/forum/ … 7019&extra=page%3D1，原来在斑竹的指导下，运行很正常，后面经过修改如下：
int bky=blockIdx.y;
int bkx=blockIdx.x;
int tdy=threadIdx.y;
int tdx=threadIdx.x;
int y=SubMatrix[bky*m_nRow+bkx+m_nRow*m_nCol]; //原图block的中心点
int x=SubMatrix[bky*m_nRow+bkx];
int i,j,tempv;
float XPos,YPos;
__shared__ float result[6],piArray[6];
__shared__ float corr,iter,fm,gm,fsum,gsum,fgsum;
__shared__ int M2;
__shared__ float shared6[5808];
if (tdy*blockDim.x+tdx==0){
for (i=0;i<6;i++){ //初始参数
result[i]=0;
piArray[i]=InitDeformPr[(bky*gridDim.x+bkx)*8+i];
}
}
M=M/2;
i=tdy-M;
j=tdx-M; //以块的中心为坐标
M2=(2*M+1)*(2*M+1); //(2*M+1)*(2*M+1)常用
/************************************************************************/
/*求以(x,y)为中心的,(2M+1)*(2M+1)为模板的的修正值(即解cdpi,cdpj的方程) */
/************************************************************************/
do
{
//////////////////////////////////////////////////////////////////////////
//提供计算数据
if(tdy*blockDim.x+tdx==0){//这里是对share memory的写入,要防止循环写入!
matrix_subGPU(piArray,6,1,result,piArray);
}
if(y+i<0||x+j<0||y+i>=Height||x+j>=Width){ //线程出界取值为0,不要return不然不能同步
shared6[0*M2+tdy*blockDim.x+tdx]=0;
}else
shared6[0*M2+tdy*blockDim.x+tdx]=m_ptrData1[(y+i)*Width+x+j]; //0~f
YPos=y+i+piArray[3]+piArray[4]*j+piArray[5]*i; //点(x,y)通过形函数对应点(g_x,g_y)
XPos=x+j+piArray[0]+piArray[1]*j+piArray[2]*i;
if(int(XPos)<0||int(YPos)<0||int(XPos)>Width-1||int(YPos)>=Height-1){
shared6[1*M2+tdy*blockDim.x+tdx]=0;
}
GetAPixelGPU(XPos,YPos,CoefMatrix16,shared6[1*M2+tdy*blockDim.x+tdx],Width); //1~g
//////////////////////////////////////////////////////////////////////////
//用规约求fm,gm
shared6[2*M2+tdy*blockDim.x+tdx]=shared6[0*M2+tdy*blockDim.x+tdx]; //数据,这里传给2、3是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=shared6[1*M2+tdy*blockDim.x+tdx];
__syncthreads();
for (tempv=M2/2;tempv>0;tempv=tempv>>1){ //奇数做规约
if (tdy*blockDim.x+tdx<tempv){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+tempv];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+tempv];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
fm=shared6[2*M2+0];
gm=shared6[3*M2+0];
fm+=shared6[2*M2+blockDim.y*blockDim.x-1];
gm+=shared6[3*M2+blockDim.y*blockDim.x-1];
fm=fm/(M2);
gm=gm/(M2);
}
//////////////////////////////////////////////////////////////////////////
//用规约求相关系数cor
shared6[2*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-fm)*(shared6[0*M2+tdy*blockDim.x+tdx]-fm); //数据,这里传给3,4,5是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=(shared6[1*M2+tdy*blockDim.x+tdx]-gm)*(shared6[1*M2+tdy*blockDim.x+tdx]-gm);
shared6[4*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-fm)*(shared6[1*M2+tdy*blockDim.x+tdx]-gm);
__syncthreads();
for (tempv=M2/2;tempv>0;tempv=tempv>>1){//规约,得到fsum,gsum,fgsum
if (tdy*blockDim.x+tdx<tempv){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+tempv];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+tempv];
shared6[4*M2+tdy*blockDim.x+tdx]+=shared6[4*M2+tdy*blockDim.x+tdx+tempv];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
fsum=shared6[2*M2+0];
gsum=shared6[3*M2+0];
fgsum=shared6[4*M2+0];
fsum+=shared6[2*M2+blockDim.y*blockDim.x-1];
gsum+=shared6[3*M2+blockDim.y*blockDim.x-1];
fgsum+=shared6[4*M2+blockDim.y*blockDim.x-1];
corr=fgsum/sqrt(fsum*gsum);
}
在这里前面的帖子是直接坐两个图像的相关，这里加入了一个对应关系：
YPos=y+i+piArray[3]+piArray[4]*j+piArray[5]*i; //点(x,y)通过形函数对应点(g_x,g_y)
XPos=x+j+piArray[0]+piArray[1]*j+piArray[2]*i;
即原图的(y+i,x+j)对应变形后的(XPos，YPos)，因为这个地方可能是小数，所以进行插值，插值函数如下：
GetAPixelGPU(XPos,YPos,CoefMatrix16,shared6[1*M2+tdy*blockDim.x+tdx],Width); //1~g
插值函数是没问题的，已经和串行的结果比对过。
现在自己发现的问题是在这里，有时候代码执行到最后时，发现
shared6[2*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-fm)*(shared6[0*M2+tdy*blockDim.x+tdx]-fm); //数据,这里传给3,4,5是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=(shared6[1*M2+tdy*blockDim.x+tdx]-gm)*(shared6[1*M2+tdy*blockDim.x+tdx]-gm);
shared6[4*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-fm)*(shared6[1*M2+tdy*blockDim.x+tdx]-gm);
这里面根本没有执行，数据都为xxxe+10，xxxe+11这种巨大的数据，有时候又是对的，有时候每次求出的均值有不同，这个到底是怎么回事啊？一直怀疑是同步的时候有问题，但是一直不知道到底是不是。。万能的斑竹。。求指导。。