有没有办法并行化下三角矩阵求解器？答案

【问题标题】：Is there a way to parallelize a lower triangle matrix solver?有没有办法并行化下三角矩阵求解器？
【发布时间】：2021-06-22 10:48:27
【问题描述】：

我们的目标是将 OpenMP 并行化添加到 for (i = 0; i < n; i++)，以用于形式 Ax=b 的下三角求解器。预期结果与没有添加到 for (i = 0; i < n; i++) 的并行化时的结果完全相同。 vector<vector<double>> 表示二维矩阵。 makeMatrix(int m, int n) 初始化大小为 mxn 的全零的 vector<vector<double>>。

cmets 中留下了两个最突出的尝试。

vector<vector<double>> lowerTriangleSolver(vector<vector<double>> A, vector<vector<double>> b)
{
    vector<vector<double>> x = makeMatrix(A.size(), 1);
    int i, j;
    int n = A.size();
    double s;

     //#pragma omp parallel for reduction(+: s)
     //#pragma omp parallel for shared(s)
     for (i = 0; i < n; i++)
     {
         s = 0.0;
         #pragma omp parallel for
         for (j = 0; j < i; j++)
         {
             s = s + A[i][j] * x[j][0];
         }
         x[i][0] = (b[i][0] - s) / A[i][i];
     }

    return x;
}

【问题讨论】：

标签： c++ multithreading performance parallel-processing openmp

【解决方案1】：

您可以尝试在线程之间分配外部循环迭代，而不是内部循环。这样，你增加了并行任务的粒度，避免了's'变量的减少。

 #pragma omp parallel for
 for (int i = 0; i < n; i++){
      double s = 0.0;
      for (int j = 0; j < i; j++){
           s = s + A[i][j] * x[j][0];
      }
      x[i][0] = (b[i][0] - s) / A[i][i];
 }

不幸的是，这是不可能的，因为s = s + A[i][j] * x[j][0]; 和x[i][0] = (b[i][0] - s) / A[i][i]; 之间存在依赖关系，更准确地说，x[j][0] 依赖于x[i][0]。

所以你可以尝试两种方法：

 for (int i = 0; i < n; i++){
      double s = 0.0;
      #pragma omp parallel for reduction(+:s)
      for (int j = 0; j < i; j++){
           s = s + A[i][j] * x[j][0];
      }
      x[i][0] = (b[i][0] - s) / A[i][i];
 }

或使用 SIMD：

 for (int i = 0; i < n; i++){
      double s = 0.0;
      #pragma omp simd reduction(+:s)
      for (int j = 0; j < i; j++){
           s = s + A[i][j] * x[j][0];
      }
      x[i][0] = (b[i][0] - s) / A[i][i];
 }

【讨论】：