Portable multithreaded matrix multiplication in C

Question

I essentially took this progam, and made it compile and run on Windows too. Also, I incorporated some nice points made by Toby Speight.

matrix.h

#ifndef NET_CODERODDE_LINEAR_ALGEBRA_MATRIX_H
#define NET_CODERODDE_LINEAR_ALGEBRA_MATRIX_H

#include <stdlib.h>

typedef struct matrix_t {
    size_t  m_rows;
    size_t  m_cols;
    double* m_data;
} matrix_t;

matrix_t* matrix_t_alloc (const size_t rows, const size_t cols);
void      matrix_t_init  (matrix_t* matrix, size_t rows, size_t cols);
void      matrix_t_clear (matrix_t* matrix);
void      matrix_t_free  (matrix_t* matrix);
double    matrix_t_get   (matrix_t* matrix, size_t x, size_t y);
void      matrix_t_set   (matrix_t* matrix, size_t x, size_t y, double value);
matrix_t* matrix_t_copy  (matrix_t* matrix);

#endif /* NET_CODERODDE_LINEAR_ALGEBRA_MATRIX_H */#pragma once

matrix_algorithm.h

#ifndef NET_CODERODDE_LINEAR_ALGEBRA_MATRIX_ALGORITHM_H
#define NET_CODERODDE_LINEAR_ALGEBRA_MATRIX_ALGORITHM_H

#include "matrix.h"

matrix_t* matrix_t_multiply(matrix_t* matrix1, matrix_t* matrix2);
matrix_t* matrix_t_multiply_parallel(matrix_t* matrix1, matrix_t* matrix2);
void      matrix_t_print(matrix_t* matrix); 

size_t get_number_of_processors();

#endif /* NET_CODERODDE_LINEAR_ALGEBRA_MATRIX_ALGORITHM_H */

matrix.c

#include "matrix.h"
#include <stdlib.h>
#include <string.h>

static size_t data_index(matrix_t* matrix, size_t x, size_t y)
{
    return y * matrix->m_cols + x;
}

matrix_t* matrix_t_alloc(size_t rows, size_t cols)
{
    matrix_t* m = malloc(sizeof *m);

    if (!m)
    {
        return m;
    }

    matrix_t_init(m, rows, cols);

    if (!m->m_data)
    {
        free(m);
        return NULL;
    }

    return m;
}

void matrix_t_init(matrix_t* matrix, size_t rows, size_t cols)
{
    matrix->m_data = malloc(sizeof *matrix->m_data * rows * cols);

    if (matrix->m_data)
    {
        matrix->m_rows = rows;
        matrix->m_cols = cols;
    }
    else
    {
        matrix->m_rows = 0;
        matrix->m_cols = 0;
    }
}

void matrix_t_clear(matrix_t* matrix)
{
    free(matrix->m_data);
    matrix->m_rows = 0;
    matrix->m_cols = 0;
}

void matrix_t_free(matrix_t* matrix)
{
    matrix_t_clear(matrix);
    free(matrix);
}

double matrix_t_get(matrix_t* matrix, size_t x, size_t y)
{
    return matrix->m_data[data_index(matrix, x, y)];
}

void matrix_t_set(matrix_t* matrix, size_t x, size_t y, double value)
{
    matrix->m_data[data_index(matrix, x, y)] = value;
}

matrix_t* matrix_t_copy(matrix_t* matrix)
{
    size_t data_len = sizeof(double) * matrix->m_rows * matrix->m_cols;
    matrix_t* copy = malloc(sizeof(*copy));
    copy->m_rows = matrix->m_rows;
    copy->m_cols = matrix->m_cols;
    copy->m_data = malloc(data_len);
    memcpy(copy->m_data, matrix->m_data, data_len);
    return copy;
}

matrix_algorithm.c

#include "matrix.h"
#include "matrix_algorithm.h"
#include <stdio.h>
#ifdef _WIN32
#include <windows.h>
#else
#include <pthread.h>
#include <unistd.h>
#endif

static size_t max_size_t(size_t x, size_t y)
{
    return x > y ? x : y;
}
static size_t min_size_t(size_t x, size_t y)
{
    return x < y ? x : y;
}

