The following MWE shows a failure when trying to write a function that makes the product between two objects. Specifically, it fails when making a product between a matrix.ldlt().solve() return value. If one uses the decltype instead, it works. I just wanted to make sure it is a bug you guys are aware of, since I don't have a bug tracker login.

Code: Select all

#include <Eigen/Dense>
#include <typeinfo>
#include <iostream>

namespace{
    template<class Lhs, class Rhs>
    auto mkprod(const Lhs& lhs, const Rhs& rhs)
    -> Eigen::ProductReturnType<Lhs,Rhs>::Type
    //-> decltype(lhs*rhs)
    { return lhs*rhs; }
}

int main(int,char*[]){
    Eigen::MatrixXd B(Eigen::MatrixXd::Random(5,5));
    Eigen::MatrixXd Id(Eigen::MatrixXd::Identity(5,5));
    std::cout << mkprod(B,Id.ldlt().solve(B)) << std::endl;
    return 0;
}

The same remark applies with the other operators ( matrix + and -, Scalar * etc.)

This is due to internal details, the object returned by solve() gets evaluated when casted to a MatrixBase<>, this make the actual return type of lhs*rhs more complicated. Fix:

Code: Select all

namespace{
    template<class Lhs, class Rhs>
    typename Eigen::ProductReturnType<Lhs,Rhs>::Type
    mkprod(const Eigen::MatrixBase<Lhs>& lhs, const Eigen::MatrixBase<Rhs>& rhs)
    { return lhs.derived() * rhs.derived(); }
}