如何解决与 OpenMP 和 Eigen 进行数据竞争，但数据是只读的？
我试图在我的代码中找到数据竞争，但我似乎无法理解为什么会发生这种情况。线程中的数据是只读的，唯一写入的变量受临界区保护。
我尝试使用 Intel Inspector，但我正在使用 g++ 9.3.0 进行编译，显然即使是 2021 版本也无法处理它的 OpenMP 实现。发行说明没有像旧版本那样明确将其声明为例外，但有一个关于误报的警告，因为它不受支持。它也总是显示 pragma 语句的数据竞争，这根本没有帮助。
我目前的怀疑是 Eigen 或我使用了对 std::vector 的引用。 Eigen 本身我用 EIGEN_DONT_ParaLLELIZE 编译以避免与嵌套并行性混淆，尽管我认为我不会使用任何会使用它的东西。
编辑： 不确定这是否真的是“数据竞争”（或错误的内存访问？），但该示例以相同输入的结果不同的形式产生非确定性输出。如果发生这种情况，则主要中断中的循环。对于多个线程，这种情况会提前发生（通常在 5-12 次迭代之后）。如果我只用一个线程运行它或者不使用 OpenMP 编译，我必须手动结束示例程序。
下面的最小（非）工作示例。
#include <Eigen/Dense>
#include <vector>
#include <iostream>
#ifdef _OPENMP
#include <omp.h>
#else
#define omp_set_num_threads(number)
#endif
typedef Eigen::Matrix<double,9,1> Vector9d;
typedef std::vector<Vector9d,Eigen::aligned_allocator<Vector9d>> Vector9dList;
Vector9d derivPath(const Vector9dList& pathPositions,int index){
int n = pathPositions.size()-1;
if(index >= 0 && index < n+1){
// path is one point,no derivative possible
if(n == 0){
return Vector9d::Zero();
}
else if(index == n){
return Vector9d::Zero();
}
// path is a line,derivative is in the direction of start to end
else {
return n * (pathPositions[index+1] - pathPositions[index]);
}
}
else{
return Vector9d::Zero();
}
}
// ********************************
// data race occurs here somewhere
double errorFunc(const Vector9dList& pathPositions){
int n = pathPositions.size()-1;
double err = 0.0;
#pragma omp parallel default(none) shared(pathPositions,err,n)
{
double err_private = 0;
#pragma omp for schedule(static)
for(int i = 0; i < n+1; ++i){
Vector9d derivX_i = derivPath(pathPositions,i);
// when I replace this with pathPositions[i][0] the loop in the main doesn't break
// (or at least I always had to manually end the program)
// but it does break if I use derivX_i[0];
double err_i = derivX_i.norm();
err_private = err_private + err_i;
}
#pragma omp critical
{
err += err_private;
}
}
err = err / static_cast<double>(n);
return err;
}
// ***************************************
int main(int argc,char **argv){
// setup data
int n = 100;
Vector9dList pathPositions;
pathPositions.reserve(n+1);
double a = 5.0;
double b = 1.0;
double c = 1.0;
Eigen::Vector3d f,u;
f << 0,-1;//-p;
u << 0,1,0;
for(int i = 0; i<n+1; ++i){
double t = static_cast<double>(i)/static_cast<double>(n);
Eigen::Vector3d p;
double x = 2*t*a - a;
double z = -b/(a*a) * x*x + b + c;
p << x,z;
Vector9d cam;
cam << p,f,u;
pathPositions.push_back(cam);
}
omp_set_num_threads(8);
//reference value
double pe = errorFunc(pathPositions);
int i = 0;
do{
double pe_i = errorFunc(pathPositions);
// there is a data race
if(std::abs(pe-pe_i) > std::numeric_limits<double>::epsilon()){
std::cout << "Difference detected at iteration " << i << " diff:" << std::abs(pe-pe_i);
break;
}
i++;
}
while(true);
}
多次运行示例的输出
Difference detected at iteration 13 diff:1.77636e-15
Difference detected at iteration 1 diff:1.77636e-15
Difference detected at iteration 0 diff:1.77636e-15
Difference detected at iteration 0 diff:1.77636e-15
Difference detected at iteration 0 diff:1.77636e-15
Difference detected at iteration 7 diff:1.77636e-15
Difference detected at iteration 8 diff:1.77636e-15
Difference detected at iteration 6 diff:1.77636e-15
如您所见，差异很小，但并不总是发生在同一次迭代中，这使得它具有不确定性。如果我单线程运行它，则没有输出，因为我通常在让它运行几分钟后结束程序。因此，它必须以某种方式与并行化有关。
我知道在这种情况下我可以使用归约，但是在我项目的原始代码中，我还必须计算并行区域中的其他内容，并且我希望使最小示例尽可能接近原始结构。
我在程序的其他部分也使用 OpenMP，我不确定那里是否也存在数据竞争，但结构相似（除了我使用 #pragma omp parallel for 和 collapse 语句） .我有一些变量或 vector 我写入，但它总是在关键区域或每个线程只写入它自己的向量子集。多线程使用的数据始终是只读的。只读数据始终是 std::vector、对 std::vector 的引用或数字数据类型，如 int 或 double。向量始终包含 Eigen 类型或 double。
解决方法
没有竞争条件。您正在观察截断浮点表示的非交换代数的自然结果。当 A、B 和 C 由于舍入误差而为有限精度浮点数时，(A + B) + C 并不总是与 A + (B + C) 相同。二进制中的 1.77636E-15 x 100（注释掉 err = err / static_cast<double>(n); 时的绝对误差）是：
0 | 01010101 | 00000000000000000001100
S exponent mantissa
如您所见，错误位于尾数的最低有效位，暗示它是舍入错误累积的结果。
问题出现在这里：
#pragma omp parallel default(none) shared(pathPositions,err,n)
{
...
#pragma omp critical
{
err += err_private;
}
}
err 的最终值取决于不同线程到达临界区的顺序以及它们的贡献被添加，这就是为什么有时您会立即看到差异，有时需要几次迭代。
为了证明这不是 OpenMP 问题本身，只需将函数修改为：
double errorFunc(const Vector9dList& pathPositions){
int n = pathPositions.size()-1;
double err = 0.0;
std::vector<double> errs(n+1);
#pragma omp parallel default(none) shared(pathPositions,errs,n)
{
#pragma omp for schedule(static)
for(int i = 0; i < n+1; ++i){
Vector9d derivX_i = derivPath(pathPositions,i);
errs[i] = derivX_i.norm();
}
}
for (int i = 0; i < n+1; ++i)
err += errs[i];
err = err / static_cast<double>(n);
return err;
}
这消除了对子和如何计算和相加的依赖，无论 OpenMP 线程的数量如何，返回值都将始终相同。
另一个版本只修复了 err_private 减少为 err 的顺序：
double errorFunc(const Vector9dList& pathPositions){
int n = pathPositions.size()-1;
double err = 0.0;
std::vector<double> errs(omp_get_max_threads());
int nthreads;
#pragma omp parallel default(none) shared(pathPositions,n,nthreads)
{
#pragma omp master
nthreads = omp_get_num_threads();
double err_private = 0;
#pragma omp for schedule(static)
for(int i = 0; i < n+1; ++i){
Vector9d derivX_i = derivPath(pathPositions,i);
double err_i = derivX_i.norm();
err_private = err_private + err_i;
}
errs[omp_get_thread_num()] = err_private;
}
for (int i = 0; i < nthreads; i++)
err += errs[i];
err = err / static_cast<double>(n);
return err;
}
同样，只要线程数保持不变，这段代码每次都会产生相同的结果。该值可能会因线程数不同而略有不同（以 LSB 为单位）。
您无法轻松绕过这种差异，只能学会接受它并采取预防措施以尽量减少其对其余计算的影响。事实上，你真的很幸运能在 2021 年偶然发现它，这是后 x87 时代的一年，那时几乎所有商品 FPU 都使用 64 位 IEEE 754 操作数，而不是在 1990 年代 x87 FPU 使用 80 位操作数和重复累加的结果取决于您是一直将值保存在 FPU 寄存器中，还是定期将其存储在其中，然后从内存中加载回来，这将 80 位表示四舍五入为 64 位。
与此同时，mandatory reading 对于任何在数字计算机上处理数学问题的人来说都是如此。
附言虽然现在是 2021 年，我们已经在后 x87 时代生活了 21 年（从 Pentium 4 在 2000 年引入 SSE2 指令集开始），如果你的 CPU 是 x86 的，你仍然可以参与 x87 的疯狂.只需使用 -mfpmath=387 编译您的代码 :)
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
