I have an application where I'd like to overload the basic operations, for example multiplication, of Eigen matrices. The difficulty is that I need to do a bunch of computations in addition to the actual computation of the matrix product, so the NumTraits approach discussed in the doc isn't applicable.

Now I can define a "multiply" function that works great,

Code: Select all

 
template<typename DerivedA, typename DerivedB>
  // enable_if magic to ensure that DerivedA and DerivedB are of the necessary type
  multiply(const Eigen::MatrixBase<DerivedA>& M1,
           const Eigen::MatrixBase<DerivedB>& M2) {
   
    const eigen_idx_t N = M1.cols();
    Eigen::Matrix<typename Eigen::MatrixBase<DerivedA>::Scalar, Eigen::Dynamic, Eigen::Dynamic> M(N, N);
    // Fill M
    return M;
   
  }

But if I try to overload operator* directly,

Code: Select all

  template<typename DerivedA, typename DerivedB>
  // enable_if magic to ensure that DerivedA and DerivedB are of the necessary type
  operator*(const Eigen::MatrixBase<DerivedA>& M1,
           const Eigen::MatrixBase<DerivedB>& M2) {
   
    const eigen_idx_t N = M1.cols();
    Eigen::Matrix<typename Eigen::MatrixBase<DerivedA>::Scalar, Eigen::Dynamic, Eigen::Dynamic> M(N, N);
    // Fill M
    return M;
   
  }


the compiler ignores it and appeals to the usual Eigen operator*.

Is it possible to do something like this with some additional template magic? Thanks!

-Mike

In the evaluator-branch (https://bitbucket.org/ggael/eigen-evaluators) which I'm going to merge very soon, you can specialize the evaluator of operator* with respect to many aspects including the scalar types of each operand.

Meanwhile, you might specialize our custom GEMM routine as we do to fallback to BLAS, see: https://bitbucket.org/eigen/eigen/src/a ... at=default

general_matrix_matrix_product::run() is automatically called by expressions looking to res += alpha * lhs * rhs, which are large enough. In particular, this does not concern inner, outer, or matrix-vector products.


Depending on what you really wanna do, there might be other options too, like creating a new field.

So for large enough dimensions MatrixBase * MatrixBase eventually resolves to calling general_matrix_matrix_product::run(),
and implementing a custom general_matrix_matrix_product class is the easiest means of implementing a custom operator*
until the evaluators are in?

Any guidance on how the templates actually resolve to general_matrix_matrix_product::run()? For example, why wouldn't
a custom

Code: Select all

template<typename DerivedA, typename DerivedB>
  // enable_if magic to ensure that DerivedA and DerivedB are of the necessary type
operator*(const Eigen::MatrixBase<DerivedA>& M1,
           const Eigen::MatrixBase<DerivedB>& M2)


be picked up over the default implementation or at least cause a resolution ambiguity? Is there a definition like

Code: Select all

template<typename Scalar>
Eigen::Matrix<typename Scalar, Eigen::Dynamic, Eigen::Dynamic>
operator*(Eigen::Matrix<typename Scalar, Eigen::Dynamic, Eigen::Dynamic>& lhs,
         Eigen::Matrix<typename Scalar, Eigen::Dynamic, Eigen::Dynamic>& rhs)


that ends up taking precedence in the resolution?

Thanks!

operator* is a member of MatrixBase<>, see: https://bitbucket.org/eigen/eigen/src/a ... ult#cl-189