为何我的AVX2二维卷积实现与标量实现速度相当？

问题概述 #

使用AVX2手动优化3x3二维卷积函数后，性能与标量实现完全一致：在8192×8192输入矩阵、基准测试迭代次数iterations=10的场景下，两者平均耗时均为750ms。编译参数为g++ -O3与-march=native，Intel VTune分析显示SIMD实现的耗时集中在_mm256_add_pd、_mm256_mul_pd指令及line_x_idx循环内，而非内存加载/存储操作。

标量实现代码 #

void kernel(int iterations, int N, double in_mat[N][N], double out_mat[N][N], double conv_mask[3][3])
{  
    for (int i = 0; i < iterations; i++)
    {
        for (int mat_y_idx = 1; mat_y_idx < N - 1; ++mat_y_idx)
        {
            for (int mat_x_idx = 1; mat_x_idx < N - 1; ++mat_x_idx)
            {
                double convolve_value_0_0 = in_mat[mat_y_idx - 1][mat_x_idx - 1] * conv_mask[0][0];
                double convolve_value_0_1 = in_mat[mat_y_idx - 1][mat_x_idx + 0] * conv_mask[0][1];
                double convolve_value_0_2 = in_mat[mat_y_idx - 1][mat_x_idx + 1] * conv_mask[0][2];

                double convolve_value_1_0 = in_mat[mat_y_idx + 0][mat_x_idx - 1] * conv_mask[1][0];
                double convolve_value_1_1 = in_mat[mat_y_idx + 0][mat_x_idx + 0] * conv_mask[1][1];
                double convolve_value_1_2 = in_mat[mat_y_idx + 0][mat_x_idx + 1] * conv_mask[1][2];

                double convolve_value_2_0 = in_mat[mat_y_idx + 1][mat_x_idx - 1] * conv_mask[2][0];
                double convolve_value_2_1 = in_mat[mat_y_idx + 1][mat_x_idx + 0] * conv_mask[2][1];
                double convolve_value_2_2 = in_mat[mat_y_idx + 1][mat_x_idx + 1] * conv_mask[2][2];

                double sum = 0.0;

                sum += convolve_value_0_0;
                sum += convolve_value_0_1;
                sum += convolve_value_0_2;

                sum += convolve_value_1_0;
                sum += convolve_value_1_1;
                sum += convolve_value_1_2;

                sum += convolve_value_2_0;
                sum += convolve_value_2_1;
                sum += convolve_value_2_2;

                out_mat[mat_y_idx][mat_x_idx] = sum;
            }
        }
    }
}

SIMD实现代码 #

void kernel(int iterations, int N, double in_mat[N][N], double out_mat[N][N], double conv_mask[3][3])
{
    for (int i = 0; i < iterations; i++)
    {
        for (int mat_y_idx = 1; mat_y_idx < N - 1; ++mat_y_idx)
        {
            int line_x_idx;
            for (line_x_idx = 1; line_x_idx < N - 1 - 4; line_x_idx += 4)
            {
                __m256d sum_0 = _mm256_setzero_pd();
                __m256d sum_1 = _mm256_setzero_pd();
                __m256d sum_2 = _mm256_setzero_pd();

                //Unrolled by the compiler when compiled with optimizations
                for (int kernel_index = -1; kernel_index <= 1; kernel_index++)
                {
                    __m256d in_values_0 = _mm256_loadu_pd(in_mat[mat_y_idx - 1] + line_x_idx + kernel_index);
                    __m256d in_values_1 = _mm256_loadu_pd(in_mat[mat_y_idx + 0] + line_x_idx + kernel_index);
                    __m256d in_values_2 = _mm256_loadu_pd(in_mat[mat_y_idx + 1] + line_x_idx + kernel_index);

                    sum_0 = _mm256_add_pd(sum_0, _mm256_mul_pd(in_values_0, _mm256_set1_pd(conv_mask[0][kernel_index + 1])));
                    sum_1 = _mm256_add_pd(sum_1, _mm256_mul_pd(in_values_1, _mm256_set1_pd(conv_mask[1][kernel_index + 1])));
                    sum_2 = _mm256_add_pd(sum_2, _mm256_mul_pd(in_values_2, _mm256_set1_pd(conv_mask[2][kernel_index + 1])));
                }

                _mm256_storeu_pd(out_mat[mat_y_idx] + line_x_idx, _mm256_add_pd(sum_0, _mm256_add_pd(sum_1, sum_2)));
            }

            //Leftover elements
            if (line_x_idx != N - 1)
            {
                for (; line_x_idx < N - 1; line_x_idx++)
                {
                    double sum = 0;

                    for (int kernel_index = -1; kernel_index <= 1; kernel_index++)
                    {
                        sum += in_mat[mat_y_idx - 1][line_x_idx + kernel_index] * conv_mask_premult[0][kernel_index + 1];
                        sum += in_mat[mat_y_idx + 0][line_x_idx + kernel_index] * conv_mask_premult[1][kernel_index + 1];
                        sum += in_mat[mat_y_idx + 1][line_x_idx + kernel_index] * conv_mask_premult[2][kernel_index + 1];
                    }

                    out_mat[mat_y_idx][line_x_idx] = sum;
                }
            }
        }
    }
}

