c++eigenstable-sort

`std::stable_sort` gives wrong results when `std::execution::par`

I wrote simple algorithm for sorting rows in Eigen matrix. This should do the same as Matlab's sortrows function:

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;
}

If I set std::execution::par i get wrong result like in the picture: enter image description here

With std::execution::seq and the same data the graph non-decreasingly grows in steps (correct result).

What should I know about execution policy to avoid such situations?

EDIT: my implementation for sort_rows that now works with std::execution::par and doesn't use recursion anymore:

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;
}
Solution

  • Here is my implementation of rowsort. I find the documentation of rowsort somewhat confusing. I work under the assumption that it is just a lexicographical sort.

    Note that your code can probably be fixed just by making a col variable local to your lambda instead of having it as a shared reference.

    template<class Derived>
void rowsort(Eigen::MatrixBase<Derived>& mat)
{
    using PermutationMatrix = 
          Eigen::PermutationMatrix<Derived::RowsAtCompileTime>;
    PermutationMatrix permut;
    permut.setIdentity(mat.rows());
    auto& indices = permut.indices();
    std::stable_sort(std::execution::par, indices.begin(), indices.end(),
          [&mat](Eigen::Index left, Eigen::Index right) noexcept -> bool          
          {
              const auto& leftrow = mat.row(left);
              const auto& rightrow = mat.row(right);
              for(Eigen::Index col = 0, cols = mat.cols();
                    col < cols; ++col) {
                  const auto& leftval = leftrow[col];
                  const auto& rightval = rightrow[col];
                  if(leftval < rightval)
                      return true;
                  if(leftval > rightval)
                      break;
              }
              return false;
          });
    mat = permut.inverse() * mat;
}

    Notes: