Search code examples
c++eigeneigen3

Write a function with transpose as argument with the Eigen 3 library


I want to write a function that can accept:

  1. a dense array/matrix, or
  2. the transpose of a dense array/matrix.

Is it possible to avoid using perfect forwarding?

I have tried to use DenseBase type argument. But that can't accept the transpose of a matrix.

I don't like to use perfect forwarding because implementing type checking with sfinae would be tedious.

Current solution:

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

using namespace Eigen;

template <typename U>
auto f(U&& x) {
    auto x2 = std::forward<U>(x);
    auto max_x = x2.colwise().maxCoeff().eval();
    x2 = x2.rowwise() + max_x;
    return max_x;
}

int main() {
    Array<float, 3, 3> M1;
    M1 << 1, 2, 3, 
          4, 5, 6, 
          7, 8, 9;
    std::cout << M1 << "\n";
    // auto here might cause problem later ...
    // see eigen.tuxfamily.org/dox/TopicPitfalls.html
    auto max_x = f(M1.transpose());
    std::cout << M1 << "\n";
    std::cout << max_x << "\n";
}

Result:

// original
1 2 3
4 5 6
7 8 9
// Increase each row by max of the row.
 4  5  6
10 11 12
16 17 18
// Max of each row (not a column vector).
3 6 9

I tried EigenBase with the following lines:

template <typename U>
auto f(EigenBase<U>& x) {
...

Compiler error:

test4.cpp:20:32: error: cannot bind non-const lvalue reference of type ‘Eigen::EigenBase<Eigen::Transpose<Eigen::Array<float, 3, 3> > >&’ to an rvalue of type ‘Eigen::EigenBase<Eigen::Transpose<Eigen::Array<float, 3, 3> > >’
     auto max_x = f(M1.transpose());
                    ~~~~~~~~~~~~^~

Solution

  • Use auto type to hold a transpose expression before calling the function.

    #include <Eigen/Eigen>
    #include <iostream>
    
    using namespace Eigen;
    
    template <typename Derived>
    auto f(DenseBase<Derived>& x) {
        auto max_x = x.colwise().maxCoeff().eval();
        x = x.rowwise() + max_x;
        return max_x;
    }
    
    int main() {
        Array<float, 3, 3> M1, M2;
        M1 << 1, 2, 3, 
              4, 5, 6, 
              7, 8, 9;
        M2 = M1;
        std::cout << M1 << "\n";
    
        std::cout << "no transpose\n";
        Array<float, 3, 1> max_x = f(M1);
        std::cout << M1 << "\n";
        std::cout << max_x << "\n";
    
        std::cout << "transpose\n";
        auto m2_t = M2.transpose();
        Array<float, 1, 3> max_x2 = f(m2_t);
        std::cout << M2 << "\n";
        std::cout << max_x2 << "\n";
    }
    

    Result:

    1 2 3
    4 5 6
    7 8 9
    no transpose
     8 10 12
    11 13 15
    14 16 18
    7
    8
    9
    transpose
     4  5  6
    10 11 12
    16 17 18
    3 6 9