主机和GPU上进行CUDA加法的结果不同

Question

我有一个函数，它接收一张彩色图片并返回它的灰度版本。 如果我在主机上运行顺序代码，一切都完美。但是如果我在设备上运行它，则结果会略有不同（每1000个像素中有一个与正确值相比要么是+1，要么是-1）。

我认为这与转换有关，但我不能确定。这是我使用的代码：

    __global__ void rgb2gray_d (unsigned char *deviceImage, unsigned char *deviceResult, const int height, const int width){
    /* calculate the global thread id*/
    int threadsPerBlock  = blockDim.x * blockDim.y;
    int threadNumInBlock = threadIdx.x + blockDim.x * threadIdx.y;
    int blockNumInGrid   = blockIdx.x  + gridDim.x  * blockIdx.y;

    int globalThreadNum = blockNumInGrid * threadsPerBlock + threadNumInBlock;
    int i = globalThreadNum;

    float grayPix = 0.0f;
    float r = static_cast< float >(deviceImage[i]);
    float g = static_cast< float >(deviceImage[(width * height) + i]);
    float b = static_cast< float >(deviceImage[(2 * width * height) + i]);
    grayPix = (0.3f * r) + (0.59f * g) + (0.11f * b);

    deviceResult[i] = static_cast< unsigned char > (grayPix);
}

void rgb2gray(unsigned char *inputImage, unsigned char *grayImage, const int width, const int height, NSTimer &timer) {

    unsigned char *deviceImage;
    unsigned char *deviceResult;

    int initialBytes = width * height * 3;  
    int endBytes =  width * height * sizeof(unsigned char);

    unsigned char grayImageSeq[endBytes];

    cudaMalloc((void**) &deviceImage, initialBytes);
    cudaMalloc((void**) &deviceResult, endBytes);
    cudaMemset(deviceResult, 0, endBytes);
    cudaMemset(deviceImage, 0, initialBytes);

    cudaError_t err = cudaMemcpy(deviceImage, inputImage, initialBytes, cudaMemcpyHostToDevice);    

    // Convert the input image to grayscale 
    rgb2gray_d<<<width * height / 256, 256>>>(deviceImage, deviceResult, height, width);
    cudaDeviceSynchronize();

    cudaMemcpy(grayImage, deviceResult, endBytes, cudaMemcpyDeviceToHost);

    ////// Sequential
    for ( int y = 0; y < height; y++ ) {
             for ( int x = 0; x < width; x++ ) {
                   float grayPix = 0.0f;
                   float r = static_cast< float >(inputImage[(y * width) + x]);
                   float g = static_cast< float >(inputImage[(width * height) + (y * width) + x]);
                   float b = static_cast< float >(inputImage[(2 * width * height) + (y * width) + x]);

                   grayPix = (0.3f * r) + (0.59f * g) + (0.11f * b);
                   grayImageSeq[(y * width) + x] = static_cast< unsigned char > (grayPix);
              }
        }

    //compare sequential and cuda and print pixels that are wrong
    for (int i = 0; i < endBytes; i++)
    {
        if (grayImage[i] != grayImageSeq[i])
        cout << i << "-" << static_cast< unsigned int >(grayImage[i]) <<
                 " should be " << static_cast< unsigned int >(grayImageSeq[i]) << endl;
        }

    cudaFree(deviceImage);
    cudaFree(deviceResult);
}

我提到我为初始图像宽度*高度*3分配内存，因为初始图像是CImg。

我使用GeForce GTX 480进行工作。

- andreea.sandu

---

“稍微不同”是什么意思？ - talonmies

这意味着大约每1000个像素中有一个与预期值相比要么是+1，要么是-1。 - andreea.sandu

Answer 1

最终我找到了答案。CUDA在单精度和双精度中都自动执行融合乘加操作。通过使用下面的文档1，第4.4节，我成功地解决了它。与其进行...

grayPix = (0.3f * r) + (0.59f * g) + (0.11f * b);

我现在正在做。

grayPix = __fadd_rn(__fadd_rn(__fmul_rn(0.3f, r),__fmul_rn(0.59f, g)), __fmul_rn(0.11f, b));

这禁用了将多个乘法和加法合并成融合乘加指令的功能。

参考链接：NVIDIA GPU 的浮点数和 IEEE 754 兼容性

- andreea.sandu

---

非常感谢您的回答，现在修饰符已经帮助我解决了CPU和GPU结果不同的问题。 - teejay

Answer 2

浮点数运算在设备代码和主机代码中可能会产生略微不同的结果。

这种情况有多种可能性。您必须考虑到这两个函数由两个不同的编译器编译成两个不同的二进制程序，在两个不同的浮点硬件实现上运行。

例如，如果浮点计算以不同的顺序执行，舍入误差可能会导致不同的结果。

此外，当在x86架构CPU上使用32位（float）或64位（double）浮点表示运行浮点计算时，浮点数学是由FPU单元完成的，该单元在内部使用80位精度，然后将结果截断回32位浮点数据类型或64位双精度数据类型。

GPU的ALU在浮点数学中使用32位精度（假设您使用float数据类型）。

关于浮点表示和算术的优秀文章可以在这里找到。

- RoBiK

---

@andreea.sandu 在图像处理中，这样微小的东西不应该有影响。如果你真的需要“正确”的结果，在临时值上使用双精度（无论是在 CPU 还是 GPU 上）。 - Pavan Yalamanchili

期望答案完全准确并不是正确的。如果将x86 FPU与SSE生成答案进行比较，您可能会遇到相同的问题。有几件事情可以尝试。(1)在GPU上禁用FMAs。(2)在CPU端生成3个256条目查找表，并将乘法替换为查找(~每个查找8种分歧)。(3)在主机上生成一个2^24条目的查找表，并在CPU上使用它。 - Greg Smith

我同意Greg的观点，除非你能够使用相同的精度、确保使用相同的浮点标准和相同的操作顺序，否则你不能期望总是得到完全相同的结果。 - RoBiK

我也在想为什么您要使用浮点数...... 您可以使用整数运算计算直方图，结果保证在GPU和CPU之间相同。您可以使用整数变量甚至int16。按照以下方式计算：grayPix = ((30 * r) + (59 * g) + (11 * b) + 50)/100; - RoBiK

最可能导致差异的原因是主机上的“grayPix”计算使用高于单精度的精度进行。为了防止这种情况，您需要在主机代码中声明一些“volatile float”变量，然后将整个计算分解为一系列具有两个输入的单独操作，并将每个操作的结果发送到一个volatile float目标。例如：{ volatile float t0,t1,t2; t0=0.3fr; t1=0.59fg; t2=0.11f*b; t0=t0+t1; t0=t0+t2; grayPix=t0; } 当我设计旨在在GPU上运行的代码时，我经常使用这种技术。效果非常好。 - njuffa