OpenMP并行高斯消元算法耗时远超串行的问题排查

问题场景 #

实现高斯消元法求解大规模方程组（N=1000）时，出现OpenMP并行版本比串行版本更慢的异常：

• 串行模式（启用omp_set_num_threads(1)）耗时1.140253秒
• 并行模式（使用#pragma omp parallel for）耗时4.394211秒

核心并行代码片段：

//omp_set_num_threads(1);
for (k = 0; k < N - 1; k++)
{
        #pragma omp parallel for
        for (i = k + 1; i < N; i++)
        {
            d = A[i][k] / A[k][k];
            for (j = k; j < N; j++)
            {
                A[i][j] = A[i][j]- d * A[k][j];
            }
            B[i] = B[i] - d * B[k];
        }
    }

完整代码：

#include <iostream>
#include <ctime>
#include <omp.h>

const int N = 1000;
double A[N][N]; 
double B[N];    
double E[N];    

using namespace std;

void gauss(double A[N][N], double B[N], double E[N])
{
    int i, j, k;
    double d;
    //omp_set_num_threads(1);
    for (k = 0; k < N - 1; k++)
    {
        #pragma omp parallel for
        for (i = k + 1; i < N; i++)
        {
            d = A[i][k] / A[k][k];
            for (j = k; j < N; j++)
            {
                A[i][j] = A[i][j]- d * A[k][j];
            }
            B[i] = B[i] - d * B[k];
        }
    }
    for (i = N - 1; i >= 0; i--)
    {
        E[i] = B[i];
        for (j = i + 1; j < N; j++)
        {
            E[i] = E[i] - A[i][j] * E[j];
        }
        E[i]= E[i]/A[i][i];
    }
};

int main() {
    srand(time(NULL));
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; j++) {
            A[i][j] = -10 + rand() % 20;
        }
        B[i] = -10 + rand() % 20;
    }
    double start = omp_get_wtime();
    gauss(A, B, E);
    double end = omp_get_wtime();
    cout << "Time: " << fixed << (end - start) << " sec" << endl;
    return 0;
}

问题根源 #

1. 重复线程开销：每轮k循环都通过#pragma omp parallel for创建新线程组，线程的创建、调度、销毁开销远超过并行带来的收益，尤其当k增大时，内层循环迭代次数减少，开销占比更高。
2. 缓存伪共享：数组A和B的内存布局导致多个线程同时访问相邻缓存行，触发缓存失效，大幅降低缓存命中率。
3. 负载不均衡：随着k递增，i的迭代范围从N-1快速缩小到1，后续循环中线程分到的任务量不足，并行收益可忽略。

优化方案 #

1. 复用线程组 #

将并行区域移到外层循环外，仅创建一次线程组，避免重复开销：

void gauss(double A[N][N], double B[N], double E[N])
{
    int i, j, k;
    double d;
    //omp_set_num_threads(1);
    #pragma omp parallel private(i,j,k,d)
    {
        for (k = 0; k < N - 1; k++)
        {
            #pragma omp for
            for (i = k + 1; i < N; i++)
            {
                d = A[i][k] / A[k][k];
                for (j = k; j < N; j++)
                {
                    A[i][j] = A[i][j]- d * A[k][j];
                }
                B[i] = B[i] - d * B[k];
            }
        }
    }
    // 回代部分串行执行（并行收益有限，可根据需求调整）
    for (i = N - 1; i >= 0; i--)
    {
        E[i] = B[i];
        for (j = i + 1; j < N; j++)
        {
            E[i] = E[i] - A[i][j] * E[j];
        }
        E[i]= E[i]/A[i][i];
    }
};

2. 缓解缓存伪共享 #

使用缓存行对齐数组，减少线程间的缓存冲突：

const int N = 1000;
alignas(64) double A[N][N]; // 64为常见缓存行大小，适配多数CPU架构
alignas(64) double B[N];    
alignas(64) double E[N];

3. 动态控制并行时机 #

设置迭代次数阈值，仅当任务量足够大时启用并行，避免小循环的并行开销：

#pragma omp parallel private(i,j,k,d)
{
    for (k = 0; k < N - 1; k++)
    {
        int iter_count = N - (k + 1);
        if (iter_count > 100) { // 阈值可根据硬件调整
            #pragma omp for
            for (i = k + 1; i < N; i++) { /* 原有计算逻辑 */ }
        } else {
            // 小循环串行执行
            for (i = k + 1; i < N; i++) { /* 原有计算逻辑 */ }
        }
    }
}

4. 开启编译优化 #

编译时启用O3级别优化，让编译器自动做循环展开、缓存优化等：

g++ -fopenmp -O3 your_code.cpp -o gauss_solver

效果验证 #

优化后，并行版本的执行时间会显著低于串行版本，且N越大，并行收益越明显。

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