昨天在该版询问了cuda改写cpu版本计算程序的问题，十分感谢ice、横扫千军两位斑竹和lianuosun站友的热情及时高水准的回复，帮我解决了久思不得其解的问题。
今天又遇到一个问题难以解决：
问题详情：
和之前问题一样，程序的物理意义是建立64512256的三维空间格点，我将数据分为444块，每块包含1612864个格点，一个物理上的格点的数据由一个thread计算。
现在希望实现空间所有格点的数据求和以及寻找最大值。
以数据求和为例：
核函数中参数的含义和值：
三个维度的格点数目 NScale_R 64,512,256
三个维度的块数目 Block_num_R 4,4,4
每一块三个维度的SM数目 Block_Width_R 16,128,64
总格电数目 Nsize = 64512256
x_g，y_g，z_g表示了物理上该格点的坐标，pos_g表示了格点数据在所有数据组成的数组中的位置。
数据求和
核函数功能：建立一个64512256的double共享数组tempD。每个thread保存一个tempD值，求所有tempD的和并存入EpsilonBar_cal；
函数调用：
…
EpsilonBar_cal=0.0；
dim3 dimGrid(Block_num_R[0],Block_num_R[1],Block_num_R[2]);
dim3 dimBlock(Block_Width_R[0],Block_Width_R[1],Block_Width_R[2]);
ComputeEpsilon_1_GPU<<<dimGrid, dimBlock>>>(… EpsilonBar_cal, Block_Width_R, NScale_R, Nsize);
…
核函数内容
global void ComputeEpsilon_1_GPU(… double *t_EpsilonBar_cal, int *Block_Width_g, int *t_NScale, int t_Nsize)
{
…//一些与求和无关的参数设定
int pos_g, pos_g1;
extern shared double tempD;
int x_g=blockIdx.xBlock_Width_g[0]+threadIdx.x;
int y_g=blockIdx.yBlock_Width_g[1]+threadIdx.y;
int z_g=blockIdx.z*Block_Width_g[2]+threadIdx.z;
//第一步：根据输入的数据计算每个格点的tempD值
if ((x_g<t_NScale[0])&&(y_g<t_NScale[1])&&(z_g<t_NScale[2]))
{
pos_g=(x_g*t_NScale[1]+y_g)*t_NScale[2]+z_g;
tempD[pos_g]=。。。。。。；
}
//第二步：单个block内的tempD求和
//总体思路：
第0步：将y_g=y0和z_g=z0的点的值求和，存入点(0,y0,z0)，这样只需求和平面x_g=0上的点；
第1步：将平面x_g=0上，z_g=z0的点的值求和，存入点(0,0,z0)，这样只需要求和直线x_g=0，y_g=0上的点；
第2步：将直线x_g=0，y_g=0上点的值求和，存入点(0,0,0)；
这样就只需要对直线上的点求和，采取的方案是反复两两求和得到结果。
unsigned int t_stride[3];
int b_x[3]={x_g,y_g,z_g};
int b_Idx[3]={blockIdx.x,blockIdx.y,blockIdx.z};
int t_Idx[3]={threadIdx.x,threadIdx.y,threadIdx.z};
int flag=0;
for (int i=0;i<3;i++)
{
flag=0;
for (t_stride[i]=1;t_stride[i]<Block_Width_g[i];t_stride[i]=2)
{
__syncthreads();
flag=1;
if (t_Idx[i]%(2t_stride[i])!=0) flag=0;
if ((x_g>=t_NScale[0])&&(y_g>=t_NScale[1])&&(z_g>=t_NScale[2])) flag=0;
if ((i==1)&&(t_Idx[0]!=0)) flag=0;
if ((i==2)&&((t_Idx[0]!=0)||(t_Idx[1]!=0))) flag=0;
int tempI=b_x[i]+t_stride[i];
if (tempI<((b_Idx[i]+1)*Block_Width_g[i])) flag=0;
if (tempI<t_NScale[i]) flag=0;
if (flag==1)
{
pos_g=(b_x[0]*t_NScale[1]+b_x[1])t_NScale[2]+b_x[2];
if (i==0) pos_g1=(tempIt_NScale[1]+b_x[1])*t_NScale[2]+b_x[2];
if (i==1) pos_g1=(b_x[0]*t_NScale[1]+tempI)*t_NScale[2]+b_x[2];
if (i==2) pos_g1=(b_x[0]*t_NScale[1]+b_x[1])*t_NScale[2]+tempI;
tempD[pos_g]+=tempD[pos_g1];
}
}
__syncthreads();
}
__syncthreads();
//第三步：将每个block中点(0,0,0)的数据加在t_EpsilonBar_cal[0]上
if ((threadIdx.x==0)&&(threadIdx.y==0)&&(threadIdx.z==0))
{
pos_g=(x_g*t_NScale[1]+y_g)*t_NScale[2]+z_g;
t_EpsilonBar_cal[0]+=tempD[pos_g]/t_Nsize;
}
}
在与CPU版本的对照中，我确认这个求和结果是存在问题的，那么，问题出在什么地方？