如何优化子矩阵-向量乘法（gesubmv）的C语言实现？

我实现了两个用于计算**y(rinds) = A(rinds, cinds)*x(cinds) + y(rinds)**的C函数：

• 函数gesubmv直接针对指定的行索引rinds和列索引cinds（均无需有序）计算，其中rinds对应rindsA_data(rindsA_1:rindsA_m)，cinds对应cindsA_data(cindsA_1:cindsA_n)。
• 函数gesubmv2采用间接方案：先计算y(:) = A(:, cinds)*x(cinds) + y(:)，再提取最终需要的y(rinds)。实际测试发现后者效率更高，推测可能和指令执行方式有关，也不确定是否有其他疏漏，希望得到优化建议。

原始实现代码 #

gesubmv（直接计算指定行） #

void gesubmv(const double A_data[], const int A_size[2],
             const int rindsA_data[], int rindsA_1, int rindsA_m,
             const int cindsA_data[], int cindsA_1, int cindsA_n,
             const double x_data[], double y_data[]) {
    int i;
    int j;
    for (i = rindsA_1; i <= rindsA_m; i++) {
        int k;
        double tmp_sum;
        tmp_sum = 0.0;
        for (j = cindsA_1; j <= cindsA_n; j++) {
            k = cindsA_data[j - 1] - 1;
            tmp_sum +=
                x_data[k] * A_data[k + (A_size[1] * (rindsA_data[i - 1] - 1))];
        }
        k = rindsA_data[i - 1] - 1;
        y_data[k] += tmp_sum;
    }
}

gesubmv2（先计算全量行再提取目标行） #

void gesubmv2(const double A_data[], const int A_size[2],
             const int cindsA_data[], int cindsA_1, int cindsA_n,
             const double x_data[], double y_data[]) {
    int i;
    int j;
    for (i = 0; i < A_size[1]; i++) {
        int k;
        double tmp_sum;
        tmp_sum = 0.0;
        for (j = cindsA_1; j <= cindsA_n; j++) {
            k = cindsA_data[j - 1] - 1;
            tmp_sum += x_data[k] * A_data[k + (A_size[1] * i)];
        }
        y_data[i] += tmp_sum;
    }
}

效率差异的核心原因 #

• 内存访问局部性差异：
  • gesubmv中，矩阵A的访问是跳着取不同行的指定列元素，属于非连续内存访问，会频繁触发CPU缓存失效，大幅降低执行效率。
  • gesubmv2中，A按连续行遍历，每次取该行的指定列元素，内存访问连续性更好，缓存命中率高，CPU能高效预加载数据。
• 循环与指令执行效率：
  • gesubmv外层循环遍历离散的rinds索引，每次都要做索引转换和数组随机访问，额外开销大，且CPU分支预测难以发挥作用。
  • gesubmv2外层循环是连续的行索引，指令流水线更顺畅，分支预测准确率高，同时更容易被编译器自动向量化，利用SIMD指令并行计算。

针对性优化建议 #

1. 优化gesubmv的内存访问模式：
   • 若rinds可提前排序，让行索引尽量连续，能显著提升缓存命中率；若无法排序，可考虑将矩阵A改为列优先存储，使指定行的列元素访问变为连续内存操作。
2. 减少索引计算冗余：
   • 提前计算并缓存循环中的重复索引值，避免重复计算和数组访问。优化后的gesubmv片段如下：

     for (i = rindsA_1; i <= rindsA_m; i++) {
         int target_row = rindsA_data[i - 1] - 1;
         int row_offset = A_size[1] * target_row;
         double tmp_sum = 0.0;
         for (j = cindsA_1; j <= cindsA_n; j++) {
             int col_idx = cindsA_data[j - 1] - 1;
             tmp_sum += x_data[col_idx] * A_data[col_idx + row_offset];
         }
         y_data[target_row] += tmp_sum;
     }
     

3. 开启编译器优化或手动SIMD加速：
   • 编译时开启高等级优化（如-O3）和SIMD指令集支持（如-mavx2），让编译器自动进行向量化；也可手动使用编译器内置函数（如__m256d）实现向量乘法累加。
4. 根据场景选择合适实现：
   • 若rinds覆盖行数远小于矩阵总行数，优化后的gesubmv更高效；若rinds接近全量行，gesubmv2本身效率足够，后续提取目标行的开销可忽略。

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