# When should I use C++ expression templates in computational science, and when should I *not* use them?

Suppose that I'm working on a scientific code in C++. In a recent discussion with a colleague, it was argued that expression templates could be a really bad thing, potentially making software compilable only on certain versions of gcc. Supposedly, this problem has affected a few scientific codes, as alluded to in the subtitles of this parody of Downfall. (These are the only examples I know of, hence the link.)

However, other people have argued that expression templates are useful because they can yield performance gains, as in this paper in SIAM Journal of Scientific Computing, by avoiding storage of intermediate results in temporary variables.

I don't know a whole lot about template metaprogramming in C++, but I do know that it is one approach used in automatic differentiation and in interval arithmetic, which is how I got into a discussion about expression templates. Given both the potential advantages in performance and the potential disadvantages in maintenance (if that's even the right word), when should I use C++ expression templates in computational science, and when should I avoid them?

• Ah, the video is too funny. I didn't know it existed. Who made it, do you know? Mar 25 '12 at 17:34
• No idea; a couple of the PETSc people sent me links at one point. I think a FEniCS developer made it. Mar 25 '12 at 18:13
• The video link is broken and I'm dying from curiosity. New link? Oct 24 '15 at 6:54
• Oh drat, nevermind, I see youtube has come for our Hitler videos. Oct 24 '15 at 7:06

My problem with expression templates is that they are a very leaky abstraction. You spend a lot of work writing very complicated code to do a simple task with nicer syntax. But if you want to change the algorithm, you have to mess with the dirty code and if you slip up with types or syntax, you get completely unintelligible error messages. If your application maps perfectly to a library based on expression templates, then it might be worth considering, but if you aren't sure, I would recommend just writing normal code. Sure, the high level code is less pretty, but you can just do what needs to be done. As a benefit, compilation time and binary sizes will go way down and you won't have to cope with huge variance in performance due to compiler and compilation flag choice.

• Yeah, I've seen some of the lengthy error messages firsthand when I had to port code from gcc 2.95 to gcc 4.x, and the compiler was throwing all sorts of errors about templates. A lab-mate of mine is developing a templated library for interval arithmetic in C++ (adding new features that aren't in Boost::Interval in order to accomplish more research), and I don't want to see the code become a nightmare to compile. Mar 15 '12 at 11:48

Others have commented on the issue of how difficult it is to write ET programs as well as the complexity of understanding error messages. Let me comment on the issue of compilers: It is true that a while back one of the big issues was finding a compiler that's compliant enough with the C++ standard to make everything work and make it work portably. As a consequence, we found lots of bugs -- I have 2-300 bug reports in my name, distributed over gcc, Intel icc, IBM xlC and Portland's pgicc. Consequently, the deal.II configuration script is a repository of a large number of compiler bug tests, primarily in the area of templates, friend declarations, namespaces, etc.

But, it turns out that the compiler makers have really gotten their act together: today, gcc and icc today pass all our tests and it's easy to write code that's portable between the two of them. I would say PGI isn't far behind but it has a number of quirks that don't seem to go away over the years. xlC, on the other hand, is a whole different story -- they fix a bug every 6 months, but despite filing bug reports with them for years, progress is extremely slow and xlC has never been able to compile deal.II successfully.

What this all means is this: if you stick with the two big compilers, you can expect that they just work today. Since most computers and OSs today typically have at least one of them, that's enough. The only platform where things are more difficult is the BlueGene, where the system compiler typically is xlC, with all its bugs.

• Just out of curiosity, have you tried compiling against the new xlc compilers on /Q? Sep 29 '12 at 13:39
• No. I will admit that I've given up on xlC. Sep 29 '12 at 16:05

I have experimented a little bit with ET's a long time ago when, as you mentioned, compilers were still struggling with them. I've used the blitz library for linear algebra in some code of mine. The problem then was getting the good compiler and as I'm not a perfect C++ programmer, interpreting the compiler error messages. The latter was simply unmanageable. The compiler would, on average, generate about 1000 lines of error messages. No way I was able to quickly find my programming error.

You can find more information on the oonumerics webpage (there are the proceedings of two ET workshops).

But I would stay far away from them....

