基于AVX2 SIMD的矩阵乘法性能偏低，如何优化？

问题背景 #

实现了基于AVX2的矩阵乘法SIMD版本，对B做转置以优化列访问，结果性能反而远低于朴素版：

• 朴素版（单线程ijk循环）：n=1000时，7.421 Gflop/s
• SIMD版（编译选项CXXFLAGS=-ggdb -Wall -pedantic -mavx2 -O3，ijk循环）：n=1000时，2.221 Gflop/s

实现代码如下：

void multMat1(int n, float *A, float *B, float *C) {
    float *Bt = (float*)malloc(sizeof(float)*n*n);
    // transposing B to Bt
    for(int col=0; col< n; col+=BLOCK_SIZE){
        for(int row=0; row< n; row+=BLOCK_SIZE){
            for (int x = row; x < row+BLOCK_SIZE && x< n; x++) {
                for (int y = col; y < col+BLOCK_SIZE && y< n; y++) {
                    Bt[y + x * n] = B[x + y * n];
                }
            }
        }
    }
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            __m256 C_temp_vect = _mm256_setzero_ps();
            for (int k = 0; k < n; k += 8) {
                // loading 8 elements from A and Bt(transposed B)
                __m256 A_vect = _mm256_loadu_ps(&A[i * n + k]);
                __m256 B_scalar = _mm256_loadu_ps(&Bt[j * n + k]);
                C_temp_vect = _mm256_add_ps(C_temp_vect,_mm256_mul_ps(A_vect, B_scalar));
            }
            float temp[8];
            float sum = 0.0f;
            _mm256_store_ps(temp, C_temp_vect);
            // summation of all the elements of the vector temp.
            for(int index=0;index<8;index++)
                sum += temp[index];
            C[i * n + j] = sum;
        }
    }
    free(Bt);
}

性能低下的核心原因 #

1. 循环顺序导致缓存命中率极低
   当前循环顺序是i->j->k，对于每个j需要加载Bt的完整一行（Bt[j*n +k]，k遍历全范围）。当n=1000时，每行占4KB内存，而L1缓存通常仅32KB，最多容纳8行。遍历1000个j时，绝大多数行无法命中缓存，每次加载都要从内存读取，带来巨大延迟，完全抵消了SIMD的计算优势。

2. 水平求和效率低下
   将向量结果存入临时数组再通过标量循环累加，额外增加了内存读写操作，没有利用AVX2原生的水平加法指令，浪费了SIMD单元的计算能力。

3. 未对齐内存访问
   使用_mm256_loadu_ps处理未对齐内存，AVX2对未对齐内存的加载会引入额外处理开销，相比对齐内存的_mm256_load_ps效率更低。

4. 转置分块可能不合理
   如果BLOCK_SIZE设置过大（超过缓存容纳范围），转置过程中也会出现大量缓存miss，增加额外预处理开销。

优化方案 #

1. 调整循环顺序+分块优化 #

将循环改为分块式的i块->j块->k块，在块内使用i->k->j的顺序，确保Bt的访问在缓存内连续：

• 选择合适的块大小（如32x32，刚好能放入L1缓存），让每个块的所有数据都能被缓存覆盖。
• 块内循环顺序改为i->k->j：对于固定的i和k，遍历所有j，此时C[i][j]的访问是连续的，Bt[j][k]的访问变为按列遍历转置后的Bt（即原B的行），结合分块可大幅提升缓存命中率。

示例伪代码结构：

const int BLOCK = 32;
for (int i_block = 0; i_block < n; i_block += BLOCK) {
    for (int j_block = 0; j_block < n; j_block += BLOCK) {
        for (int k_block = 0; k_block < n; k_block += BLOCK) {
            // 块内计算：i从i_block到min(i_block+BLOCK, n)
            for (int i = i_block; i < std::min(i_block+BLOCK, n); i++) {
                for (int k = k_block; k < std::min(k_block+BLOCK, n); k += 8) {
                    __m256 a_vect = _mm256_load_ps(&A[i*n + k]);
                    for (int j = j_block; j < std::min(j_block+BLOCK, n); j++) {
                        __m256 b_vect = _mm256_load_ps(&Bt[j*n + k]);
                        __m256 c_vect = _mm256_load_ps(&C[i*n + j]);
                        c_vect = _mm256_add_ps(c_vect, _mm256_mul_ps(a_vect, b_vect));
                        _mm256_store_ps(&C[i*n + j], c_vect);
                    }
                }
            }
        }
    }
}

2. 使用AVX2水平求和指令优化 #

替换标量循环累加为AVX2原生指令，避免内存读写：

__m256 sum_vec = _mm256_hadd_ps(C_temp_vect, C_temp_vect);
sum_vec = _mm256_hadd_ps(sum_vec, sum_vec);
// 提取最终求和结果
float sum = _mm_cvtss_f32(_mm256_extractf128_ps(sum_vec, 1));
sum += _mm_cvtss_f32(_mm256_castps256_ps128(sum_vec));

或者使用GCC 8+支持的简化指令：

float sum = _mm256_reduce_add_ps(C_temp_vect);

3. 内存对齐优化 #

使用对齐内存分配函数替代malloc，确保AVX2向量操作的内存对齐：

// 32字节对齐（AVX2向量长度）
float *Bt = (float*)_aligned_alloc(32, sizeof(float)*n*n);
// 使用后释放
_aligned_free(Bt);

同时将_mm256_loadu_ps替换为_mm256_load_ps，保持对齐存储。

4. 优化转置分块大小 #

将BLOCK_SIZE设置为与缓存匹配的值（如32或64），确保转置时每个块都能完全放入L1缓存，减少转置过程中的缓存miss。

5. 减少冗余计算 #

提前计算循环中的偏移量，避免重复乘法：

int i_offset = i * n;
__m256 A_vect = _mm256_load_ps(&A[i_offset + k]);

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