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

Question

这是this 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.66s 才能运行。

所以我想也许是数组访问时间太长了，所以我把那行改成了

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 ++;
            }
        }
    } 

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

所以两者的区别 5.15 秒和 0.853 秒不仅代表加法，而且实际上代表您的函数所做的所有有趣计算。

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

currentPixelValue = 8 * kernelCenterX * kernelCenterY * kernelCenterZ;

但它可能没那么聪明。