具有数据依赖性的for循环矢量化

我有基于BiCCG(共轭梯度)的矩阵求解器的实现,它也考虑了周期性。 碰巧的情况是,实现是计算密集型的,并且由于依赖性问题,循环不会自动向量化。 我进行了一些探索,似乎红黑高斯seidel算法比vanilla版本(也有类似的依赖性问题)更有效并行化。

可以更改此循环/算法,以便可以有效地进行矢量化吗?

// FORWARD #pragma omp for schedule(static, nx/NTt) for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) { dummy = res_sparse_s[i][j][k]; dummy -= COEFF[i][j][k][7] * RLL[i-1][j][k]; if (PeriodicBoundaryX && i==nx)dummy -= COEFF[i][j][k][8] * RLL[1 ][j][k]; dummy -= COEFF[i][j][k][2] * RLL[i][j-1][k]; if (PeriodicBoundaryY && j==ny) dummy -= COEFF[i][j][k][3] * RLL[i][1 ][k]; dummy -= COEFF[i][j][k][4] * RLL[i][j][k-1]; if (PeriodicBoundaryZ && k==nz) dummy -= COEFF[i][j][k][5] * RLL[i][j][1 ]; RLL[i][j][k] = dummy / h_sparse_s[i][j][k]; } 

对于任何迭代i,j,k的PS RLL在i-1,j-1和k-1处通过变量dummy合并了更新的“RLL”。 此外,现在循环仅使用指令schedule(static, nx/NTt)在x方向上分解schedule(static, nx/NTt)其中NTt只是可用线程数的宏。 可以使用指令collapse向各个方向分解吗?

——-主要编辑————————–以下Ajay的回答是一个最小的工作示例

 #include #include #include #include typedef double lr; #define nx 4 #define ny 4 #define nz 4 void print3dmatrix(double a[nx+2][ny+2][nz+2]) { for(int i=1; i<= nx; i++) { for(int j=1; j<= ny; j++) { for(int k=1; k<= nz; k++) { printf("%f ", a[i][j][k]); } printf("\n"); } printf("\n"); } } int main() { double a[nx+2][ny+2][nz+2]; double b[nx+2][ny+2][nz+2]; srand(3461833726); // matrix filling // b is just a copy of a for(int i=0; i< nx+2; i++) for(int j=0; j< ny+2; j++) for(int k=0; k< nz+2; k++) { a[i][j][k] = rand() % 5; b[i][j][k] = a[i][j][k]; } // loop 1 //#pragma omp parallel for num_threads(1) for(int i=1; i<= nx; i++) for(int j=1; j<= ny; j++) for(int k=1; k<= nz; k++) { a[i][j][k] = -1*a[i-1][j][k] - 1*a[i][j-1][k] -1 * a[i][j][k-1] + 4 * a[i][j][k]; } print3dmatrix(a); printf("******************************\n"); // loop 2 //#pragma omp parallel for num_threads(1) for(int i=1; i<= nx; i++) for(int j=1; j<= ny; j++) // #pragma omp simd for(int m=j+1; m<= j+nz; m++) { b[i][j][mj] = -1*b[i-1][j][mj] - 1*b[i][j-1][mj] -1 * b[i][j][mj-1] + 4 * b[i][j][mj]; } print3dmatrix(b); printf("=========================\n"); return 0; } 

主要观察 –

  1. 矩阵a用0到5之间的随机数填充,循环1是非变换原始循环,而循环2是变换循环
  2. 变换后的循环已经偏斜,以便消除依赖性。
  3. 如果在没有openmp并行化的情况下运行,则操作矩阵a和b相同
  4. 如果部署了open mp,则答案会发生变化(因为可能是竞争条件)[无论是否放置了pragma,循环都不可解析]
  5. 如果#pragma omp simd用于强制最内层循环的矢量化,则失败。

这是循环携带依赖的经典问题。 你的每次迭代都依赖于其他一些迭代(已经完成),因此可以连续地调度它的唯一方式。

但这只是因为你的循环是如何编写的。

你提到R[i][j][k]取决于R[i-1][j][k]R[i][j-1][k]R[i][j][k-1] 。 我看到三个依赖 –

  1. [1,0,0]
  2. [0,1,0]
  3. [0,0,1]

我希望这种表述很直观。

对于你现在的场景,依赖1)和2)不是问题,因为k0i / j1 ,这意味着迭代不依赖于k的先前迭代来完成这两个依赖。

问题是因为3)。 由于k1 ,因此每次迭代都取决于它之前的迭代。 如果以某种方式我们能够在i / j带一个>0的数字,我们就会完成。 循环偏斜变换让我们完全相同。

3D示例有点难以理解。 那么让我们看看ij 2D示例。

假设 – R[i][j]取决于R[i-1][j]R[i][j-1] 。 我们有同样的问题。

如果我们必须在图片中表示它,它看起来像这样 –

 . <- . <- . | | vv . <- . <- . | | vv . . . 

在该图中,每个点表示和迭代(i,j)和源自每个点的箭头指向它所依赖的迭代。 很明显,为什么我们不能在这里并行化最内层的循环。

但是假设我们做了偏差 -

  . /| / | . . /| /| / | / | . . . /| / | . . . 

如果你绘制与上图相同的箭头(我无法在ASCII艺术中绘制对角线箭头)。

您将看到所有箭头都指向下方,即它们至少继续向下迭代,这意味着您可以并行化水平循环。

现在说你的新循环尺寸是y (外循环)和x (内循环),

你的原始变量ij将是

j = xi = x - y

你的循环体因此变成 -

 for ( y = 0; y < j_max + i_max; y++) for ( x = 0; x < j_max; x++) R_dash[y][x] = R_dash[y-1][x-1] + R_dash[y-1][x]; 

其中R_dash是偏斜域并且具有到R一对一映射

您将看到R_dash[y-1][x-1]R_dash[y-1][x]都将在y的某个先前迭代中计算。 因此,您可以完全并行化x循环。

这里应用的转换是

i -> i, j -> i + j

您可以类似地将其用于3个维度。

为了进一步了解仿射变换的工作原理以及它们如何用于引入并行性,您可以看到这些讲义 。