c++ - 有没有办法并行化下三角矩阵求解器？

Question

目标是将 OpenMP 并行化添加到 for (i = 0; i < n; i++)用于形式 Ax=b 的下三角求解器.预期结果与未向 for (i = 0; i < n; i++) 添加并行化时的结果完全相同.vector<vector<double>>表示一个二维矩阵。 makeMatrix(int m, int n)初始化 vector<vector<double>>所有大小为零的 mxn .
评论中留下了两个最突出的尝试。

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

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