Pascal's triangle with libgmp and pthread in C

Question

I have written a program to generate a pascal triangle of predefined size.

Code

#include <stdio.h>
#include <assert.h>

#include <pthread.h>
#include <gmp.h>

#define N 50
#define SIZE(x) (((x) >> 1) + 1)
#define NTHREADS 8

int min(int a, int b) {
    return (a < b) ? a : b;
}

mpz_t row[2][SIZE(N)];


typedef struct _tls_t {
    int
        id,
        begin,
        end;
} tls_t;

void *calc_elem(void *args) {
    tls_t *local = (tls_t *)(args);
    for(int j = local->begin; j < local->end; ++j) {
        mpz_add(row[1][j], row[0][j - 1], row[0][j]);
    }
    return NULL;
}


void calc_row(pthread_t threads[NTHREADS], tls_t thrtls[NTHREADS], int i) {
    if((i & 1) == 0)
        mpz_set(row[0][SIZE(i) - 1], row[0][SIZE(i) - 2]);

    const int t_step = SIZE(i) / NTHREADS + 1;
    for(int t = 0; t < min(NTHREADS, SIZE(i) / t_step + 1); ++t) {
        thrtls[t].id = t;
        thrtls[t].begin = t_step * t + 1;
        thrtls[t].end = min(t_step * (t + 1) + 1, SIZE(i));
        assert(!pthread_create(threads + t, NULL, calc_elem, (void *)(thrtls + t)));
        assert(!pthread_join(threads[t], NULL));
    }
}


void init(), clear(), print_row(int), mirror_rows(int);
void pascal_triangle() {
    pthread_t threads[NTHREADS];
    tls_t thrtls[NTHREADS];
    init();
    for(int i = 0; i < N; ++i) {
        calc_row(threads, thrtls, i);
        print_row(i);
        mirror_rows(i);
    }
    clear();
}


void init() {
    for(int i = 0; i < SIZE(N); ++i) {
        mpz_init(row[0][i]);
        mpz_init(row[1][i]);
        mpz_set_si(row[0][i], 0);
    }
    mpz_set_si(row[0][0], 1);
    mpz_set_si(row[1][0], 1);
}

void print_row(int rsize) {
    int idx = -1;
    for(int i = 0; i <= rsize; ++i) {
        if(i < SIZE(rsize))
            ++idx;
        else if(i > SIZE(rsize) || !(rsize & 1))
            --idx;

        assert(idx >= 0 && idx < SIZE(rsize));
        gmp_printf("%Zd ", row[1][idx]);
    }
    gmp_printf("\n");
}

void mirror_rows(int i) {
    for(int j = 0; j < SIZE(i); ++j) {
        mpz_set(row[0][j], row[1][j]);
    }
}

void clear() {
    for(int i = 0; i < SIZE(N); ++i) {
        mpz_clear(row[0][i]);
        mpz_clear(row[1][i]);
    }
    pthread_exit(NULL);
}


main() {
    pascal_triangle();
}

The program only has to remember two arrays of mpz_t, each of size (N >> 1) + 1 (left half) being respectively, previous and current row. The i-th row is calculated from (i-1)-th.

Although with this strategy the algorithm is not very efficient to parallel, it should consume \$\mathcal{O}(N)\$ space and \$\mathcal{O}(N \cdot \mathcal{O}(\texttt{mpz_add}))\$ operations.

The whole parallelization block only accelerates calculations of a single row.

N=25 output

1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
1 7 21 35 35 21 7 1
1 8 28 56 70 56 28 8 1
1 9 36 84 126 126 84 36 9 1
1 10 45 120 210 252 210 120 45 10 1
1 11 55 165 330 462 462 330 165 55 11 1
1 12 66 220 495 792 924 792 495 220 66 12 1
1 13 78 286 715 1287 1716 1716 1287 715 286 78 13 1
1 14 91 364 1001 2002 3003 3432 3003 2002 1001 364 91 14 1
1 15 105 455 1365 3003 5005 6435 6435 5005 3003 1365 455 105 15 1
1 16 120 560 1820 4368 8008 11440 12870 11440 8008 4368 1820 560 120 16 1
1 17 136 680 2380 6188 12376 19448 24310 24310 19448 12376 6188 2380 680 136 17 1
1 18 153 816 3060 8568 18564 31824 43758 48620 43758 31824 18564 8568 3060 816 153 18 1
1 19 171 969 3876 11628 27132 50388 75582 92378 92378 75582 50388 27132 11628 3876 969 171 19 1
1 20 190 1140 4845 15504 38760 77520 125970 167960 184756 167960 125970 77520 38760 15504 4845 1140 190 20 1
1 21 210 1330 5985 20349 54264 116280 203490 293930 352716 352716 293930 203490 116280 54264 20349 5985 1330 210 21 1
1 22 231 1540 7315 26334 74613 170544 319770 497420 646646 705432 646646 497420 319770 170544 74613 26334 7315 1540 231 22 1
1 23 253 1771 8855 33649 100947 245157 490314 817190 1144066 1352078 1352078 1144066 817190 490314 245157 100947 33649 8855 1771 253 23 1
1 24 276 2024 10626 42504 134596 346104 735471 1307504 1961256 2496144 2704156 2496144 1961256 1307504 735471 346104 134596 42504 10626 2024 276 24 1

Could you suggest how to improve the code and/or algorithm?

Answer 1

Work not split evenly

Your computation to split the work isn't correct. It is more apparent at low row numbers, where often only the first thread does all the work and the other threads don't do anything. The problem is that you are performing a rounding up operation by just adding one. I would change this code:

const int t_step = SIZE(i) / NTHREADS + 1;
for(int t = 0; t < min(NTHREADS, SIZE(i) / t_step + 1); ++t) {

to:

const int t_step = (SIZE(i) + NTHREADS - 1) / NTHREADS;
for(int t = 0; t < (SIZE(i) + t_step - 1) / t_step; ++t) {

Notice that to round x/y up to the nearest integer I use the expression (x+y-1)/y instead of your current expression x/y + 1. This gets rid of the need for the min as well.

Additionally, you are splitting the wrong amount of work because you are splitting SIZE(i) operations but you only need to do SIZE(i)-1 operations. This is a minor problem but it occasionally results in a slightly uneven split of the work.

Answer 2

pthread_join is blocking. It means that the loop

for(int t = 0; t < min(NTHREADS, SIZE(i) / t_step + 1); ++t) {
    ....
    assert(!pthread_create(threads + t, NULL, calc_elem, (void *)(thrtls + t)));
    assert(!pthread_join(threads[t], NULL));
}

waits for a thread to terminate before starting next one. At no time threads are running in parallel. You need to start and join threads in different loops:

for(int t = 0; t < min(NTHREADS, SIZE(i) / t_step + 1); ++t) {
    ....
    assert(!pthread_create(threads + t, NULL, calc_elem, (void *)(thrtls + t)));
}

for(int t = 0; t < min(NTHREADS, SIZE(i) / t_step + 1); ++t) {
    assert(!pthread_join(threads[t], NULL));
}

---

Never put a production code inside assert. Id NDEBUG is defined, assert is optimized away, along with the predicate.