3

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?

janou195
  • 1,175
  • 2
  • 10
  • 25
  • 1
    If you are sure that it references an `Eigen::MatrixXd`, why not pass `Eigen::MatrixXd & M`? – chtz Sep 21 '17 at 09:02
  • 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 Sep 21 '17 at 09:05
  • I don't want a template function, because I want it to be virtual... – janou195 Sep 21 '17 at 09:11
  • 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 Sep 21 '17 at 09:19
  • No, M do not always have the same maximum size. Actually I can even have Eigen::Matrix 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 Sep 21 '17 at 09:35
  • As chtz said, `Ref<>` does not have at hand the necessary information to resize the references objects, and anyway, resizing a MatrixXd or your statically allocated matrix M is not the same. So the only workaround I see would be to write your own little wrapper with a virtual resize member, and a virtual conversion to Ref. – ggael Sep 21 '17 at 10:08
  • Thanks for your comment ggael, but I am not sure to get your workaround: do you suggest a wrapper over the Eigen::Ref class? What do you mean by "virtual conversion to Ref"? – janou195 Sep 21 '17 at 10:19

1 Answers1

2

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
chtz
  • 17,329
  • 4
  • 26
  • 56