Outperforming LAPACK with C++ metaprogramming

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.

As an anecdote of just how powerful this can be, I once wrote a prototype of a simulation (written in C) that used a lot of integer divisions and needed to be run many times with different parameters.

As a quick and dirty hack that took about half hour to implement, I decided to just have the program as a template string and the parameters be substituted in by a python script that called GCC. To the amazement of a very sceptical friend, the result was a 5x speedup.

A lot of this speedup was consistent with replacing the integer divisions with multiplicative inverses, but also merging constants in a way that would have made things quite unreadable. I was quite impressed with how far GCC had gone.

Considering each instance of this simulation ran for tens of minutes, it was well worth compilation time of a few seconds. It kept the code really simple, too.

I think that's why it uses `volatile double` in stead of `double`?

> With maximal inlining, there's a risk of constant propagation resulting in a lot of computation happening at compile time rather than run time.

Without even reading TFA, this was my first thought, remembering back to implementing Fibonnaci's sequence in template metaprogramming from "Thinking in C++, Volume 2": http://web.mit.edu/merolish/ticpp/TicV2.html#_Toc53985733

It's cool and useful to be sure, but as costrouc pointed out, LAPACK has been optimized for a very long time by some very smart people. Most people who have not studied the field for a bit are unlikely to come across a more optimal method.

Poking at the various compilers the machine instructions generated from the article appear to produce an unrolled version of the operation (as claimed by the article): https://godbolt.org/g/Wmj8ar

Also, it seems the author solved a triangular system. LAPACK has special routines for that. Were they used?

LAPACK is typically optimized for larger systems, not 5 unknowns. Also, 5 is not a great number for vectorized operations - it might even be beneficial to zero/one pad the matrix.

Optimized LAPACK is often 5-10x faster than «basic LAPACK».

BLAS/LAPACK was originally written in the 70s/80s and sometimes unroll loops to the tune of 5 or 7 - it made sense then. Not so much these days.

I didn’t read the article in detail, but there appear to be a lot of holes in it.

Thanks for sharing the godbolt site, can see myself using it for debugging.

What bothers me more is the fact that the article seems to compare LAPACKs function for solving any-sized systems of equations with a specialized function for solving 5x5 systems.

I mean, it's neat trick and all, but the comparison with LAPACK does not really stand on equal ground.

> instead of the bug elevated to feature that metaprogramming is in C++.

That inflammatory wording doesn’t add anything to your comment. If you wanted to compare D with C++ it’d be a lot more useful to link to something which does that rather than just the documentation index.

>and metaprogramming in D is much nicer, because it was designed from the ground up instead of the bug elevated to feature that metaprogramming is in C++

Yet neither of them are as elegant as LISP macros.

Or even Dlang .

These kinds of tricks certainly are more pleasant to do in languages with real metaprogramming support.

For numerical applications, one interesting project that comes to mind is the Terra programming language. It uses a flexible high level language (Lua) to metaprogram a fast and predictable C-like language (Terra). The advantage compared to the examples you mentioned is that each language can focus on its strength: the low level language being similar to C allows the programmer to better reason about performance and having a high level language for metaprogramming keep things simple.

https://terralang.org

There are also systems out there specifically designed for numeric applications that let you specify an algorithm in a high level DSL and generate highly performing C code, some of which might even tesr various versions of the code to find one that runs better on your machine (for example, the best way to iterate the data might be affected by cache sizes). I don't remember the names of these off the top of my head but the Terra PhD thesis has lots of references to them.

It was playing around with exactly these sorts of tricks in c++ (with late 90's compilers) that drove me to learn common lisp. It's a much nicer environment to do this sort of thing in, but difficult to compete head to head with a good optimizing c++ compiler on code generation in general. Fwiw, in common lisp particularly having access to to compiler macros allows for interesting approaches.

I was doing primarily research code at the time though, and the overall improvement in productivity more than made up for the fractional difference in run times for me. Of course, ymmv.

It does lead to performance gains. This kind of stuff is done all the time in Julia's numerical codes. Flexible macros definitely make this easier, and multiple dispatch + automatic function specialization on types + inlining + interprocedural optimizations makes a lot of what's discussed here automatic. If it's not automatic, then libraries like StaticArrays.jl make functions that use small arrays do these kinds of tricks automatically (it overloads and unrolls operations specifically to match the size of the array which is compile time information). I am sure other languages can do this too, but having it semi-automatic is pretty cool if what you do all day is write numerical codes.

Just be wary of performance gains claims. Per other comments this implementation is highly optimized for a specific use and uses a numerically unstable algorithm.

It's not nearly that simple. You can use template meta-programming to make expressions both easy to read and efficient, by avoiding creation of temporary objects, fusion, lifting, etc. etc.

In some domains (e.g. linear algebra) this can be very powerful because you don't have write non idiomatic expressions just to feed the compiler things it can understand.

There needs to be a j in the exit condition of the inner loop of the clobber function.

I've read your comment several times and can't follow it. I think you're saying that even if the author continued to use Cramer's rule, he would lose to LAPACK on big problems.