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 这是一个四核 Xeon E5-4620 (Sandy-EP)（总共 32 个内核）服务器，具有 480 GB 的 DDR3。不，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。还假设每个线程分别有一个私有变量ta和tb。

Ta这个动作其实是三个操作：

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

Tb 也是如此：

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

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

但是，如果任何操作被穿插，我们可能会得到 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];
                }
            }
        }
    }
}