与 OpenMP 和 Eigen 进行数据竞争但数据是只读的？
Data race with OpenMP and Eigen but data is read-only?
我试图在我的代码中找到数据竞争，但我似乎无法理解为什么会发生。线程中的数据以只读方式使用，唯一写入的变量受临界区保护。
我尝试使用 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, 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, 0, 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 线程的数量如何，return 值将始终相同。
另一个版本只修正了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, errs, 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 年偶然发现它，那是 post-x87 时代的一年，当时几乎所有商品 FPU 使用 64 位 IEEE 754 操作数，而不是在 1990 年代 x87 FPU 使用 80-位操作数和重复累加的结果将取决于您是将值一直保存在 FPU 寄存器中，还是定期将其存储在内存中，然后从内存中将其加载回来，这会将 80 位表示形式四舍五入为 64 位表示形式.
同时，mandatory reading 任何在数字计算机上处理数学的人。
P.S。虽然现在是 2021 年，我们已经在 post-x87 时代生活了 21 年（从 2000 年 Pentium 4 引入 SSE2 指令集开始），如果您的 CPU 是 x86 时代，你仍然可以参与 x87 的疯狂。只需使用 -mfpmath=387 编译您的代码 :)
我试图在我的代码中找到数据竞争，但我似乎无法理解为什么会发生。线程中的数据以只读方式使用，唯一写入的变量受临界区保护。
我尝试使用 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, 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, 0, 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 线程的数量如何，return 值将始终相同。
另一个版本只修正了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, errs, 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 年偶然发现它，那是 post-x87 时代的一年，当时几乎所有商品 FPU 使用 64 位 IEEE 754 操作数，而不是在 1990 年代 x87 FPU 使用 80-位操作数和重复累加的结果将取决于您是将值一直保存在 FPU 寄存器中，还是定期将其存储在内存中，然后从内存中将其加载回来，这会将 80 位表示形式四舍五入为 64 位表示形式.
同时，mandatory reading 任何在数字计算机上处理数学的人。
P.S。虽然现在是 2021 年，我们已经在 post-x87 时代生活了 21 年（从 2000 年 Pentium 4 引入 SSE2 指令集开始），如果您的 CPU 是 x86 时代，你仍然可以参与 x87 的疯狂。只需使用 -mfpmath=387 编译您的代码 :)