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），它花费了 ~7s，我认为这远低于典型 CPU.

如果可以的话，能不能帮我优化一下这个套路？我乐于接受任何可以提供帮助的事情，如果您可以考虑以下几点：

我正在处理的问题可能很大（例如，大小为 200 x 200 x 200 的图像，内核大小为 50 x 50 x 50 甚至更大）。我知道优化这个的一种方法是将这个问题转换为矩阵乘法问题并使用 blas GEMM 例程，但恐怕内存无法容纳这么大的矩阵
由于问题的性质，我更喜欢直接卷积而不是FFTConvolve，因为我的模型是在考虑直接卷积的情况下开发的，而我对FFT卷积的印象是它给出的略有不同结果比直接卷积特别是对于快速变化的图像，这是我试图避免的差异。也就是说，我绝不是这方面的专家。所以如果你有一个很好的基于 FFTconvolve 的实现 and/or 我对 FFT convolve 的印象是完全有偏见的，如果你能帮助我，我将不胜感激。
假定输入图像是周期性的，因此需要周期性填充
我知道利用 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;
}

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

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

用这样的东西装饰你的最外层循环：#pragma omp parallel for schedule(guided)
在该循环内移动您的函数级变量。这也意味着对于每次迭代，您必须根据 k 计算初始 imageIndex。

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

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

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

Optimization of 3D Direct Convolution Implementation in C

c

optimization

simd

convolution

blas