Is there a way to parallelize a lower triangle matrix solver?

Question

The goal is to add OpenMP parallelization to for (i = 0; i < n; i++) for the lower triangle solver for the form Ax=b. Expected result is exactly same as the result when there is NO parallelization added to for (i = 0; i < n; i++). vector<vector<double>> represents a 2-D matrix. makeMatrix(int m, int n) initializes a vector<vector<double>> of all zeroes of size mxn.

Two of the most prominent tries have been left in comments.

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

dreamcrash · Accepted Answer · 2021-03-25T23:37:59.433

You could try to assign the outer loop iterations among threads, instead of the inner loop. In this way, you increase the granularity of the parallel tasks and avoid the reduction of the 's' variable.

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

Unfortunately, that is not possible because there is a dependency between s = s + A[i][j] * x[j][0]; and x[i][0] = (b[i][0] - s) / A[i][i];, more precisely x[j][0] depends upon the x[i][0].

So you can try two approaches:

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

or using 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];
 }

Is there a way to parallelize a lower triangle matrix solver?

1 Answers1