size_t get_number_of_processors()
{
#ifdef _WIN32
    SYSTEM_INFO info;
    GetSystemInfo(&info);
    return (size_t)info.dwNumberOfProcessors;
#else
    return (size_t)sysconf(_SC_NPROCESSORS_ONLN);
#endif
}

matrix_t* matrix_t_multiply(matrix_t* matrix1, matrix_t* matrix2)
{
    matrix_t* result;
    size_t x;
    size_t y;
    size_t i;
    double sum;

    if (!matrix1 || !matrix2)
    {
        return NULL;
    }

    if (!matrix1->m_data || !matrix2->m_data)
    {
        return NULL;
    }

    if (matrix1->m_cols != matrix2->m_rows)
    {
        return NULL;
    }

    result = malloc(sizeof(*result));
    matrix_t_init(result, matrix1->m_rows, matrix2->m_cols);

    for (y = 0; y != matrix1->m_rows; ++y)
    {
        for (x = 0; x != matrix2->m_cols; ++x)
        {
            sum = 0.0;

            for (i = 0; i != matrix1->m_cols; ++i)
            {
                sum += matrix_t_get(matrix1, i, y) *
                    matrix_t_get(matrix2, x, i);
            }

            matrix_t_set(result, x, y, sum);
        }
    }

    return result;
}

static const size_t MINIMUM_THREAD_LOAD = 1; /*10 * 1000;*/

typedef struct thread_info {
    matrix_t* left_matrix;
    matrix_t* right_matrix;
    matrix_t* result_matrix;
    size_t    start_row;
    size_t    rows;
} thread_info;

#ifdef _WIN32
DWORD WINAPI ThreadProc(LPVOID arg)
#else
static void* thread_func(void* arg)
#endif
{
    size_t i;
    size_t x;
    size_t y;
    double sum;
    thread_info* info = (thread_info*) arg;

    for (y = info->start_row; y < info->start_row + info->rows; ++y)
    {
        for (x = 0; x < info->right_matrix->m_cols; ++x)
        {
            sum = 0.0;

            for (i = 0; i < info->left_matrix->m_cols; ++i)
            {
                sum += matrix_t_get(info->left_matrix, i, y) *
                    matrix_t_get(info->right_matrix, x, i);
            }

            matrix_t_set(info->result_matrix, x, y, sum);
        }
    }

    return NULL;
}

matrix_t* matrix_t_multiply_parallel(matrix_t* left_matrix,
    matrix_t* right_matrix)
{
    size_t i;
    size_t rows_reserved;
    size_t work_load;
    size_t num_threads;
    size_t basic_rows_per_thread;
    thread_info* thread_info_structs;
    matrix_t* result_matrix;

    work_load = left_matrix->m_rows *
        right_matrix->m_cols *
        right_matrix->m_rows;

    num_threads = get_number_of_processors();
    num_threads = min_size_t(num_threads, work_load / MINIMUM_THREAD_LOAD);
    num_threads = min_size_t(num_threads, left_matrix->m_rows);
    num_threads = max_size_t(num_threads, 1);

    if (num_threads == 1)
    {
        return matrix_t_multiply(left_matrix, right_matrix);
    }

    result_matrix = matrix_t_alloc(left_matrix->m_rows, right_matrix->m_cols);

    if (!result_matrix)
    {
        return NULL;
    }

    basic_rows_per_thread = left_matrix->m_rows / num_threads;
    thread_info_structs = calloc(num_threads, sizeof(thread_info));

    if (!thread_info_structs)
    {
        matrix_t_free(result_matrix);
        return NULL;
    }

    for (i = 0, rows_reserved = 0;
        i != num_threads;
        i++, rows_reserved += basic_rows_per_thread)
    {
        thread_info_structs[i].left_matrix = left_matrix;
        thread_info_structs[i].right_matrix = right_matrix;
        thread_info_structs[i].result_matrix = result_matrix;
        thread_info_structs[i].rows = basic_rows_per_thread;
        thread_info_structs[i].start_row = rows_reserved;
    }

    /* Last thread function will be run in current thread: */
    thread_info_structs[num_threads - 1].rows +=
        left_matrix->m_rows % basic_rows_per_thread;

#ifdef _WIN32
    HANDLE* threads = calloc(num_threads, sizeof(HANDLE));
#else
    pthread_t* threads = calloc(num_threads, sizeof(pthread_t));
#endif

    for (i = 0; i < num_threads - 1; ++i)
    {
#ifdef _WIN32
        threads[i] = CreateThread(NULL,
                                  1000000,
                                  ThreadProc,
                                  (void*) &thread_info_structs[i],
                                  0, 
                                  NULL);
#else
        pthread_create(&threads[i],
            NULL,
            thread_func,
            (void*)&thread_info_structs[i]);
#endif
    }

#ifdef _WIN32
    ThreadProc((void*) &thread_info_structs[num_threads - 1]);
#else
    thread_func((void*)&thread_info_structs[num_threads - 1]);
#endif

#ifdef _WIN32
    WaitForMultipleObjects(num_threads - 1, threads, TRUE, INFINITE);
#else 
    for (i = 0; i < num_threads - 1; ++i)
    {
        pthread_join(threads[i], NULL);
    }
#endif
    return result_matrix;
}

