在 CUDA 中找到线程大小的总和减少问题

Question

在之前的帖子here 中，我询问了如何通过归约计算数组的总和。现在我有一个新问题，图像更大，我的结果不正确，每次运行它都会改变。 我用96*96的图片大小测试array sample

第一次结果：28169.046875

第二次结果：28169.048828

预期结果：28169.031250

这是我的代码：

#include <stdio.h>
#include <cuda.h>

__global__  void calculate_threshold_kernel(float * input, float * output)
{

   int idx = blockIdx.x * blockDim.x + threadIdx.x;
   int t = threadIdx.x;
   __shared__ float partialSum[256];
   partialSum[t] = input[idx];
   __syncthreads();

   for (int stride = 1; stride < blockDim.x; stride *= 2)
   {
       if (t % (2 * stride) == 0)
           partialSum[t] += partialSum[t + stride];

       __syncthreads();
   }

   if (t == 0)
   {
       atomicAdd(output,partialSum[0]);
   }

}
int main( void )
{

    float *d_array, *d_output,*h_input, *h_output;
    int img_height = 96;
    int img_width = 96;
    int input_elements = img_height * img_width;

    h_input = (float*) malloc(sizeof(float) * input_elements);
    cudaMalloc((void**)&d_output, sizeof(float));
    cudaMemset(d_output, 0, sizeof(float));
    h_output = (float*)malloc(sizeof(float));
    cudaMalloc((void**)&d_array, input_elements*sizeof(float));

    float array[] = {[array sample]};
    for (int i = 0; i < input_elements; i++)
    {
        h_input[i] = array[i];
    }

    cudaMemcpy(d_array, h_input, input_elements*sizeof(float), cudaMemcpyHostToDevice);

        dim3 blocksize(256);
        dim3 gridsize(input_elements/blocksize.x);

        calculate_threshold_kernel<<<gridsize,blocksize>>>(d_array, d_output);

        cudaMemcpy(h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);

        printf("Sum from GPU = %f\n", *h_output);

    return 0;
}

Answer 1

虽然康世印关于浮点精度和浮点算术不可交换的回答是正确的，但他对每次运行结果不同的原因并不正确。

• 浮点运算是不可交换的，这意味着以不同顺序执行的操作可以返回不同的结果。例如，对于 a、b、c 和 d 的某些值，(((a+b)+c)+d) 可能与 ((a+b)+(c+d)) 略有不同。但是这两个结果都不应该因运行而异。

• 您的结果在不同的运行之间会有所不同，因为atomicAdd 导致添加顺序不同。使用 double 也不能保证不同运行之间的结果相同。

• 有一些方法可以在没有 atomicAdd 作为最后一步的情况下实现并行缩减（例如：使用第二个内核启动来添加第一次启动的部分总和），这可以提供一致的（但与 CPU 略有不同）的结果。

Answer 2

float 的精度有限，最高可达 7 个小数位，如此处所述。

https://en.wikipedia.org/wiki/Floating_point#Accuracy_problems

结果会发生变化，因为float 上的操作是不可交换的，并且您正在使用并行归约。

结果会发生变化，因为float 上的操作是不可交换的，并且您使用的是atomicAdd()，它无法保持加法的顺序。

如果您想要更准确的结果，您可以改用double。