I will assume that you've already made the changes suggested by L. F., especially the data layout and passing the Matrix
objects by reference. You asked for a performance review too, so here it is.
Baseline benchmark
For timing, I used std::chrono
, and took the average of 10 runs, while testing the multiplication of two equal-size square matrixes with a power-of-two side-length. For example a 64x64 matrix times an other 64x64 matrix. The matrixes were just filled with ones. The actual value is not important, but leaving memory zeroed out can have funny effects sometimes. The times in the table are the time per element of the result. Since normal matrix multiplication is an O(n³) time algorithm with O(n²) output elements, a reasonable hypothesis could be that those times increase linearly with the size.
Results (in microseconds per output element):
size baseline
64 0.89
128 0.37
256 0.21
512 0.27
1024 1.51
2048 ????
4096 ????
These times very much do not increase linearly with size, it's more of a "bath tub curve". For large matrixes, my entire computer lagged so hard that the cursor stopped moving.
#1, too many threads
The cause of the lag, and also the source of a lot of the overhead for small matrixes (where each thread has comparatively less work to do, even though there are fewer of them), is that there are too many threads. For example 8x8 matrix multiplication is a trivial calculation which should not have any threads created for it, and on the other end of the spectrum, a 1024x1024 matrix multiplication would create 1024 threads which is extremely excessive.
How many threads should there be? Maybe (to keep things simple) std::thread::hardware_concurrency()
, or otherwise limited some reasonable number based on the expected number of cores. These threads are (intended to be) compute-bound, so having more threads than cores won't do anything. Whether "cores" in this story is physical cores or logical cores depends on how good the code is: good code saturates a whole physical core, making SMT pointless or even harmful.
Anyway, fewer threads. Distributing the work over at most std::thread::hardware_concurrency()
threads (which is 8 on my system, for 4 physical cores with 2-way SMT) but also ensuring that small matrixes didn't get too many threads assigned to them by giving each thread at least 64 rows to handle, now the timing table looks like this (still in microseconds per output element):
size baseline #1
64 0.89 0.11
128 0.37 0.07
256 0.21 0.08
512 0.27 0.20
1024 1.51 1.49
2048 ???? ????
4096 ???? ????
In order to have one thread handle multiple rows, an extra parameter is given and an extra loop is put around the two nested loops, so that the thread itself can loop over a range of rows.
That didn't solve the memory access pattern, so bix matrixes are still slow. Let's solve that next.
Mind the cache
Memory is slow, cache is fast. Memory often looks fast, but that's because the caches work hard to hide the slowness of memory. Memory stops looking fast if you access it in a way that makes the caches ineffective. One of those ways is accessing lots of tiny chunks of memory spaced far apart - just like accessing the column of a matrix.
There are broadly two techniques that mitigate that problem, and they should both be used:
- Chopping the matrix into parts, small enough to fit in cache entirely. But which cache? L3 is still quite slow so probably not that one, but L1 or L2? opinions are divided but IME sizing the blocks for L2 works best, L1 leaves the blocks too small: a 64x64 block is already big for L1.
- Re-using data when we've loaded it anyway. For example, if we load an extry from a column, we may as well multiply it by a couple of partial rows. That needs to happen at some point anyway, and we have the data now, we may as well use it now - instead of waiting until later and having to load the data again.
I'll divide this into two chapters..
#2, chopping the matrix into parts
In order to do this, three new loops will be introduced around the three existing nested loops. It looks a bit intimidating. Feel free to factor the inner part into its own function. For example:
const size_t rowblock = 128;
const size_t cblock = 128;
const size_t iblock = 64;
// loop over tiles of the matrixes
for (size_t rr = rowStart; rr < rowEnd; rr += rowblock) {
size_t rlim = std::min(rr + rowblock, rowEnd);
for (size_t cc = 0; cc < m2.cols; cc += cblock) {
size_t clim = std::min(cc + cblock, m2.cols);
for (size_t ii = 0; ii < m1.cols; ii += iblock) {
size_t ilim = std::min(ii + iblock, m1.cols);
// multiply tile by tile
for (size_t row = rr; row < rlim; row++) {
for (size_t c = cc; c < clim; ++c) {
int t = result(row, c);
for (size_t i = ii; i < ilim; ++i)
t += m1(row, i) * m2(i, c);
result(row, c) = t;
}
}
}
}
}
Note that t
is initialized with an element of the output, because it an element of the output is not computed all in one go this way. The computation is done in parts. That means the product is added to the result, of course the result can be zeroed out first to avoid the effect of that, and it's implicitly zeroed out (when creating the result matrix) in this program anyway.
Intimidating or not, it clearly helps:
size baseline #1 #2
64 0.89 0.11 0.11
128 0.37 0.07 0.07
256 0.21 0.08 0.08
512 0.27 0.20 0.12
1024 1.51 1.49 0.30
2048 ???? ???? 0.78
4096 ???? ???? ????
Below 512 there is essentially no difference, which is expected. Actually, for even smaller matrixes we should expect a slowdown, since there is more overhead in the code now and small matrixes get "chopped" into one piece anyway.
The improvement for matrixes of size 512x512 is already almost a factor of 2, and for 1024x1024 it's even better. 2048 was also within my "patience limit" with this speedup, sadly there is no earlier score to compare it with.
This could be tweaked further if you want. For example, above a certain size it becomes interesting to repack tiles (making them more linear in memory, rather than spread out in a big matrix) to reduce TLB misses. When done well, the time per output element should scale approximately proportionally to the size of the matrix, rather than the slightly worse scaling seen here.
#3, unrolling the row
and c
loops
The other aspect of data reuse can be accompished by unrolling the row
and c
loops somewhat. Interestingly, unrolling the i
loop is not very interesting, despite being the innermost loop.
Unrolling two loops by 2 and 3 respectively means that in the inner loop, there will be 2 elements to load from m1
and 3 elements to load from m2
, and using these 5 elements there will be 6 different products to compute. As you can see, this significantly reduces the ratio of loads-to-arithmetic.
The inner loop might look like this, there are some other changes to the code but they're boring and the code is quite large.
for (size_t i = ii; i < ilim; ++i) {
int m1A = m1(row, i);
int m1B = m1(row + 1, i);
int m2A = m2(i, c);
int m2B = m2(i, c + 1);
int m2C = m2(i, c + 2);
t0 += m1A * m2A;
t1 += m1A * m2B;
t2 += m1A * m2C;
t3 += m1B * m2A;
t4 += m1B * m2B;
t5 += m1B * m2C;
}
The block sizes can be tweaked again (the unrolling slightly changes what the best sizes are) to get the times down to the ones shown in column #3B (the result for 512 actually gets worse though). Optimal block sizes depend on the chache sizes of your CPU though.
const size_t rowblock = 64;
const size_t cblock = 128;
const size_t iblock = 32;
Times:
size baseline #1 #2 #3 #3B
64 0.89 0.11 0.11 0.12 0.11
128 0.37 0.07 0.07 0.07 0.06
256 0.21 0.08 0.08 0.06 0.05
512 0.27 0.20 0.12 0.09 0.09
1024 1.51 1.49 0.30 0.20 0.18
2048 ???? ???? 0.78 0.42 0.41
4096 ???? ???? ???? 0.95 0.90
#4, repacking
Accessing a tile of a matrix still has an issue, even if it fits into say the L2 cache, it may be "spread out" across many 4KB pages. Each 4KB page (or 2MB page, if large pages are used, but they're difficult to use) costs an TLB entry, and the TLB is not that big. There is a tile that's accessed by column, which would hit all the TLB entries one by one, and then start the cycle over again: not good. This can be fixed by loading the tile into contiguous memory once, and then do various loops over it. By the way, in order to do this I switched the order of the outer 3 loops a bit, to ensure that this tile repacking is done as little as possible (nice as it is, it does represent an overhead).
Times:
size baseline #1 #2 #3 #3B #4
64 0.89 0.11 0.11 0.12 0.11 0.11
128 0.37 0.07 0.07 0.07 0.06 0.06
256 0.21 0.08 0.08 0.06 0.05 0.05
512 0.27 0.20 0.12 0.09 0.09 0.07
1024 1.51 1.49 0.30 0.20 0.18 0.15
2048 ???? ???? 0.78 0.42 0.41 0.28
4096 ???? ???? ???? 0.95 0.90 0.52
As you can see, above a certain size it suddenly starts to matter a lot.
Is that all?
No. 0.52 µs per entry for 4096² entries is about 8.7 seconds. I estimate that, even staying with scalar code, it could be done twice as fast on my PC, based on:
- 4096 multiplications per element
- 1 multiplication per cycle per core
- 4 cores
- 4GHz
- (4096 / 4 / 4GHz) * 4096² = 4.295 seconds
Why is this code twice as bad as it should be? First of all the codegen by MSVC 2019 is not as good at it could be. Probably some of that is my fault for not writing the code in the way the compiler prefers it, but I could not find any way that made it do the right thing. Let's take a look
MultiplyRow3+240h:
mov edx,dword ptr [rcx+r12]
mov r11d,dword ptr [r12]
lea r12,[r12+4]
mov r8d,dword ptr [r15]
mov ecx,r8d
mov r9d,dword ptr [r15+4]
mov r10d,dword ptr [r15+8]
add r15,0Ch
imul ecx,edx
imul r8d,r11d
add dword ptr [t0],ecx
mov ecx,r9d
add dword ptr [t3],r8d
imul ecx,edx
imul r9d,r11d
add dword ptr [t1],ecx
mov ecx,r10d
add dword ptr [t4],r9d
imul ecx,edx
imul r10d,r11d
add dword ptr [t2],ecx
add dword ptr [t5],r10d
mov rcx,qword ptr [rsp+50h]
sub r14,1
jne MultiplyRow3+240h
I disagree with the decision to put all t
-variables on the stack. While the register pressure is high, it is made higher here by unfortunate instruction scheduling, and some general purpose registers are still unused within the loop body: rax, rbx, rsi, rdi, r13
. I don't think that should make the loop fully twice as slow as intended, but LLVM-LCA suggests that it this contributes to making this loop take over 9 cycles per iteration instead of approximately 6 (which it takes if I replace the memory destinations by registers). That's not quite a 2x difference, but it's about 1.5x so the hypothetical remaining difference is only a factor of 1.3, but that's uncertain in practice, and one inefficiency can hide another.
Some of the remaining difference can be found by tweaking the block sizes, but that has the annoying property of usually not resulting in an across-the-board improvement: often the high sizes get better at the expense of mid-sizes, etc.
Further, the limit actually isn't 4.295 seconds, but even less. Basically all modern processors have SIMD, which enables them to do more calculations with fewer instructions. I've used that in the answer to this other integer matrix multiplication question. The code there got quite close to the theoretical single-threaded performance limit of my PC (perhaps partly because the compiler has an easier job: most of the variables are vectors, which use a separate set of registers, so there's less of a fight over scarce general purpose registers). However, SIMD is a niche subject even for experts, so I expect it is a bit advanced for this question.