3
#include <iostream>

#include <Eigen/Core>

namespace Eigen {

// float op double -> double
template <typename BinaryOp>
struct ScalarBinaryOpTraits<float, double, BinaryOp> {
  enum { Defined = 1 };
  typedef double ReturnType;
};

// double op float -> double
template <typename BinaryOp>
struct ScalarBinaryOpTraits<double, float, BinaryOp> {
  enum { Defined = 1 };
  typedef double ReturnType;
};

}

int main() {
    Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> m1(2, 2);
    m1 << 1, 2, 3, 4;

    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> m2(2, 2);
    m2 << 1, 2, 3, 4;

    std::cerr << m1 * m2 <<std::endl;  // <- boom!!
}

I'd like to know why the above code does not compile. Here is the full error messages. Please note that if I define m1 and m2 to have fixed sizes, it works fine.

I'm using Eigen3.3.1. It's tested on a Mac running OSX-10.12 with Apple's clang-800.0.42.1.

Soonho Kong
  • 365
  • 1
  • 9

1 Answers1

4

This is because the general matrix-matrix product is highly optimized with aggressive manual vectorization, pipelining, multi-level caching, etc. This part does not support mixing float and double. You can bypass this heavily optimized implementation with m1.lazyProduct(m2) that corresponds to the implementations used fro small fixed-size matrices, but there is only disadvantages of doing so: the ALUs does not support mixing float and double, so float values have to be promoted to double anyway and you will loose vectorization. Better cast the float to double explicitly:

m1.cast<double>() * m2
ggael
  • 28,425
  • 2
  • 65
  • 71
  • 1
    Actually, my real problem is to define custom scalar types `T1` and `T2`. I've found that using `ScalarBinaryOpTraits` works mostly but not for the multiplication of two dynamic-sized Eigen matrices of `T1` and `T2`. I used `double` and `float` to make the example simpler. – Soonho Kong Jan 05 '17 at 22:09
  • I'm wondering if one of the followings looks a good solution to you: 1. Provide operator overloading `operator*` for Eigen matrices of `T1` and `T2` explicitly (probably using cast operation in it) 2. Provide a specialization of `Eigen::internal::general_matrix_matrix_product` for `T1` and `T2`. If this is the way to go, is there a good example of doing this? – Soonho Kong Jan 05 '17 at 22:12
  • 3
    We should really bypass this heavy code for custom scalar types. This will be the case for 3.4 and probably backported to some 3.3.x release. Meanwhile, you can indeed specialize `general_matrix_matrix_product` as [there](https://bitbucket.org/eigen/eigen/src/8dc0c20fada6963cdc42e97d4468b067b24b1be2/Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h?at=default&fileviewer=file-view-default), or also follow this [piece of code](https://bitbucket.org/eigen/eigen/src/8dc0c20fada6963cdc42e97d4468b067b24b1be2/unsupported/Eigen/MPRealSupport?at=default&fileviewer=file-view-default#MPRealSupport-140). – ggael Jan 06 '17 at 09:27
  • @ggael the developers of Stan (and I) [are discussing this matrix multiplication issue](https://discourse.mc-stan.org/t/help-testing-release-candidate-for-std-complex/12314/28) on their forum in the context of their implemention of `std::complex` math on their automatic differentiation types. Would appreciate any thoughts you might have (either there or here). – Chris Chiasson Dec 17 '19 at 20:44
  • @ggael not sure if I properly tagged you in the last comment, so I'm just making sure I tab completed it correctly. Sorry if it double notifies you. – Chris Chiasson Dec 17 '19 at 20:45