[stephencanon]

There's no mystery here, author is only solving a 5x5 system. LAPACK will have barely finished checking function call parameters by the time a fully specialized and inlined solver is finished. Also, from a cursory inspection, author does no pivoting [edit: worse, it just uses Cramer's rule, which is insane], so their solver will be much less numerically robust (but this is admittedly less of a problem for the 5x5 case).

The performance is completely expected--in fact, I would expect an optimized 5x5 solver to be faster than this.

[electricslpnsld]

A better comparison might be a C++ library like Eigen that is intended for use on these small systems, that should optimize out function calls in a similar way, and that uses intrinsics.

[stochastic_monk]

If your matrices are of arbitrary size, Blaze [0] would be my pick. It has the best performance on their benchmarks, Baptiste Wicht's benchmarks, and my own.

Though if the matrices are guaranteed to be small (like in this example), libxsmm [1] is specialized and highly optimized for this use case and beats its competitors above.

And yes, it's absolutely essential to make sure you didn't just move computations to compile-time.

[0]: https://bitbucket.org/blaze-lib/blaze

[1]: https://github.com/hfp/libxsmm

[stabbles]

Cramer's rule is unstable even for 2x2 problems

[mattnewton]

I think gp is missing a “less” before “numerically robust”

[stephencanon]

Yeah, thanks. Fixed.