C 中 3D 直接卷积实现的优化

Question

对于我的项目，我编写了直接 3D 卷积的简单 C 实现，并在输入上定期填充。不幸的是，由于我是 C 新手，所以性能不是很好……这是代码：

int mod(int a, int b)
{
    // calculate mod to get the correct index with periodic padding
    int r = a % b;
    return r < 0 ? r + b : r;
}
void convolve3D(const double *image, const double *kernel, const int imageDimX, const int imageDimY, const int imageDimZ, const int stencilDimX, const int stencilDimY, const int stencilDimZ, 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 imageIndex = 0, kernelIndex = 0;
    
    // Loop through each voxel
    for (k = 0; k < imageDimZ; k++){
        for ( j = 0; j < imageDimY; j++) {
            for ( i = 0; i < imageDimX; i++) {
                stencilIndex = 0;
                // for each voxel, loop through each kernel coefficient
                for (n = 0; n < kernelDimZ; n++){
                    for ( m = 0; m < kernelDimY; m++) {
                        for ( l = 0; l < kernelDimX; l++) {
                            // find the index of the corresponding voxel in the output image
                            xShift = l - kernelCenterX;
                            yShift = m - kernelCenterY;
                            zShift = n - kernelCenterZ;

                            outI = mod ((i - xShift), imageDimX);
                            outJ = mod ((j - yShift), imageDimY);
                            outK = mod ((k - zShift), imageDimZ);
                            
                            outIndex = outK * imageDimX * imageDimY + outJ * imageDimX + outI;

                            // calculate and add
                            result[outIndex] += stencil[stencilIndex]* image[imageIndex];
                            stencilIndex++;
                        }
                    }
                } 
                imageIndex ++;
            }
        }
    } 
}

• 按照惯例，所有矩阵（图像、内核、结果）都以列优先方式存储，这就是为什么我以这种方式循环它们，以便它们在内存中更接近（听说这会有所帮助）。

我知道实现非常幼稚，但由于它是用 C 编写的，我希望性能会很好，但结果有点令人失望。我用大小为 100^3 的图像和大小为 10^3 的内核对其进行了测试（如果只计算乘法和加法，总共 ~1GFLOPS），它花了 ~7 秒，我认为这远低于典型 CPU 的能力。

如果可能的话，你们能帮我优化这个程序吗？ 我愿意接受任何可以提供帮助的东西，如果您可以考虑的话，请提供以下几点：

1. 我正在处理的问题可能很大（例如，大小为 200 x 200 x 200 的图像，内核大小为 50 x 50 x 50 甚至更大）。我知道优化此问题的一种方法是将这个问题转换为矩阵乘法问题并使用 blas GEMM 例程，但恐怕内存无法容纳如此大的矩阵

2. 由于问题的性质，我更喜欢直接卷积而不是 FFTConvolve，因为我的模型是在考虑直接卷积的情况下开发的，我对 FFT 卷积的印象是它给出的结果与直接卷积略有不同，特别是对于快速改变形象，我试图避免的差异。 也就是说，我绝不是这方面的专家。因此，如果您有一个基于 FFTconvolve 的出色实现和/或我对 FFT convolve 的印象完全有偏见，如果您能帮助我，我将不胜感激。

3. 假设输入图像是周期性的，因此需要周期性填充

4. 我知道使用 blas/SIMD 或其他较低级别的方法在这里肯定会有很大帮助。但由于我是这里的新手，我真的不知道从哪里开始......如果你有这些图书馆的经验，如果你能帮助我指出正确的方向，我将不胜感激，

非常感谢您的帮助，如果您需要有关问题性质的更多信息，请告诉我

Answer 1

第一步，将您的 mod ((i - xShift), imageDimX) 替换为以下内容：

inline int clamp( int x, int size )
{
    if( x < 0 ) return x + size;
    if( x >= size ) return x - size;
    return x;
}

这些分支是非常可预测的，因为它们对大量连续元素产生相同的结果。整数取模比较慢。

现在，下一步（按成本/利润排序）将是并行化。如果您有任何现代 C++ 编译器，只需在项目设置中的某处启用 OpenMP。之后，您需要进行 2 次更改。

1. 用这样的东西装饰你的最外层循环：#pragma omp parallel for schedule(guided)
2. 在该循环中移​​动您的函数级变量。这也意味着每次迭代都必须从 k 计算初始 imageIndex。

下一个选项，重新编写您的代码，以便您只编写每个输出值一次。在最里面的 3 个循环中计算最终值，从图像和内核的随机位置读取，并且只写入一次结果。当您在内部循环中有result[outIndex] += 时，CPU 会停止等待来自内存的数据。当您在一个寄存器而不是内存的变量中累积时，就没有访问延迟。

SIMD 是最复杂的优化。但简而言之，您需要硬件具有的 FMA 的最大宽度（如果您有 AVX 并且需要双精度，则宽度为 4），并且您还需要多个独立的累加器用于 3 个最内层循环，以避免撞到延迟而不是使吞吐量饱和。这里以my answer to much easier problem 为例说明我的意思。