使用std::execution::par时std::stable_sort执行Eigen矩阵行排序返回错误结果的问题

我写了一个和Matlab sortrows 功能一致的Eigen矩阵行排序算法，初始实现代码如下：

template <typename D>
void _sort(
    const D &M,
    Eigen::VectorX<ptrdiff_t>& idx,
    std::function<bool(ptrdiff_t, ptrdiff_t)> cmp_fun)
{
  // initialize original index locations
  idx = Eigen::ArrayX<ptrdiff_t>::LinSpaced(
        M.rows(), 0, M.rows()-1);

  std::stable_sort(std::execution::par, idx.begin(), idx.end(), cmp_fun);
}


/// \brief sort_rows  sorts the rows of a matrix in ascending order
/// based on the elements in the first column. When the first column
/// contains repeated elements, sortrows sorts according to the values
/// in the next column and repeats this behavior for succeeding equal values.
/// M_sorted = M(ind, Eigen::all)
/// \param M
/// \return ind
template <typename D>
Eigen::VectorX<ptrdiff_t> sort_rows(const Eigen::DenseBase<D> &M){

  // initialize original index locations
  Eigen::VectorX<ptrdiff_t> idx;

  ptrdiff_t col = 0;

  std::function<bool(ptrdiff_t, ptrdiff_t)> cmp_fun;
  cmp_fun = [&M, &col, &cmp_fun](
      const ptrdiff_t& row1,
      const ptrdiff_t& row2)->bool
  {
    if (M(row1, col) < M(row2, col)){
      col = 0;
      return true;
    }

    if (M(row1, col) > M(row2, col)){
      col = 0;
      return false;
    }

    // only 'M(row1, col) == M(row2, col)'  option is left
    // it will return 'true' only if this is the last column
    // i.e. all other columns at these rows are also equal
    if (col == M.cols()-1){
      if (M(row1, col) == M(row2, col)){
        col = 0;
        return true;
      }

      col = 0;
      return false;
    }

    col++;
    return cmp_fun(row1, row2);
  };

  _sort(M.derived(), idx, cmp_fun);

  return idx;
}

当指定std::execution::par并行执行策略时，得到了错误的结果；而使用std::execution::seq串行执行策略处理相同数据时，排序结果正确（呈现非递减的阶梯式增长）。

为避免此类情况，需要了解执行策略的哪些要点？

编辑：现在我实现的sort_rows可正常工作于std::execution::par策略下，且不再使用递归，代码如下：

template <typename D>
void _sort(
    const D &M,
    Eigen::VectorX<ptrdiff_t>& idx,
    std::function<bool(ptrdiff_t, ptrdiff_t)> cmp_fun)
{
  // initialize original index locations
  idx = Eigen::ArrayX<ptrdiff_t>::LinSpaced(
        M.rows(), 0, M.rows()-1);

  std::stable_sort(std::execution::par, idx.begin(), idx.end(), cmp_fun);
}


/// \brief sort_rows  sorts the rows of a matrix in ascending order
/// based on the elements in the first column. When the first column
/// contains repeated elements, sortrows sorts according to the values
/// in the next column and repeats this behavior for succeeding equal values.
/// M_sorted = M(ind, Eigen::all)
/// \param M
/// \return ind
template <typename D>
Eigen::VectorX<ptrdiff_t> sort_rows(const Eigen::DenseBase<D> &M){

  // initialize original index locations
  Eigen::VectorX<ptrdiff_t> idx;

  std::function<bool(ptrdiff_t, ptrdiff_t)> cmp_fun;
  cmp_fun = [&M](
      const ptrdiff_t& row1,
      const ptrdiff_t& row2)->bool
  {
    ptrdiff_t N = M.cols()-1;
    for (ptrdiff_t col = 0; col < N; col++){
      if (M(row1, col) < M(row2, col))
        return true;

      if (M(row1, col) > M(row2, col))
        return false;
    }

    // notice the operator is '<=' as it is the last column check
    // i.e. when all other columns are equal at these rows
    if (M(row1, Eigen::last) <= M(row2, Eigen::last))
      return true;

    return false;
  };

  _sort(M.derived(), idx, cmp_fun);

  return idx;
}

需要了解的并行执行策略核心要点 #

• 比较函数必须无状态且线程安全：并行排序时，多个线程会同时调用比较函数。原始代码里的col是共享变量，多线程并发读写会直接打乱比较逻辑。并行场景下，比较函数的所有逻辑必须仅依赖输入参数，不能修改或依赖任何外部共享状态。
• 严格遵循严格弱序规则：不管串行还是并行，排序算法要求比较函数满足严格弱序：
  • 无自反性：cmp(a,a)必须返回false
  • 非对称性：如果cmp(a,b)为true，那么cmp(b,a)必须为false
  • 传递性：如果cmp(a,b)和cmp(b,c)都为true，则cmp(a,c)也必须为true
    原始递归实现不仅有状态问题，递归逻辑也可能在并行环境下破坏这种顺序一致性。
• 限制函数副作用：使用std::execution::par时，标准库会假设用户提供的函数没有数据竞争，也不会产生影响其他线程的副作用。任何共享资源的访问都必须同步，或者直接避免使用共享资源。
• 递归在并行场景的风险：原始代码用递归实现比较逻辑，并行执行时多个线程会同时触发递归调用，既加剧了共享变量的竞争问题，还可能导致调用栈异常或逻辑混乱。改用迭代实现（修复后的代码用循环）更适配并行环境。

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