在 OpenMP 中使用临界区

Question

我有一段代码想使用 OpenMP 进行并行处理。该程序演示了执行多项式插值的函数。首先，我想到在代码的众多函数之一中并行化一个 for 循环。这个 for 的串行版本是这样的：

void polint (double xa[], double ya[], int n, double x, double *y, double *dy) {
    int i, m, ns=1;
    double den,dif,dift,ho,hp,w;
    double *c, *d;
    dif = fabs(x-xa[1]); c = vector(1,n);
    d = vector(1,n); for (i=1; i<= n; i++) { 
        dift = fabs (x - xa[i]); 
        if (dift<dif) {
            ns = i;
            dif = dift;
        }
        c[i] = ya[i];
        d[i] = ya[i]; 
    }
    *y = ya[ns--];
    for (m = 1; m < n; m++) { 
        for (i = 1; i<= n-m; i++) {
            ho = xa[i] - x;
            hp = xa[i+m] - x;
            w = c[i+1] - d[i]; den = ho - hp;
            den = w / den; d[i] = hp * den; c[i] = ho * den;
        }
    *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--])); 
    }
    free_vector (d, 1, n); free_vector (c, 1, n); 
}

现在，我使用 OpenMP 构建了这段代码：

void polint (double xa[], double ya[], int n, double x, double *y, double *dy) {
    int i, m, ns=1;
    double den,dif,dift,ho,hp,w;
    double *c, *d;
    int nth;
    dif = fabs(x-xa[1]); c = vector(1,n);
    d = vector(1,n); for (i=1; i<= n; i++) { 
        dift = fabs (x - xa[i]); 
        if (dift<dif) {
            ns = i;
            dif = dift;
        }
        c[i] = ya[i];
        d[i] = ya[i]; 
    }
    #pragma omp parallel
    {
    *y = ya[ns--];    
        nth = omp_get_num_threads();
        for (m = 1; m < n; m++) {
            #pragma omp for 
            for (i = 1; i<= n-m; i++) {
                #pragma omp critical 
                {
                    ho = xa[i] - x;
                    hp = xa[i+m] - x;
                    w = c[i+1] - d[i]; den = ho - hp;
                    den = w / den; d[i] = hp * den; c[i] = ho * den;
                }
            }
            #pragma omp critical
            *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--])); 
        }
        free_vector (d, 1, n); free_vector (c, 1, n);
    } 
}

程序编译正确，但是当我运行它时问题就来了。我收到以下错误消息：

我该如何解决这个问题？我一开始以为数组被越界了，所以我添加了临界区，这样多个进程就不会同时访问数组了。然而，这并没有解决问题。提前致谢！

Answer 1

除了导致您遇到崩溃的 c 和 d 指针的明显释放之外，您的代码中还有很多问题。为了更好地指出它们，我用更合理的变量声明和正确的缩进重写了你的顺序函数。这是它给我的：

void polint( double xa[], double ya[], int n, double x, double *y, double *dy ) {
    double *c = vector( 1, n );
    double *d = vector( 1, n );

    double dif = fabs( x - xa[1] );
    int ns = 1;
    for ( int i = 1; i <= n; i++ ) { 
        double dift = fabs( x - xa[i] ); 
        if ( dift < dif ) {
            ns = i;
            dif = dift;
        }
        c[i] = ya[i];
        d[i] = ya[i]; 
    }

    *y = ya[ns--];
    for ( int m = 1; m < n; m++ ) { 
        for ( int i = 1; i <= n - m; i++ ) {
            double ho = xa[i] - x;
            double hp = xa[i + m] - x;
            double den = ( c[i + 1] - d[i] ) / ( ho - hp );
            d[i] = hp * den;
            c[i] = ho * den;
        }
        *dy = 2 * ns < ( n - m ) ? c[ns + 1] : d[ns--];
        *y += *dy; 
    }

    free_vector( d, 1, n );
    free_vector( c, 1, n ); 
}

现在，阅读代码变得更简单了（尽管它仍然很模糊）并且一些问题变得明显：

您声明数组的索引范围是从 1 到 n 而不是从 0 到 n-1？如果是这样，那是不寻常的（并且容易出错）。如果不是，您的循环中可能存在一些问题（特别是第一个 i 循环）。
dy 的值在 m 循环的每次迭代中被重写。知道这是一个输出参数，这意味着它的最终值将只对应于 m 的最后一个值。这不是一个问题，但看起来很可疑......你确定它是正确的吗？
现在，关于可能的并行化，我们可以从den的定义和c的后续更新中看出，c[i]依赖于c[i+1]。从并行化的角度来看，这是一个大问题。这意味着代码现在的编写方式无法并行化。这并没有完全关闭并行化的大门，但这只是意味着除了放置一些 OpenMP 指令来获得它之外，还有更多的工作要做。

现在，除了这些问题，如果我们看一下您的并行化版本，还有很多其他问题：

已经提到，并行区域内矢量的释放：它绝对应该移到外面。
巨大的问题，您没有声明 private 应该是的数据：至少，m、ho、hp、w 和den 应该声明为 private。此外，我强烈怀疑 *y 应该声明 reduction(+) 并且 *dy 应该声明 lastprivate （尽管正确地执行此操作的确切方法需要一些更精细的调整）。
为了弥补缺少 private 声明，您将并行循环的整个主体封装到 critical 区域中，实质上是重新序列化它。不幸的是，它既具有超级感染力又是错误的，因为正如我们所见，迭代的顺序很重要，而这并不能保留它...

总而言之，我强烈建议您做几件事：

这非常重要，去混淆并清理你的代码：正确缩进，为你的变量赋予有意义的名称，每行一个语句，尽可能晚地在代码中声明你的变量并立即初始化它们...在这里和那里的一两个明智的评论可能会派上用场。
打破 c 的 i 依赖关系。这可以很容易地通过保留它的两个副本来完成（一个对应于它在 m 循环的前一次迭代中的值，一个对应于它的当前值）。
添加 OpenMP 指令时，请小心将变量正确声明为 private 等。好消息是，如果你按照我的建议去做并尽可能地延迟它们的声明，你可能很少需要照顾它们，因为大多数声明将在 parallel 区域内，从而使相应的变量自动 private 因为它们应该是。

这样，您应该能够获得正确且（希望）有效的代码并行版本。

在 OpenMP 中使用临界区

Using critical region with OpenMP

c

parallel-processing

openmp