1

以我之前回答的问题为例：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，它不支持加倍。而且我使用无标志/无选项进行编译。

以 Np = 1k 为例。当我编译然后运行时，结果是 dE = 6.557993。当我运行第二、第三、第四次时，无论何时，它都是 dE = -0.3515406。有谁知道这是从哪里来的？

PS：我忘了提一下，可以在My Previous Question中找到的同一个内核 AvgDistance在 DeltaE 之前被调用。我不知道这是否有任何关系，但我认为值得一提。

PS2：nn 是任意选择的点（粒子）。

4

1 回答 1

2

正如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;

    }

}
于 2013-03-16T10:56:57.423 回答