void matrix_t_print(matrix_t* matrix)
{
    for (size_t y = 0; y < matrix->m_rows; ++y)
    {
        for (size_t x = 0; x < matrix->m_cols; ++x)
        {
            printf("%f ", matrix_t_get(matrix, x, y));
        }

        puts("");
    }
}

main.c

#include "matrix.h"
#include "matrix_algorithm.h"
#include <stdio.h>
#include <stdlib.h>

#ifndef _WIN32
#include <sys/time.h>
#include <unistd.h>
#else
#include <windows.h>
#endif

static matrix_t* create_random_matrix(const size_t rows, const size_t cols)
{
    size_t x;
    size_t y;
    matrix_t* m = matrix_t_alloc(rows, cols);

    if (!m)
    {
        return NULL;
    }

    for (x = 0; x < cols; ++x)
    {
        for (y = 0; y < rows; ++y)
        {
            matrix_t_set(m, x, y, ((double)rand()) / RAND_MAX);
        }
    }

    return m;
}

static size_t get_milliseconds()
{
#ifdef _WIN32
    return (size_t) GetTickCount64();
#else
    struct timeval time;
    gettimeofday(&time, NULL);
    return time.tv_sec * 1000 + time.tv_usec / 1000;
#endif
}

static int matrix_equals(const matrix_t const* a, const matrix_t const* b)
{
    size_t x;
    size_t y;

    if (a->m_cols != b->m_cols || a->m_rows != b->m_rows)
    {
        return 0;
    }

    for (y = 0; y < a->m_rows; ++y)
    {
        for (x = 0; x < a->m_cols; ++x)
        {
            if (matrix_t_get(a, x, y) != matrix_t_get(b, x, y))
            {
                return 0;
            }
        }
    }

    return 1;
}

int main() {
    size_t start_time;
    size_t end_time;
    matrix_t* a;
    matrix_t* b;
    matrix_t* ab1;
    matrix_t* ab2;

    srand((unsigned int)time(NULL));

    a = create_random_matrix(500, 500);
    b = matrix_t_copy(a);

    start_time = get_milliseconds();
    ab1 = matrix_t_multiply(a, b);
    end_time = get_milliseconds();

    printf("Single-threaded multiplication in %zu milliseconds.\n",
        end_time - start_time);

    start_time = get_milliseconds();
    ab2 = matrix_t_multiply_parallel(a, b);
    end_time = get_milliseconds();

    printf("%zu-threaded multiplication in %zu milliseconds.\n",
        get_number_of_processors(),
        end_time - start_time);

    printf("Algorithms agree: %d\n", matrix_equals(ab1, ab2));

    matrix_t_free(a);
    matrix_t_free(b);
    matrix_t_free(ab1);
    matrix_t_free(ab2);

#ifdef _WIN32
    Sleep(5000);
#endif

    return 0;
}

The above compiles in Visual C++ 2017 and Xcode 8. Please criticize my code.

Answer 1

speedup