• The compiler error messages are indeed one of my concerns. With some of the templated C++ code I compile in order to build libraries for my projects, the compiler may generate hundreds of lines of warning messages. However, it's not my code, I don't understand it, and generally speaking, it works, so I leave it alone. Lengthy, cryptic error messages don't bode well for debugging. Mar 15 '12 at 11:43

The problem already starts with the term 'expression templates (ET)'. I don't know if there is a precise definition for it. But in its common usage it somehow couples 'how you code linear algebra expressions' and 'how it gets computed'. For example:

You code the vector operation

v = 2*x + 3*y + 4*z;                    // (1)


And it gets computed by a loop

for (int i=0; i<n; ++i)                 // (2)
v(i) = 2*x(i) + 3*y(i) + 4*z(i);


In my opinion this are two different things and need to be decoupled: (1) is a interface and (2) one possible implementation. I mean this is common practice in programming. Sure (2) might be a good default implementation, but in general I want to be able to utilize a specialized, dedicated implementation. For instance, I want that a function like

myGreatVecSum(alpha, x, beta, y, gamma, z, result);    // (3)


get called when I am coding (1). Maybe (3) just uses internally a loop like in (2). But depending on the vector size other implementations might be more efficient. Anyway, some expert in high performance can implement and tune (3) as much as possible. So if (1) can not be mapped to a call of (3) then I rather avoid the syntactic sugar of (1) and directly call (3) right away.

What I describe is nothing new. On the contrary, it's the idea behind BLAS/LPACK:

• All performance critical operations in LAPACK are done by calling BLAS functions.
• BLAS just defines a interface for those linear algebra expressions that are commonly needed.
• For BLAS different optimized implementations exist.

If the scope of BLAS is not sufficient (e.g. it does not provide a function like (3)) then one can extend the scope of BLAS. So this dinosaur from the 60s and 70s realizes with its stone age tool a clean and orthogonal separation of interface and implementation. It's kind of funny that (most) numerical C++ libraries do not achieve this level of software quality. Although the programming language itself is so much more sophisticated. So it is no surprise that BLAS/LAPACK is still alive and actively developed.

So in my opinion ETs are not evil per se. But how they are commonly used in numerical C++ libraries earned them a very bad reputation in the scientific computing circles.

• Michael, I think you are missing one of the points of expression templates. Your code example (1) does not actually map to any optimized BLAS calls. In fact, even when a BLAS routine exists, the overhead of a BLAS function call makes it fairly terrible for small vectors and matrices. Sophisticated expression template libraries like Blaze and Eigen can use deferred expression evaluation to avoid the use of temporaries, but I'm convinced that almost nothing short of a domain specific language is going to be able to beat hand-rolled linear algebra. Sep 29 '12 at 13:31
• No, I think you are missing the point. You have to distinguish between (a) BLAS as a specification of some frequently needed linear algebra operation (b) an implementation of BLAS like ATLAS, GotoBLAS, etc. BTW that how it works in FLENS: By default an expression like (1) would be evaluated by calling axpy from BLAS three times. But without modifying (1) I also could evaluate it like in (2). So what happens logically is the following: if an operation like in (1) is important then the set of specified BLAS operations (a) can be extended. Sep 29 '12 at 13:46
• So the key point is: Notation like 'v=x+y+z' and how it finally gets finally computed should be separated. Eigen, MTL, BLITZ, blaze-lib completely fail in this respect. Sep 29 '12 at 13:54
• Right, but the number of frequently needed linear algebra operations is combinatoric. If you are going to use a language like C++, you have the choice of either implementing-as-needed using expression templates (this is the Eigen/Blaze approach) by combining sub-blocks and algorithms intelligently using deferred evaluation, or implementing a massive library of every possible routine. I don't advocate either approach, as recent work in Numba and Cython shows we can get similar or better performance working from high-level scripting languages like Python. Sep 29 '12 at 14:31
• But again, what I complain about is the fact that such sophisticated (in the sense of complicated but inflexible) libraries like Eigen tightly couple notation and evaluation-mechanism and even think it is a good thing. If I use a tool like Matlab I just want to code things and rely that Matlab is doing the best possible thing. If I use a language like C++ then I want to be in control. So appreciate if a default evaluation-mechanism exists but it must be possible to change it. Otherwise I go back and call functions in C++ directly. Sep 29 '12 at 23:49