使用OpenMP并行化下三角矩阵向量乘法的更高效实现方案

核心优化方案 #

• 首先修复原始代码的笔误：你当前内层循环的迭代变量误写为k++，实际应为j++，该问题会导致程序逻辑错误、死循环或内存越界，是首先要解决的问题。
• 解决负载不均问题：下三角矩阵的行计算量随i增大线性递增，默认schedule(static)是将连续大块行分配给单个线程，会出现前序线程早跑完、后序线程负载过重的问题，导致整体效率低下。可以改用小chunk的静态交错调度：schedule(static, 256)，将256行作为一个小分块交替分配给各个线程，既保证了静态调度的低开销，又实现了负载均衡。也可以测试schedule(guided, 128)，调度器会自动从大到小调整分块大小，适配负载递增的场景，开销远低于dynamic调度。
• 消除伪共享开销：result数组是共享连续存储，多线程同时修改相邻下标的result[i]会触发缓存行一致性同步，带来额外开销。可以为每一行的计算引入私有临时变量，计算完成后再一次性写回result数组。
• 开启SIMD向量化：内层循环的乘加操作符合向量化条件，可以添加#pragma omp simd制导指令，配合编译器的-mavx2（或-march=native）编译选项，利用CPU的向量指令集并行提升单线程计算效率，注意要保证matrix、vector、result数组按64字节地址对齐，否则向量化效率会打折扣。

优化后的参考代码：

// 提前将result数组初始化为0
#pragma omp parallel shared(matrix, result, vector) private(i, j)
{
    double tmp; // 私有临时变量，避免伪共享
    #pragma omp for schedule(static, 256)
    for (i = 0; i < n; i++)
    {
        tmp = 0.0;
        #pragma omp simd reduction(+:tmp)
        for (j = 0; j <= i; j++) // 下三角特性，j不可能超过i，无需额外判断j<n
        {
            tmp += matrix[i * n + j] * vector[j];
        }
        result[i] = tmp;
    }
}

针对30000阶的下三角矩阵，上述优化可以在原有基础上再获得1.5~3倍的性能提升，具体数值取决于你的CPU核心数和内存带宽

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