c - 具有数据依赖性的 for 循环向量化

Question

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

是否可以更改此循环/算法以便对其进行有效矢量化？

 // 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 的 RLL 通过变量 dummy 在 i-1 、 j-1 和 k-1 合并更新的“RLL”。此外，现在仅使用指令 schedule(static, nx/NTt) 在 x 方向上分解循环，其中 NTt 只是可用线程数的宏。是否可以使用指令 collapse 在所有方向上分解它？

-------- 主要编辑------------------------以下是 Ajay 的回答，这是一个最低限度的工作示例

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<omp.h>

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][m-j] = -1*b[i-1][j][m-j] - 1*b[i][j-1][m-j] -1 * b[i][j][m-j-1] + 4 * b[i][j][m-j];
            }

    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 用于强制执行最内层循环的矢量化，则会失败。

Answer 1

这是循环携带依赖的经典问题。您的每个迭代都依赖于其他一些迭代(完成)，因此唯一可以安排的方式是串行。

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

你提到 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) 不是问题，因为 k 中有一个 0 并且有 1在i/j中，这意味着对于这两个依赖，本次迭代并不依赖k的前几次迭代来完成。

问题是因为3)。由于 k 中有一个 1，因此每次迭代都取决于它的前一次迭代。如果我们能够以某种方式在 i/j 中输入一个数字 >0，我们就完成了。循环倾斜变换让我们可以做完全相同的事情。

3D 例子有点难理解。因此，让我们看一下带有 i 和 j 的二维示例。

假设 - R[i][j] 依赖于 R[i-1][j] 和 R[i][j-1] 。我们有同样的问题。

如果我们必须用图片表示它，它看起来像这样 -

. <- . <- .
     |    |
     v    v
. <- . <- .
     |    |
     v    v
.    .    .

在这张图中，每个点代表迭代(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 个维度。

要进一步了解仿射变换的工作原理以及它们如何用于引入并行性，您可以查看这些 lecture notes .