问题描述
我正在尝试并行化我使用的代码，因为它没有并行化，并且在相对较大的 C++ 项目（在 C++11 上编译）中没有问题。代码严重依赖 Eigen 包（此处使用 3.3.9 版）。 这是代码的简约版本，它必须并行化并且我会崩溃（我希望在压缩事物时不会引入错误......）：
需要并行化的两个for循环的main函数：
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,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,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],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 屏障但没有成功。但我对 omp 很陌生，我可能误解了屏障的工作原理。
谁能告诉我为什么我在这里崩溃以及如何解决它？对于一些在 omp 方面具有良好经验的人来说，解决方案可能是显而易见的。
解决方法
nu_p、nu_g、sols_iter 和 test 应该是私有的。由于这些变量被声明为共享的，多个线程可能会以非线程安全的方式在同一内存区域中写入。这可能是您的问题。