Search code examples
c++templatesmatrixeigen3

Template overloaded operator in Eigen::Matrix


I got this basic code working:

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

using T12 = Eigen::Matrix<double, 2, 1>;
using T22 = Eigen::Matrix<double, 2, 2>;

template <typename Der1, typename Der2> T22
operator^ (const Eigen::MatrixBase<Der1>& d1, const Eigen::MatrixBase<Der2>& d2) {
  return d1 * d2.transpose();
}

void testT12DyadT12 () {
  T12 t12_1, t12_2;
  t12_1 << 2.0, 3.0;  t12_2 << 1.5, 2.5;

  const double c1{2.0}, c2{0.5};
  T22 t22_1 = c1*t12_1 ^ t12_2*c2;
  std::cout << t22_1 << "\n"; 
}

int main () {
  testT12DyadT12();
  //testT1nDyadT1n ();
  return 0;
}

This works and produces the correct 2 x 2 matrix output. Now I want to templatize the number of rows and columns in the matrix as

template<int nc>
using T1n = Eigen::Matrix<double, nc, 1>;
template<int nc>
using T2n = Eigen::Matrix<double, nc, nc>;

template <int nc, typename Der1, typename Der2> T2n<nc>
operator& (const Eigen::MatrixBase<Der1>& d1, const Eigen::MatrixBase<Der2>& d2) {
  return d1 * d2.transpose();
}

void testT1nDyadT1n () {
  T1n<2> t12_1, t12_2;
  t12_1 << 2.0, 3.0;  t12_2 << 1.5, 2.5;

  T2n<2> t22_1 = t12_1 & t12_2;
  std::cout << t22_1 << "\n"; 
}

But it gives this compiler error:

error: invalid operands to binary expression ('T1n<2>' (aka 'Eigen::Matrix<double, 2, 1, 0, 2, 1>') and 'T1n<2>')
T2n<2> t22_1 = t12_1 & t12_2;
               ~~~~~ ^ ~~~~~
note: candidate template ignored: couldn't infer template argument 'nc'
operator& (const Eigen::MatrixBase<Der1>& d1, const Eigen::MatrixBase<Der1>& d2) {
^

I understand the error message (that the compiler is not able to infer the template arg). How do I fix this? I tried passing the args to the operator as const Der1&, const Der2&, to no avail. I am sure it's something silly ... that I have not understood or paid attention to. Eventually my goal is to templatize Eigen scalar type as well: template< typename T, int nc > using T1n = Eigen::Matrix<T, nc, 1>; etc., and overload operator^.


Solution

  • You can derive the size of the return type from Der1::RowsAtCompileTime:

    template <typename Der1, typename Der2> 
    T2n<Der1::RowsAtCompileTime>
    operator& (const Eigen::MatrixBase<Der1>& d1, const Eigen::MatrixBase<Der2>& d2) {
      return d1 * d2.transpose();
    }
    

    Or fully generic (assuming Der1 and Der2 have the same scalar type):

    template <typename Der1, typename Der2> 
    Eigen::Matrix<typename Der1::Scalar, Der1::RowsAtCompileTime, Der2::RowsAtCompileTime>
    operator& (const Eigen::MatrixBase<Der1>& d1, const Eigen::MatrixBase<Der2>& d2) {
      return d1 * d2.transpose();
    }
    

    With C++14 you could of course just return auto (this will return a lazy expression, which may or may not be what you expect):

    template <typename Der1, typename Der2> 
    auto
    operator& (const Eigen::MatrixBase<Der1>& d1, const Eigen::MatrixBase<Der2>& d2) {
      return d1 * d2.transpose();
    }