# Code to multiply matrices together

I've got a short algorithm used in my code to multiply two matrices and put them into a destination Matrix. I decided to use the temporary array t_elements, which allows the user to make one of the input matrices the output matrix as well.

Usually I like to save optimization for later, but this code is going to be called often, so I was wondering if anyone could suggest ways to make it faster.

#include <vector>

struct Matrix{
std::vector<std::vector<float>> elements;
int height;
int width;

Matrix(int w, int h)
{
width = w;
height = h;

for (int i = 0; i < height; i++) {
elements.push_back(std::vector<float>());

for (int j = 0; j < width; j++) {
elements[i].emplace_back(0.0f);
}
}
}

float get_element(int x, int y) const { return elements[y][x]; }

void set_element(int x, int y, float value) { elements[y][x] = value; }
};

void dot(Matrix const& left, Matrix const& right, Matrix& dest)
{
float t_elements[right.width][left.height];

for (int i = 0; i < left.height; i++) {
for (int j = 0; j < right.width; j++) {
float elementVal = 0;

for (int k = 0; k < left.width; k++)
elementVal += (left.get_element(k, i) * right.get_element(j, k));

t_elements[j][i] = elementVal;
}
}

for (int i = 0; i < dest.width; i++)
for (int j = 0; j < dest.height; j++)
dest.set_element(i, j, t_elements[i][j]);
}

int main()
{
Matrix m1 = Matrix(3,3);
Matrix m2 = Matrix(3,3);
Matrix m3 = Matrix(3,3);

dot(m1, m2, m3);

return 0;
}

• What sorts of Matrices are you going to be working with? The ideal approaches for matrices that are 3x3 is completely different from the best approaches when your matrices are 3000x3000, and different again when one is 3x3 and the other is 3x3000. – Josiah May 4 '18 at 22:52
• Do you really need to actually multiply them together? This often requires both memory allocation and lots of accesses. – mathreadler May 6 '18 at 8:13
• If you are using the Intel compiler, check to see that it isn't already detecting a matrix-multiply loop and optimizing in its own BLAS. If this is anything more than a toy project to learn how to optimize matrix multiplication for its own sake, you should be using a vendor-optimized implementation. MKL, ATLAS, cuBLAS, etc. – Alex Reinking May 6 '18 at 19:48
• In addition to Josiah's comment: Did you consider if the matrix is sparse (has a lot of 0 entries)? If so, you can save space and time by choosing different storage types (coordinated storage or compressed column are keywords) – ab.o2c May 8 '18 at 4:58

Both of the other answers make some very good points, especially when it comes to replacing vector. Two dimensional vectors are particularly iffy because the whole lump of memory will not be contiguous. If you access elements[y][x] you first have to look in elements[y] for the address of the inner vector, and then have to look into that inner vector. By contrast if you have an array that is say, 4x4, the compiler would typically internally have a single 1D array of length 16. Then to get elements[2][2] it might just grab the entry in offset 11 (i.e. 2*4+3)

There are (at least) three distinct performance advantages associated with this when you think about low level hardware.

1. You don't waste time chasing pointers to find which bit of memory you need.
2. Having the values that you're working with all be close together gives you nice cache locality.
3. This simplifies the code within the key loops, giving the compiler a better chance of automatically vectorizing your code.

Snowbody mentions templating your Matrix so that you could have different primitive types in it (A matrix of ints or of doubles or whatever). I would further suggest templating your matrix to specify the sizes. Comparing with how std::array is done, where you might have std::array<int, 3>;, you could have matrix<float, 3, 3>. (In fact it would be necessary if you wanted your internal data storage to be a std::array as vnp suggests.)

This again gives the compiler more guarantees about your program ahead of time, which again means it has a better chance of doing helpful optimisations. For example a classic optimisation that compilers might do is "Loop unrolling", whereby they actually rewrite

for (int i = 0; i < 3; ++i) {
doSomething();
}


as

doSomething();
doSomething();
doSomething();


They can do this because they know ahead of time how many iterations of the loop to expect, and calculate whether duplicating the machine code is justified by getting rid of the loop overhead.

