Eigen Matrix-Vector multiplication w/ noalias

user_1745861 0
Viewed 544

I've been going through the Eigen documentation for using noalias to avoid unnecessary temporary allocations when doing Matrix-Matrix products, but I was wondering if it's possible to use noalias in the following situation:

Eigen::VectorXf x(3);
x << 1, 2, 3;

Eigen::MatrixXf A(3, 3);
A << 1, 2, 3, 4, 5, 6, 7, 8, 9;

// Is this valid? Is it valid under certain size assumptions for A and x but not others?
x.noalias() = A * x;

Naively, it seems like noalias might be valid in this case because you really only need to access the Vector elements once per column in the Matrix.

On the other hand, x clearly appears on both sides of the expression, and matrix multiplication involves all sorts of low level black magic that makes situations like this hard to reason about.

eigen c++
1 Answers

After testing with Clang trunk and Eigen trunk, with optimizations turned on, the case with noalias produces incorrect results. GCC trunk has the same behaviour.

The outputs:

Program returned: 0
.noalias
0
0
0
alias
14
32
50

Code:

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

auto no_alias() -> Eigen::VectorXf
{
    Eigen::VectorXf x(3);
    x << 1, 2, 3;

    Eigen::MatrixXf A(3, 3);
    A << 1, 2, 3, 4, 5, 6, 7, 8, 9;

    x.noalias() = A * x;
    return x;
}

auto alias() -> Eigen::VectorXf
{
    Eigen::VectorXf x(3);
    x << 1, 2, 3;

    Eigen::MatrixXf A(3, 3);
    A << 1, 2, 3, 4, 5, 6, 7, 8, 9;

    x = A * x;
    return x;
}

int main () 
{
    std::cout << ".noalias" << std::endl;
    for (const auto& no : no_alias())
    {
        std::cout << no << std::endl;
    }

    std::cout << "alias" << std::endl;
    for (const auto& no : alias())
    {
        std::cout << no << std::endl;
    }
}
user_34769540
Related