如何在 openmp 中并行化 while 循环 - 共轭梯度

How to parallelize a while loop in openmp - conjugate gradient

我有一个 while 循环并希望使用 OpenMP 在 2 个线程上将其并行化。循环内的变量不依赖于它们在前一次迭代中的值,所以我认为必须有某种方法来并行化它。我有 2 个线程,因此每次都可能同时发生 2 次 while 循环迭代,每个循环执行自己的计算。这个循环的目标是找到 alfa 的值,这是共轭梯度法寻找最优点的步长。

我想我必须以某种方式利用 alfaalfaIter 变量和 OpenMP 语句来使这个并行循环工作,但不知道如何。

    #pragma omp parallel num_threads(2)
    {
    while (alfaSet == false) {
        alfaIter++;
        alfa = pow(gamma, alfaIter - 1);

        b = 0;

        for (i = 0; i < dim; i++) {
            testX[i] = x[i] + alfa * d[i];
        }
        for (i = 0; i < dim; i++) {
            b += d[i] * g[i];
        }
        if (shanno(testX, dim) - shanno(x, dim) <= delta * alfa * b) {
            alfaIter = 0;
            alfaSet = true;
        }
    }
    }

编辑 1:这个实现似乎没问题:

    #pragma omp parallel num_threads(alfaThreads)
    {
    int alfaIter = omp_get_num_threads();
    int step = omp_get_num_threads();
    double localAlfa = alfa;
    double *testX = (double *) malloc(dim * sizeof(double));
    while (!alfaSet) {
        #pragma omp barrier
        alfaIter += step;
        localAlfa = pow(gamma, alfaIter - 1);
        for (i = 0; i < dim; i++) {
            testX[i] = x[i] + localAlfa * d[i];
        }
        if (func(testX, dim) - delta * localAlfa * b <= oldFunc) {
            #pragma omp critical
            {
                if (!alfaSet) {
                    alfaSet = true;
                    alfaIter = 0;
                    alfa = localAlfa;
                }
            }
        } 
    }
    free(testX);
    }

所以在玩了这段代码一段时间后,我发现没有任何同步,所以线程不会互相等待,它们以不可预测的方式到达部分代码。 OpenMP 屏障现在同步它们,我总是获得相同的迭代次数和性能增益。但是,现在有时程序会崩溃。死锁?如何检查导致崩溃的原因以及如何防止崩溃?

下面是算法的完整实现:​​https://gist.github.com/mazury/394adc82ab51ce36acfae297ae0555ce

#pragma omp parallel 运行s 以下代码在多个线程上并行执行。因此,多个循环将同时 运行。所有这些版本都会获取全局变量并或多或少同时更新它们,您不能简单地控制发生的事情。

例如,很可能以不受控制的方式修改 alfaIter 导致未定义的行为。

这是处理器如何执行代码的第一行

1 read alfaIter in local var (register)
2 var++
3 write register var in alfaIter
4 fetch alfaIter to call pow and put it in stack or register
5 call pow(...)

假设这些指令在线程 A 中是 1a 2a 3a 4a 5a,在线程 B 中是 1b 2b 3b 4b 5b。

现在实际的执行顺序是什么?

假设是

1a 2a 3a 4a 5a 1b 2b 3b 4b 5b. 

该行为将符合预期。 Pow 在线程 A 中以 alfaIter=1 调用,在线程 B 中以 alfaIter=2

调用

但其他顺序可能会导致不同的行为

例如

1a 1b (both local regs in thrd A and B have initial value of 0)
2a 3a (alfaIter=1 written back to memory by thead A)
2b 3b (alfaIter=1 written back to memory by thead B)
4a 4b 5a 5c (pow called by both threads with the same value of alfaIter=1)

由于 任何 排序都是可能的,因此循环的行为是不可预测的。

一种使其可预测的解决方案是通过原子操作。 在这种情况下,您可以确保对内存的访问是顺序的,并且 while 循环的行为符合预期。

但这有一个主要缺点。原子操作 非常 长,在现代处理器上通常需要大约 100 个周期。这将显着降低您的代码速度,并使其 比顺序版本慢 很多。

一般来说最好是用for循环,不过你好像不行

我的建议是将大多数 var 呈现为局部变量,运行 将 alfaIter 递增 2(或线程数)的并行线程,并仅对终止条件使用全局操作。

示例代码:

#pragma omp parallel num_threads(2)
{
  int alfaIter=omp_get_thread_num();
  int step=omp_get_num_threads();
  float alfa;
  float testX[dim],b; 
      // and maybe d[] and g[] but I do not understand what they do
  while (alfaSet == false) { // no problem to read a global var
    alfaIter+=step;
    alfa = pow(gamma, alfaIter - 1);
    b = 0;
    for (i = 0; i < dim; i++) {
        testX[i] = x[i] + alfa * d[i];
    }
    for (i = 0; i < dim; i++) {
        b += d[i] * g[i];
    }
    if (shanno(testX, dim) - shanno(x, dim) <= delta * alfa * b) {
     #pragma omp critical
      if (! alfaSet) { // you can do safe operations here
        alfaIter = 0;
        alfaSet = true;
      }
    }
  }
} 

未经测试,但可以作为起点。