我正试图在一个相对较大的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_allVectorXd上存在并发写入。我尝试使用#pragma-omp屏障，但没有成功。但我对omp还很陌生，我可能误解了屏障的工作原理。
有人能告诉我为什么我在这里发生车祸以及如何解决吗？对于一些在omp方面经验丰富的人来说，解决方案可能是显而易见的。
nu_p、nu_g、sols_iter和test应该是私有的。由于这些变量被声明为共享，多个线程可能会以非线程安全的方式在同一内存区域中进行写入。这可能是你的问题。