Multithreaded Matrix Multiplication in C++, improving efficiency?

Question

Here is a very simple Matrix Multiplication console application using multiple threads (one per row), I'm just looking for feedback to keep improving :) thank you in advance

    #include <iostream>
    #include <fstream>
    #include <sstream>
    #include <vector>
    #include <string>
    #include <thread>

using namespace std;

struct Matrix {
    unsigned int cols = 0;
    unsigned rows = 0;
    int **data;

    Matrix() {};
    Matrix(int r, int c) {
        cols = c;
        rows = r;
        data = new int* [rows];
        for (int i = 0; i < rows; ++i) {
            data[i] = new int[cols];
            for (int j = 0; j < cols; ++j) 
                data[i][j] = 0;            
        }
    }
};

Matrix ReadMatrix(string fileName) {    
    fstream file(fileName, fstream::in);
    Matrix matrix;
    try {
        if (file.is_open()) {            
            int rows, cols;
            file >> rows >> cols;
            matrix = Matrix(rows, cols);

            for (int r = 0; r < matrix.rows; ++r) 
                for (int c = 0; c < matrix.cols; ++c) 
                    file >> matrix.data[r][c];                
            
            file.close();
            cout << "File: " << fileName << " read.\n";
        }
        else cout << "could not open file/matrix: " << fileName << "\n";             
        return matrix;
    }
    catch (exception ex) {
        if (file.is_open())file.close();
        throw ex;
    }
}

void MultiplyRow(Matrix m1, Matrix m2, Matrix result, int row) {    
    for (int c = 0; c < m2.cols; ++c)
        for (int i = 0; i < m1.cols; ++i) 
            result.data[row][c] += (m1.data[row][i] * m2.data[i][c]);       
    string str = ("Row: " + to_string(row) + " complete.\n");
    cout << str;
}


Matrix MultiplyMatrices(Matrix m1, Matrix m2) {
    try {
        Matrix matrix;
        if (m1.cols == m2.rows) {        

            vector<thread> threads;            
            matrix = Matrix(m1.rows, m2.cols);            
            
            for (int i = 0; i < m1.rows; ++i)
                threads.emplace_back(MultiplyRow, m1, m2, matrix, i);

            for (int i = 0; i < threads.size(); ++i)
                threads[i].join();

            cout << "Matrices multiplied.\n";
        }
        else cout << "cannot multiply matrices\n";
        return matrix;
    }
    catch (exception ex) {
        throw ex;
    }
}

void WriteMatrixToFile(int **data, int rows, int cols, string fileName) {
    try {
        ofstream file(fileName);
        char SPACE = ' ';
        if (file.is_open()) {
            file << rows << SPACE << cols << "\n";
            for (int r = 0; r < rows; ++r) {
                for (int c = 0; c < cols; ++c) {
                    file << data[r][c] << SPACE;
                }
                file << "\n";
            }
        }
        cout << "Matrix written to: " << fileName << ".\n";
        file.close();
    }
    catch (exception ex) {
        throw ex;
    }
}

int main()
{
    //read matrices
    Matrix matrix1 = ReadMatrix("matrix1.txt"); 
    Matrix matrix2 = ReadMatrix("matrix2.txt");

    //multiply matrices
    Matrix result = MultiplyMatrices(matrix1, matrix2);

    //put result into results file
    WriteMatrixToFile(result.data, result.rows, result.cols, "result.txt");
}

Answer 1

Welcome to Code Review!

Algorithm

Currently, a thread is started for every row of the result matrix. This works best if the matrix has a large number of columns and the number of rows is close to the number of native threads.

A more efficient strategy in general might be to partition the rows to several chunks according to the native thread count instead. Without careful benchmarking, it's hard to say for sure.

Encapsulation

Matrix has a class invariant — for a non-empty matrix, data points to an array of rows pointers that point to arrays of cols elements each. To ensure synchronization, the data members can be kept private.

Memory issues

Matrix fails to observe the rule of zero:

Classes that have custom destructors, copy/move constructors or copy/move assignment operators should deal exclusively with ownership (which follows from the Single Responsibility Principle). Other classes should not have custom destructors, copy/move constructors or copy/move assignment operators.

This rule also appears in the C++ Core Guidelines as C.20: If you can avoid defining default operations, do.

To fix this issue, simply change the pointers to std::vector.

Storage of matrices

Currently, you are storing a pointer to a dynamic array of pointers, which, in turn, point to dynamic arrays, incurring double indirection overhead and causing memory segmentation. A more efficient method is to store all elements in a flattened array and provide a helper function to access the element at a given place:

class Matrix {
    std::vector<int> elements;
    std::size_t n_rows;
    std::size_t n_cols;
public:
    Matrix()
        : Matrix(0, 0)
    {
    }
    Matrix(std::size_t n_rows, std::size_t n_cols)
        : elements(n_rows * n_cols), n_rows(n_rows), n_cols(n_cols)
    {
    }
    
    std::size_t row_count() const noexcept {
        return n_rows;
    }
    std::size_t col_count() const noexcept {
        return n_cols;
    }
    
    int& operator()(std::size_t i, std::size_t j) {
        return elements[i * n_cols + j];
    }
    const int& operator()(std::size_t i, std::size_t j) const {
        return elements[i * n_cols + j];
    }
};

Note that this does not preserve the aliasing semantics of your Matrix and performs deep copying instead, so the by-value parameters should be changed to references.

Miscellaneous

using namespace std; is considered bad practice in C++.

If you turn on high warning levels, you will see warnings about mixing signed and unsigned integers in comparisons, which is because you are mixing int and unsigned int. Instead of them, use std::size_t to store sizes, which is what the language functionalities and the standard library use.

In ReadMatrix, file streams close automatically when they go out of scope, so there's no need to do that manually. std::ifstream can be used to simplify the construction of the file stream.

Instead of printing messages, throw exceptions when errors occur.

Answer 2

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:

1. 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.
2. 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.