[eigen] Avoiding memory allocations with triangular solves

[Wenzel Jakob <wenzel@xxxxxxxxxx>]

Hi,

I've run into performance issues when using Eigen to solve linear
systems with PartialPivLU.

Despite having pre-allocated all memory, the PartialPivLU::solve()
function appears to perform dynamic memory allocation, which creates a
bottleneck in my application. The code is structured as follows:

/// Initialization
Eigen::PartialPivLU<Eigen::MatrixXf>  lu(n);
lu.compute(...);

while (...) { // main computation done here
     ...
     output.block(i, j, n, 1) = lu.solve(input.block(k, l, n, 1));
     ...
}

When profiling this code, most time is spent in malloc() and free()
calls inside Eigen::internal::triangular_solve_matrix.

Hence, I would like to ask how it is possible to avoid dynamic memory
allocation when solving linear systems using Eigen's decompositions.

Thanks!
Wenzel