OpenMP与类及Eigen结合时结果异常的原因分析

我正在编写结合OpenMP、自定义类与Eigen库的代码，简化后仅实现矩阵求逆与乘法功能。

未启用OpenMP时的正常输出 #

solution in for loop at 0th iteration is... 1 1 1
 solution in class at 0th iteration is... 1 1 1
 solution in for loop at 1th iteration is... 4 4 4
 solution in class at 1th iteration is... 4 4 4
 solution in for loop at 2th iteration is... 9 9 9
 solution in class at 2th iteration is... 9 9 9
 solution in for loop at 3th iteration is... 16 16 16
 solution in class at 3th iteration is... 16 16 16

启用OpenMP后的混乱输出 #

9   11 91  
 1
16 16
 solution in class at 3th iteration is... 116   16 161
 1 1
 solution in class at 2th iteration is... 1
9 9 9
 solution in for loop at 1th iteration is...  solution in for loop at 9 
13th iteration is... 1
 solution in class at 0th iteration is...  4 4  4 16 solution in for loop at 1 
 solution in for loop at  01 11 
1 solution in for loop at  1 solution in for loop at  41 3th iteration is... th iteration is... 
16 9  solution in for loop at th iteration is... 1 solution in class at 33th iteration is...  th iteration is... 161 
 16  solution in for loop at  016 16 2th iteration is... 1616   1616   16 161616th iteration is... 
 solution in for loop at 19  solution in class at 3th iteration is... 9 9
16 

最初仅用#pragma omp parallel for时出现LDLT未初始化错误，改用#pragma omp parallel for private(j,k) firstprivate(foo)后结果仍混乱，请问原因是什么？

自定义类代码如下：

class Foo {
 private:
  std::vector<Eigen::Matrix<double, 3, 3>> a;
  Eigen::LDLT<Eigen::MatrixXd> Minv; // Minv
  std::vector<Eigen::Matrix<double, 3, 1>> b;
  std::vector<Eigen::Matrix<double, 3, 1>> solution;
 public:
  Foo() {};
  ~Foo() {};
  void Initialization() {
    a.resize(4);
    b.resize(4);
    solution.resize(4);
    for (int i = 0; i < 4; i++) {
      a.at(i).setZero();
      b.at(i).setZero();
      solution.at(i).setZero();
    }
  }

  void SetInv(int idx) {
    a.at(idx) = (1.0/(idx+1.0))*Eigen::Matrix3d::Identity();
    b.at(idx) = (idx+1)*Eigen::Vector3d::Ones();
    Minv.compute(a.at(idx));
  }

  void Calculation(int idx) {
    solution.at(idx) = Minv.solve(b.at(idx));
  }

  Eigen::Matrix<double, 3, 1> &ReturnSolution(int idx) {
    return solution.at(idx);
  }

  void ShowSolution(int idx) {
    std::cout << " solution in class at " << idx << "th iteration is..." << solution.at(idx).transpose() << std::endl;
  }

};

1. Minv成员的线程竞争问题 #

你的Foo类中，Minv是全局类成员变量。当用firstprivate(foo)时，虽然每个线程会复制foo实例，但并行循环中多个线程会同时调用SetInv（修改Minv）和Calculation（读取Minv），而Eigen的LDLT类不是线程安全的——它的compute和solve会修改内部分解数据结构，并发操作会直接破坏Minv的内部状态，导致计算结果错误、未初始化报错。

2. std::cout的输出原子性问题 #

输出混乱的直接原因是多个线程同时调用std::cout。std::cout的单个operator<<是线程安全的，但它不会保证整行输出的原子性，不同线程的输出会被打断、交叉拼接，最终形成混乱的文本。

修复Minv竞争问题 #

• 将Minv改为方法局部变量，让每个迭代拥有独立的分解实例：

  void Calculation(int idx) {
    Eigen::LDLT<Eigen::MatrixXd> Minv_local;
    Minv_local.compute(a.at(idx));
    solution.at(idx) = Minv_local.solve(b.at(idx));
  }
  

  这样每个线程的计算都用自己的LDLT实例，完全避免竞争。

修复输出混乱问题 #

• 在输出时添加OpenMP临界区，保证整行输出的完整性：

  void ShowSolution(int idx) {
    #pragma omp critical
    {
      std::cout << " solution in class at " << idx << "th iteration is..." << solution.at(idx).transpose() << std::endl;
    }
  }
  

内容的提问来源于stack exchange，提问作者Junie