Introduction into Decomposition Methods of Linear Equation Systems in C++

The math parts are good. C++ parts are less so.

Arrays can be on the stack just as vectors. The main cause of performance difference between them is cache misses, stack is almost always well cached, newly allocated heap memory almost never.

It's rarely a good idea to write matrix libraries these days. Eigen switched to developer-friendly MPL license in 2012, I've been using it since then and it's awesome. It has substantially more internal complexity, but when used correctly it's faster that naïve solutions by an order of magnitude. Template metaprogramming pays off for their use case.

Quite often you know sizes of your matrices at compile time. If that size is 3x3 or 27x27 or something else reasonably small, switching to compile time sized matrices can improve performance by a lot. These fixed-size matrices and vectors are essentially value types, they have no pointers to chase and they don't use dynamic memory.

When you have complex formulae operating on vectors and matrices, if you use auto for intermediate variables, Eigen code often compiles into a single loop over the source data. These intermediate things of human-unfriendly types are essentially compile-time futures.