CUDA C/C++: Same executable gives different results from the first run

Question

Take my previous answered question: My Previous Question, which, by the way, was properly answered by Robert Crovella.

I came up with another kernel that computes a random step to a point (by using the same RNG from my previous question) and calculates the difference of energy of that point with respect to its previous position (coordinates). This is the kernel:

__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;


    }
}

Then I do the final sum on the 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];
}

My kernel is launched with the following configuration:

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);

Where Np is the number of points (I used 1k - 4k). I have a GeForce 9500 GT, which has no support to double. And I compile using no flag/no option.

Take Np = 1k, for example. When I compile and then run, the result is dE = 6.557993. When I run a second, third, fourth, whatever time, it is dE = -0.3515406. Does anyone know where this come from?

P.S.: I forgot to mention, the same kernel AvgDistance that can be found at My Previous Question is called right before DeltaE. I don't know if this has anything to do, but I thought it was worth to mention.

P.S.2: nn is any chosen point(particle).

Answer 1

As was pointed out by Robert Crovella via comment above, what probably was happening is that while tIdx = n was computing Xn[n] and Yn[n], other threads was using this value and it may not be computed yet. In that case, the only reason to the other runs (other then the first one) to get ways the same (correct) value is that the memory pointed to by Xn and Yn is already occupied with the right value, and even with that synchronization problem the application returns the right value.

In any case, I avoided the synchronization problem by splitting the kernel into two, just as I was advised by Robert Crovella via comment:

__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;

    }

}