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

【问题标题】：OPENMP - Parallelize Schwarz algorithm with preconditionsOPENMP - 使用先决条件并行化 Schwarz 算法
【发布时间】：2017-11-23 07:58:54
【问题描述】：

我需要在下面并行化 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），但没有成功...

非常感谢！

【问题讨论】：

标签： c parallel-processing mpi openmp

【解决方案1】：

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

您可以使用具有依赖关系的任务，这将是最简单的方法。 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
    }
}

【讨论】：