cuda 共享内存 - 结果不一致

我正在尝试并行缩减以对 CUDA 中的数组求和。目前我传递一个数组来存储每个块中元素的总和。这是我的代码：

#include <cstdlib>
#include <iostream>
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <helper_cuda.h>
#include <host_config.h>
#define THREADS_PER_BLOCK 256
#define CUDA_ERROR_CHECK(ans) { gpuAssert((ans), __FILE__, __LINE__); }

using namespace std;

inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
   if (code != cudaSuccess) 
   {
      fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}

struct double3c {
    double x; 
    double y;
    double z;
    __host__ __device__ double3c() : x(0), y(0), z(0) {}
    __host__ __device__ double3c(int x_, int y_, int z_) : x(x_), y(y_), z(z_) {}
    __host__ __device__ double3c& operator+=(const double3c& rhs) { x += rhs.x; y += rhs.y; z += rhs.z;}
    __host__ __device__ double3c& operator/=(const double& rhs) { x /= rhs; y /= rhs; z /= rhs;}

};

class VectorField {
public:
    double3c *data;
    int size_x, size_y, size_z;
    bool is_copy;  

    __host__ VectorField () {}

    __host__ VectorField (int x, int y, int z) {
        size_x = x; size_y = y; size_z = z;
        is_copy = false;
        CUDA_ERROR_CHECK (cudaMalloc(&data, x * y * z * sizeof(double3c))); 
    }

    __host__ VectorField (const VectorField& other) {
        size_x = other.size_x; size_y = other.size_y; size_z = other.size_z;
        this->data = other.data;
        is_copy = true;
    }

    __host__ ~VectorField() {     
        if (!is_copy) CUDA_ERROR_CHECK (cudaFree(data));
    }
};

__global__ void KernelCalculateMeanFieldBlock (VectorField m, double3c* result) {
    __shared__ double3c blockmean[THREADS_PER_BLOCK];    
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if (index < m.size_x * m.size_y * m.size_z) blockmean[threadIdx.x] = m.data[index] = double3c(0, 1, 0);
    else blockmean[threadIdx.x] = double3c(0,0,0);
    __syncthreads();
    for(int s = THREADS_PER_BLOCK / 2; s > 0; s /= 2) {
        if (threadIdx.x < s) blockmean[threadIdx.x] += blockmean[threadIdx.x + s];
        __syncthreads();
    }


    if(threadIdx.x == 0) result[blockIdx.x] = blockmean[0];   
}

double3c CalculateMeanField (VectorField& m) { 
    int blocknum = (m.size_x * m.size_y * m.size_z - 1) / THREADS_PER_BLOCK + 1;
    double3c *mean = new double3c[blocknum]();
    double3c *cu_mean;
    CUDA_ERROR_CHECK (cudaMalloc(&cu_mean, sizeof(double3c) * blocknum));
    CUDA_ERROR_CHECK (cudaMemset (cu_mean, 0, sizeof(double3c) * blocknum));

        KernelCalculateMeanFieldBlock <<<blocknum, THREADS_PER_BLOCK>>> (m, cu_mean);
        CUDA_ERROR_CHECK (cudaPeekAtLastError());
        CUDA_ERROR_CHECK (cudaDeviceSynchronize());
        CUDA_ERROR_CHECK (cudaMemcpy(mean, cu_mean, sizeof(double3c) * blocknum, cudaMemcpyDeviceToHost));

    CUDA_ERROR_CHECK (cudaFree(cu_mean));
    for (int i = 1; i < blocknum; i++) {mean[0] += mean[i];}
    mean[0] /= m.size_x * m.size_y * m.size_z;
    double3c aux = mean[0];
    delete[] mean;
    return aux;
}



int main() {
    VectorField m(100,100,100);
    double3c sum = CalculateMeanField (m);
    cout <<  sum.x << '\t' << sum.y << '\t' <<sum.z;  


    return 0;
}

EDIT

贴出功能代码。构建一个VectorField使用 10x10x10 元素可以正常工作并给出平均值 1，但是使用 100x100x100 元素构建它可以得到平均值 ~0.97（每次运行都会有所不同）。这是进行并行减少的正确方法，还是应该坚持每个块启动一个内核？

当我在 Linux 上编译现在的代码时，我收到以下警告：

t614.cu(55): warning: __shared__ memory variable with non-empty constructor or destructor (potential race between threads)

不应忽视此类警告。它与这行代码相关：

__shared__ double3c blockmean[THREADS_PER_BLOCK]; 

由于存储在共享内存中的这些对象的初始化（通过构造函数）将以某种任意顺序发生，并且在该初始化与也将设置这些值的后续代码之间没有障碍，不可预测的事情（*） 可以发生。

如果我插入一个__syncthreads()在将构造函数活动与后续代码隔离的代码中，我得到了预期的结果：

__shared__ double3c blockmean[THREADS_PER_BLOCK];    
int index = threadIdx.x + blockIdx.x * blockDim.x;
__syncthreads();  // add this line
if (index < m.size_x * m.size_y * m.size_z) blockmean[threadIdx.x] = m.data[index] = double3c(0, 1, 0);
else blockmean[threadIdx.x] = double3c(0,0,0);
__syncthreads();

然而，这仍然给我们留下了警告。解决此问题并使警告消失的修改是分配必要的__shared__动态大小。将共享内存声明更改为：

extern __shared__ double3c blockmean[];

并修改你的内核调用：

KernelCalculateMeanFieldBlock <<<blocknum, THREADS_PER_BLOCK, THREADS_PER_BLOCK*sizeof(double3c)>>> (m, cu_mean);

这将消除警告，产生正确的结果，并避免共享内存变量上不必要的构造函数流量。 （以及额外的__syncthreads()上面的描述就不再需要了。）

*关于“不可预测的事情”，如果您通过检查生成的 SASS (库对象转储 http://docs.nvidia.com/cuda/cuda-binary-utilities/index.html#abstract-萨斯...）或 PTX(**) (nvcc -ptx ...)，你会看到每个线程初始化entire __shared__对象数组为零（默认构造函数的行为）。因此，一些线程（即扭曲）可以提前竞争并开始根据以下行填充共享内存区域：

if (index < m.size_x * m.size_y * m.size_z) blockmean[threadIdx.x] = m.data[index] = double3c(0, 1, 0);

然后，当其他扭曲开始执行时，这些线程将清除entire再次共享内存阵列。这种赛车行为会导致不可预测的结果。

** 我通常不建议通过检查 PTX 来判断代码行为，但在本例中它同样具有指导意义。最后的编译阶段不会优化构造函数的行为。