If this code still runs in half the single-thread time, with four cores, then we at least have an aspirational goal to get to 3x or even close to 4x. But stopping short of actual speedup, we should first at least document the source of the overhead. I don't see comments discussing hyperthread, cache effects, or similar.

If there is a particular BLAS routine this is competing with, it would be helpful to name it and to say that we benchmark at e.g. 60% of its speed, and then describe our advantages, such as clarity, portability, don't need SSE, whatever.

public API

I'm thinking about how easy it is for callers to use your API correctly. I'm not crazy about the return here:

matrix_t* m = malloc(sizeof *m);
if (!m)
{
    return m;
}

It is correct, of course, but now you're burdening the caller to check for a rare event. For any malloc fail, my API preference would be something fatal like throw std::bad_alloc. It would be a global decision - you consistently return NULL, and I'm advocating for a consistent throw instead.

Breaking apart alloc & init is an interesting API design choice, and possibly a good one, possibly there's some motivating use case I'm not seeing. (And this will sound weird coming from me, since usually I'm advocating in the other direction, breaking up large functions by extracting a helper.) Consider simplifying the API by removing init from it. If we view malloc fail as fatal, we probably don't need:

    matrix->m_rows = 0;
    matrix->m_cols = 0;

Or maybe we have some minor concern about debugability in case there are dangling pointers, and so calloc() of that very small struct would suffice.

In a similar vein, looking at free & clear, maybe we don't want to publicly expose clear. Also, having free'd something, maybe you then want to explicitly write a zero into the pointer?

I'm looking at set & get. You did not publicly expose data_index, so that's good. Are we confident data_index will get inlined? Would we like to offer the compiler a hint? Would preprocessor macro expansion possibly be more appropriate?

Your triply nested loop does a lot of data_index integer multiplies. Perhaps you'd like to bench a competing approach that turns "multiply by m_cols" into "pointer addition"?

fatal malloc fail

Back to this topic. I notice that copy will segfault when assigning rows if 1st malloc failed, and will give memcpy a NULL destination if 2nd malloc failed.

cores

get_number_of_processors() is very nice. The caller cannot choose to use a subset of available cores, but you've made a decision to simplify the public API and that seems fair to me. You might possibly want to add benchmark() to the public API, and have it privately test with 1 up to max number of cores, for evaluating speedup. I assume more cores always translates to lower elapsed time, that cache thrashing never gets ugly. Embedded or mobile (battery constrained) apps might value "energy consumed" over "elapsed time", and so might want to tune to a subset of cores, but now I'm wandering off into the weeds with creeping featuritis - better to wait for the use case to arise. On the topic of feature requests: I wouldn't mind seeing you post a .PNG graph of speedup vs. num_threads.

Your signature is perfectly nice, but perhaps slightly verbose. I would have been happy with just A, B:

matrix_t* matrix_t_multiply(matrix_t* a, matrix_t* b)

The first two return NULL error checks impress me as fatal errors that are worth a throw. Checking if A cols == B rows perhaps would throw a different error.

In

result = malloc(sizeof(*result));
matrix_t_init(result, matrix1->m_rows, matrix2->m_cols);

you ignore a NULL malloc result, and also you don't check whether init succeeded. This is why I've been advocating for throw - it results in an API that is easier to use correctly. In this case the exception would simply keep bubbling up the stack. If init fails, the Y loop body never executes (zero rows), and the caller receives a hard-to-interpret result. I predict caller would typically segfault. An API that requires caller to be very very careful seems like a complex API to me.

row major order

I'm looking at the for (i = 0; i != matrix1->m_cols; ++i) inner loop. Here's a crazy idea, I'll just throw it out there. There's 3 nested loops. We would prefer to be prefetch- and cache-friendly, to advance through sequential locations, which clearly doesn't happen for matrix2. Would it make any sense to run 2 nested loops beforehand, which create a temporary copy of matrix2 rotated 90-degrees, so it's in column major order? Might make for an interesting benchmark. If this is a win, then maybe your datastructure winds up being more complex, featuring a flag that specifies the matrix is stored row- or column-major at the moment. One can do a fancy swap approach that rotates without extra storage, but I'd just go for the naive approach using malloc.

DRY

Your matrix_t_multiply() is simple and good. It does allocation with error checking and then is mostly about the 3 nested loops. Your thread_func has sort of a vague name, and it is all about those 3 nested loops. I suggest you could extract the 3 loops into a private helper function which both of them call.

As a separate matter, maybe the multiply API you offer to a caller is a good one, where caller understands your layer shall allocate, and caller shall be responsible for freeing. Another perspective on API design asks that a library provide separate "allocate" and "compute" functions, to clarify object lifetime and caller's responsibilities. I'm not advocating for API changes, I'm just inviting you to make design choices explicit. (Personally I find your current multiply API to be convenient and good.)

thread_func should probably throw if we fail a NULL check, or a NULL check, or if cols != rows, or if result is of the wrong shape. As it stands I imagine we segfault if there is such an error.

autotune

Consider moving this:

static const size_t MINIMUM_THREAD_LOAD = 1; /*10 * 1000;*/

down to just above the multiply parallel function. You probably want the 1e4 value. Or, a more principled way to autotune it would be to bench the cost of thread creation overhead on the target system, and bench 1e4 or 1e5 memory reads, and set minimum thread load based on those two measurements.

API

This is a perfectly lovely signature:

matrix_t* matrix_t_multiply_parallel(matrix_t* left_matrix, matrix_t* right_matrix)

But it's not consistent with the matrix1 & matrix2 you used before. I don't really care how you spell it (I suggested A & B), I'm just encouraging consistency in the API.

Usual grumble, when looking at if (!result_matrix) and if (!thread_info_structs): consider a throw. Of course one would want to be consistent throughout, it's an all-or-nothing API design decision.

Naming it rows_reserved seemed a little opaque to me - your start_row name seems a perfectly good one.

thread overhead, and even workload

This logic appears after the loop, and it is correct:

thread_info_structs[num_threads - 1].rows +=
    left_matrix->m_rows % basic_rows_per_thread;

However, it can give an interesting amount of extra work to that last thread. Putting the round-off logic within the loop could arrange for all threads to be responsible for no more than a single extra row relative to siblings.

It looks like you're allocating a megabyte for each thread's stack. Seems pretty arbitrary. Maybe we don't care, but maybe there's a more principled way to do it? Might become more interesting for high core counts.

The /* Last thread function will be run in current thread: */ comment is helpful for understanding the N-1 creates and joins. Just before #ifdef ... thread_func would you please add a repeat comment, explaining that this is where we do the work. I think you believe thread creation consumes an interesting number of cycles, so complexifying the code in this way is a useful optimization. Ok, fine, I accept that you may be right. It would be helpful to justify that design decision by adding a comment describing some benchmark result, showing a cost of 100 usec or 100 msec or whatever. If the cost is "small" I would be inclined to go with N creates and joins and no special case.

Partitioning by sequential rows makes perfect sense, so each thread has a row range. Rows are not perfect multiples of cache line size. Let me float an alternate proposal. What if we partitioned by a row stride (of num_threads) instead? Then each thread i starts on row i, then each starts on row i + row_stride, and so on. There are two effects I'm going for. One is a tendency for a thread that suffers an L2 miss to wait and therefore synchronize with sibling threads. The other is allowing end-of-row (or column) cached values in a cache line that was dragged in by one thread to soon be used by another thread, so our cache friendly code enjoys more cache hits. Are you interested enough to bench it and compare?

Consider adding matrix_equals() to your public API. Consider replacing the two nested loops with a single call to memcmp(). This is similar to the use in copy of memcpy().

I'm looking at

srand((unsigned int)time(NULL));

I feel you should print out the seed value, so that in the unlikely event that a benchmark run reports a mismatch, you could reproduce the run. So I am seconding Toby Speight's remark. I also agree with him that multiplying square matrices (500 x 500) is too easy of a test, that, say, 300 x 800 would be more rigorous. Additionally, multiplying skinny matrices, say 2 x 120000 or 4 x 60000, may reveal that for performance you want two different code variants selected according to whether rows or cols is bigger.

I already touched on benching with num_threads < num_cores. It might be interesting to benchmark with extra threads, from 1 up to say 2 * num_cores. We anticipate worse performance due to extra thread overhead, and the point of the experiment would be to measure the magnitude of the overhead, both for scheduling effects and cache effects. It is interesting because the caller may already have threads running, again leading to a non-ideal number of concurrent threads even when you use exactly num_cores threads. Also, you never voluntarily relinquish, which is fine, but if the caller has a latency-sensitive multimedia thread in the background, possibly he'll be annoyed by high dispatch latency. If benching reveals acceptably low overhead, then deliberately creating very large number of tiny tasks (in the limit, single-row tasks) could be a useful way to give the scheduler more opportunities to run a high-priority audio playback thread.

Ok, that's it. Overall it looks like solid code, quite clear, that does what it claims to do.

Answer 2

Accessors not inlined

I ran your program and found it that it was quite slow. After some investigation, I determined that one of the biggest factors causing the slow speed was the fact that your accessor functions matrix_t_get() and matrix_t_set() were not inlined.

Inlining in this case helps in two ways:

1. Eliminates function call overhead for an operation that is performed \$O(n^3)\$ times.
2. Provides the opportunity for the compiler to make additional optimizations to the loop.

I moved matrix_t_get(), matrix_t_set(), and data_index() to matrix.h and changed them to static inline functions. Here are the results:

Before inlining (500x500)
-------------------------
Single-threaded multiplication in 767 milliseconds.
4-threaded multiplication in 211 milliseconds.

After inlining (500x500)
------------------------
Single-threaded multiplication in 356 milliseconds.
4-threaded multiplication in 100 milliseconds.

Before inlining (1000x1000)
---------------------------
Single-threaded multiplication in 9221 milliseconds.
4-threaded multiplication in 2389 milliseconds.

After inlining (1000x1000)
--------------------------
Single-threaded multiplication in 7031 milliseconds.
4-threaded multiplication in 1900 milliseconds.

As you can see, the 500x500 case ran twice as fast using the inlined functions. However, the 1000x1000 case did not improve as much.

Cache friendly access order

To further improve the speed of your matrix multiplication, you should consider the access order of your elements. Accessing items sequentially in the same row is much better for the cache than accessing items sequentially in the same column. There are plenty of articles online which can show you how to rearrange your loops and multiplications to accomplish this. I went ahead and rewrote your function like this:

matrix_t* matrix_t_multiply(matrix_t* matrix1, matrix_t* matrix2)
{
    matrix_t* result;
    double sum;

    if (!matrix1 || !matrix2)
    {
        return NULL;
    }

    if (!matrix1->m_data || !matrix2->m_data)
    {
        return NULL;
    }

    if (matrix1->m_cols != matrix2->m_rows)
    {
        return NULL;
    }

    result = malloc(sizeof(*result));
    matrix_t_init(result, matrix1->m_rows, matrix2->m_cols);

    for (size_t i = 0; i < matrix1->m_rows; i++)
    {
        for (int k = 0; k < matrix1->m_cols; k++)
        {
            double t = matrix_t_get(matrix1, k, i);
            for (int j = 0; j < matrix1->m_cols; j++ )
            {
                double res = matrix_t_get(result, j, i);
                res += t * matrix_t_get(matrix2, j, k);
                matrix_t_set(result, j, i, res);
            }
        }
    }

    return result;
}

And here are the results (note: all times are with inlined functions. I did not modify the threaded version):

Before reorder (500x500)
------------------------
Single-threaded multiplication in 356 milliseconds.

After reorder (500x500)
-----------------------
Single-threaded multiplication in 106 milliseconds.

Before reorder (1000x1000)
--------------------------
Single-threaded multiplication in 7031 milliseconds.

After reorder (1000x1000)
-------------------------
Single-threaded multiplication in 972 milliseconds.

The 500x500 case got 3 times faster, and the 1000x1000 case got 7 times faster. Overall, including both inlining and reordering, the single threaded case got about 8-9x faster, which made the final single threaded case even faster than the original multi threaded case.

Answer 3

You are leaking the memory for the "thread_info_struct" and "threads" array.