3
votes

I want to use Eigen::Ref to have non-template functions using Eigen::Matrix arguments. My problem is that in these functions, I may have to resize the matrices referenced by the Eigen::Ref. I understand that for generality an Eigen::Ref should not be resized because it can map to an expression or a matrix block, but In my case, I am sure that what is behind my Eigen::Ref is an Eigen::Matrix.

To illustrate this:

#include "Eigen/Dense"

void add(Eigen::Ref<Eigen::MatrixXd> M, const Eigen::Ref<const Eigen::MatrixXd> &A, const Eigen::Ref<const Eigen::MatrixXd> &B) {
  M=A+B;
}

int main() {
  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor, 32, 32> M(2,3);
  Eigen::Matrix<double, 2, 2> A;
  Eigen::Matrix<double, 2, 2> B;

  add(M,A,B);
}

gives at runtime:

void Eigen::DenseBase<Derived>::resize(Eigen::Index, Eigen::Index) [with Derived = Eigen::Ref<Eigen::Matrix<double, -1, -1> >; Eigen::Index = long int]: Assertion `rows == this->rows() && cols == this->cols() && "DenseBase::resize() does not actually allow to resize."' failed.

I tried to cheat it:

void add(Eigen::Ref<Eigen::MatrixXd> M, const Eigen::Ref<const Eigen::MatrixXd> &A, const Eigen::Ref<const Eigen::MatrixXd> &B) {
  Eigen::Ref<Eigen::Matrix<double,2,2>> MM(M);
  MM=A+B;
}

but I got at runtime:

Eigen::internal::variable_if_dynamic<T, Value>::variable_if_dynamic(T) [with T = long int; int Value = 2]: Assertion `v == T(Value)' failed.

So, what could I do to handle this? In the Eigen documentation, the problem of resizing is addressed for template functions using MatrixBase as arguments, but for Eigen::Ref?

1
If you are sure that it references an Eigen::MatrixXd, why not pass Eigen::MatrixXd & M?chtz
Sorry, did not completely read your question. Assuming the 32 is not always the same, you could pass that as a template parameter. Otherwise, I don't see a safe way to implement this.chtz
I don't want a template function, because I want it to be virtual...janou195
The MaxRowsAtCompileTime and MaxColsAtCompileTime of M in main will not be stored inside the Ref object. And the actual number of rows and cols will be copied into the Ref. So there is no way to safely resize the Ref object. Will the M object always have a fixed maximum size at compile time?chtz
No, M do not always have the same maximum size. Actually I can even have Eigen::Matrix<double, NR, NC> as M, but in this case I am sure that no resizing will occur in my equivalent of the add function. I understand that there is no "safe" way to resize the Eigen::Ref, but I would like to tell Eigen to trust me!janou195

1 Answers

2
votes

Here is a hacky solution using member function pointers and brutal casting:

#include <iostream>
#include <Eigen/Core>
template<class MatrixType>
struct ResizableRef
{
    typedef typename MatrixType::Scalar Scalar;
    class MatrixDummy;
    typedef void (MatrixDummy::*ResizeFunctor)(Eigen::Index rows, Eigen::Index Cols);
    typedef Scalar* (MatrixDummy::*DataGetter)();

    MatrixDummy *m;
    const ResizeFunctor resizer;
    const DataGetter getData;

    template<class Derived>
    ResizableRef(Eigen::MatrixBase<Derived>& M)
      : m(reinterpret_cast<MatrixDummy*>(&M))
      , resizer(reinterpret_cast<ResizeFunctor>((void (Derived::*)(Eigen::Index, Eigen::Index)) &Derived::resize))
      , getData(reinterpret_cast<DataGetter>((Scalar* (Derived::*)()) &Derived::data))
    { }

    template<class Derived>
    ResizableRef& operator=(const Eigen::EigenBase<Derived>& other)
    {
        (m->*resizer)(other.rows(), other.cols());
        MatrixType::Map((m->*getData)(), other.rows(), other.cols()) = other;
    }
};

void foo(ResizableRef<Eigen::MatrixXd> A)
{
    A = Eigen::Matrix2d::Identity();
}

int main(int argc, char *argv[])
{
    using namespace Eigen;
    MatrixXd A;
    Matrix<double, Dynamic, Dynamic, Eigen::ColMajor, 20, 12> B;
    Matrix<double, 2, Dynamic, Eigen::ColMajor, 2, 4> C;
    Matrix2d D;
    foo(A);
    foo(B);
    foo(C);
    foo(D);

    std::cout << A << "\n\n" << B << "\n\n" << C << "\n\n" << D << '\n';
}

This likely breaks strict aliasing rules and I would generally advise to re-think your design. It should however work without unnecessary run-time allocations, and it is safe against some wrong usages:

    MatrixXf fail1;
    Matrix3d fail2;
    foo(fail1); // fails at compile time
    foo(fail2); // fails at runtime