Minimize temporaries when adding std::arrays

[u/nomadwolf0]

I have a bunch of code of the form (sometimes in more convoluted fashion)

using aVec = std::array<double, 32>;
aVec a, b, c, d, e;      // some are constexpr, others runtime values
double x, y;
for ( int i = 0; i < 32; ++i)
  a[i] = x * b[i] + y * c[i] + y*y*d[i] + e[i];

I'd like to rewrite it such that a = x * b + y * c + y*y*d + e; to generally be easier to read intent, but I don't want all of those operations to create & destroy a bunch of temporary std::arrays.

Is there any straightforward way to achieve this? These generally lives in the inner (or mid-level) loops.

Only thing I could think of is to have the addition & scalar multiplications operators return a proxy type that is essentially a fixed-length std::vector that implicitly converts to std::array. It'll be a bit slower than the code I'm trying to replace, but the move semantics should reduce that impact.

[u/the_poope]

There's no straightforward way, but there is a way and it's called expression templates.

Instead of rolling your own you could use one of the existing libraries that use this technique, such as Eigen, blaze, xtensor or armadillo. But maybe your usecase is so small that it doesn't warrant including large libraries (they are mostly header-only being template libraries, but they do add quite a bit to compile times).

[u/KazDragon]

Have you used e.g. Compiler Explorer to check that your optimizer isn't doing it for you anyway?

[u/aiusepsi]

Personally, I wouldn't worry about temporaries here, using std::array. All allocating an array of 32 doubles does is cause the stack pointer register to be decremented by 256 bytes. It is, effectively, free. Returning a std::array from a function doesn't cause any overhead either; because of copy elision, the return result is written directly into the caller's stack, no copy required.

Even better, in an optimised build, the compiler can just optimise away use of the stack if it's not necessary because your data will fit into registers. I implemented an wrapper around std::array, implemented the * and + operators on it, and the expression a = x * b + y * c + y*y*d + e optimises into a load of the data from b, c, etc. into vector registers, does the maths using the registers, and then spits out the result into a store to a. See it here: https://godbolt.org/z/EMvchf94G

[u/nomadwolf0 → u/aiusepsi]

Not sure why the wrapper is necessary, but point taken.

As noted by u/fortsnek274, MSVC is less than great ([https://godbolt.org/z/M33W6saWd\](https://godbolt.org/z/M33W6saWd)). Not even sure why it's doing most of the operations on QWORDs instead of XMMWORDs.

When I reduced the array size to 4, it did inline it all within the function: ([https://godbolt.org/z/q9nnoeYcM\](https://godbolt.org/z/q9nnoeYcM)).

[u/StaticCoder → u/aiusepsi]

Allocating an array on the stack requires more than just increasing the stack pointer. It also requires accessing that memory, which is going to affect CPU caches.

[u/aiusepsi → u/StaticCoder]

A std::array for non-class types (like double) starts uninitialised, which can be surprising. Because of that, just allocating one doesn’t touch memory at all.

[u/StaticCoder → u/aiusepsi]

Sure, but we're not talking about unused temporaries here. They would be used. gcc notably also tends to be quite bad about stack reuse.

[u/n1ghtyunso]

so essentially you want a different type than std::array that lets you write vector-based code instead.
And when you do this with regular operator overloading, you'd need to create intermediary results for each operator return value.
And the common way around this is expression templates.

Ideally you could use a math library to do this for you.

[u/fortsnek274]

Nice exercise in seeing how bad MSVC is: https://godbolt.org/z/abxxr6sWv

[u/nomadwolf0 → u/fortsnek274]

Fantastic, because I'm on MSVC.

[u/fortsnek274 → u/nomadwolf0]

It's not all lost. You can still use clang-cl for the performance build!

[u/Independent_Art_6676]

you can certainly do it without a temporary:
a = e;
a += yyd:
a += yc;
a += xb;

this stinks because you now did 4 loops instead of 1. But it didn't use a temporary! However if the inputs are DISPOSABLE you can do the above into pieces in parallel, eg c*=y and d*= ysquared and e+= xb all threaded then put it together? It would take a rather large size to justify threading like this and it still cost you the time to iterate twice (once for the partials and once for the final assembly). Also you could look at some of the inputs... perhaps where C is computed you can go ahead and add in a multiply by y so its ready to use? Sometimes you can get away with that sort of thing even if its a little convoluted. Nothing here really strikes me as incredibly useful ways to approach it other than if you can do the scalar multiplications ahead of time.

move is slick, but as often as not you can rig the code so that isn't necessary. That is, if you want to move a into b, sure, ok, but... just keep using a and not have b is slightly better and with this kind of math code, you can often write it that way instead. if you can't avoid the move, its cheap enough compared to everything else.

For linear algebra I often have a pool of global, always existing temporaries that I reuse for intermediates. Its frequently impossible to avoid them, but you can at least avoid create/destroy repeated costs. Normally about 5 is sufficient for me. That may be the only practical thing I can offer tonight. Once you realize that avoiding the temps is more or less impossible any time you have terms like this... minimize the damage is all you can do (or go with the loop you started with, which is exactly how a library or expression template would do it). What you started with seems close to ideal.

** When I say you can't avoid the temporary, I mean without rolling it all out into an equation specific loop like what you started with. If you wrote a little matrix class that could return scalar * vector or matrix + matrix etc, it would crank out temporaries like mad. Obviously they can be avoided if you have a way to boil it all down to a single loop. A pool of temps is only useful if you have something that can't be rolled out into one loop, and its been a while since I had one of those. Not even sure it can crop up with just vectors, but with full matrix, it can.

[u/mc_pm]

In Python I'd use numpy for that, it permits just that sort of construction. It looks like maybe Eigen for C++?

[u/Total-Box-5169]

You can do it with alias declaration: one intermediary that alias the array element and another that alias the intermediary and shadows the array variable. To avoid having to write lots of alias declarations you need a MACRO that takes the expression and the array variables in the expression. To make it work with up to N arrays variables you need several macros that allows to do that trick.

FOREACH(a = x * b + y * c + y * y * d + e,
    a, b, c, d, e);

https://godbolt.org/z/5dGa5Eq4x