C语言中三维直接卷积实现的优化

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的能力。

如果可能的话，你们能帮我优化这个程序吗？我愿意接受任何有帮助的东西，如果你能考虑的话，只有几件事：

>

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

由于问题的性质，我更喜欢直接卷积，而不是FFTConvolve，因为我的模型是在考虑直接卷积的情况下开发的，我对FFT卷积的印象是，它给出的结果与直接卷积略有不同，尤其是对于快速变化的图像，这是我试图避免的差异。也就是说，我根本不是这方面的专家。因此，如果您有一个基于FFTconvolve的很棒的实现，并且/或者我对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;
}

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

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

1. 用这样的东西装饰你的外环：\pragma omp parallel for schedule（guided）

下一个选项是重新编写代码，这样每个输出值只写一次。在最里面的3个循环中计算最终值，从映像和内核的随机位置读取，只写一次结果。当内部循环中有这个结果时，CPU会暂停等待内存中的数据。当您在寄存器而非内存的变量中累积时，就不会有访问延迟。

SIMD是这方面最复杂的优化。但简而言之，您需要硬件具有的FMA的最大宽度（如果您有AVX并且需要双精度，则宽度为4），并且您还需要为您的3个最里面的循环使用多个独立的累加器，以避免达到延迟，而不是使吞吐量饱和。这是我对更简单问题的回答，作为我的意思的一个例子。