如何使用OpenMP归约加速矩阵元素Fisher_M[FX][FY]的求和计算

原有OpenMP性能不及预期的核心原因 #

• 多个线程同时对共享变量Fisher_M[FX][FY]做累加，存在严重的竞态条件与缓存伪共享问题，不仅计算结果可能出错，缓存反复失效也会极大拖慢性能
• 你的双层循环每轮迭代计算量完全一致，没有负载不均的问题，使用dynamic调度会额外增加调度开销，反而拖累性能

方案1：OpenMP归约优化 #

你要实现的本质是两个矩阵的Hadamard乘积求和，完全可以用OpenMP归约实现，性能会有明显提升，修改后代码如下：

#define num_threads 8

double tmp_sum = 0.0;
// 加入reduction子句自动处理多线程累加合并，调度改为static消除调度开销
#pragma omp parallel for schedule(static) num_threads(num_threads) reduction(+:tmp_sum)
for(int i=0; i<CO_CL_WL.size(); i++){
    for(int j=0; j<CO_CL_WL.size(); j++){
        // 可省略if判断：只要两个元素有一个为0，乘积就是0，累加不影响结果，还能避免分支预测失败开销
        // 若矩阵稀疏度极高（99%以上元素为0），可保留if判断自行测试性能取舍
        tmp_sum += CO_CL_WL[i][j] * CO_CL_WL_D[i][j];
    }
}
Fisher_M[FX][FY] = tmp_sum;

方案2：MAGMA/CUDA实现 #

如果你的矩阵尺寸非常大（比如边长大于4096），可以用GPU方案加速：

• 该计算等价于两个矩阵拉平为一维向量后的点积运算，直接调用cuBLAS的cublasDdot（双精度）/cublasSdot（单精度）接口即可得到结果，不需要手动写CUDA核函数，官方库已经做了极致优化
• 用MAGMA的话也有对应的矩阵点积接口，还支持混合精度计算，对精度要求不高的场景可以进一步提速
• 注意GPU存在显存拷贝的固定开销，如果是单次计算小矩阵，OpenMP方案反而更快，只有矩阵尺寸足够大、计算量远大于拷贝开销时，GPU才有明显收益

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