如何加速5000×5000稠密矩阵乘法？一维/二维数组哪种更优？

我希望加速标准矩阵乘法运算（即 ( C[i][j] += \sum_{k=0}^{size-1} A[i][k] \times B[k][j] )），以下是需要优化的原始代码：

void multiply(int size, int **matA, int **matB, int ** matC) {
    for(int i=0;i<size;i++) {
        for(int j=0;j<size;j++) {
            int t = matC[i][j];
            for(int k=0;k<size;k++) {
                t += matA[i][k] * matB[k][j];
            }
            matC[i][j] += t;
        }
    }
}

我使用的输入矩阵和结果矩阵尺寸均为5000×5000，如此大的矩阵可能无法完全载入缓存，循环中是否会产生大量页错误？我想了解如何加速该乘法运算，以及该如何组织数据（使用一维数组还是二维数组？）

我尝试了如下优化代码：用一维数组模拟二维数组（每个矩阵调用一次new []分配连续内存），同时使用临时矩阵存储matB的转置以避免循环中的非连续内存访问，并加入AVX2指令集提升性能。但不确定一维数组是否比直接使用二维数组更快，也想知道还有哪些其他优化思路。

尝试的优化代码：

int** allocate(int rows, int cols) {
    int ** mat;
    mat = new int*[rows];
    int *temp = new int[rows*cols];
    for(int i = 0; i<rows; i++) {
        mat[i] = temp + i * cols;
    }
    return mat;
}


void multiply(int size, int **matA, int **matB, int ** matC) {
int i;
int n = size*size; // total size


// column-major order
int **transMatB = allocate(size, size);
int *transArrB = transMatB[0];

//copy transposed data, maybe many page faults here.
#pragma omp parallel for
    for(i = 0; i < n; i++) {
        transMatB[i/size][i%size] = matB[i%size][i/size];
    }

#pragma omp parallel for
    for(i = 0; i < n; i ++) {
        int *row = matA[i / size];
        int *col = transMatB[i % size];
        int temp;
#ifdef __AVX2__
        temp = multiplyAndSumArrays(row, col, size);
#else
        temp = 0;
        for (int k = 0; k < size; k ++) {
            temp += row[k] * col[k];
        }
#endif
        matC[i / size][i % size] += temp;
    }
    // remove temp transposed mastrix
    delete[] transMatB[0];
    delete[] transMatB; transMatB = nullptr;
}

一、一维数组 vs 二维数组：一维数组更优 #

对于5000×5000的大矩阵，用一维数组模拟二维数组（连续内存分配）的性能远优于传统指针数组式二维数组，原因如下：

• 传统int**型二维数组是指针数组，每一行的内存可能离散分布，会导致频繁的缓存失效和页错误——CPU缓存按连续内存块加载，离散的行内存无法被缓存高效利用。
• 你实现的allocate函数用单块连续内存分配，再通过指针偏移模拟二维访问，既保留了二维数组的语法便利性，又保证了内存连续性，能最大化缓存命中率，减少页错误。

如果追求极致性能，甚至可以直接用一维数组（int*），通过index = i*size + j计算下标，避免指针数组的间接寻址开销，代价是代码可读性稍降。

二、核心优化技巧 #

1. 转置优化的细节调整 #

你的转置思路是对的——原始代码中访问matB[k][j]是按列访问，属于非连续内存访问，缓存效率极低。转置后按行访问能匹配缓存加载模式，但当前实现可以优化：

• 避免用i/size和i%size这类除法/取模运算，换成嵌套循环，减少运算开销：

  #pragma omp parallel for collapse(2)
  for(int i=0; i<size; i++){
      for(int j=0; j<size; j++){
          transMatB[j][i] = matB[i][j];
      }
  }
  

  用collapse(2)能让OpenMP更好地并行化嵌套循环，提升转置效率。
• 转置时按块处理（比如64×64块），利用缓存加载连续内存，减少转置过程中的缓存失效。

2. 分块优化（Blocking/Tiling）——解决大矩阵缓存放不下的问题 #

5000×5000的矩阵（单矩阵约95MB，三个矩阵近285MB）远超普通CPU的L3缓存容量（16-64MB），分块优化能让计算过程中只把当前需要的小块数据载入缓存，大幅提升命中率：

• 选择匹配缓存容量的块大小（比如64×64），将矩阵拆分成多个小块按块计算：

  const int BLOCK_SIZE = 64;
  #pragma omp parallel for collapse(2)
  for(int i_block=0; i_block<size; i_block+=BLOCK_SIZE){
      for(int j_block=0; j_block<size; j_block+=BLOCK_SIZE){
          for(int i=i_block; i<min(i_block+BLOCK_SIZE, size); i++){
              for(int j=j_block; j<min(j_block+BLOCK_SIZE, size); j++){
                  int temp = 0;
                  for(int k_block=0; k_block<size; k_block+=BLOCK_SIZE){
                      for(int k=k_block; k<min(k_block+BLOCK_SIZE, size); k++){
                          temp += matA[i][k] * transMatB[j][k];
                      }
                  }
                  matC[i][j] += temp;
              }
          }
      }
  }
  

  分块后，每个块的数据能被缓存反复利用，避免频繁从内存加载数据。

3. AVX2指令集的高效实现 #

multiplyAndSumArrays需要针对性实现才能发挥AVX2的性能，AVX2支持256位寄存器，一次可处理8个int：

int multiplyAndSumArrays(const int* a, const int* b, int size){
    __m256i sum = _mm256_setzero_si256();
    int i;
    // 处理能被8整除的部分
    for(i=0; i<=size-8; i+=8){
        __m256i vec_a = _mm256_loadu_si256((const __m256i*)(a+i));
        __m256i vec_b = _mm256_loadu_si256((const __m256i*)(b+i));
        __m256i mul = _mm256_mullo_epi32(vec_a, vec_b);
        sum = _mm256_add_epi32(sum, mul);
    }
    // 处理剩余不足8个的元素
    int total = 0;
    int temp[8];
    _mm256_storeu_si256((__m256i*)temp, sum);
    for(int j=0; j<8; j++) total += temp[j];
    for(; i<size; i++) total += a[i]*b[i];
    return total;
}

如果能保证内存对齐（比如用_mm_malloc分配对齐内存），可以把_mm256_loadu_si256换成_mm256_load_si256，进一步提升性能。

4. OpenMP并行优化细节 #

• 用collapse(2)并行化最外层的块循环，最大化利用多核CPU，避免负载不均。
• 可以通过omp_set_num_threads()设置线程数（一般等于CPU核心数），或通过OMP_NUM_THREADS环境变量指定。
• 避免在并行区域内进行内存分配/释放，减少线程同步开销。

5. 页错误与内存预取 #

• 大矩阵第一次访问时必然产生页错误，可以提前进行一次遍历写入（比如初始化0），触发操作系统提前分配物理页，避免计算过程中频繁缺页。
• 用__builtin_prefetch（GCC）或_mm_prefetch（MSVC）预取下一块数据，减少后续计算的等待时间：

  __builtin_prefetch(&matA[i+1][k], 0, 3);
  

6. 其他细节优化 #

• 确保matC在调用multiply前被初始化为0，避免不确定的结果。
• 如果场景允许，直接在matB的转置版本上进行后续计算，省去重复转置的内存拷贝开销。

内容的提问来源于stack exchange，提问作者Z.Lun