没有线程本地副本的 OpenMP 减少数组

Question

我希望 OpenMP 将 大 数组缩减为较小的动态数组。例如，在哪里

large[] = {1, 1, 1, 2, 2, 2, 3, 3, 3};

// OpenMP reduce length-3 sublists of large, to produce:

small[] = {3, 6, 9};

我的要求类似于this question，但我有重要的附加限制：

• 必须支持 OpenMP 3.1（所以我不能像 this answer 那样使用 OpenMP 4.5 的数组缩减）
• 不能拥有small 的线程私有副本（如this answer），因为small 可以和large 一样大，因此线程私有副本可能会导致堆栈溢出
• large 数组必须按顺序迭代，以获得良好的缓存（因为它可能非常大）

其他一些细节：

• large 的求和元素不一定相邻
• 这种减少发生在函数中。 small（输出）数组由用户预先分配，并通过引用传递。理想情况下，该函数不应分配 small 的新本地副本
• small 和large 的长度都是2 的幂（例如2、4、8、16 ...）。 small 的每个元素都从相同数量的 large 元素减少（2 的另一个幂）。

为了清楚起见，这里是示例串行伪代码：

void myfunc(double* small, double* large, int lLen, // other args) {

    for (int i=0; i<lLen; i++) {
    
        int sInd = // determined by other args
        
        small[sInd] += large[i];
    }
}

这是一个示例实现，由于它利用了 OpenMP 4.5 的数组缩减，因此不合格。此外，它不希望地使用了small 的本地副本（根据使用reduction 的要求）。

void myfunc(double* small, int sLen, double* large, int lLen, // other args) {

    double smallLocal[sLen];
    int i, sInd;

    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        smallLocal[i] = 0;
    
    #pragma omp parallel private (i, sInd) reduction (+:smallLocal) for
    for (i=0; i<largeLen; i++) {

        sInd = // determined by other args
        
        smallLocal[sInd] += large[i];
    }

    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        small[i] = smallLocal[i];
}

我希望我可以在 OpenMP 3.1 中实现相同的目标，甚至可以通过自己管理数组元素减少来避免 smallLocal 的本地副本。我该怎么做？

Answer 1

您可以使用 OpenMPs atomic update 构造，它已经存在于 OpenMP 3.1 中：

void myfunc(double* small, double* large, int lLen, // other args) {

    #pragma omp parallel for
    for (int i=0; i<lLen; i++) {
    
        int sInd = // determined by other args
        #pragma omp atomic update
        small[sInd] += large[i];
    }
}

这应该比使用锁要快。

Answer 2

这是一个使用锁的有点不令人满意的解决方案，它需要一个与 small 一样大的本地 locks 数组（可能需要进行 malloc 处理）。

void myfunc(double* small, int sLen, double* large, int lLen, // other args) {

    int i, sInd;

    // clear small
    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        small[i] = 0;

    // initialise locks
#ifdef _OPENMP
    omp_lock_t locks[sLen];
    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        omp_init_lock(&locks[i]);
#endif 
    
    // main loop
    #pragma omp parallel private (i, sInd) for
    for (i=0; i<largeLen; i++) {

        sInd = // determined by other args
        
        // update small with locks
#ifdef _OPENMP
        omp_set_lock(&locks[sInd]);
        small[sInd] += large[i];
        omp_unset_lock(&locks[sInd]);
#else 
        small[sInd] += large[i];
#endif
    }

    // destroy locks
#ifdef _OPENMP
    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        omp_destroy_lock(&locks[i]);
#endif 
}

正如预期的那样，这比使用 OpenMP 4.5 的数组缩减要慢：

• 对于sLen = 8、lLen = 1048576，使用锁26x慢了
• 对于sLen = 32768、lLen = 1048576，使用锁1.5x慢

需要说明的是，OpenMP 4.5 的数组缩减只兼容本地栈数组，不兼容动态内存。这将缩减数组的大小限制为 ~1 MiB，或 2^17 双倍 (lLen = 131072)。相比之下，这个锁解决方案可以处理任何大小，因为locks 可以是 malloc 的。

Answer 3

你可以做这种事情，这在科学代码中很常见：

#pragma omp parallel for
for(int i = 0; i < N/3; i++)
{
    for (int j = 0; j < 3; j++)
    {
        small[i] += large[3*i + j];
    }
}