由于并发写入，在嵌套循环崩溃中使用 omp parallel？

Question

我正在尝试并行化我使用了一段时间的代码，而没有并行化，并且在一个相对较大的 C++ 项目（在 C++11 上进行的编译）中没有问题。该代码严重依赖 Eigen 包（此处使用 3.3.9 版）。这是代码的简约版本，必须并行化并且我会为此崩溃（我希望在压缩东西时没有引入错误......）：

主函数有两个需要并行化的for循环：

Data_eigensols solve(const VectorXd& nu_l0_in, const int el, const long double delta0l, 
    const long double DPl, const long double alpha, const long double q, const long double sigma_p, 
        const long double resol){

const int Nmmax=5000;

int np, ng, s, s0m;
double nu_p, nu_g;
VectorXi test;
VectorXd nu_p_all, nu_g_all, nu_m_all(Nmmax);
Data_coresolver sols_iter;
Data_eigensols nu_sols;

//----
// Some code that initialize several variables (that I do not show here for clarity).
// This includes initialising freq_max, freq_min, tol, nu_p_all, nu_g_all, deriv_p and deriv_g structures
//----

// Problematic loops that crash with omp parallel but works fine when not using omp
s0m=0;
#pragma omp parallel for collapse(2) num_threads(4) default(shared) private(np, ng)
        for (np=0; np<nu_p_all.size(); np++)
        {
                for (ng=0; ng<nu_g_all.size();ng++)
                {
                        nu_p=nu_p_all[np];
                        nu_g=nu_g_all[ng];
   
                        sols_iter=solver_mm(nu_p, nu_g); // Depends on several other function but for clarity, I do not show them here (all those 'const long double' or 'const int')
                        //printf("np = %d, ng= %d, threadId = %d \n", np, ng, omp_get_thread_num());
                        for (s=0;s<sols_iter.nu_m.size();s++)
                        {
                                // Cleaning doubles: Assuming exact matches or within a tolerance range
                                if ((sols_iter.nu_m[s] >= freq_min) && (sols_iter.nu_m[s] <= freq_max))
                                {
                                        test=where_dbl(nu_m_all, sols_iter.nu_m[s], tol, 0, s0m); // This function returns -1 if there is no match for the condition (here means that sols_iter.nu_m[s] is not found in nu_m_all)
                                        if (test[0] == -1) // Keep the solution only if it was not pre-existing in nu_m_all
                                        {
                                                nu_m_all[s0m]=sols_iter.nu_m[s];
                                                s0m=s0m+1;
                                        }
                                }
                        }              
                }
        }
    nu_m_all.conservativeResize(s0m); // Reduced the size of nu_m_all to its final size

    nu_sols.nu_m=nu_m_all;
    nu_sols.nu_p=nu_p_all;
    nu_sols.nu_g=nu_g_all;
    nu_sols.dnup=deriv_p.deriv; 
    nu_sols.dPg=deriv_g.deriv;

return, nu_sols;
}

类型Data_coresolver和Data_eigensols定义为：

struct Data_coresolver{
    VectorXd nu_m, ysol, nu,pnu, gnu; 
};

struct Data_eigensols{
    VectorXd nu_p, nu_g, nu_m, dnup, dPg;
};

where_dbl()如下：

VectorXi where_dbl(const VectorXd& vec, double value, const double tolerance){
/*
 * Gives the indexes of values of an array that match the value.
 * A tolerance parameter allows you to control how close the match
 * is considered as acceptable. The tolerance is in the same unit
 * as the value
 *
*/
   int cpt;
   VectorXi index_out;
  
   index_out.resize(vec.size());
    
    cpt=0;
    for(int i=0; i<vec.size(); i++){
        if(vec[i] > value - tolerance && vec[i] < value + tolerance){
            index_out[cpt]=i;
            cpt=cpt+1;
        }       
    }
    if(cpt >=1){
        index_out.conservativeResize(cpt);
    } else{
        index_out.resize(1);
        index_out[0]=-1;
    }
    return index_out;
}

关于 solver_mm(): 我没有详细说明这个函数，因为它调用的子程序很少，而且可能太长，无法在此处显示，而且我认为它与此处无关。它基本上是一个搜索求解隐式方程的函数。

主要功能应该做什么： 主函数 solve() 迭代调用 solver_mm() 以求解不同条件下的隐式方程，其中唯一的变量是 nu_p 和 nu_g。有时一对 (nu_p(i),nu_g(j)) 的解决方案与另一对 (nu_p(k), nu_g(l)) 导致重复的解决方案.这就是为什么有一个部分调用 where_dbl() 来检测那些重复的解决方案并抛出它们，只保留唯一的解决方案。

问题是什么： 没有#pragma 调用，代码工作正常。但是它在随机执行时失败了。经过几次测试，罪魁祸首似乎与删除重复解决方案的部分有关。我的猜测是 nu_m_all VectorXd 上有并发写入。我尝试使用 #pragma omp barrier 但没有成功。但我对 omp 很陌生，我可能误解了屏障的工作原理。

谁能告诉我为什么我会在这里崩溃以及如何解决？对于一些在 omp 方面有丰富经验的人来说，解决方案可能是显而易见的。

Answer 1

nu_p、nu_g、sols_iter 和 test 应该是私有的。由于这些变量被声明为共享变量，多个线程可能以非线程安全的方式写入同一内存区域。这可能是你的问题。

由于并发写入，在嵌套循环崩溃中使用 omp parallel？

Using omp parallel in a nested loop crash due to concurrent writing?

c++

c++11

openmp

eigen3