CUDA C/C++：相同的可執行文件給出了第一次運行的不同結果

Question

以我之前的回答問題：My Previous Question，順便說一句，正確回答了Robert Crovella。

我想出了另一個計算隨機步驟到一個點的內核（通過使用與我以前的問題相同的RNG）並計算該點相對於其先前位置（座標）的能量差異。這是內核：

__global__ 
void DeltaE(float *X, float *Y,float *Xn, float *Yn, float *step,float *DELTA, curandState *state, const int N,const int n){ 

    int tIdx = blockIdx.x*blockDim.x + threadIdx.x; 
    int bIdx = blockIdx.x; 
    //int sIdx = blockIdx.x*blockDim.x; 

    float x , y; 
    float rdmn1, rdmn2; 

    float dIntE = 0.0e0f, dConfE = 0.0e0f, dTotE = 0.0e0f; 
    if(tIdx < N){ 

     if(tIdx == n){ 
      step[tIdx] = 0.2; 
      rdmn1 = curand_uniform(&state[tIdx]); 
      rdmn2 = curand_uniform(&state[tIdx]); 

      Xn[tIdx] = X[tIdx] + step[tIdx]*(2.0e0f*rdmn1 - 1.0e0f); 
      Yn[tIdx] = Y[tIdx] + step[tIdx]*(2.0e0f*rdmn2 - 1.0e0f); 
      dConfE = - (X[tIdx]*X[tIdx] + Y[tIdx]*Y[tIdx]); 
      dConfE += Xn[tIdx]*Xn[tIdx] + Yn[tIdx]*Yn[tIdx]; 

     } 
     else{ 
      x = X[tIdx] - X[n]; 
      y = Y[tIdx] - Y[n]; 

      dIntE += -1.0e0f/sqrt(x*x + y*y); 
     } 
     __syncthreads(); 
     if(tIdx != n){ 
      x = X[tIdx] - Xn[n]; 
      y = Y[tIdx] - Yn[n]; 

      dIntE += 1.0e0f/sqrt(x*x + y*y); 
     }  
     dTotE = dConfE + dIntE; 
     dTotE = ReduceSum2(dTotE); 
     if(threadIdx.x == 0)DELTA[bIdx] = dTotE; 


    } 
} 

然後我做CPU上最後一筆：

cudaMemcpy(&delta[0],&d_delta[0],blocks.x*sizeof(float),cudaMemcpyDeviceToHost); 
float dE = 0; 
for(int i = 0; i < blocks.x; i++){ 
    dE += delta[i]; 
} 

我的內核具有以下配置推出：

dim3 threads(BlockSize,BlockSize); 
dim3 blocks(ceil(Np/threads.x),ceil(Np/threads.y)); 
DeltaE<<<blocks.x,threads.x,threads.x*sizeof(float)>>>(d_rx,d_ry,d_rxn,d_ryn,d_step,d_delta,d_state,Np,nn); 

中，Np是多少點（我用1k - 4k）。我有一個GeForce 9500 GT，它不支持雙倍。我編譯時不使用flag/no選項。例如，取Np = 1k。當我編譯然後運行時，結果是dE = 6.557993。當我跑第二，第三，第四，無論什麼時候，它都是dE = -0.3515406。有誰知道這是從哪裏來的？

P.S .:我忘了提及，在DeltaE之前就可以找到My Previous Question的相同內核AvgDistance。我不知道這是否有什麼可做的事，但我認爲值得一提。 P.2：nn是任何選擇的點（粒子）。

Answer 1

正如通過評論指出，通過上述Robert Crovella，什麼可能發生的事情是，雖然tIdx = n被計算Xn[n]和Yn[n]，其他線程使用這個值，它可能還沒有被計算。在這種情況下，其他運行（除第一個之外）唯一的原因是獲得相同（正確）值的方式是由Xn和Yn指向的內存已被正確的值佔用，即使使用該同步問題應用程序返回正確的值。

在任何情況下，我通過拆分內核分爲二避免了同步的問題，就像我被Robert Crovella通過評論建議：

__global__ 
void DeltaE1(float *X, float *Y,float *Xn, float *Yn, float *step,float *DELTA, curandState *state, const int N,const int n){ 

    int tIdx = blockIdx.x*blockDim.x + threadIdx.x; 
    float x , y; 
    float rdmn1, rdmn2; 

    if(tIdx < N){ 
     DELTA[tIdx] = 0.0e0f; 
     if(tIdx == n){ 
      step[tIdx] = 0.2e0f; 
      rdmn1 = curand_uniform(&state[tIdx]); 
      rdmn2 = curand_uniform(&state[tIdx]); 

      Xn[tIdx] = X[tIdx] + step[tIdx]*(2.0e0f*rdmn1 - 1.0e0f); 
      Yn[tIdx] = Y[tIdx] + step[tIdx]*(2.0e0f*rdmn2 - 1.0e0f); 
      DELTA[tIdx] = - (X[tIdx]*X[tIdx] + Y[tIdx]*Y[tIdx]); 
      DELTA[tIdx] += Xn[tIdx]*Xn[tIdx] + Yn[tIdx]*Yn[tIdx]; 

     } 
     else{ 
      x = X[tIdx] - X[n]; 
      y = Y[tIdx] - Y[n]; 

      DELTA[tIdx] += -1.0e0f/sqrt(x*x + y*y); 
     }   
    } 
} 

__global__ 
void DeltaE2(float *X, float *Y,float *Xn, float *Yn,float *DELTA,const int N,const int n){ 

    int tIdx = blockIdx.x*blockDim.x + threadIdx.x; 
    int bIdx = blockIdx.x; 

    float x , y; 
    float dTotE = 0.0e0f; 
    if(tIdx < N){ 
     if(tIdx != n){ 
      x = X[tIdx] - Xn[n]; 
      y = Y[tIdx] - Yn[n]; 

      DELTA[tIdx] += 1.0e0f/sqrt(x*x + y*y); 

     } 
     dTotE = DELTA[tIdx]; 
     dTotE = ReduceSum2(dTotE); 
     if(threadIdx.x == 0)DELTA[bIdx] = dTotE; 

    } 

}