为什么 C 中的 += 花费这么多时间，以及如何提高性能？

Question

这是这个问题的后续，我试图在其中实现具有周期性边界条件的 3D 直接卷积。

在社区的帮助下，我能够将性能提高约 50%（非常感谢社区）：

int mod(int a, int b)
{
    if (a<0)
        return a + b;
    else if (a >= b)
        return a - b;
    else
        return a;
}

void convolve(const double *image, const double *kernel, const int imageDimX, const int imageDimY, 
const int imageDimZ, const int kernelDimX, const int kernelDimY, const int kernelDimZ, double *result)
{
    int imageSize = imageDimX * imageDimY * imageDimZ;
    int kernelSize = kernelDimX * kernelDimY * kernelDimZ;

    int i, j, k, l, m, n;
    int kernelCenterX = (kernelDimX - 1) / 2;
    int kernelCenterY = (kernelDimY - 1) / 2;
    int kernelCenterZ = (kernelDimZ - 1) / 2;
    int xShift,yShift,zShift;
    int outIndex, outI, outJ, outK;
    int outputIndex = 0, imageIndex = 0, kernelIndex = 0;
    double currentPixelValue;

    for (k = 0; k < imageDimZ; k++){
        for ( j = 0; j < imageDimY; j++) {
            for ( i = 0; i < imageDimX; i++) {
                currentPixelValue = 0.0;

                kernelIndex = 0;
                for (n = 0, zShift = - kernelCenterZ; n < kernelDimZ; n++, zShift++){
                    outK = mod ((k - zShift), imageDimZ);
                    for ( m = 0, yShift = - kernelCenterY; m < kernelDimY; m++, yShift++) {
                        outJ = mod ((j - yShift), imageDimY);
                        for ( l = 0, xShift = - kernelCenterX; l < kernelDimX; l++, xShift++) {
                            outI = mod ((i - xShift), imageDimX);
                            
                            imageIndex = outK * imageDimX * imageDimY + outJ * imageDimX + outI;
                            
                            // This mysterious line
                            currentPixelValue += kernel[kernelIndex]* image[imageIndex]; 

                            kernelIndex++;
                        }
                    }
                }

                result[outputIndex] = currentPixelValue;
                outputIndex ++;
            }
        }
    } 
}

作为参考，对于我的测试用例，这需要大约 5.65 秒才能运行。

出于好奇，我试图准确指出哪个操作是瓶颈，结果发现它是最内层循环中的神秘线。

通过删除该行，运行需要0.66 秒。

所以我想也许是访问数组的时间太长了，所以我将该行更改为

currentPixelValue += 1.0;

但运行时性能仅提高到5.185s，这绝对出乎我的意料，

所以我试图将 += 更改为 = 就像测试一样：

currentPixelValue = 1.0;

性能大幅提升至0.853s

很明显，瓶颈是+=操作，这对我来说再次非常有趣。

仅仅访问一个变量的值并向它添加一个常量是如何成为算法的瓶颈？你们能否帮助我提供一些见解，并希望进一步提高性能？

编辑：

作为另一个比较案例，我尝试将行更改为

currentPixelValue = stencil[stencilIndex]* image[imageIndex];

花了~5.15s运行

我正在努力理解这一点，我认为测试表明任何类型的值访问都会成为算法的瓶颈。但是，它正上方的行，也在最内层的循环中，也有值访问，似乎没有瓶颈......

这对我来说非常神秘和有趣，哈哈

编译信息：

CC=mpicc
CFLAGS = -O3 -Wall -g -std=gnu99

Answer 1

通过替换currentPixelValue += ...为currentPixelValue = 1您使大部分其余函数变得不必要，从而允许编译器对其进行优化。

必须阅读生成的程序集才能确定，但​​我们可以这样推理：

• 假设您替换currentPixelValue += ...为currentPixelValue = 1.

• 现在imageIndex从不使用上一行计算的值。所以编译器会删除那个计算。

• 现在outI, outJ, outK也不再使用，并且由于编译器可以内联mod函数，它知道它没有副作用，因此也可以优化这些计算。

• kernelIndex 现在也从未使用过，因此也将其删除。

你的i循环现在看起来像这样：

                currentPixelValue = 0.0;

                kernelIndex = 0;
                for (n = 0, zShift = - kernelCenterZ; n < kernelDimZ; n++, zShift++){
                    for ( m = 0, yShift = - kernelCenterY; m < kernelDimY; m++, yShift++) {
                        for ( l = 0, xShift = - kernelCenterX; l < kernelDimX; l++, xShift++) {
                            currentPixelValue = 1;
                        }
                    }
                }

                result[outputIndex] = currentPixelValue;
                outputIndex ++;

换句话说，在n,m,l循环什么也不做，除了设置currentPixelValue来1一遍又一遍。编译器知道这是没有意义的，所以它只需要做一次。这也使得初始化变得currentPixelValue = 0.0不必要，并且kernelIndex根本不使用。因此，我们只剩下：

    for (k = 0; k < imageDimZ; k++){
        for ( j = 0; j < imageDimY; j++) {
            for ( i = 0; i < imageDimX; i++) {
                currentPixelValue = 1;
                result[outputIndex] = currentPixelValue;
                outputIndex ++;
            }
        }
    } 

也就是说，该函数现在除了result用 1填充矩阵之外什么都不做，并且编译器可以以非常有效的方式做到这一点。

因此，5.15 秒和 0.853 秒之间的差异不仅代表加法，而且代表您的函数所做的几乎所有有趣的计算。

（如果您改为将行更改为 ，currentPixelValue += 1那么编译器仍然必须运行所有n,m,l循环以将其增加正确的次数。如果它足够聪明，它可以将循环替换为

currentPixelValue = 8 * kernelCenterX * kernelCenterY * kernelCenterZ;

但它可能没有那么聪明。