核心问题分析 #

1. 编译器自动向量化：在-O3优化级别下，g++会自动对标量代码的内层循环进行SIMD向量化与循环展开。由于标量代码的内存访问模式连续、计算逻辑规整，编译器能生成与手写AVX2代码效率相当的指令序列，导致两者性能无差异。
2. 手写SIMD的额外开销：当前SIMD实现中，_mm256_set1_pd在kernel_index循环内重复调用，即使编译器会展开该循环，仍可能产生不必要的常量加载指令；而标量代码中编译器会提前将卷积掩码常量优化到寄存器，避免重复操作。
3. 内存访问对齐问题：手写代码使用_mm256_loadu_pd/_mm256_storeu_pd处理非对齐内存，相比对齐的_mm256_load_pd/_mm256_store_pd会有微小开销；但标量代码的内存访问同样是非对齐的，编译器自动向量化后也会生成类似指令，因此两者在内存开销上持平。
4. 计算强度匹配：3x3卷积每个输出元素需要9次乘加操作，计算强度与内存带宽需求的比例刚好让SIMD单元被充分利用，无论是手写还是编译器自动生成的代码，都无法进一步提升计算吞吐量。

优化建议 #

1. 验证自动向量化效果：添加编译选项-fopt-info-vec-all，查看标量代码的向量化日志，确认编译器是否已生成AVX2指令。若已完全向量化，手写SIMD的收益空间极小。
2. 预加载卷积掩码：将卷积掩码的所有元素预加载为__m256d常量并存储到寄存器，避免循环内重复调用_mm256_set1_pd：

   // 在最外层循环外预加载掩码
   __m256d mask0_0 = _mm256_set1_pd(conv_mask[0][0]);
   __m256d mask0_1 = _mm256_set1_pd(conv_mask[0][1]);
   __m256d mask0_2 = _mm256_set1_pd(conv_mask[0][2]);
   __m256d mask1_0 = _mm256_set1_pd(conv_mask[1][0]);
   __m256d mask1_1 = _mm256_set1_pd(conv_mask[1][1]);
   __m256d mask1_2 = _mm256_set1_pd(conv_mask[1][2]);
   __m256d mask2_0 = _mm256_set1_pd(conv_mask[2][0]);
   __m256d mask2_1 = _mm256_set1_pd(conv_mask[2][1]);
   __m256d mask2_2 = _mm256_set1_pd(conv_mask[2][2]);
   

3. 使用融合乘加指令：替换分开的乘法与加法为_mm256_fmadd_pd指令，减少指令数并提升流水线效率，例如：

   sum_0 = _mm256_fmadd_pd(in_values_0, mask0_k, sum_0);
   

4. 内存对齐优化：确保输入/输出矩阵的内存是32字节对齐的，使用aligned_alloc(32, size)或编译器属性__attribute__((aligned(32)))分配内存，将_mm256_loadu_pd/_mm256_storeu_pd替换为对齐版本，消除非对齐访问的额外开销。
5. 循环分块优化：针对8192×8192的大矩阵，采用循环分块（tiling）策略，将矩阵划分为适配L3缓存的子块，提升数据缓存命中率，减少内存访问延迟。例如将y方向划分为64行的块，x方向划分为64列的块，每次处理一个子块内的卷积计算。

内容的提问来源于stack exchange，提问作者Tom Clabault