Search code examples
c++vectorizationeigen

Loop vectorising of vector3d products using Eigen

Hi I'm using Eigen to roll or vectorise a number of loop operations in a particle filter.

In essence I have two matricies

Eigen::Matrix<N, 3> A;
Eigen::Matrix<3, N> B;

Where N is a large number.

And I would like a one line which does the equivalent of:

Eigen::Matrix<N, 1> D;
for (size_t i=0; i< N; i++)
{
   D.row(i) = A.row(i)*B.col(i);
}

I had been trying to use D =A.rowwise()*B.colwise() but these broadcasting methods do not define an operator*() between them.

Solution

  • Here is a version that is pretty much optimal for small vectors like these.

    Eigen::MatrixX3d A;
Eigen::Matrix3Xd B;
Eigen::VectorXd D = (A.array() * B.array().transpose()).rowwise().sum();

    Fine-tuning this expression for larger vector size, e.g. if these were square matrices, is a bit of a challenge but for 3 rows/columns it works fine.

    If you can chose the shape of your matrices, consider changing to this:

    Eigen::Matrix3Xd A, B;
Eigen::VectorXd D = (A.array() * B.array()).colwise().sum();

    or this:

    Eigen::MatrixX3d A, B;
Eigen::VectorXd D = (A.array() * B.array()).rowwise().sum();

    Both produce assembly that looks very good and makes full use of vector instructions. Tested with Eigen-3.4.0, GCC-11.3 -O3 -DNDEBUG. The second one is better with low dimensions (MatrixX2d up to MatrixX4d). The first will scale better to large sizes.