为何OpenMP任务中无法使用reduction？附曼德博集合代码问题

问题描述 #

我有一段计算曼德博集合面积数值近似的代码，在声明OpenMP并行区域时，尝试用reduction(+ : numoutside)替代shared(numoutside)会得到错误结果，想知道这背后的原因。

代码示例 #

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

#define NPOINTS 2000
#define MAXITER 2000

struct complex
{
    double real;
    double imag;
};

int main()
{
    int i, j, iter, numoutside = 0;
    double area, error, ztemp;
    struct complex z, c;

    double time1, time2, elapsed;

    /*
     *
     *
     *     Outer loops run over npoints, initialise z=c
     *
     *     Inner loop has the iteration z=z*z+c, and threshold test
     */

    time1 = omp_get_wtime();

#pragma omp parallel default(none) private(i, j, iter, c, z, ztemp) \
    shared(numoutside) // why not reduction(+ : numoutside) ???
    {
#pragma omp single
        {
            for (i = 0; i < NPOINTS; i++)
            {
#pragma omp task
                for (j = 0; j < NPOINTS; j++)
                {
                    // #pragma omp task - LOSE NE TREBA OVDE - NAPRAVICE SE PREVISE
                    // TASKOVA (NPOINTS * NPOINTS) I USPORICEMO SISTEM UMESTO DA UBRZAMO
                    {
                        c.real = -2.0 + 2.5 * (double)(i) / (double)(NPOINTS) + 1.0e-7;
                        c.imag = 1.125 * (double)(j) / (double)(NPOINTS) + 1.0e-7;
                        z = c;
                        for (iter = 0; iter < MAXITER; iter++)
                        {
                            ztemp = (z.real * z.real) - (z.imag * z.imag) + c.real;
                            z.imag = z.real * z.imag * 2 + c.imag;
                            z.real = ztemp;
                            if ((z.real * z.real + z.imag * z.imag) > 4.0e0)
                            {
#pragma omp atomic
                                numoutside++;
                                break;
                            }
                        }
                    }
                }
            }
        }
    }

    time2 = omp_get_wtime();
    elapsed = time2 - time1;
    printf("elapsed %f\n", elapsed);

    /*
     *  Calculate area and error and output the results
     */

    area = 2.0 * 2.5 * 1.125 * (double)(NPOINTS * NPOINTS - numoutside) / (double)(NPOINTS * NPOINTS);
    error = area / (double)NPOINTS;

    printf("Area of Mandlebrot set = %12.8f +/- %12.8f\n", area, error);
}

原因解释 #

1. Reduction的生命周期与Task异步执行冲突 #

OpenMP的reduction(+ : var)机制逻辑是：并行区域启动时，为每个参与执行的线程创建var的私有副本；并行区域结束时，自动将所有私有副本的值累加合并到原始变量。

但你的代码里，#pragma omp single块创建的是异步任务（task）——这些任务会被放入OpenMP的任务队列，由线程池中的线程调度执行，执行时间不确定，甚至可能在parallel区域结束后才完成。当parallel区域退出时，reduction的合并操作会立即执行，此时大部分task还没跑完，它们对私有副本的修改根本没被计入合并，最终numoutside的值不完整，自然得到错误结果。

2. Shared+Atomic的适配逻辑 #

用shared(numoutside)配合#pragma omp atomic时，所有task直接操作同一个共享变量，原子操作保证了numoutside++自增动作的完整性，不会出现多线程同时修改导致的竞争条件。不管task什么时候执行，只要触发自增就会正确累加，直到所有task完成，numoutside的值就是准确的。

3. 若要使用Reduction的正确姿势 #

如果一定要用reduction，不能直接在parallel区域声明。可以改为每个task内部统计自己的局部计数，最后再将局部值合并到全局变量：

// 修改task内部逻辑
#pragma omp task private(local_count)
{
    int local_count = 0;
    for (j = 0; j < NPOINTS; j++)
    {
        // ... 原有计算逻辑 ...
        if (/* 条件满足 */)
        {
            local_count++;
            break;
        }
    }
    #pragma omp atomic
    numoutside += local_count;
}

这样既利用局部变量避免了频繁原子操作的开销，又保证了最终结果的正确性，本质还是通过原子操作完成合并，和直接用shared+atomic的核心逻辑一致。

内容的提问来源于stack exchange，提问作者Pavle Šarenac