CODEWITHSUNDEEP
c++matrixsliceeigen
Kaiyakha
Kaiyakha

Reputation: 1913

Pointer to sliced Eigen matrix

I have a function which takes a pointer to an Eigen matrix. I need to slice this matrix to several pieces and pass a portion of the matrix to the function. I want to avoid copying data, so I want to pass it with a pointer.

I slice a matrix, for example like this:

MatrixXd A;
// some operations with the matrix
// now I want something like this:
MatrixXd *B = &A(seq(2, 4), all);

The last line is invalid. How can I get a pointer to a piece of matrix so that it can still be treated as an Eigen matrix? Are there any workarounds? Again, I'd like to avoid creating new matrices as in my real program they are quite big.

Upvotes: 1

Views: 1312

Answers (1)

Homer512
Homer512

Reputation: 13310

That's what Eigen::Ref was built for. See the section "How to write generic, but non-templated function?"

Something like this:


void consume_matrix(const Eigen::Ref<const Eigen::MatrixXd>& in);
void produce_matrix(Eigen::Ref<Eigen::MatrixXd> out);


void foo()
{
    Eigen::MatrixXd matrix = ...;
    produce_matrix(matrix.middleRows(start, len));
    consume_matrix(matrix.leftCols(len));
}

Storing block references

Eigen doesn't come with a type that is well-suited to keep a "borrowed" block expression of a matrix. I wouldn't recommend Ref for this case because it is far too easy to make it point at temporary allocations which can easily lead to dangling pointers.

Ultimately, storing borrowed views is a fundamentally unsafe operation. But if you absolutely have to do it, I would use an Eigen::Map. At least that makes it clear you aren't owning the memory to which you refer. Something like this:

struct MatrixReference
{
  using map_type = Eigen::Map<
    Eigen::MatrixXd, Eigen::Unaligned, Eigen::OuterStride<> >;

  map_type matrix;

  template<class Derived>
  static map_type to_map(const Eigen::DenseBase<Derived>& matrix) noexcept
  {
    assert(matrix.innerStride() == 1);
    return map_type(const_cast<double*>(matrix.derived().data()),
                    matrix.rows(), matrix.cols(),
                    Eigen::OuterStride<>(matrix.outerStride()));
  }

  MatrixReference() noexcept
  : matrix(nullptr, 0, 0, Eigen::OuterStride<>(1))
  {}
  template<class Derived>
  explicit MatrixReference(const Eigen::DenseBase<Derived>& matrix) noexcept
    : matrix(to_map(matrix))
  {}
  MatrixReference(const MatrixReference&) = default;
  MatrixReference& operator=(const MatrixReference& o) noexcept
  {
    // placement-new because assignment would overwrite the content
    new (&matrix) map_type(o.matrix);
    return *this;
  }
};
int main()
{
  Eigen::MatrixXd mat = Eigen::MatrixXd::Random(4, 5);
  MatrixReference ref(mat.topRows(2));
  std::cout << mat << "\n\n" << ref.matrix << "\n\n";
}

A slightly safer option is to keep a shared_ptr to the matrix and store the information what slice you use separately. Maybe like this:

struct CountedBlockRef
{
  std::shared_ptr<Eigen::MatrixXd> matrix;
  Eigen::Index first_row, first_col, rows, cols;

  using block_type = Eigen::Block<
    Eigen::MatrixXd, Eigen::Dynamic, Eigen::Dynamic>;
  using const_block_type = Eigen::Block<
    const Eigen::MatrixXd, Eigen::Dynamic, Eigen::Dynamic>;

  CountedBlockRef() noexcept
  : matrix(), first_row(), first_col(), rows(), cols()
  {}
  CountedBlockRef(std::shared_ptr<Eigen::MatrixXd> matrix,
                  Eigen::Index first_row, Eigen::Index first_col,
                  Eigen::Index rows, Eigen::Index cols) noexcept
    : matrix(std::move(matrix)),
      first_row(first_row),
      first_col(first_col),
      rows(rows),
      cols(cols)
  {}
  block_type block() noexcept
  { return matrix->block(first_row, first_col, rows, cols); }

  const_block_type block() const noexcept
  {
    const Eigen::MatrixXd& matrix = *this->matrix;
    return matrix.block(first_row, first_col, rows, cols);
  }
};
int main()
{
  std::shared_ptr<Eigen::MatrixXd> matptr =
    std::make_shared<Eigen::MatrixXd>(Eigen::MatrixXd::Random(4, 5));
  CountedBlockRef aref(matptr, 0, 0, 2, 5);
  std::cout << *matptr << "\n\n" << aref.block() << '\n';
}

That way, the reference will remain valid even after resizing the matrix and dangling pointers are avoided. Only shrinking remains a threat but that is caught by assertions, assuming you compile with them.

Upvotes: 4

Related Questions