c - 多个线程的 openmp for loop 错误结果

Question

double compute_inner_contraction_0( double* X, double* T, int n, int f, int m, int e, int No, int Nv )
{
    double contraction_value = 0.0;
    int l_begin = 0;
    int l_end   = m;
    for ( int l = l_begin; l < l_end; l++ )
    {
        int k_begin = 0;
        int k_end   = l;
        for ( int k = k_begin; k < k_end; k++ )
        {
            contraction_value += (-1)*X[n * Nv * No * No * No * No + e * No * No * No * No + m * No * No * No + l * No * No + l * No + k]*T[k * Nv + f];
        }
    }
    return contraction_value;
}

void compute_contraction_0( double* X, double* T, int No, int Nv, double* ExEx_block , int nthd )
{
    omp_set_dynamic(0);
    omp_set_num_threads(nthd);
    #pragma omp parallel for
    for ( int n = 0; n < No; n++ )
    {
        int f_begin = 0;
        int f_end   = Nv;
        for ( int f = f_begin; f < f_end; f++ )
        {
            int m_begin = 0;
            int m_end   = n;
            for ( int m = m_begin; m < m_end; m++ )
            {
                int e_begin = 0;
                int e_end   = f;
                for ( int e = e_begin; e < e_end; e++ )
                {
                    ExEx_block[n * Nv * No * Nv + f * No * Nv + m * Nv + e] += compute_inner_contraction_0( X, T, n, f, m, e, No, Nv);
                }
            }
        }
    }
}

大家好，抱歉函数太长了。请不要担心内存泄漏的索引，因为它们已经过广泛测试。基本上它是将数组 X 和 T 中的数字相乘，然后将结果添加到另一个数组 ExEx_block。此代码将在多线程上运行，但与在一个线程上运行的结果(正确)相比，会给出非常错误的计算结果。能否请你帮忙？它被构建为一个共享库(具有许多其他类似的功能)并从 Python 包装器加载。操作系统为CentOS 6，编译器为gcc-4.4.7它是一个具有 480 GB DDR3 的四路至强 E5-4620 (Sandy-EP)(总共 32 个内核)服务器。不，Nv 只是计算中需要的整数。 nthd(线程数)通常只有 10 或 20。这些数组的大小可以达到 17 GB。非常感谢您的宝贵时间。

Answer 1

我已经分析了您的代码，您在更新 ExEx_block 时遇到了竞争条件。没有其他可能的方法可以得到错误的结果 [因为其他一切都只是读取数据]。

为了让我们达成共识，我已经解释了我所说的竞争条件。

我还为您的程序制作了一个经过清理和简化的版本，没有任何修复。但是，它为固定版本奠定了基础。

并且，第二个版本包含一些示例 omp 指令，我认为您需要(例如 omp private 和 omp atomic)

请注意，我没有尝试编译这些示例，因此仅将它们视为引用。

---

在您的程序中，如果两个不同的线程尝试同时更新 ExEx_block[i]，就会出现竞争条件。为了简化，让我们使用 int 标量代替(例如)target。原理是一样的。

假设 target 的值为 5。我们有两个线程 Ta 和 Tb。 Ta 想要执行 target += 3 而 Tb 想要执行 target += 7。还假设每个线程分别有一个私有(private)变量ta和tb。

Ta Action 实际上是三个操作:

ta = target;
ta += 3
target = ta;

同样适用于 Tb:

tb = target;
tb += 7;
target = tb;

只要它们不相交，一切都很好，我们将 [最终] 获得 15 的 target 值。

但是，如果任何操作被穿插，我们可能会得到 8 或 12，而不是正确的值 15:

ta = target;
tb = target;

ta += 3;
target = ta;

tb += 7;
target = tb;

---

程序的简化版本。

请注意，我用括号对事物进行了分组，并使用了诸如 int No4 = No * No * No * No 之类的事物来简化方程式并使其更具可读性。不是严格要求，但它帮助我发现问题和解决方案。我试着小心点，但我很容易搞砸了一个分组。

double
compute_inner_contraction_0(double *X,double *T,int n,int f,int m,int e,
    int No,int Nv)
{
    double contraction_value = 0.0;
    int l_begin = 0;
    int l_end = m;
    int No2 = No * No;
    int No3 = No2 * No;
    int No4 = No3 * No;

    for (int l = l_begin; l < l_end; l++) {
        int k_begin = 0;
        int k_end = l;

        for (int k = k_begin; k < k_end; k++) {
            contraction_value += (-1) *
                X[(n * Nv * No4) + (e * No4) + (m * No3) + (l * No2) +
                (l * No) + k] * T[(k * Nv) + f];
        }
    }

    return contraction_value;
}

void
compute_contraction_0(double *X,double *T,int No,int Nv,double *ExEx_block,
    int nthd)
{
    omp_set_dynamic(0);
    omp_set_num_threads(nthd);
    int NoNv = Nv * No;
    int NoNv2 = NoNv * Nv;

    #pragma omp parallel for
    for (int n = 0; n < No; n++) {
        int f_begin = 0;
        int f_end = Nv;
        int nx = n * NoNv2;

        for (int f = f_begin; f < f_end; f++) {
            int m_begin = 0;
            int m_end = n;
            int fx = f * NoNv;

            for (int m = m_begin; m < m_end; m++) {
                int e_begin = 0;
                int e_end = f;
                int mx = m * Nv;
                int ax = nx + fx + mx;

                for (int e = e_begin; e < e_end; e++) {
                    ExEx_block[ax + e] +=
                        compute_inner_contraction_0(X,T,n,f,m,e,No,Nv);
                }
            }
        }
    }
}

---

您需要采取哪些措施来防止竞争条件的示例。

为了提高效率，我将 e 循环拆分为两个循环。一种执行繁重的计算，另一种更新全局数组。这有助于最大限度地减少线程停顿，在繁重的计算阶段相互等待。 “锁定”更新阶段是一个简单快速的循环。

您需要一个使用 omp 的 private 指令的每线程本地数组，并且需要将 omp atomic 添加到更新 ExEx_block 的循环中.

警告:虽然我已经使用 pthreads 进行了大量多线程编程并且知道如何修复竞争条件等，但我对特定的等效 omp 指令。

因此，根据我对 omp 文档的阅读，以下是我对您需要什么的最佳猜测。所以，请......彻底检查文档，了解我在这里使用的指令。

double
compute_inner_contraction_0(double *X,double *T,int n,int f,int m,int e,
    int No,int Nv)
{
    double contraction_value = 0.0;
    int l_begin = 0;
    int l_end = m;
    int No2 = No * No;
    int No3 = No2 * No;
    int No4 = No3 * No;

    for (int l = l_begin; l < l_end; l++) {
        int k_begin = 0;
        int k_end = l;

        for (int k = k_begin; k < k_end; k++) {
            contraction_value += (-1) *
                X[(n * Nv * No4) + (e * No4) + (m * No3) + (l * No2) +
                (l * No) + k] * T[(k * Nv) + f];
        }
    }

    return contraction_value;
}

void
compute_contraction_0(double *X,double *T,int No,int Nv,double *ExEx_block,
    int nthd)
{
    omp_set_dynamic(0);
    omp_set_num_threads(nthd);
    int NoNv = Nv * No;
    int NoNv2 = NoNv * Nv;

    // ExEx_local must be at least Nv:
    // NOTE: placement here may be wrong
    double ExEx_local[Nv];

    #pragma omp parallel for private(ExEx_local)
    for (int n = 0; n < No; n++) {
        int f_begin = 0;
        int f_end = Nv;
        int nx = n * NoNv2;

        for (int f = f_begin; f < f_end; f++) {
            int m_begin = 0;
            int m_end = n;
            int fx = f * NoNv;

            for (int m = m_begin; m < m_end; m++) {
                int e_begin = 0;
                int e_end = f;
                int mx = m * Nv;
                int ax = nx + fx + mx;

                // separate computation from global update
                for (int e = e_begin; e < e_end; e++) {
                    ExEx_local[e] =
                        compute_inner_contraction_0(X,T,n,f,m,e,No,Nv);
                }

                // now do global update
                for (int e = e_begin; e < e_end; e++) {
                    #pragma omp atomic update
                    ExEx_block[ax + e] += ExEx_local[e];
                }
            }
        }
    }
}