I'm a student, and I'm trying to make the product of two square matrix with some threads in the soup. I've already worked on some fast single-threaded product (with some cache-friendly tricks), and I'm interested now on multithreading.The following code works, but I don't fully understand my outputs in terms of performance.I'd like to improve my code a bit and of course, improve my implementation (and in a perfect world, achieving a multithreaded version faster than a naive singlethreaded version).Makefile:CC=gccCCO=# -O3 -march=nativeCCF=-std=gnu11 -Wall -Wextra -pedantic -Wno-unusedLDF=-pthreadEXTRA=-DNB_TH=8all: modular B_bloc K_bloc benchmodular:    $(CC) -Dmodular $(EXTRA) $(LDF) $(CCF) $(CCO) ./mat_mult.c -o modularB_bloc:    $(CC) -DB_bloc $(EXTRA) $(LDF) $(CCF) $(CCO) ./mat_mult.c -o B_blocK_bloc:    $(CC) -DK_bloc $(EXTRA) $(LDF) $(CCF) $(CCO) ./mat_mult.c -o K_blocmr_proper:    rm -f modular B_bloc K_blocbench:    ./modular 128    ./B_bloc 128    ./K_bloc 128    ./modular 208    ./B_bloc 208    ./K_bloc 208    ./modular 400    ./B_bloc 400    ./K_bloc 400    ./modular 1024    ./B_bloc 1024    ./K_bloc 1024mat_mult.c:#include <stdio.h>#include <assert.h>#include <stdlib.h>#include <unistd.h>#include <pthread.h>// Pick your favorite flavor// #define modular// #define B_bloc// #define K_bloc#ifdef modular    #define FOO_TH modular_    #define MALLOC_P_TH default_alloc_#elif defined(B_bloc)    #define FOO_TH B_bloc_    #define MALLOC_P_TH default_alloc_#elif defined(K_bloc)    #define FOO_TH K_bloc_    #define MALLOC_P_TH K_bloc_alloc_    #ifndef BLK        #define BLK 2    #endif#endif#ifndef BLK    #define BLK 16#endif#ifndef NB_TH    #define NB_TH 8#endif#define THREAD_WAIT 1#define THREAD_OVER 0// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// TOOLS// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// Tool to build a matrixvoid build_mat_(int*** target, int width, int height){    void* origin = malloc(sizeof(int*) * width + sizeof(int) * width * height);    int** header = (int**) origin;    int* cursor = (int*) (header + width);    for(int i = 0; i < width; ++i)    {        header[i] = cursor;        cursor += width;    }    *target = (int**) origin;}// Tool to allocate a struct-like void* thread parameter (default)void* default_alloc_(){    return malloc(0            + sizeof(int)   * 3            + sizeof(int*)  * 1            + sizeof(int**) * 3            );}// Tool to allocate a struct-like void* thread parameter (K_bloc only)void* K_bloc_alloc_(){    return malloc(0            + sizeof(int)               * 3            + sizeof(int*)              * 1            + sizeof(int**)             * 3            + sizeof(pthread_mutex_t*)  * 1            );}// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// THREADS FOO_TH// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// Split the mult' by picking one-out-of NB_TH column.void* modular_(void* data){    // Data retrieving    int* int_p = (int*) data;    int width = int_p[0];    int height = int_p[1];    int id_th = int_p[2];    int** vec_p = (int**) (int_p + 3);    int* thread_status = vec_p[0];    int*** mat_p = (int***) (vec_p + 1);    int** A = mat_p[0];    int** B = mat_p[1];    int** C = mat_p[2];    // Iterating on columns    for(int i = 0; i < width; ++i)        if(i % NB_TH == id_th)            for(int j = 0; j < height; ++j)                for(int k = 0; k < width; ++k)                    C[i][j] += A[i][k] * B[k][j];#ifdef UNJOIN    for(;thread_status[-1];sleep(0))        thread_status[id_th] = THREAD_OVER;#else    thread_status[id_th] = THREAD_OVER;#endif    pthread_exit(NULL);}// Split the mult' in many blocs on j axis. Cache-friendly on B.void* B_bloc_(void* data){    // Data retrieving    int* int_p = (int*) data;    int width = int_p[0];    int height = int_p[1];    int id_th = int_p[2];    int** vec_p = (int**) (int_p + 3);    int* thread_status = vec_p[0];    int*** mat_p = (int***) (vec_p + 1);    int** A = mat_p[0];    int** B = mat_p[1];    int** C = mat_p[2];    // Calculating blocs size    int blk_width = BLK;    int cursor = id_th * blk_width;    // Iterating on blocs    while(cursor < width)    {        for(int i = 0; i < width; ++i)            for(int j = cursor; j < cursor + blk_width; ++j)                for(int k = 0; k < width; ++k)                    C[i][j] += A[i][k] * B[k][j];        cursor += NB_TH * blk_width;    }    // Various stuff and exit.#ifdef UNJOIN    for(;thread_status[-1];sleep(0))        thread_status[id_th] = THREAD_OVER;#else    thread_status[id_th] = THREAD_OVER;#endif    pthread_exit(NULL);}// Split the mult' in blocs on i, j, and k axis. Ensure local data.void* K_bloc_(void* data){    // Data retrieving    int* int_p = (int*) data;    int width = int_p[0];    int height = int_p[1];    int id_th = int_p[2];    int** vec_p = (int**) (int_p + 3);    int* thread_status = vec_p[0];    int*** mat_p = (int***) (vec_p + 1);    int** A = mat_p[0];    int** B = mat_p[1];    int** C = mat_p[2];    pthread_mutex_t** mut_p = (pthread_mutex_t**) (mat_p + 3);    pthread_mutex_t* mutex = mut_p[0];    // Calculating bloc size    int blk_width = width / BLK;    int blk_i = (id_th % BLK) * blk_width;    int blk_j = ((id_th / BLK) % BLK) * blk_width;    int blk_k = ((id_th / (BLK * BLK)) % BLK) * blk_width;    int blk_id = id_th % (BLK * BLK);    // Creating local matrix for saving local results    int** C_local;    build_mat_(&C_local, blk_width, blk_width);    // Initializing matrix to 0    for(int i = 0; i < blk_width; ++i)        for(int j = 0; j < blk_width; ++j)            C_local[i][j] = 0;    // Iterating on the bloc    for(int i = blk_i, i_l = 0; i < blk_i + blk_width; ++i, ++i_l)        for(int j = blk_j, j_l = 0; j < blk_j + blk_width; ++j, ++j_l)            for(int k = blk_k; k < blk_k + blk_width; ++k)                C_local[i_l][j_l] += A[i][k] * B[k][j];    // Locking the bloc on targeted matrix to write local results    pthread_mutex_lock(mutex + blk_id);    for(int i = blk_i, i_l = 0; i < blk_i + blk_width; ++i, ++i_l)        for(int j = blk_j, j_l = 0; j < blk_j + blk_width; ++j, ++j_l)            C[i][j] += C_local[i_l][j_l];    pthread_mutex_unlock(mutex + blk_id);    // Free memory    free(C_local);    // Various stuff and exit.#ifdef UNJOIN    for(;thread_status[-1];sleep(0))        thread_status[id_th] = THREAD_OVER;#else    thread_status[id_th] = THREAD_OVER;#endif    pthread_exit(NULL);}// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// MAIN// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /int main(int argc, char** argv){#ifdef TM    int width = TM;    int height = width;#else    int width = atoi(argv[1]);    int height = width;#endif// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// CREATINGMATRIX// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /    // Contiguous columns to help the cache    int **A, **B, **C, **D;    build_mat_(&A, width, height);    build_mat_(&B, width, height);    build_mat_(&C, width, height);    build_mat_(&D, width, height);    // Initialize them to various values    for(int i = 0; i < width; ++i)        for(int j = 0; j < height; ++j)        {            A[i][j] = (i * i) % (j + 1) + i;            B[i][j] = (i * j) % (i + 1) + j;            C[i][j] = 0;            D[i][j] = 0;        }// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// CALCULATEREFERENCE RESULT (NAIVE)// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /    // Naive calculus of the resulting matrix    clock_t naive_start = clock();    for(int i = 0; i < width; ++i)        for(int j = 0; j < height; ++j)            for(int k = 0; k < width; ++k)                D[i][j] += A[i][k] * B[k][j];    clock_t naive_end = clock();// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// PREPARING THREADS// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /    pthread_t thread_id[NB_TH];    void* thread_p[NB_TH];    // Thread status vector, mostly used with -DUNJOIN    int* thread_status = (int*) malloc(sizeof(int) * (NB_TH + 1));    for(int i = 0; i < (NB_TH + 1); ++i)        thread_status[i] = THREAD_WAIT;    ++thread_status;    // Checking asserts and execute method-related initializations.    assert(width == height);#ifdef B_bloc    assert(width % BLK == 0);#elif defined(K_bloc)    assert(NB_TH == BLK * BLK * BLK);    assert(width % BLK == 0);    pthread_mutex_t* K_bloc_mutex =        malloc(sizeof(pthread_mutex_t) * (NB_TH / BLK));    for(int i = 0; i < NB_TH / BLK; ++i)        pthread_mutex_init(K_bloc_mutex + i, NULL);#endif    // Initialize thread's parameters.    for(int i = 0; i < NB_TH; ++i)    {        thread_p[i] = MALLOC_P_TH();        int* int_p = (int*) thread_p[i];        int_p[0] = width;        int_p[1] = height;        int_p[2] = i;        int** vec_p = (int**) (int_p + 3);        vec_p[0] = thread_status;        int*** mat_p = (int***) (vec_p + 1);        mat_p[0] = A;        mat_p[1] = B;        mat_p[2] = C;#ifdef K_bloc        pthread_mutex_t** mut_p = (pthread_mutex_t**) (mat_p + 3);        mut_p[0] = K_bloc_mutex;#endif    }// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// STARTING THREADS// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /    // Start the thread's timer    clock_t thread_start = clock();    for(int i = 0; i < NB_TH; ++i)        pthread_create(&thread_id[i], NULL, FOO_TH, thread_p[i]);// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// WAITINGTHREADS// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /#ifndef UNJOIN    for(int i = 0; i < NB_TH; ++i)        pthread_join(thread_id[i], NULL);#else    for(int i, status = THREAD_WAIT; status == THREAD_WAIT;)        for(i = 0, status = THREAD_OVER; i < NB_TH; ++i, sleep(0))            status = status || thread_status[i];#endif    // End the thread's timer    clock_t thread_end = clock();    for(int i = 0; i < NB_TH; ++i)        assert(thread_status[i] == THREAD_OVER);    thread_status[-1] = THREAD_OVER;// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// CHECKING THERESULT// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /    for(int i = 0; i < width; ++i)        for(int j = 0; j < height; ++j)            if(C[i][j] != D[i][j])                return 1;// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /// FREEING MEMORY// / // / // / // / // / // / // / // / // / // / // / // / // / // / // / // /#ifdef K_bloc    free(K_bloc_mutex);#endif    free(--thread_status);    for(int i = 0; i < NB_TH; ++i)        free(thread_p[i]);    free(A);    free(B);    free(C);    free(D);    float naive_time = ((float) naive_end - naive_start) / CLOCKS_PER_SEC;    float thread_time = ((float) thread_end - thread_start) / CLOCKS_PER_SEC;    printf(%f\n, naive_time);    printf(%f\n, thread_time);    printf(%f\n, thread_time / naive_time);    return 0;}Sample output:gcc -Dmodular -DNB_TH=8 -pthread -std=gnu11 -Wall -Wextra -pedantic -Wno-unused  ./mat_mult.c -o modulargcc -DB_bloc -DNB_TH=8 -pthread -std=gnu11 -Wall -Wextra -pedantic -Wno-unused  ./mat_mult.c -o B_blocgcc -DK_bloc -DNB_TH=8 -pthread -std=gnu11 -Wall -Wextra -pedantic -Wno-unused  ./mat_mult.c -o K_bloc./modular 1280.0099410.0091690.922342./B_bloc 1280.0100690.0084710.841295./K_bloc 1280.0097710.0103461.058848./modular 2080.0386600.0822872.128479./B_bloc 2080.0381760.0483731.267105./K_bloc 2080.0379380.0518021.365438./modular 4000.3297820.5023671.523331./B_bloc 4000.3271450.5018081.533901./K_bloc 4000.3249830.6389291.966038./modular 10249.32289116.2631611.744433./B_bloc 10249.13641416.8749431.846999./K_bloc 10248.65123315.3782311.777577
Multiplying square matrix with C pthreads (POSIX threads)
performance;c;multithreading;matrix;pthreads
Wrong loop orderIn your naive matrix multiplication, you have your loops in a suboptimal ordering.  I compared your implementation with this one, where I swapped the second and third loops:for(int i = 0; i < width; ++i)    for(int k = 0; k < width; ++k)        for(int j = 0; j < height; ++j)            D[i][j] += A[i][k] * B[k][j];When I ran B_bloc 1024, I got these numbers (no threads involved):Original loops: 9.89 secSwapped loops : 1.07 secSo clearly you need to think about cache usage a little bit more.  The version I listed here is more cache friendly because A[i][k] is constant in the inner loop while the other two matrices are being iterated in memory order.  With the original loops where the inner loop iterated on k, the B[k][j] term was jumping cache lines on every loop iteration.Use a struct!The code you use to pass information between the main thread and the pthreads is not good at all.  You are reinventing the wheel by constructing your own handcoded struct:    // This is what MALLOC_P_TH() calls:    void* default_alloc_()    {        return malloc(0                + sizeof(int)   * 3                + sizeof(int*)  * 1                + sizeof(int**) * 3                );    }    // In main():    thread_p[i] = MALLOC_P_TH();    int* int_p = (int*) thread_p[i];    int_p[0] = width;    int_p[1] = height;    int_p[2] = i;    int** vec_p = (int**) (int_p + 3);    vec_p[0] = thread_status;    int*** mat_p = (int***) (vec_p + 1);    mat_p[0] = A;    mat_p[1] = B;    mat_p[2] = C;You should be using an actual struct, like this:typedef struct ThreadArgs {    int width;    int height;    int threadNum;    int *pStatus;    int **matrixA;    int **matrixB;    int **matrixC;} ThreadArgs;// In main():thread_p[i]            = malloc(sizeof(ThreadArgs));thread_p[i]->width     = width;thread_p[i]->height    = height;thread_p[i]->threadNum = i;thread_p[i]->pStatus   = thread_status;thread_p[i]->matrixA   = A;thread_p[i]->matrixB   = B;thread_p[i]->matrixC   = C;As you can see, this improves the code in many ways:You don't have to have a custom function just to compute the size of the struct.Your fields are actually named, so you can refer to each field by name rather than number.If you need to add or remove fields, you simply make modifications to the struct rather than fiddling with the handcoded arrays.Using the wrong clock functionWhen I ran your program, I was surprised to see that none of your multithreaded versions ran faster than the single threaded version.  However, I discovered that you were using the wrong function to time your program.  The clock() function measures cpu cycles used across all threads.  So even if your program ran 8 times faster, you would still get the same reading from clock() as for a single threaded program.  What you want is a wall clock measurement.  I would suggest using something like clock_gettime() or gettimeofday().  After I modified your program to use clock_gettime() and fixed the loop orders, I ran B_bloc and got the following output, which makes a lot more sense (I have a 4 core system, so 0.25 on the third line is expected):$ B_bloc 10241.1029120.2767520.250928