并行化 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 私下讨论后，我们得出了他在cmets 中提出的解决方案，为每个线程创建了一对向量x 和y，它们将在单个列上独立工作。

Answer 1

在与 OP 就 cmets（后来被删除）进行讨论后，我们都得出以下结论：

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

#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];
}

而不是有一个parallel region 每个 j 迭代。您可以提取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 数组，以便有在这些数组的更新/读取期间没有竞态条件。

#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);
}