fix transposed matrix product bug

Reference issue

Fixes #2266 (closed)

What does this implement/fix?

The transpose / adjoint of matrix products is being evaluated correctly, but incurs unnecessary allocations even if noalias is called. For example, the below snippet will allocate a useless temporary matrix:

#define EIGEN_RUNTIME_NO_MALLOC
#include <Eigen/Core>
int main()
{
    int n = 10;
    MatrixXd A = MatrixXd::Random(n, n);
    MatrixXd B = MatrixXd::Random(n, n);
    MatrixXd C = MatrixXd::Zero(n, n);

    Eigen::internal::set_is_malloc_allowed(false);
    C.noalias() = (A * B).transpose(); // triggers malloc despite noalias()
    Eigen::internal::set_is_malloc_allowed(true);
}   

This is because evaluator<Transpose<Product<Lhs,Rhs>>> will first create the evaluator for Product<Lhs,Rhs>, which by default creates a temporary. An easy fix is to write out the expanded matrix product C.noalias() = B.transpose() * A.transpose(); which is handled correctly. This patch specializes transpose() and adjoint() to return an expanded expression of the transposed product. This is arguably syntactic sugar but the issue is so subtle it could be causing unnecessary allocations in many codebases.

Additional information

Edited by Charles Schlosser

Merge request reports

Loading