CUDA C/C++：计算每点距离的平均值(相互作用能量,也许？)

我能够重现你的结果.

首先,我将sinf,cosf和sqrtf的所有实例切换为相应的双版本.这对结果没有影响.

接下来我包含了一个typedef,所以我可以很容易地将精度从float切换到double和back,用mytype取代我的typedef中代码中每个相关的float实例.

当我使用float的typedef运行代码并且数据大小为4096时,我得到以下结果：

GPU average = 108.294922
CPU average = 109.925285

当我运行typedef为double且数据大小为4096的代码时,我得到以下结果：

GPU average = 108.294903
CPU average = 108.294903

当我使用float的typedef运行代码并且数据大小为8192时,我得到以下结果：

GPU average = 153.447327
CPU average = 131.473526

当我运行typedef为double且数据大小为8192的代码时,我得到以下结果：

GPU average = 153.447380
CPU average = 153.447380

至少有2个观察结果：

> GPU结果在float和double之间不变,除了小数点后5位
>在float和double之间CPU结果变化1-20％左右,但是当选择double时,它们与GPU结果完全对齐(到第6个小数位).

基于此,我相信CPU正在提供可变的,可疑的行为.

这是我的代码供参考：

#include <stdio.h>
#include <curand.h>
#include <curand_kernel.h>
#define DSIZE 8192
#define BlockSize 32
#define pi 3.14159f


#define cudaCheckErrors(msg) 
    do { 
        cudaError_t __err = cudaGetLastError(); 
        if (__err != cudaSuccess) { 
            fprintf(stderr,"Fatal error: %s (%s at %s:%d)n",
                msg,cudaGetErrorString(__err),
                __FILE__,__LINE__); 
            fprintf(stderr,"*** FAILED - ABORTINGn"); 
            exit(1); 
        } 
    } while (0)

typedef double mytype;

__global__ void InitRNG ( curandState * state,&state[tIdx]);
}

__global__
void SortPoints(mytype* X,mytype* Y,curandState *state){

    mytype rdmn1,rdmn2;

    unsigned int tIdx = blockIdx.x*blockDim.x + threadIdx.x;
    mytype range;
    if(tIdx < N){
        rdmn1 = curand_uniform(&state[tIdx]);
        rdmn2 = curand_uniform(&state[tIdx]);
        range = sqrt(0.25e0f*N*rdmn1);
            X[tIdx] = range*cos(2.0e0f*pi*rdmn2);
            Y[tIdx] = range*sin(2.0e0f*pi*rdmn2);
    }
}

__device__
mytype ReduceSum2(mytype In){

    __shared__ mytype data[BlockSize];

    unsigned int tIdx = threadIdx.x;

    data[tIdx] = In;
    __syncthreads();

    for(unsigned int i = blockDim.x/2; i > 0; i >>= 1){

        if(tIdx < i){

            data[tIdx] += data[tIdx + i];
        }

        __syncthreads();
    }

    return data[0];

}

__global__
void AvgDistance(mytype *X,mytype *Y,mytype *Avg,const int N){

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

    mytype x,y;
    mytype d = 0.0f;
    if(tIdx < N){

        for(int i = tIdx + 1; i < N ; i++){

            x = X[tIdx] - X[i];
            y = Y[tIdx] - Y[i];

            d += 1.0e0f/(sqrt(x*x + y*y));
        }
        __syncthreads();
        Avg[bIdx] = ReduceSum2(d);

    }
}

mytype cpu_avg(const mytype *rx,const mytype *ry,const int size){
  mytype average = 0.0f;
  for(int i = 0; i < size; i++){
    for(int j = i + 1; j < size; j++){
        average += 1.0e0f/sqrt((rx[i]-rx[j])*(rx[i]-rx[j]) + (ry[i]-ry[j])*(ry[i]-ry[j]));
    }
  }
  average = average/(mytype)size;
  return average;
}

int main() {

  int Np = DSIZE;
  mytype *rx,*ry,*d_rx,*d_ry,*d_Avg,*Avg;
  curandState *d_state;
  int seed = 1;

  dim3 threads(BlockSize,BlockSize);
  dim3 blocks((int)ceilf(Np/(float)threads.x),(int)ceilf(Np/(float)threads.y));
  printf("number of blocks = %dn",blocks.x);
  printf("number of threads= %dn",threads.x);
  rx = (mytype *)malloc(DSIZE*sizeof(mytype));
  if (rx == 0) {printf("malloc failn"); return 1;}
  ry = (mytype *)malloc(DSIZE*sizeof(mytype));
  if (ry == 0) {printf("malloc failn"); return 1;}

  cudaMalloc((void**)&d_rx,DSIZE * sizeof(mytype));
  cudaMalloc((void**)&d_ry,DSIZE * sizeof(mytype));
  cudaMalloc((void**)&d_Avg,blocks.x * sizeof(mytype));
  cudaMalloc((void**)&d_state,DSIZE * sizeof(curandState));
  cudaCheckErrors("cudamalloc");



  InitRNG<<<blocks.x,seed);
  SortPoints<<<blocks.x,d_state);
  AvgDistance<<<blocks.x,threads.x*sizeof(mytype)>>>(d_rx,Np);
  cudaCheckErrors("kernels");


  Avg = new mytype[blocks.x];

  cudaMemcpy(Avg,blocks.x*sizeof(mytype),cudaMemcpyDeviceToHost);
  cudaMemcpy(rx,d_rx,DSIZE*sizeof(mytype),cudaMemcpyDeviceToHost);
  cudaMemcpy(ry,cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudamemcpy");
  mytype average = 0;

  for(int i = 0; i < blocks.x; i++){
    average += Avg[i];
  }
  average = average/(mytype)Np;

  printf("GPU average = %fn",average);
  average = cpu_avg(rx,ry,DSIZE);
  printf("CPU average = %fn",average);

  return 0;
}

我在RHEL 5.5,CUDA 5.0,Intel Xeon X5560上运行

编译：

nvcc -O3 -arch=sm_20 -lcurand -lm -o t93 t93.cu

编辑：
在观察到CPU的可变性之后,我发现通过修改CPU平均代码可以消除大部分CPU的可变性,如下所示：

mytype cpu_avg(const mytype *rx,const int size){
  mytype average = 0.0f;
  mytype temp = 0.0f;
  for(int i = 0; i < size; i++){
    for(int j = i + 1; j < size; j++){
        temp += 1.0e0f/sqrt((rx[i]-rx[j])*(rx[i]-rx[j]) + (ry[i]-ry[j])*(ry[i]-ry[j]));
    }
    average += temp/(mytype)size;
    temp = 0.0f;
  }
  return average;
}

所以我想说CPU端的中间结果有问题.有趣的是,它没有显示在GPU结果上.我怀疑其原因是GPU平均值的最终总和是在CPU上完成的(因此每个单独的GPU块结果按大小缩小,例如8192),并且这些可能具有足以存活的中间精度.最后的师.如果您内联CPU平均值计算,您可能会再次观察到不同的内容.