Laugh at your own risk: A couple of threads walk into a bar, order some drinks, and have a great time. Later, they realise that they have run out of cache and won't be able to remember this core memory.
I have been building my very own linear algebra library for the past few weeks. At first, no effort was initially made to make the operations performant, just that they did what they were supposed to. After finishing the initial version ( -(tests + documentation + a few functions) ), I was ready to dip my toes into optimising the functions.
Before I get into that, I must confess that I have been so stressed out about my time management lately, that I have started binge watching C++ tutorials while I eat my meals. However bad that sounds, it's probably worse. This gluttonous consumption did have a silver lining - I learnt about multithreading in C++ (it's about time honestly).
Nick aka CoffeeBeforeArch is my favourite programming educator on YouTube. He only teaches C++, and the way he does it makes C++ seem like a language simpler than the girthy Python itself. From him I first learnt about std::thread, and it will be my honour to pass on his knowledge to you from the perspective of the challenges I overcame using it.
For executing a basic program, certain resources are assigned. To be particular, a thread. A hardware thread (HT from now) is a small part of the CPU core (there are typically 2 threads per core) that can be assigned tasks to execute while leaving the rest of the CPU resources (the other threads) free for other tasks. Now, depending on the system, these resources can generally be in high demand (consistently high CPU usage), or be left idle for most of the time. Somehow there is no in between. My situation fell on the latter end of the spectrum.
The threads can be thought of as workers. Initially, I had been writing programs (including this library) that was just using a single worker to execute tasks. Nick taught me that we can intentionally harness the power of the extra threads available to us to divide the work between these "workers".
That is literally all multithreading is about.
If you have heard of intimidating terms like parallelization and high-performance computing (HPC), this is literally that. Divide the task between many workers that execute simultaneously.
Coming to std::thread, that is how C++ is able to spawn OS threads (OST from now). OSTs bridge the gap between hardware threads and software. They are a resource which needs its own stack memory, scheduling priorities among other stuff. The operating system is responsible for distributing OSTs among the available HTs. This means that there can be more OSTs than HTs, in which case the OS must use time-slicing to allocate hardware fairly among all OSTs. Later, you will see why crossing the hardware limit isn't always the right option.
The way std::thread works is very simple. It takes 1 + n arguments - a callable, and then the callable's arguments. An OST is then allocated which executes the callable separately from the main thread (which is responsible for the overall main program).
#include <iostream>
#include <thread>
int main() {
int r1{ 0 }, r2{ 1 };
auto func = [](int& accumulator) {
for (int i = 0; i < 1000; i++) {
accumulator += i;
}
};
std::thread t1(func, std::ref(r1));
std::thread t2(func, std::ref(r2));
t1.join(); t2.join();
std::cout << "r1 = " << r1 << "\n";
std::cout << "r2 = " << r2 << "\n";
int total = r1 + r2;
std::cout << "total = " << total << std::endl;
return 0;
}
r1 = 499500
r2 = 499501 // r2 was initialized to 1!!
total = 999001
Let's go through this example:
- I initialize
r1andr2, the two objects to be modified. funcis a pretty straightforward lambda. The sum of all whole numbers in the range[0, 999]is added to the accumulator. Notice that the accumulator is passed by reference, so the changes naturally reflect in the outer scope.- The two threads are initialized with the function name,
r1, andr2are passed in the thread constructor. These threads begin executing their instructions immediately upon initialization.- There is a real reason why
r1andr2are wrapped instd::ref(). I learnt about this while practicing usingstd::thread. Unlike normal functions where arguments are automatically passed as references if defined accordingly,std::threadscopy lvalues liker1andr2to have safe-to-use copies. This circumvents situations where the lvalues go out of scope before the thread actually starts executing. This makes it necessary to explicitly wrap them instd::ref()to force them to be processed as references.
- There is a real reason why
- The thing about using
std::threadis that while it executes tasks on the side, the primary program thread continues executing what is next. Which is why situations can arise where the main thread reaches the end of the program before a thread is done with its task.- In order to avoid this, waiting for the threads to complete their execution is necessary. Setting up this waiting point at the latest possible point in the program (before the results are actually needed) is best for performance, so that the main program mustn't wait prematurely. For this we use
t1.join()andt2.join(), which basically tells the main thread to wait for the parallel threads and their results, which will effectively "join" them back to the main thread.
- In order to avoid this, waiting for the threads to complete their execution is necessary. Setting up this waiting point at the latest possible point in the program (before the results are actually needed) is best for performance, so that the main program mustn't wait prematurely. For this we use
Just to hammer down the value proposition of multi-threaded execution, here is what would happen in a sequential, non-multithreaded program - The two loops will be executed one after the other, taking almost double the time that is taken by the multithreaded version. I underlined 'almost' because it is doing a lot of heavy lifting here. More on performance later.
Whew, we can finally end this detour and focus on the problem initially posed. Matrix multiplication is a fundamental operation for any linear algebra library. So getting it as fast as possible is a matter of priority, especially as matrices scale to 1000s of dimensions when solving real deep learning problems. LLMs run on NNs. And NNs are nothing without matrix multiplication.
The interesting thing is that I use the exact naive approach almost everyone first learns:
A (n x m) * B (m x p) = C (n x p)
for row_idx from 1 to n
for col_idx from 1 to p
sum <- 0
for inner_idx from 1 to m
sum <- sum + A[row_idx][inner_idx] * B[inner_idx][col_idx]
C[row_idx][col_idx] <- sum
Here was my first attempt at the multithreaded algorithm:
// This function is part of the Matrix class template
template <typename T>
Matrix<T> fast_mult(const Matrix<T>& m1, const Matrix<T>& m2) {
if (m1.cols != m2.rows)
throw std::invalid_argument("Matrix dimensions incompatible.");
Matrix<T> result(m1.rows, m2.cols);
auto l = [&m1, &m2, &result](size_t row, size_t col) {
auto res = T{};
for (size_t i = 0; i < m1.cols; i++)
res += m1(row, i) * m2(i, col);
result(row, col) = res;
};
std::vector<std::jthread> jobs;
for (size_t r = 0; r < m1.rows; r++)
for (size_t c = 0; c < m2.cols; c++)
jobs.emplace_back(l, r, c);
return result;
}
This gets really stupid really fast:
- As I explained before, a thread/worker needs a task, and we need to define a blueprint for that task in the form of a callable. My callable is designed such that it works on a single row-column dot product on each call, as illustrated by its parameters. Therefore, according to this design, each new thread will handle a single dot product.
- I then create a vector that will hold my threads. I make use of
std::jthreadinstead of vanillastd::thread, as that allows me to not write.join()'s. It automatically takes care of waiting and merging. By the way, try not to usestd::thread. RAII for life! - In the loop, for each column in each row, I
emplace_back()astd::jthreadin which the callable, row and column are passed. - For a 3x3 times 3x3, I will have spawned 9 threads for 9 dot products. That is just a tiny bit of addition and multiplication per thread. Clearly, I wrote about this program before even understanding the cost of
std::thread.
See, operations like addition are executed extremely efficiently by the CPU. That is literally one of its main jobs. But the task given to the thread isn't the entire work itself. Spawning the OST itself is also a mighty task. The request for an additional worker must be given to the OS and upon availability, an additional thread becomes ready to execute tasks within microseconds (10⁻⁶). Compare that to the nanoseconds (10⁻⁹) it takes for the CPU to speed through sums and products. This makes std::thread a double edged sword. It can cut down work done per unit time significantly, but getting one ready in the first place is quite expensive.
In this case, this means that the CPU actually spent most of the execution time spawning and getting the threads ready, not calculating the actual dot products. Through this example, I hope you can yourself reach to the conclusion that std::thread is best used when the task significantly outweighs the spawning costs.
But wait, there is more - HTs (hardware threads, if you were skimming) are a limited resource. I've got an CPU with 4 p-cores (2 threads/core) and 8 e-cores (1 thread/core) on my laptop, so that is a maximum of 16 threads. Your mileage may vary. This is why the target system is also an important variable to consider before integrating multithreading into your program. Usually, the more workers you can get to complete complex calculations, the better. But the final performance will naturally vary due to the hardware constraints.
My code did not care about the number of threads available to me. I just pushed more and more threads depending on the dimensions of the matrices. For a 4x4 times 4x4, that is 16 threads, just within my personal hardware limits. But any bigger than that, and two things would drag me down - The first is the resource cost of spawning a thread (stack memory, spawn time etc.). Secondly, each HT will have to context switch between multiple OSTs in order to maintain the illusion of simultaneous execution, which introduces a lot of overhead.
At this point, before reading any further, I would like to invite you to think of ways this code could actually work and properly use the resources (threads) available.
If you thought of chunking/grouping work based on the matrix dimensions or something along those lines, you get a C++ grade!
Here is that approach:
// This function is part of the Matrix class template
template <typename U>
Matrix<T> operator*(const Matrix<U>& other) const {
if (this->cols != other.rows)
throw std::invalid_argument("Matrix dimensions incompatible");
auto thread_count = std::thread::hardware_concurrency();
float row_thread_ratio = static_cast<float>(this->rows + thread_count - 1) / static_cast<float>(thread_count);
int chunks = row_thread_ratio > 4 ? static_cast<int>(row_thread_ratio) : 4;
Matrix<T> result(this->rows, other.cols);
std::vector<std::jthread> workers;
for (unsigned int t = 0; t < thread_count; t++){
size_t start_row = t * chunks;
size_t end_row = std::min(start_row + chunks, this->rows);
if (start_row < end_row) {
workers.emplace_back([&, start_row, end_row]() {
for (size_t r = start_row; r < end_row; r++) {
for (size_t c = 0; c < other.cols; c++) {
T res{};
for (size_t k = 0; k < this->cols; k++) {
res += (*this)(r, k) * other(k, c);
}
result(r, c) = res;
}
}
});
}
}
return result;
}
And here is the breakdown:
- Firstly, in order to check the number of threads available to me, I use
std::thread::hardware_concurrency(), which returns an integer amount. - I ceiling divide the number of rows in the left matrix by this thread count to get the
row_thread_ratio. This ratio will determine the amount of work to be handled by each thread. - Notice that for small matrices, this ratio can be quite small. 8 threads for a 4x4 times 4x4 matrix multiplication equates to a ratio of half a row per thread. That is just stupid, because it shows that (a) we don't need that many threads anyways and (b) processing half a row at a time is plain inefficient from the CPU's perspective.
- That is why I make use of a minimum chunk size of 4. For example, 30 rows and 8 threads mean a chunk size less than 4. Instead of having each thread handle way too small chunks, I set a minimum chunk size so that each thread can actually save on the execution time.
- Again, I create an
std::vectorto storestd::jthread's. - Then, a loop goes through the number of available threads and allocates the correctly-sized chunk of work to each thread. As you can see, this time, I have written the lambda directly inside the thread creation statement. Also notice that I haven't passed any extra arguments inside
.emplace_back()because the required values are captured by the lambda instead. - Apart from the different looping ranges, the multiplication algorithm is basically the same.
I feel proud to share the performance comparison with you. Both versions were run with the highest compiler optimisation settings, vectorization enabled, and my CPU (i5-12500H Laptop) on turbo -
Multiplying two 1180x1180 matrices together is 8 times faster with the multi-threaded version. Again, this result is system dependent but the takeaway is clear - There is massive potential for performance improvements by implementing concurrency.
The volatility you see in the red curve by the end is due to my CPU throttling. I was pushing it to its limits, with temperatures at a consistent (and quite harmful) 93-95 degrees Celsius. But I would do it again. For science :) The alignment with the theoretical curve is still clear. It is not discernible in this graph, but before n=160, the naive approach is consistently faster, with much higher performance in multiplying smaller matrices. But now you know why that is. You do remember, right?
With this result, I was ready to upgrade my matrix multiplication function to use multithreading. But a question remained - the two matrices don't necessarily need to be square, so what would be the threshold to cross for me to opt for the multithreaded version? At first, I thought of using 160 x 160 = 25600 as my threshold, because that would be the number of dot products (160 rows of Matrix A times 160 columns of Matrix B) that would be carried out.
But that is only part of the full picture. How many pairs would I be multiplying and adding together in that dot product? For a 160x160 matrix, that's 160 multiply-add operations per dot product and over 4 million in total.
But in the case of 160x2 times 2x160, that is just 160 x 2 x 160 = 51200 operations! But if I had used my older threshold, multithreading would have triggered for a matrix that doesn't need it at all, leading to poor performance. This means that using the total number of operations, rather than number of dot products, is the right way to choose between standard and multithreaded approaches.
For an $M \times N$ matrix multiplied by an $N \times P$ matrix, the total operation count is $M \times N \times P$. And since I know the threshold for my system was 160^3 operations, that is the number I used.
For transparency, here is the script I used for benchmarking and graphing the results, in C++ and Python respectively -
int main() {
Timer t, g;
for (int i = 5; i < 1200; i += 25) {
auto m1 = mat_random_int_range(i, i, 1, 9);
auto m2 = mat_random_int_range(i, i, 1, 9);
float sum_naive = 0.0f;
for (int repeats = 0; repeats < 10; repeats++) {
t.reset();
auto slow = m1 * m2;
sum_naive += t.elapsed();
}
std::cout << i << ", " << sum_naive / 10.0f << std::endl;
}
for (int i = 5; i < 1200; i += 25) {
auto m1 = mat_random_int_range(i, i, 1, 9);
auto m2 = mat_random_int_range(i, i, 1, 9);
float sum_multi = 0.0f;
for (int repeats = 0; repeats < 10; repeats++) {
g.reset();
auto fast = m1.fast_mult(m2);
sum_multi += g.elapsed();
}
std::cout << i << ", " << sum_multi / 10.0f << std::endl;
}
}
- If you are wondering why I am using two different loops that are almost doing the same thing, I wanted to separate the two modes of execution (single thread vs multithread) so that each result isn't influenced by the other.
- Also, in this benchmark,
fast_mult()uses the multithreaded version and the*operatoruses the standard one.
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-deep')
import pandas as pd
col_sizes = "n_x"
col_naive = "time_x"
col_multi = "time_y"
data = pd.read_csv("data.csv")
data.columns = data.columns.str.strip()
print(data.head())
sizes = data[col_sizes]
naive = data[col_naive]
mt = data[col_multi]
constant = naive.iloc[-1] / (sizes.iloc[-1] ** 3)
ideal_times = constant * (sizes ** 3)
plt.plot(sizes, ideal_times, '--c', label='O(n^3) reference', linewidth=1)
plt.plot(sizes, naive, 'r', label='Naive', linewidth=2)
plt.plot(sizes, mt, 'b', label='Multi-threaded', linewidth=2)
plt.vlines(160, 0, 4, 'g', linestyles=":", label="Diverging point (n≈160)")
plt.xlabel("nXn Matrix Multiplication")
plt.ylabel("Execution Time (s)")
plt.legend(frameon=True, loc='upper left')
plt.title("Matrix Multiplication: Naive vs Multi-threaded")
plt.grid()
plt.tight_layout()
plt.show()