As well as probably helping performance, this opens the door to writing compile time checks for your program. For example you could get the compiler to stop you from accidentally adding a 4x4 matrix to a 2x2 matrix, just like it would stop you from adding a double and a string.

Finally, because this information is known at compile time it does not need to be stored or accessed during the running of the program. Instead of i < dest.width the compiler would be writing machine code for i < 3. (But, importantly, you would still not have to write i < 3 which is inflexible and harder to understand than width.)

Modern optimisers have their own cost models to decide when to inline methods, but sensible hints still wouldn't go amiss. I'd suggest that get_element and set_element are particularly good candidates for inlining.

I actually disagree with vnp about directly accessing elements, although I do agree that the current getter and setter aren't any better. The advantages of using such methods are that they give you opportunity to

1. add precondition and postcondition checks, such as requiring that the indices are within bounds.

2. mess with the internal representation, such as the aforementioned storing your data in a 1D array and directly indexing it as if it was a 2D one. (If for some reason you didn't trust the compiler to do that indexing for you.)

Note that

float t_elements[right.width][left.height];


Is not legal c++, because c++ requires that array size (for stack arrays) must be a constant expression. You do have the option of using new, in which case you must remember to delete[]

I would also suggest, though this strays away from performance questions, that doing some error checking would be sensible. If you don't want to include sizes in the templates as above, it's still worth a basic check that you're not multiplying a 2x5 matrix by a 3x4 one.

Please note that all of the above is based on some theoretical understanding of the sorts of machine code that compilers can write. It is not based on empirical profiling of your code with a representative set of test data, which should be your first move when looking to optimise anything. (Closely followed of course by questioning whether there might be a suitable third party library that has already done all the hard work for you.)

• Good recommendations, except for adding inline hints. As I understand it, these are totally ignored by compilers. The keyword is only useful when defining functions in a header file to prevent multiple definitions. On balance this answer is a +1 though. – Cris Luengo May 5 '18 at 1:17
• I have heard such things, which is why the most I'd say is that the hint wouldn't hurt. However, this random guy on the internet decided that wasn't enough, so he went digging through the clang/llvm source. The summary is that clang may inline what you don't hint or may refrain when you do hint, but does seem to adjust its threshold based on the presence of the hint. blog.tartanllama.xyz/inline-hints – Josiah May 5 '18 at 17:57
• Of course those getters and setters are such obvious inline candidates, it is unlikely that the compiler would miss them anyway. I probably would be amenable to the line that because a modern compiler (presuming we can use a modern compiler) is more likely to get it right than the programmer, it's better not to meddle unless you need the multiple definitions semantics. – Josiah May 5 '18 at 18:03
• std::vector doesn't seem to be the best choice. The real strength of std::vector comes from its ability to dynamically change size. I don't think the Matrix class ever benefits from it. Consider std::array.

• Matrix::height and Matrix::width look redundant. This information is readily available as elements.size() and elements[0].size().

• I don't see any benefit in unrestricted getters and setters, especially with elements being public. Particularly the class method has no reason for not subscripting elements directly.

• Logically the dot method should be Matrix::operator*, and return the result, rather than have an in-out parameter. This may also avoid unnecessary copying.

If you insist of copying, prefer std::copy to a manual loop.

• Re performance, there is not much to be done with a naive implementation. It'd remain $O(n^3)$ no matter what. You may want to invest into Strassen algorithm.

• Using std::array would require compile time knowledge of the matrix size. Although there is a lot to be said for templating the size of the matrix (such as compile time checking that the sizes are compatible), that isn't what's been done here. – Josiah May 4 '18 at 22:55
• At least historically, giving matrix operations an out parameter was massively more performant than allocating new memory for results; I haven't done heavy matrix math since move constructors came so don't know now. Matrices also have other multiplication operations, so you wouldn't always know that * meant 'dot' rather than 'Hadamard' without context. – Pete Kirkham May 5 '18 at 8:57
• There is actually a ton of stuff to be done with a cubic time algorithm – harold May 5 '18 at 13:13
• The constants are only negligible when n is sufficiently large. And 3 is far from that. – Deduplicator May 5 '18 at 19:45
for (int j = 0; j < width; j++) {
elements[i].emplace_back(0.0f);


Instead of writing a loop, use the constructor form to create a vector of the desired length; all elements will be initialized to zero automatically.

for (int i = 0; i < dest.width; i++)
for (int j = 0; j < dest.height; j++)
dest.set_element(i, j, t_elements[i][j]);
}


Calling set_element on every element to copy the results — Copy the elements in the order in which they appear in the row, and use an iterator to place them efficiently.

for (int k = 0; k < left.width; k++)
elementVal += (left.get_element(k, i) * right.get_element(j, k));


Likewise for the get_element function, called repeatedly for every output value. For the inner dimension, just use iterator increment to advance to the next value.

Here is the way we did raster graphics in the old days: Don’t make a vector of vectors. Make a single vector, and an access function that multiplies the row by the row size and adds the column to produce a single index.

Now, given a pointer to any cell, you can efficiently move to the next in any direction. Going right, increment by one. Going down, increment by the row size. So trace through the source matrices this way, one going right, one going down.

If the compiler grasps this and auto-vectorizes the code, you are golden! So you might find out what coding idioms are understood by the compiler. Having the loop structured to traverse both inputs with constant (though different) strides is probably key, though. Getting the compiler to see you are doing a dot-product (on each row) is the most significant speed-up you can do.

Of course, you can invest in something like Intel IPP library or find some free code that uses the AVX2 instructions. A dot-product with two inputs and an output taking different strides is a very general function you can find to reuse.

Oh, if you’re not compiling in 64-bit, do so: you get more vector registers.

Getting good performance out of matrix multiplication is tricky. Many concerns have been addressed already, but there are a couple of big ones that I have not seen mentioned.

The main issue in medium-large matrix multiplication (for tiny matrices it's 100% micro-optimization) has historically been, still is, and will always be: reuse of data at several levels. The most obvious is at the cache level, but even at the register level you must economize.

For example, on a modern Intel processor (similar principles are fairly universal, but that's where the exact numbers will come from) it is possible to execute two vector FMAs (fused multiply-add) in a cycle, and it is possible to load two vectors from memory provided the loads hit L1 cache. Two FMAs together have 6 inputs, but you can load only 2 things, so in order to get anywhere near a good performance there has to be significant reuse. In a dot-product, you can reuse the summation variable(s), but that is not enough, that still requires twice as many loads as can possibly be executed.

So a perhaps surprising conclusion is that you cannot do just one dot-product at the time. There would have to be several in progress side by side, to enable sufficient reuse of data. Also you might think about different computational structures than dot products, such as small outer products (of rectangular, yet not too narrow, tiles).

An other (but related) issue is that while many FMAs can be executed per unit time, any individual FMA is actually quite slow in terms of latency. Of course it is not surprising that it should take at least one cycle, and even that is already too long to be executing two FMAs every cycle with a naive dot-product. The problem there is that the sum variable takes a while to be updated, and this will delay the next FMA if the code is written the obvious way. Actually a latency of 4 or 5 cycles and a throughput of 2/cycle is common enough (Haswell through Skylake, Ryzen). That means a simple dot-product, even if vectorized (scalar is not even worth considering), can be 8 to 10 times as slow as a proper implementation, though "fast math" flags can alleviate this. Anyway this problem also disappears when performing a sufficient number of simultaneous dot-products - every individual one of them will still be about as slow as a naive implementation, but the total throughput can be increased to near the theoretical maximum.

Once you have written a fast kernel like that (for reference, "fast" means close to 32 sp fp-ops per cycle on Haswell-Skylake and close to 16 on Ryzen), there are still more levels of data reuse that must be taken care of to get good results for large matrices. Loop blocking is often mentioned, but what's less commonly noted (and actually super important) is repacking the current tile from the matrix into a contiguous block. Without repacking, the tile may indeed fit in the cache, but may span more pages than the TLB can simultaneously cache address translations for. What the right tile size and shape ratio are, is also an interesting question that I frankly don't know the answer to, except to experiment.

I don't really have great ready-made code to show you for this, it really is difficult and takes a lot of engineering and experimentation to get right, taking care of more (and weirder) details than I addressed. The good implementations that you can find elsewhere have had a lot of work put into them. I did give it a try, but the results were not as close to what they should be as I had hoped, reaching 28.5 sp fl-ops per cycle at best on as Haswell, and that was after writing a highly unrolled loop in assembly using all the vector registers and with various dirty tricks such as loading some data an iteration in advance and pre-computing addresses to avoid some complex memory operands with associated µop delamination, and despite my best efforts with tiling/repacking I still saw a performance drop-off for larger matrices.

I think the most important thing for a programmer is being efficient, and that means mostly not reinventing the wheel. I would understand if you had done this code just for fun, but you say it will be used a lot. So my advice is to use some existing matrix library like BLAS instead. It will be faster than whatever you might write since some parts are likely written in assembly.

The constructor could be much shorter:

Matrix(int w, int h) : height(h), width(w)
{
for (int i = 0; i < height; i++) {
elements.push_back(std::vector<float>(width, 0.0f));
}
}


And since height and width are now in the initialization list, they can be declared as const members. I'm using the vector constructor which can set all values to 0.0f.

I would rename get_element to just get, and same with setter. Maybe operator overloading with [] would be possible instead.

Others already made valid points, especially regarding using vector<vector<>> instead of a simpler double[][] or such.

• When you are constructing something and pushing it into a vector, it is often better to use emplace_back to avoid temporaries: elements.emplace_back(width, 0.0f); – Josiah May 6 '18 at 7:33
• There’s no need for any loop. Just initialise elements directly: : elements(std::vector<float>(w, 0.0f), h) — and get rid of width and height. – Konrad Rudolph May 6 '18 at 14:17
• +1 Unless the project is to learn how to optimize MM, use a real implementation. – Alex Reinking May 6 '18 at 19:52
• Why are you using the float type? Almost all math processors today are natively 64-bit (double), so using a 32-bit float doesn't save any speed. In addition, if you're running on a 64-bit machine, the float will be padded with 32 bits in order to keep things word aligned, so you won't save any space.
• Why aren't you using templating for the Matrix so that the caller can choose which type they want inside the matrix?
• Why aren't you specifying the size of the std::vector when you construct it? The vector constructor takes a preferred size and default value, so you can also preload it with zeroes.
• But why are you even using std::vector at all? The array size is fixed, you could just use std::array.
• dot() should be either a static member of Matrix or a friend function. This is necessary in order so that...
• elements, height, and width should be private: so that they can only be touched by class methods. You don't want other users of the class to mess with them.
• Why not have a Matrix constructor that takes the 2d array?
• float padding is implementation defined, as is the performance of float vs double. On many 64 bit systems the padding is still 4 bytes, so using floats can save space (as no padding will be used), particularly when used with SIMD instructions. – 1201ProgramAlarm May 4 '18 at 21:01
• I have never worked with a system where floats are padded to 8 bytes. The same as you would never expect char to be padded to 8 bytes. Maybe such padding is allows by the standard, but Intel-compatible processors certainly don’t require it. – Cris Luengo May 5 '18 at 1:12
• The array size might be fixed, but it’s fixed at runtime. Sure, in OP’s code the width and height are constant but who knows what the real use-case is? Anyway, the variable size isn’t the only advantage of a std::vector over a std::array. The other advantage is the location of the storage. For slightly larger matrices, the automatic storage of std::array is unsuitable, and you’d resort to using a std::unique_ptr<T[N]>. Might as well use a std::vector. – Konrad Rudolph May 6 '18 at 14:12

Try using openMP. It's compiler-directive-based and you can get huge performance gains by simply using a directive like:

#pragma omp parallel for


Large matrices may call for GPU processing or even just using SIMD operations on your CPU. Matrix multiplication is also one of the often illustrated examples for parallel processing libraries so you won't even have to dive very deep into the libraries to get the performance gains.