c - OPENMP - 使用先决条件并行化 Schwarz 算法

我需要在下面并行化 Schwarz 算法，但我不知道如何处理前提条件以及存在嵌套循环的事实。
我必须使用 OpenMP 或 MPI。

void ssor_forward_sweep(int n, int i1, int i2, int j1, int j2, int k1, int k2, double* restrict Ax, double w)
{
#define AX(i,j,k) (Ax[((k)*n+(j))*n+(i)])

int i, j, k;
double xx, xn, xe, xu;

for (k = k1; k < k2; ++k) {
   for (j = j1; j < j2; ++j) {
      for (i = i1; i < i2; ++i) {
         xx = AX(i,j,k);
         xn = (i > 0)   ? AX(i-1,j,k) : 0;
         xe = (j > 0)   ? AX(i,j-1,k) : 0;
         xu = (k > 0)   ? AX(i,j,k-1) : 0;
         AX(i,j,k) = (xx+xn+xe+xu)/6*w;
      }
   }
}
#undef AX
}

考虑到每个循环都使用之前循环中的值，如何并行化此函数以获得最佳时间。

我已经尝试将循环两两并行化或通过分块(如 Stencil Jacobi 3D)来并行化，但没有成功......

非常感谢你 !

最佳答案

不幸的是，循环间数据依赖性限制了您可以在嵌套循环中获得的并行量。

您可以使用具有依赖性的任务，这将是最简单的方法。 OpenMP 运行时库将负责调度，您只关注您的算法。另一个好的方面是在任何循环结束时都没有同步，而只是在代码的相关部分之间。

#pragma omp parallel
#pragma omp single
for (int k = 0; k < k2; k += BLOCK_SIZE) {
  for (int j = 0; j < j2; j += BLOCK_SIZE) {
    for (int i = 0; i < i2; i += BLOCK_SIZE) {
      #pragma omp task depend (in: AX(i-1,j,k),AX(i,j-1,k),AX(i,j,k-1)) \
                       depend (out: AX(i,j,k))
      // your code here
    }
  }
}

任务有时比并行循环更昂贵(取决于工作负载和同步粒度)，因此另一种选择是 wavefront parallelization 模式，它基本上转换迭代空间，使内循环中的元素相互独立(因此您可以在那里使用并行)。



如果您想采用任何一种方法，我强烈建议您将算法转换为阻塞算法:展开您的 3 嵌套循环以分两个阶段进行计算:

在固定大小的块/立方体之间迭代(让我们称你的新归纳变量为 ii 、 jj 和 kk )。

对于每个块，调用循环的原始串行版本。

阻塞的目的是增加并行部分的粒度，使并行化开销不那么明显。

这是阻塞部分的一些伪代码:

#define min(a,b) ((a)<(b)?(a):(b))
// Inter block iterations
for (int kk = 0; kk < k2; kk += BLOCK_SIZE) {
  for (int jj = 0; jj < j2; jj += BLOCK_SIZE) {
    for (int ii = 0; ii < i2; ii += BLOCK_SIZE) {

      // Intra block iterations
      for (int k = kk; k < min(k2,k+BLOCK_SIZE); k++) {
        for (int j = jj; j < min(j2,j+BLOCK_SIZE); j++) {
          for (int i = ii; i < min(i2,i+BLOCK_SIZE); i++) {
              // Your code goes here
          }
        }
      }

    }
  }
}

在波前并行化的情况下，最后一步是将外部循环(块间迭代)转换为波前，以便您迭代彼此之间不相关的元素。在 3D 迭代空间中，它基本上是一个从 (0,0,0) 前进到 (i2,j2,k2) 的对角平面。类似于下图中以红色突出显示的那个。



我将举一个二维波前的例子，因为它更容易理解。

#define min(a,b) ((a)<(b)?(a):(b))
#pragma omp parallel
for (int d = 1; d < i2+j2; d++ ) {
    int i = min(d,i2) - 1;
    int j = 0;
    // Iterations in the inner loop are independent
    // Implicit thread barrier (synchronization) at the end of the loop
    #pragma omp for
    for ( ; i >= 0 && j < min(d,j2); i--, j++) {
      // your code here
    }
}

关于c - OPENMP - 使用先决条件并行化 Schwarz 算法，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/47450267/