具有数据依赖性的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; }
主要观察 –
- 矩阵a用0到5之间的随机数填充,循环1是非变换原始循环,而循环2是变换循环
- 变换后的循环已经偏斜,以便消除依赖性。
- 如果在没有openmp并行化的情况下运行,则操作矩阵a和b相同
- 如果部署了open mp,则答案会发生变化(因为可能是竞争条件)[无论是否放置了pragma,循环都不可解析]
- 如果
#pragma omp simd
用于强制最内层循环的矢量化,则失败。
这是循环携带依赖的经典问题。 你的每次迭代都依赖于其他一些迭代(已经完成),因此可以连续地调度它的唯一方式。
但这只是因为你的循环是如何编写的。
你提到R[i][j][k]
取决于R[i-1][j][k]
, R[i][j-1][k]
, R[i][j][k-1]
。 我看到三个依赖 –
- [1,0,0]
- [0,1,0]
- [0,0,1]
我希望这种表述很直观。
对于你现在的场景,依赖1)和2)不是问题,因为k
有0
, i
/ j
有1
,这意味着迭代不依赖于k的先前迭代来完成这两个依赖。
问题是因为3)。 由于k
有1
,因此每次迭代都取决于它之前的迭代。 如果以某种方式我们能够在i
/ j
带一个>0
的数字,我们就会完成。 循环偏斜变换让我们完全相同。
3D示例有点难以理解。 那么让我们看看i
和j
2D示例。
假设 – R[i][j]
取决于R[i-1][j]
和R[i][j-1]
。 我们有同样的问题。
如果我们必须在图片中表示它,它看起来像这样 –
. <- . <- . | | vv . <- . <- . | | vv . . .
在该图中,每个点表示和迭代(i,j)
和源自每个点的箭头指向它所依赖的迭代。 很明显,为什么我们不能在这里并行化最内层的循环。
但是假设我们做了偏差 -
. /| / | . . /| /| / | / | . . . /| / | . . .
如果你绘制与上图相同的箭头(我无法在ASCII艺术中绘制对角线箭头)。
您将看到所有箭头都指向下方,即它们至少继续向下迭代,这意味着您可以并行化水平循环。
现在说你的新循环尺寸是y
(外循环)和x
(内循环),
你的原始变量i
, j
将是
j = x
, i = 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个维度。
为了进一步了解仿射变换的工作原理以及它们如何用于引入并行性,您可以看到这些讲义 。