并行化 BW 和 FW 算法的最佳方法

Question

我已经实现了 BW 和 FW 算法来求解 L 和 U 三角矩阵。我以串行方式非常快速地实现运行的算法，但我不知道这是否是并行化它的最佳方法。我认为我已经考虑了所有可能的数据竞争（alpha 版），对吗？

void solveInverse (double **U, double **L, double **P, int rw, int cw) {
    double **inverseA = allocateMatrix(rw,cw);
    double* x = allocateArray(rw);
    double* y = allocateArray(rw);
    
    double alpha;
    
    //int i, j, t;
    
    // Iterate along the column , so at each iteration we generate a column of the inverse matrix
    for (int j = 0; j < rw; j++) {
        
        // Lower triangular solve Ly=P
        y[0] = P[0][j];
        #pragma omp parallel for reduction(+:alpha)
        for (int i = 1; i < rw; i++) {
            alpha = 0;
            for (int t = 0; t <= i-1; t++)
                alpha += L[i][t] * y[t];
            y[i] = P[i][j] - alpha;
        }

        // Upper triangular solve Ux=P
        x[rw-1] = y[rw-1] / U[rw-1][rw-1];
        #pragma omp parallel for reduction(+:alpha)
        for (int i = rw-2; (i < rw) && (i >= 0); i--) {
            alpha = 0;
            
            for (int t = i+1; t < rw; t++)
                alpha += U[i][t]*x[t];
            x[i] = (y[i] - alpha) / U[i][i];
        }
        
        for (int i = 0; i < rw; i++)
            inverseA[i][j] = x[i];  
        }
    freeMemory(inverseA,rw);
    free(x);
    free(y);
}

在与用户 dreamcrash 私下讨论后，我们找到了他在评论中提出的解决方案，为每个线程创建一对向量 x 和 y单列。

Answer 1

在与 OP 就评论（后来被删除）进行讨论后，我们得出的结论是：

您不需要减少 alpha 变量，因为在第一个parallel region 之外它会再次初始化为零。相反，将 alpha 变量设为 private.

#pragma omp parallel for
for (int i = 1; i < rw; i++) {
    double alpha = 0;
    for (int t = 0; t <= i-1; t++)
        alpha += L[i][t] * y[t];
    y[i] = P[i][j] - alpha;
}

同样适用于第二个parallel region。

#pragma omp parallel for
for (int i = rw-2; (i < rw) && (i >= 0); i--) {
    double alpha = 0;
    for (int t = i+1; t < rw; t++)
        alpha += U[i][t]*x[t];
    x[i] = (y[i] - alpha) / U[i][i];
}

而不是每个 j 迭代一次 parallel region 。您可以提取 parallel region 来封装整个最外层循环，并使用 #pragma omp for 而不是 #pragma omp parallel for。尽管如此，虽然我们使用这种方法将 parallel regions 的创建数量从 rw 减少到只有 1，但通过这种优化实现的加速应该不会那么显着，因为高效的 OpenMP 实现将使用线程池，其中线程在第一个 parallel region 上初始化，但在下一个 parallel regions 上重复使用。因此，节省了创建和销毁线程的开销。

#pragma omp parallel { for (int j = 0; j < rw; j++) { y[0] = P[0][j]; #pragma omp for for (int i = 1; i < rw; i++) { double alpha = 0; for (int t = 0; t <= i-1; t++) alpha += L[i][t] * y[t]; y[i] = P[i][j] - alpha; } x[rw-1] = y[rw-1] / U[rw-1][rw-1]; #pragma omp for for (int i = rw-2; (i < rw) && (i >= 0); i--) { double alpha = 0; for (int t = i+1; t < rw; t++) alpha += U[i][t]*x[t]; x[i] = (y[i] - alpha) / U[i][i]; } #pragma omp for for (int i = 0; i < rw; i++) inverseA[i][j] = x[i]; } }

我已经向您展示了此代码转换，以便您可以看到一些潜在的技巧，您可以在未来的其他并行化中使用这些技巧。不幸的是，并行化将不起作用。

为什么？

让我们看看第一个循环：

#pragma omp parallel for for (int i = 1; i < rw; i++) { double alpha = 0; for (int t = 0; t <= i-1; t++) alpha += L[i][t] * y[t]; y[i] = P[i][j] - alpha; }

在 alpha += L[i][t] * y[t]; 中读取的 y[t] 和在 y[i] = P[i][j] - alpha; 中写入的 y[i] 之间存在依赖关系。

所以你可以做的是并行化最外层的循环（即，将每一列分配给线程）并创建单独的 x 和 y 数组，以便在这些数组的 updates/reads 期间没有 race-conditions。

#pragma omp parallel { double* x = allocateArray(rw); double* y = allocateArray(rw); #pragma omp for for (int j = 0; j < rw; j++) { y[0] = P[0][j]; for (int i = 1; i < rw; i++) { double alpha = 0; for (int t = 0; t <= i-1; t++) alpha += L[i][t] * y[t]; y[i] = P[i][j] - alpha; } x[rw-1] = y[rw-1] / U[rw-1][rw-1]; for (int i = rw-2; i >= 0; i--) { double alpha = 0; for (int t = i+1; t < rw; t++) alpha += U[i][t]*x[t]; x[i] = (y[i] - alpha) / U[i][i]; } for (int i = 0; i < rw; i++) inverseA[i][j] = x[i]; } free(x); free(y); }

并行化 BW 和 FW 算法的最佳方法

Best approach to parallelize BW and FW algorithms

c

parallel-processing

performance

multithreading

openmp