From 34b3eedbb94cba9a40f4c4aa5a26240f32d54682 Mon Sep 17 00:00:00 2001 From: Weijia Wang <9713184+wegank@users.noreply.github.com> Date: Fri, 2 Oct 2026 01:06:11 +0200 Subject: [PATCH] Replace __sync builtins with C11 atomics in linear algebra The parallel linear algebra publishes new pivot rows with __sync_bool_compare_and_swap, a GCC builtin that is not available with MSVC. Moreover, the other threads read these pivot arrays with plain loads, which is a data race under the C11 memory model. The pivot arrays that threads add pivots to concurrently are now arrays of _Atomic pointers. A new pivot is published with a compare and swap using release ordering, the reducers load entries with acquire ordering, so a published row and its coefficients are fully visible to them. Sequential code after the parallel regions uses relaxed accesses. Arrays only read during parallel regions stay plain, and the compare and swap in the sequential kernel computation of the saturation is replaced by a plain check and store. The prototypes of the internal 32 bit reducers taking atomic arguments move from the installed data.h to data.c, so that the public headers stay usable from C++. configure now checks for C11 atomics. Co-Authored-By: Claude Opus 5.5 --- configure.ac | 22 +++ src/neogb/data.c | 38 +++- src/neogb/data.h | 46 +---- src/neogb/io.c | 8 +- src/neogb/la_ff_16.c | 397 +++++++++++++++++++------------------- src/neogb/la_ff_32.c | 441 ++++++++++++++++++++++--------------------- src/neogb/la_ff_8.c | 397 +++++++++++++++++++------------------- src/neogb/la_qq.c | 144 +++++++------- 8 files changed, 769 insertions(+), 724 deletions(-) diff --git a/configure.ac b/configure.ac index 2221b5a0..a71a72bf 100644 --- a/configure.ac +++ b/configure.ac @@ -56,6 +56,28 @@ fi # Checks for header files. AC_CHECK_HEADERS([inttypes.h stdint.h]) +# The parallel linear algebra needs C11 atomics. +AC_CACHE_CHECK([for C11 atomics], [msolve_cv_c11_atomics], + [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ +#ifdef __STDC_NO_ATOMICS__ +#error no C11 atomics +#endif +#include +#include +]], [[ +int x = 0; +_Atomic(int *) p; +int *expected = NULL; +atomic_init(&p, NULL); +return !atomic_compare_exchange_strong_explicit(&p, &expected, &x, + memory_order_release, memory_order_relaxed); +]])], + [msolve_cv_c11_atomics=yes], + [msolve_cv_c11_atomics=no])]) +if test "x$msolve_cv_c11_atomics" != xyes; then + AC_MSG_ERROR([C11 atomics () are required, try a C11 compiler or a newer -std= option.]) +fi + # Checks for typedefs, structures, and compiler characteristics. AC_C_INLINE AC_TYPE_INT16_T diff --git a/src/neogb/data.c b/src/neogb/data.c index e64622fd..7ae7983c 100644 --- a/src/neogb/data.c +++ b/src/neogb/data.c @@ -19,8 +19,38 @@ * Mohab Safey El Din */ +#include + #include "data.h" +/* Pivot rows are stored in arrays indexed by their leading column. During + * the parallel linear algebra several threads add new pivots to the same + * array, so these arrays are atomic: a thread publishes a fully reduced + * and normalized row with a compare and swap (release), and readers load + * the entries with acquire semantics, so a published row and its + * coefficients are completely visible to them. Outside of parallel regions + * relaxed accesses are sufficient. + * + * Allocates such an array with n entries, the first nrows of them are + * initialized with rows, all others with NULL. */ +static inline _Atomic(hm_t *) *allocate_atomic_pivots( + const len_t n, + hm_t * const * const rows, + const len_t nrows + ) +{ + len_t i; + _Atomic(hm_t *) *pivs = (_Atomic(hm_t *) *)malloc( + (uint64_t)n * sizeof(_Atomic(hm_t *))); + for (i = 0; i < nrows; ++i) { + atomic_init(&pivs[i], rows[i]); + } + for (; i < n; ++i) { + atomic_init(&pivs[i], NULL); + } + return pivs; +} + /* functions */ /* bs_t *initialize_basis( * const int32_t ngens @@ -126,7 +156,7 @@ hm_t *reduce_dense_row_by_known_pivots_sparse_ff_32( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t *const *pivs, + _Atomic(hm_t *) *pivs, const hi_t dpiv, const hm_t tmp_pos, const len_t mh, /* multiplier hash for tracing */ @@ -140,7 +170,7 @@ hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_32( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t *const *pivs, + _Atomic(hm_t *) *pivs, const hi_t dpiv, const hm_t tmp_pos, const len_t mh, @@ -154,7 +184,7 @@ cf32_t *reduce_dense_row_by_all_pivots_ff_32( const bs_t * const bs, len_t *pc, hm_t *const *pivs, - cf32_t *const *dpivs, + _Atomic(cf32_t *) *dpivs, const uint32_t fc ); @@ -162,7 +192,7 @@ cf32_t *reduce_dense_row_by_all_pivots_ff_32( cf32_t *reduce_dense_row_by_dense_new_pivots_ff_32( int64_t *dr, len_t *pc, - cf32_t * const * const pivs, + _Atomic(cf32_t *) * const pivs, const len_t ncr, const uint32_t fc ); diff --git a/src/neogb/data.h b/src/neogb/data.h index 4dc09bf2..2f5fa6cc 100644 --- a/src/neogb/data.h +++ b/src/neogb/data.h @@ -537,48 +537,8 @@ extern hm_t *sba_reduce_dense_row_by_known_pivots_sparse_ff_32( md_t *st ); -extern hm_t *reduce_dense_row_by_known_pivots_sparse_ff_32( - int64_t *dr, - mat_t *mat, - const bs_t * const bs, - hm_t *const *pivs, - const hi_t dpiv, - const hm_t tmp_pos, - const len_t mh, /* multiplier hash for tracing */ - const len_t bi, /* basis index of generating element */ - const len_t tr, /* trace data? */ - md_t *st - ); - -extern hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_32( - rba_t *rba, - int64_t *dr, - mat_t *mat, - const bs_t * const bs, - hm_t *const *pivs, - const hi_t dpiv, - const hm_t tmp_pos, - const len_t mh, - const len_t bi, - md_t *st - ); - -extern cf32_t *reduce_dense_row_by_all_pivots_ff_32( - int64_t *dr, - mat_t *mat, - const bs_t * const bs, - len_t *pc, - hm_t *const *pivs, - cf32_t *const *dpivs, - const uint32_t fc - ); - -extern cf32_t *reduce_dense_row_by_dense_new_pivots_ff_32( - int64_t *dr, - len_t *pc, - cf32_t * const * const pivs, - const len_t ncr, - const uint32_t fc - ); +/* The reducers working on pivot arrays shared between threads take C11 + * atomic arguments, they are declared in data.c to keep this installed + * header usable from C++. */ #endif diff --git a/src/neogb/io.c b/src/neogb/io.c index d2a83160..55706476 100644 --- a/src/neogb/io.c +++ b/src/neogb/io.c @@ -1096,7 +1096,7 @@ hm_t *reduce_dense_row_by_known_pivots_sparse_ff_32( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t *const *pivs, + _Atomic(hm_t *) *pivs, const hi_t dpiv, const hm_t tmp_pos, const len_t mh, @@ -1118,7 +1118,7 @@ hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_32( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t *const *pivs, + _Atomic(hm_t *) *pivs, const hi_t dpiv, const hm_t tmp_pos, const len_t mh, @@ -1140,7 +1140,7 @@ cf32_t *reduce_dense_row_by_all_pivots_ff_32( const bs_t * const bs, len_t *pc, hm_t *const *pivs, - cf32_t *const *dpivs, + _Atomic(cf32_t *) *dpivs, const uint32_t fc ) { @@ -1153,7 +1153,7 @@ cf32_t *reduce_dense_row_by_all_pivots_ff_32( cf32_t *reduce_dense_row_by_dense_new_pivots_ff_32( int64_t *dr, len_t *pc, - cf32_t * const * const pivs, + _Atomic(cf32_t *) * const pivs, const len_t ncr, const uint32_t fc ) diff --git a/src/neogb/la_ff_16.c b/src/neogb/la_ff_16.c index de7ba538..7350a052 100644 --- a/src/neogb/la_ff_16.c +++ b/src/neogb/la_ff_16.c @@ -98,7 +98,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_16( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -157,7 +157,9 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_16( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -166,7 +168,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_16( } /* found reducer row, get multiplier */ const uint32_t mul = (uint32_t)(fc - dr[i]); - dts = pivs[i]; if (i < ncl) { /* set corresponding bit of reducer in reducer bit array */ if (tr > 0) { @@ -417,7 +418,7 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_16( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -442,7 +443,9 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_16( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -451,7 +454,6 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_16( } /* found reducer row, get multiplier */ const uint32_t mul = (uint32_t)(fc - dr[i]); - dts = pivs[i]; if (i < ncl) { cfs = bs->cf_16[dts[COEFFS]]; /* set corresponding bit of reducer in reducer bit array */ @@ -504,7 +506,7 @@ static cf16_t *reduce_dense_row_by_all_pivots_ff_16( const bs_t * const bs, len_t *pc, hm_t *const * const pivs, - cf16_t *const * const dpivs, + _Atomic(cf16_t *) * const dpivs, const uint32_t fc ) { @@ -554,15 +556,14 @@ static cf16_t *reduce_dense_row_by_all_pivots_ff_16( if (dr[i] == 0) { continue; } - if (dpivs[i-ncl] == NULL) { + red = atomic_load_explicit(&dpivs[i-ncl], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } k++; continue; } - - red = dpivs[i-ncl]; const uint32_t mul = (uint32_t)(fc - dr[i]); const len_t os = (ncols - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { @@ -662,7 +663,7 @@ static cf16_t *reduce_dense_row_by_old_pivots_ff_16( static cf16_t *reduce_dense_row_by_dense_new_pivots_ff_16( int64_t *dr, len_t *pc, - cf16_t * const * const pivs, + _Atomic(cf16_t *) * const pivs, const len_t ncr, const uint32_t fc ) @@ -678,7 +679,8 @@ static cf16_t *reduce_dense_row_by_dense_new_pivots_ff_16( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + const cf16_t * const red = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } @@ -689,13 +691,13 @@ static cf16_t *reduce_dense_row_by_dense_new_pivots_ff_16( const uint32_t mul = (uint32_t)(fc - dr[i]); const len_t os = (ncr - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { - dr[j] += mul * pivs[i][l]; + dr[j] += mul * red[l]; } for (; j < ncr; l += UNROLL, j += UNROLL) { - dr[j] += mul * pivs[i][l]; - dr[j+1] += mul * pivs[i][l+1]; - dr[j+2] += mul * pivs[i][l+2]; - dr[j+3] += mul * pivs[i][l+3]; + dr[j] += mul * red[l]; + dr[j+1] += mul * red[l+1]; + dr[j+2] += mul * red[l+2]; + dr[j+3] += mul * red[l+3]; } } if (k == 0) { @@ -731,8 +733,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( const len_t ncl = mat->ncl; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); j = nrl; for (i = 0; i < mat->nru; ++i) { mat->cf_16[j] = bs->cf_16[mat->rr[i][COEFFS]]; @@ -834,7 +835,10 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( } cfs = mat->cf_16[npiv[COEFFS]]; sc = npiv[OFFSET]; - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); } while (!k); bctr++; } @@ -854,8 +858,8 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -869,15 +873,16 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( hm_t sc, cfp; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_16[pivs[k][COEFFS]]; - cfp = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_16[piv[COEFFS]]; + cfp = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -888,12 +893,13 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_16( dr, mat, bs, pivs, sc, cfp, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(pivs); @@ -924,8 +930,7 @@ static void exact_trace_sparse_reduced_echelon_form_ff_16( const int32_t nthrds = st->in_final_reduction_step == 1 ? 1 : st->nthrds; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); j = nrl; for (i = 0; i < mat->nru; ++i) { mat->cf_16[j] = bs->cf_16[mat->rr[i][COEFFS]]; @@ -984,7 +989,10 @@ static void exact_trace_sparse_reduced_echelon_form_ff_16( normalize_sparse_matrix_row_ff_16( mat->cf_16[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_16[npiv[COEFFS]]; } while (!k); } @@ -994,8 +1002,8 @@ static void exact_trace_sparse_reduced_echelon_form_ff_16( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -1008,15 +1016,16 @@ static void exact_trace_sparse_reduced_echelon_form_ff_16( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_16[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_16[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -1027,12 +1036,13 @@ static void exact_trace_sparse_reduced_echelon_form_ff_16( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_16( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(pivs); @@ -1065,12 +1075,14 @@ static void exact_sparse_reduced_echelon_form_ff_16( len_t bad_prime = 0; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = NULL; if (st->in_final_reduction_step == 0) { - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); } else { + pivs = allocate_atomic_pivots(ncols, NULL, 0); for (i = 0; i < mat->nru; ++i) { - pivs[mat->rr[i][OFFSET]] = mat->rr[i]; + atomic_store_explicit(&pivs[mat->rr[i][OFFSET]], mat->rr[i], + memory_order_relaxed); } } j = nrl; @@ -1143,7 +1155,10 @@ static void exact_sparse_reduced_echelon_form_ff_16( normalize_sparse_matrix_row_ff_16( mat->cf_16[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_16[npiv[COEFFS]]; } } while (!k); @@ -1152,8 +1167,8 @@ static void exact_sparse_reduced_echelon_form_ff_16( if (bad_prime == 1) { for (i = 0; i < ncl+ncr; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } mat->np = 0; if (st->info_level > 0) { @@ -1169,8 +1184,8 @@ static void exact_sparse_reduced_echelon_form_ff_16( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -1184,15 +1199,16 @@ static void exact_sparse_reduced_echelon_form_ff_16( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_16[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_16[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -1203,12 +1219,13 @@ static void exact_sparse_reduced_echelon_form_ff_16( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_16( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } mat->tr = realloc(mat->tr, (uint64_t)npivs * sizeof(hi_t *)); @@ -1239,8 +1256,7 @@ static int exact_application_sparse_reduced_echelon_form_ff_16( const int32_t nthrds = st->in_final_reduction_step == 1 ? 1 : st->nthrds; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); j = nrl; for (i = 0; i < mat->nru; ++i) { mat->cf_16[j] = bs->cf_16[mat->rr[i][COEFFS]]; @@ -1301,7 +1317,10 @@ static int exact_application_sparse_reduced_echelon_form_ff_16( normalize_sparse_matrix_row_ff_16( mat->cf_16[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_16[npiv[COEFFS]]; } while (!k); } @@ -1312,8 +1331,8 @@ static int exact_application_sparse_reduced_echelon_form_ff_16( } /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -1326,15 +1345,16 @@ static int exact_application_sparse_reduced_echelon_form_ff_16( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_16[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_16[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -1345,12 +1365,13 @@ static int exact_application_sparse_reduced_echelon_form_ff_16( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_16( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(pivs); @@ -1446,8 +1467,8 @@ static cf16_t **sparse_AB_CD_linear_algebra_ff_16( return drs; } -static cf16_t **interreduce_dense_matrix_ff_16( - cf16_t **dm, +static _Atomic(cf16_t *) *interreduce_dense_matrix_ff_16( + _Atomic(cf16_t *) *dm, const len_t ncr, const uint32_t fc ) @@ -1457,33 +1478,35 @@ static cf16_t **interreduce_dense_matrix_ff_16( for (i = 0; i < ncr; ++i) { k = ncr-1-i; - if (dm[k]) { + cf16_t *row = atomic_load_explicit(&dm[k], memory_order_relaxed); + if (row) { /* printf("interreduce dm[%d]\n", i); */ memset(dr, 0, (uint64_t)ncr * sizeof(int64_t)); const len_t npc = ncr - k; const len_t os = npc % UNROLL; for (j = k, l = 0; l < os; ++j, ++l) { - dr[j] = (int64_t)dm[k][l]; + dr[j] = (int64_t)row[l]; } for (; l < npc; j += UNROLL, l += UNROLL) { - dr[j] = (int64_t)dm[k][l]; - dr[j+1] = (int64_t)dm[k][l+1]; - dr[j+2] = (int64_t)dm[k][l+2]; - dr[j+3] = (int64_t)dm[k][l+3]; + dr[j] = (int64_t)row[l]; + dr[j+1] = (int64_t)row[l+1]; + dr[j+2] = (int64_t)row[l+2]; + dr[j+3] = (int64_t)row[l+3]; } - free(dm[k]); - dm[k] = NULL; + free(row); + atomic_store_explicit(&dm[k], NULL, memory_order_relaxed); /* start with previous pivot the reduction process, so keep the * pivot element as it is */ - dm[k] = reduce_dense_row_by_dense_new_pivots_ff_16( + row = reduce_dense_row_by_dense_new_pivots_ff_16( dr, &k, dm, ncr, fc); + atomic_store_explicit(&dm[k], row, memory_order_relaxed); } } free(dr); return dm; } -static cf16_t **exact_dense_linear_algebra_ff_16( +static _Atomic(cf16_t *) *exact_dense_linear_algebra_ff_16( cf16_t **dm, mat_t *mat, md_t *st @@ -1495,7 +1518,11 @@ static cf16_t **exact_dense_linear_algebra_ff_16( const len_t ncr = mat->ncr; /* rows already representing new pivots */ - cf16_t **nps = (cf16_t **)calloc((uint64_t)ncr, sizeof(cf16_t *)); + _Atomic(cf16_t *) *nps = (_Atomic(cf16_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf16_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } /* rows to be further reduced */ cf16_t **tbr = (cf16_t **)calloc((uint64_t)nrows, sizeof(cf16_t *)); int64_t *dr = (int64_t *)malloc( @@ -1511,16 +1538,16 @@ static cf16_t **exact_dense_linear_algebra_ff_16( while (dm[i][k] == 0) { ++k; } - if (nps[k] == NULL) { + if (atomic_load_explicit(&nps[k], memory_order_relaxed) == NULL) { /* we have a pivot, cut the dense row down to start * at the first nonzero entry */ /* printf("ncr - k = %d - %d\n", ncr, k); */ memmove(dm[i], dm[i]+k, (uint64_t)(ncr-k) * sizeof(cf16_t)); dm[i] = realloc(dm[i], (uint64_t)(ncr-k) * sizeof(cf16_t)); - nps[k] = dm[i]; - if (nps[k][0] != 1) { - nps[k] = normalize_dense_matrix_row_ff_16(nps[k], ncr-k, st->fc); + if (dm[i][0] != 1) { + dm[i] = normalize_dense_matrix_row_ff_16(dm[i], ncr-k, st->fc); } + atomic_store_explicit(&nps[k], dm[i], memory_order_relaxed); } else { tbr[j++] = dm[i]; } @@ -1562,29 +1589,17 @@ static cf16_t **exact_dense_linear_algebra_ff_16( if (npc == -1) { break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, npiv); + cf16_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, npiv, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); } /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -1596,7 +1611,7 @@ static cf16_t **exact_dense_linear_algebra_ff_16( return nps; } -static cf16_t **probabilistic_dense_linear_algebra_ff_16( +static _Atomic(cf16_t *) *probabilistic_dense_linear_algebra_ff_16( cf16_t **dm, mat_t *mat, md_t *st @@ -1610,7 +1625,11 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( const len_t ncr = mat->ncr; /* rows already representing new pivots */ - cf16_t **nps = (cf16_t **)calloc((uint64_t)ncr, sizeof(cf16_t *)); + _Atomic(cf16_t *) *nps = (_Atomic(cf16_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf16_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } /* rows to be further reduced */ cf16_t **tbr = (cf16_t **)calloc((uint64_t)nrows, sizeof(cf16_t *)); @@ -1624,15 +1643,15 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( while (dm[i][k] == 0) { ++k; } - if (nps[k] == NULL) { + if (atomic_load_explicit(&nps[k], memory_order_relaxed) == NULL) { /* we have a pivot, cut the dense row down to start * at the first nonzero entry */ memmove(dm[i], dm[i]+k, (uint64_t)(ncr-k) * sizeof(cf16_t)); dm[i] = realloc(dm[i], (uint64_t)(ncr-k) * sizeof(cf16_t)); - nps[k] = dm[i]; - if (nps[k][0] != 1) { - nps[k] = normalize_dense_matrix_row_ff_16(nps[k], ncr-k, st->fc); + if (dm[i][0] != 1) { + dm[i] = normalize_dense_matrix_row_ff_16(dm[i], ncr-k, st->fc); } + atomic_store_explicit(&nps[k], dm[i], memory_order_relaxed); } else { tbr[j++] = dm[i]; } @@ -1714,7 +1733,10 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( bctr = nrbl; break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, tmp); + cf16_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, tmp, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); @@ -1727,23 +1749,8 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( } } /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -1756,7 +1763,7 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( return nps; } -static cf16_t **probabilistic_sparse_dense_echelon_form_ff_16( +static _Atomic(cf16_t *) *probabilistic_sparse_dense_echelon_form_ff_16( mat_t *mat, const bs_t * const bs, md_t *st @@ -1777,7 +1784,11 @@ static cf16_t **probabilistic_sparse_dense_echelon_form_ff_16( hm_t **upivs = mat->tr; /* rows already representing new pivots */ - cf16_t **nps = (cf16_t **)calloc((uint64_t)ncr, sizeof(cf16_t *)); + _Atomic(cf16_t *) *nps = (_Atomic(cf16_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf16_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } const uint32_t fc = st->fc; const int64_t mod2 = (int64_t)fc * fc; @@ -1852,7 +1863,10 @@ static cf16_t **probabilistic_sparse_dense_echelon_form_ff_16( bctr = nrbl; break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, tmp); + cf16_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, tmp, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); @@ -1866,23 +1880,8 @@ static cf16_t **probabilistic_sparse_dense_echelon_form_ff_16( } npivs = 0; /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -1904,7 +1903,7 @@ static cf16_t **probabilistic_sparse_dense_echelon_form_ff_16( static void convert_to_sparse_matrix_rows_ff_16( mat_t *mat, - cf16_t * const * const dm + _Atomic(cf16_t *) * const dm ) { if (mat->np == 0) { @@ -1925,7 +1924,8 @@ static void convert_to_sparse_matrix_rows_ff_16( l = 0; for (i = 0; i < ncr; ++i) { m = ncr-1-i; - if (dm[m] != NULL) { + const cf16_t * const row = atomic_load_explicit(&dm[m], memory_order_relaxed); + if (row != NULL) { cfs = malloc((uint64_t)(ncr-m) * sizeof(cf16_t)); dts = malloc((uint64_t)(ncr-m+OFFSET) * sizeof(hm_t)); const hm_t len = ncr-m; @@ -1934,26 +1934,26 @@ static void convert_to_sparse_matrix_rows_ff_16( dss = dts + OFFSET; for (k = 0, j = 0; j < os; ++j) { - if (dm[m][j] != 0) { - cfs[k] = dm[m][j]; + if (row[j] != 0) { + cfs[k] = row[j]; dss[k++] = j+shift; } } for (; j < len; j += UNROLL) { - if (dm[m][j] != 0) { - cfs[k] = dm[m][j]; + if (row[j] != 0) { + cfs[k] = row[j]; dss[k++] = j+shift; } - if (dm[m][j+1] != 0) { - cfs[k] = dm[m][j+1]; + if (row[j+1] != 0) { + cfs[k] = row[j+1]; dss[k++] = j+1+shift; } - if (dm[m][j+2] != 0) { - cfs[k] = dm[m][j+2]; + if (row[j+2] != 0) { + cfs[k] = row[j+2]; dss[k++] = j+2+shift; } - if (dm[m][j+3] != 0) { - cfs[k] = dm[m][j+3]; + if (row[j+3] != 0) { + cfs[k] = row[j+3]; dss[k++] = j+3+shift; } } @@ -2127,10 +2127,10 @@ static void exact_sparse_dense_linear_algebra_ff_16( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf16_t **dm; - dm = sparse_AB_CD_linear_algebra_ff_16(mat, bs, st); + _Atomic(cf16_t *) *dm = NULL; + cf16_t **drs = sparse_AB_CD_linear_algebra_ff_16(mat, bs, st); if (mat->np > 0) { - dm = exact_dense_linear_algebra_ff_16(dm, mat, st); + dm = exact_dense_linear_algebra_ff_16(drs, mat, st); dm = interreduce_dense_matrix_ff_16(dm, ncr, st->fc); } @@ -2141,7 +2141,7 @@ static void exact_sparse_dense_linear_algebra_ff_16( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -2177,10 +2177,10 @@ static void probabilistic_sparse_dense_linear_algebra_ff_16_2( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf16_t **dm; - dm = sparse_AB_CD_linear_algebra_ff_16(mat, bs, st); + _Atomic(cf16_t *) *dm = NULL; + cf16_t **drs = sparse_AB_CD_linear_algebra_ff_16(mat, bs, st); if (mat->np > 0) { - dm = probabilistic_dense_linear_algebra_ff_16(dm, mat, st); + dm = probabilistic_dense_linear_algebra_ff_16(drs, mat, st); dm = interreduce_dense_matrix_ff_16(dm, mat->ncr, st->fc); } @@ -2191,7 +2191,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_16_2( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -2227,7 +2227,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_16( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf16_t **dm = NULL; + _Atomic(cf16_t *) *dm = NULL; mat->np = 0; dm = probabilistic_sparse_dense_echelon_form_ff_16(mat, bs, st); dm = interreduce_dense_matrix_ff_16(dm, mat->ncr, st->fc); @@ -2239,7 +2239,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_16( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -2289,12 +2289,13 @@ static void interreduce_matrix_rows_ff_16( mat->cf_16 = realloc(mat->cf_16, (uint64_t)ncols * sizeof(cf16_t *)); memset(mat->cf_16, 0, (uint64_t)ncols * sizeof(cf16_t *)); - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, NULL, 0); /* copy coefficient arrays from basis in matrix, maybe * several rows need the same coefficient arrays, but we * cannot share them here. */ for (i = 0; i < nrows; ++i) { - pivs[mat->rr[i][OFFSET]] = mat->rr[i]; + atomic_store_explicit(&pivs[mat->rr[i][OFFSET]], mat->rr[i], + memory_order_relaxed); } int64_t *dr = (int64_t *)malloc((uint64_t)ncols * sizeof(int64_t)); @@ -2305,14 +2306,15 @@ static void interreduce_matrix_rows_ff_16( k = nrows - 1; for (i = 0; i < ncols; ++i) { l = ncols-1-i; - if (pivs[l] != NULL) { + hm_t *piv = atomic_load_explicit(&pivs[l], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = bs->cf_16[pivs[l][COEFFS]]; - const len_t bi = pivs[l][BINDEX]; - const len_t mh = pivs[l][MULT]; - const len_t os = pivs[l][PRELOOP]; - const len_t len = pivs[l][LENGTH]; - const hm_t * const ds = pivs[l] + OFFSET; + cfs = bs->cf_16[piv[COEFFS]]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -2323,11 +2325,12 @@ static void interreduce_matrix_rows_ff_16( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[l]); - pivs[l] = NULL; - pivs[l] = mat->tr[k--] = + free(piv); + atomic_store_explicit(&pivs[l], NULL, memory_order_relaxed); + mat->tr[k] = reduce_dense_row_by_known_pivots_sparse_ff_16( dr, mat, bs, pivs, sc, l, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[l], mat->tr[k--], memory_order_relaxed); } } for (i = 0; i < ncols; ++i) { diff --git a/src/neogb/la_ff_32.c b/src/neogb/la_ff_32.c index c5794cf8..7a470551 100644 --- a/src/neogb/la_ff_32.c +++ b/src/neogb/la_ff_32.c @@ -309,7 +309,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_17_bit( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -353,7 +353,9 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_17_bit( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -363,7 +365,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_17_bit( /* found reducer row, get multiplier */ const int64_t mul = mod - dr[i]; - dts = pivs[i]; if (i < ncl) { /* set corresponding bit of reducer in reducer bit array */ if (tr > 0) { @@ -584,7 +585,7 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_17_bit( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -609,7 +610,9 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_17_bit( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -619,7 +622,6 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_17_bit( /* found reducer row, get multiplier */ const int64_t mul = mod - dr[i]; - dts = pivs[i]; if (i < ncl) { cfs = bs->cf_32[dts[COEFFS]]; /* set corresponding bit of reducer in reducer bit array */ @@ -968,7 +970,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_31_bit( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t *const *pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -1019,7 +1021,9 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_31_bit( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -1029,7 +1033,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_31_bit( /* found reducer row, get multiplier */ const int64_t mul = (int64_t)dr[i]; - dts = pivs[i]; if (i < ncl) { /* cfs = bs->cf_32[dts[COEFFS]]; */ /* set corresponding bit of reducer in reducer bit array */ @@ -1483,7 +1486,7 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_31_bit( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t *const *pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -1515,7 +1518,9 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_31_bit( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -1525,7 +1530,6 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_31_bit( /* found reducer row, get multiplier */ const int64_t mul = (int64_t)dr[i]; - dts = pivs[i]; if (i < ncl) { cfs = bs->cf_32[dts[COEFFS]]; /* set corresponding bit of reducer in reducer bit array */ @@ -1845,7 +1849,7 @@ static cf32_t *reduce_dense_row_by_all_pivots_17_bit( const bs_t * const bs, len_t *pc, hm_t *const * const pivs, - cf32_t *const * const dpivs, + _Atomic(cf32_t *) * const dpivs, const uint32_t fc ) { @@ -1895,15 +1899,14 @@ static cf32_t *reduce_dense_row_by_all_pivots_17_bit( if (dr[i] == 0) { continue; } - if (dpivs[i-ncl] == NULL) { + red = atomic_load_explicit(&dpivs[i-ncl], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } k++; continue; } - - red = dpivs[i-ncl]; const int64_t mul = mod - dr[i]; const len_t os = (ncols - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { @@ -1942,7 +1945,7 @@ static cf32_t *reduce_dense_row_by_all_pivots_31_bit( const bs_t * const bs, len_t *pc, hm_t * const * const pivs, - cf32_t * const * const dpivs, + _Atomic(cf32_t *) * const dpivs, const uint32_t fc ) { @@ -1998,15 +2001,14 @@ static cf32_t *reduce_dense_row_by_all_pivots_31_bit( if (dr[i] == 0) { continue; } - if (dpivs[i-ncl] == NULL) { + red = atomic_load_explicit(&dpivs[i-ncl], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } k++; continue; } - - red = dpivs[i-ncl]; const int64_t mul = (int64_t) dr[i]; const len_t os = (ncols - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { @@ -2184,7 +2186,7 @@ static cf32_t *reduce_dense_row_by_old_pivots_31_bit( static cf32_t *reduce_dense_row_by_dense_new_pivots_17_bit( int64_t *dr, len_t *pc, - cf32_t * const * const pivs, + _Atomic(cf32_t *) * const pivs, const len_t ncr, const uint32_t fc ) @@ -2200,7 +2202,8 @@ static cf32_t *reduce_dense_row_by_dense_new_pivots_17_bit( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + const cf32_t * const red = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } @@ -2211,13 +2214,13 @@ static cf32_t *reduce_dense_row_by_dense_new_pivots_17_bit( const int64_t mul = mod - dr[i]; const len_t os = (ncr - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { - dr[j] += mul * pivs[i][l]; + dr[j] += mul * red[l]; } for (; j < ncr; l += UNROLL, j += UNROLL) { - dr[j] += mul * pivs[i][l]; - dr[j+1] += mul * pivs[i][l+1]; - dr[j+2] += mul * pivs[i][l+2]; - dr[j+3] += mul * pivs[i][l+3]; + dr[j] += mul * red[l]; + dr[j+1] += mul * red[l+1]; + dr[j+2] += mul * red[l+2]; + dr[j+3] += mul * red[l+3]; } } if (k == 0) { @@ -2242,7 +2245,7 @@ static cf32_t *reduce_dense_row_by_dense_new_pivots_17_bit( static cf32_t *reduce_dense_row_by_dense_new_pivots_31_bit( int64_t *dr, len_t *pc, - cf32_t * const * const pivs, + _Atomic(cf32_t *) * const pivs, const len_t ncr, const uint32_t fc ) @@ -2259,7 +2262,8 @@ static cf32_t *reduce_dense_row_by_dense_new_pivots_31_bit( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + const cf32_t * const red = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } @@ -2270,14 +2274,14 @@ static cf32_t *reduce_dense_row_by_dense_new_pivots_31_bit( const int64_t mul = (int64_t)dr[i]; const len_t os = (ncr - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { - dr[j] -= mul * pivs[i][l]; + dr[j] -= mul * red[l]; dr[j] += (dr[j] >> 63) & mod2; } for (; j < ncr; l+=4, j += UNROLL) { - dr[j] -= mul * pivs[i][l]; - dr[j+1] -= mul * pivs[i][l+1]; - dr[j+2] -= mul * pivs[i][l+2]; - dr[j+3] -= mul * pivs[i][l+3]; + dr[j] -= mul * red[l]; + dr[j+1] -= mul * red[l+1]; + dr[j+2] -= mul * red[l+2]; + dr[j+3] -= mul * red[l+3]; dr[j] += (dr[j] >> 63) & mod2; dr[j+1] += (dr[j+1] >> 63) & mod2; dr[j+2] += (dr[j+2] >> 63) & mod2; @@ -2318,9 +2322,8 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( const len_t ncl = mat->ncl; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - if (mat->rr != NULL) - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots( + ncols, mat->rr, mat->rr != NULL ? mat->nru : 0); j = nrl; for (i = 0; i < mat->nru; ++i) { mat->cf_32[j] = bs->cf_32[mat->rr[i][COEFFS]]; @@ -2439,7 +2442,10 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( } cfs = mat->cf_32[npiv[COEFFS]]; sc = npiv[OFFSET]; - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); } while (!k); bctr++; } @@ -2454,8 +2460,8 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -2469,15 +2475,16 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( hm_t sc, cfp; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_32[pivs[k][COEFFS]]; - cfp = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_32[piv[COEFFS]]; + cfp = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -2488,12 +2495,13 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_32( dr, mat, bs, pivs, sc, cfp, mh, bi, 0, st); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(mat->rr); @@ -2616,12 +2624,14 @@ static void exact_sparse_reduced_echelon_form_ff_32( len_t bad_prime = 0; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = NULL; if (st->in_final_reduction_step == 0) { - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); } else { + pivs = allocate_atomic_pivots(ncols, NULL, 0); for (i = 0; i < mat->nru; ++i) { - pivs[mat->rr[i][OFFSET]] = mat->rr[i]; + atomic_store_explicit(&pivs[mat->rr[i][OFFSET]], mat->rr[i], + memory_order_relaxed); } } j = nrl; @@ -2694,7 +2704,10 @@ static void exact_sparse_reduced_echelon_form_ff_32( normalize_sparse_matrix_row_ff_32( mat->cf_32[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_32[npiv[COEFFS]]; } } while (!k); @@ -2703,8 +2716,8 @@ static void exact_sparse_reduced_echelon_form_ff_32( if (bad_prime == 1) { for (i = 0; i < ncl+ncr; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } mat->np = 0; if (st->info_level > 0) { @@ -2720,8 +2733,8 @@ static void exact_sparse_reduced_echelon_form_ff_32( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -2735,15 +2748,16 @@ static void exact_sparse_reduced_echelon_form_ff_32( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_32[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_32[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -2754,12 +2768,13 @@ static void exact_sparse_reduced_echelon_form_ff_32( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_32( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } mat->tr = realloc(mat->tr, (uint64_t)npivs * sizeof(hi_t *)); @@ -2971,7 +2986,11 @@ static void exact_sparse_reduced_echelon_form_sat_ff_32( pivcf[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + /* the kernel computation is sequential, no other thread adds pivots */ + k = pivs[npiv[OFFSET]] == NULL; + if (k) { + pivs[npiv[OFFSET]] = npiv; + } } while (!k); } @@ -3023,8 +3042,7 @@ static void exact_trace_sparse_reduced_echelon_form_ff_32( const int32_t nthrds = st->in_final_reduction_step == 1 ? 1 : st->nthrds; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); /* unkown pivot rows we have to reduce with the known pivots first */ hm_t **upivs = mat->tr; @@ -3076,7 +3094,10 @@ static void exact_trace_sparse_reduced_echelon_form_ff_32( mat->cf_32[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); st->trace_nr_mult += npiv[LENGTH] / 1000.0; } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_32[npiv[COEFFS]]; } while (!k); } @@ -3086,8 +3107,8 @@ static void exact_trace_sparse_reduced_echelon_form_ff_32( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -3100,15 +3121,16 @@ static void exact_trace_sparse_reduced_echelon_form_ff_32( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_32[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_32[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -3119,12 +3141,13 @@ static void exact_trace_sparse_reduced_echelon_form_ff_32( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_32( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(pivs); @@ -3153,8 +3176,7 @@ static int exact_application_sparse_reduced_echelon_form_ff_32( const int32_t nthrds = st->in_final_reduction_step == 1 ? 1 : st->nthrds; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); /* unkown pivot rows we have to reduce with the known pivots first */ hm_t **upivs = mat->tr; @@ -3210,7 +3232,10 @@ static int exact_application_sparse_reduced_echelon_form_ff_32( mat->cf_32[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); st->application_nr_mult += npiv[LENGTH] / 1000.0; } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_32[npiv[COEFFS]]; } while (!k); } @@ -3221,8 +3246,8 @@ static int exact_application_sparse_reduced_echelon_form_ff_32( } /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -3235,15 +3260,16 @@ static int exact_application_sparse_reduced_echelon_form_ff_32( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_32[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_32[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -3254,12 +3280,13 @@ static int exact_application_sparse_reduced_echelon_form_ff_32( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_32( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(pivs); @@ -3355,8 +3382,8 @@ static cf32_t **sparse_AB_CD_linear_algebra_ff_32( return drs; } -static cf32_t **interreduce_dense_matrix_ff_32( - cf32_t **dm, +static _Atomic(cf32_t *) *interreduce_dense_matrix_ff_32( + _Atomic(cf32_t *) *dm, const len_t ncr, const uint32_t fc ) @@ -3366,32 +3393,34 @@ static cf32_t **interreduce_dense_matrix_ff_32( for (i = 0; i < ncr; ++i) { k = ncr-1-i; - if (dm[k]) { + cf32_t *row = atomic_load_explicit(&dm[k], memory_order_relaxed); + if (row) { memset(dr, 0, (uint64_t)ncr * sizeof(int64_t)); const len_t npc = ncr - k; const len_t os = npc % UNROLL; for (j = k, l = 0; l < os; ++j, ++l) { - dr[j] = (int64_t)dm[k][l]; + dr[j] = (int64_t)row[l]; } for (; l < npc; j += UNROLL, l += UNROLL) { - dr[j] = (int64_t)dm[k][l]; - dr[j+1] = (int64_t)dm[k][l+1]; - dr[j+2] = (int64_t)dm[k][l+2]; - dr[j+3] = (int64_t)dm[k][l+3]; + dr[j] = (int64_t)row[l]; + dr[j+1] = (int64_t)row[l+1]; + dr[j+2] = (int64_t)row[l+2]; + dr[j+3] = (int64_t)row[l+3]; } - free(dm[k]); - dm[k] = NULL; + free(row); + atomic_store_explicit(&dm[k], NULL, memory_order_relaxed); /* start with previous pivot the reduction process, so keep the * pivot element as it is */ - dm[k] = reduce_dense_row_by_dense_new_pivots_ff_32( + row = reduce_dense_row_by_dense_new_pivots_ff_32( dr, &k, dm, ncr, fc); + atomic_store_explicit(&dm[k], row, memory_order_relaxed); } } free(dr); return dm; } -static cf32_t **exact_dense_linear_algebra_ff_32( +static _Atomic(cf32_t *) *exact_dense_linear_algebra_ff_32( cf32_t **dm, mat_t *mat, md_t *st @@ -3403,7 +3432,11 @@ static cf32_t **exact_dense_linear_algebra_ff_32( const len_t ncr = mat->ncr; /* rows already representing new pivots */ - cf32_t **nps = (cf32_t **)calloc((uint64_t)ncr, sizeof(cf32_t *)); + _Atomic(cf32_t *) *nps = (_Atomic(cf32_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf32_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } /* rows to be further reduced */ cf32_t **tbr = (cf32_t **)calloc((uint64_t)nrows, sizeof(cf32_t *)); int64_t *dr = (int64_t *)malloc( @@ -3419,15 +3452,15 @@ static cf32_t **exact_dense_linear_algebra_ff_32( while (dm[i][k] == 0) { ++k; } - if (nps[k] == NULL) { + if (atomic_load_explicit(&nps[k], memory_order_relaxed) == NULL) { /* we have a pivot, cut the dense row down to start * at the first nonzero entry */ memmove(dm[i], dm[i]+k, (uint64_t)(ncr-k) * sizeof(cf32_t)); dm[i] = realloc(dm[i], (uint64_t)(ncr-k) * sizeof(cf32_t)); - nps[k] = dm[i]; - if (nps[k][0] != 1) { - nps[k] = normalize_dense_matrix_row_ff_32(nps[k], ncr-k, st->fc); + if (dm[i][0] != 1) { + dm[i] = normalize_dense_matrix_row_ff_32(dm[i], ncr-k, st->fc); } + atomic_store_explicit(&nps[k], dm[i], memory_order_relaxed); } else { tbr[j++] = dm[i]; } @@ -3469,29 +3502,17 @@ static cf32_t **exact_dense_linear_algebra_ff_32( if (npc == -1) { break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, npiv); + cf32_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, npiv, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); } /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -3503,7 +3524,7 @@ static cf32_t **exact_dense_linear_algebra_ff_32( return nps; } -static cf32_t **probabilistic_dense_linear_algebra_ff_32( +static _Atomic(cf32_t *) *probabilistic_dense_linear_algebra_ff_32( cf32_t **dm, mat_t *mat, md_t *st @@ -3517,7 +3538,11 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( const len_t ncr = mat->ncr; /* rows already representing new pivots */ - cf32_t **nps = (cf32_t **)calloc((uint64_t)ncr, sizeof(cf32_t *)); + _Atomic(cf32_t *) *nps = (_Atomic(cf32_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf32_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } /* rows to be further reduced */ cf32_t **tbr = (cf32_t **)calloc((uint64_t)nrows, sizeof(cf32_t *)); @@ -3531,15 +3556,15 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( while (dm[i][k] == 0) { ++k; } - if (nps[k] == NULL) { + if (atomic_load_explicit(&nps[k], memory_order_relaxed) == NULL) { /* we have a pivot, cut the dense row down to start * at the first nonzero entry */ memmove(dm[i], dm[i]+k, (uint64_t)(ncr-k) * sizeof(cf32_t)); dm[i] = realloc(dm[i], (uint64_t)(ncr-k) * sizeof(cf32_t)); - nps[k] = dm[i]; - if (nps[k][0] != 1) { - nps[k] = normalize_dense_matrix_row_ff_32(nps[k], ncr-k, st->fc); + if (dm[i][0] != 1) { + dm[i] = normalize_dense_matrix_row_ff_32(dm[i], ncr-k, st->fc); } + atomic_store_explicit(&nps[k], dm[i], memory_order_relaxed); } else { tbr[j++] = dm[i]; } @@ -3640,7 +3665,10 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( bctr = nrbl; break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, tmp); + cf32_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, tmp, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); @@ -3653,23 +3681,8 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( } } /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -3682,7 +3695,7 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( return nps; } -static cf32_t **probabilistic_sparse_dense_echelon_form_ff_32( +static _Atomic(cf32_t *) *probabilistic_sparse_dense_echelon_form_ff_32( mat_t *mat, const bs_t * const bs, md_t *st @@ -3703,7 +3716,11 @@ static cf32_t **probabilistic_sparse_dense_echelon_form_ff_32( hm_t **upivs = mat->tr; /* rows already representing new pivots */ - cf32_t **nps = (cf32_t **)calloc((uint64_t)ncr, sizeof(cf32_t *)); + _Atomic(cf32_t *) *nps = (_Atomic(cf32_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf32_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } const uint32_t fc = st->fc; const int64_t mod2 = (int64_t)fc * fc; @@ -3778,7 +3795,10 @@ static cf32_t **probabilistic_sparse_dense_echelon_form_ff_32( bctr = nrbl; break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, tmp); + cf32_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, tmp, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); @@ -3792,23 +3812,8 @@ static cf32_t **probabilistic_sparse_dense_echelon_form_ff_32( } npivs = 0; /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -3830,7 +3835,7 @@ static cf32_t **probabilistic_sparse_dense_echelon_form_ff_32( static void convert_to_sparse_matrix_rows_ff_32( mat_t *mat, - cf32_t * const * const dm + _Atomic(cf32_t *) * const dm ) { if (mat->np == 0) { @@ -3851,7 +3856,8 @@ static void convert_to_sparse_matrix_rows_ff_32( l = 0; for (i = 0; i < ncr; ++i) { m = ncr-1-i; - if (dm[m] != NULL) { + const cf32_t * const row = atomic_load_explicit(&dm[m], memory_order_relaxed); + if (row != NULL) { cfs = malloc((uint64_t)(ncr-m) * sizeof(cf32_t)); dts = malloc((uint64_t)(ncr-m+OFFSET) * sizeof(hm_t)); const hm_t len = ncr-m; @@ -3860,26 +3866,26 @@ static void convert_to_sparse_matrix_rows_ff_32( dss = dts + OFFSET; for (k = 0, j = 0; j < os; ++j) { - if (dm[m][j] != 0) { - cfs[k] = dm[m][j]; + if (row[j] != 0) { + cfs[k] = row[j]; dss[k++] = j+shift; } } for (; j < len; j += UNROLL) { - if (dm[m][j] != 0) { - cfs[k] = dm[m][j]; + if (row[j] != 0) { + cfs[k] = row[j]; dss[k++] = j+shift; } - if (dm[m][j+1] != 0) { - cfs[k] = dm[m][j+1]; + if (row[j+1] != 0) { + cfs[k] = row[j+1]; dss[k++] = j+1+shift; } - if (dm[m][j+2] != 0) { - cfs[k] = dm[m][j+2]; + if (row[j+2] != 0) { + cfs[k] = row[j+2]; dss[k++] = j+2+shift; } - if (dm[m][j+3] != 0) { - cfs[k] = dm[m][j+3]; + if (row[j+3] != 0) { + cfs[k] = row[j+3]; dss[k++] = j+3+shift; } } @@ -4128,10 +4134,10 @@ static void exact_sparse_dense_linear_algebra_ff_32( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf32_t **dm; - dm = sparse_AB_CD_linear_algebra_ff_32(mat, bs, st); + _Atomic(cf32_t *) *dm = NULL; + cf32_t **drs = sparse_AB_CD_linear_algebra_ff_32(mat, bs, st); if (mat->np > 0) { - dm = exact_dense_linear_algebra_ff_32(dm, mat, st); + dm = exact_dense_linear_algebra_ff_32(drs, mat, st); dm = interreduce_dense_matrix_ff_32(dm, ncr, st->fc); } @@ -4142,7 +4148,7 @@ static void exact_sparse_dense_linear_algebra_ff_32( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -4178,10 +4184,10 @@ static void probabilistic_sparse_dense_linear_algebra_ff_32_2( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf32_t **dm; - dm = sparse_AB_CD_linear_algebra_ff_32(mat, bs, st); + _Atomic(cf32_t *) *dm = NULL; + cf32_t **drs = sparse_AB_CD_linear_algebra_ff_32(mat, bs, st); if (mat->np > 0) { - dm = probabilistic_dense_linear_algebra_ff_32(dm, mat, st); + dm = probabilistic_dense_linear_algebra_ff_32(drs, mat, st); dm = interreduce_dense_matrix_ff_32(dm, mat->ncr, st->fc); } @@ -4192,7 +4198,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_32_2( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -4228,7 +4234,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_32( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf32_t **dm = NULL; + _Atomic(cf32_t *) *dm = NULL; mat->np = 0; dm = probabilistic_sparse_dense_echelon_form_ff_32(mat, bs, st); dm = interreduce_dense_matrix_ff_32(dm, mat->ncr, st->fc); @@ -4240,7 +4246,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_32( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -4281,12 +4287,13 @@ static void interreduce_matrix_rows_ff_32( mat->cf_32 = realloc(mat->cf_32, (uint64_t)ncols * sizeof(cf32_t *)); memset(mat->cf_32, 0, (uint64_t)ncols * sizeof(cf32_t *)); - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, NULL, 0); /* copy coefficient arrays from basis in matrix, maybe * several rows need the same coefficient arrays, but we * cannot share them here. */ for (i = 0; i < nrows; ++i) { - pivs[mat->rr[i][OFFSET]] = mat->rr[i]; + atomic_store_explicit(&pivs[mat->rr[i][OFFSET]], mat->rr[i], + memory_order_relaxed); } int64_t *dr = (int64_t *)malloc((uint64_t)ncols * sizeof(int64_t)); @@ -4297,14 +4304,15 @@ static void interreduce_matrix_rows_ff_32( k = nrows - 1; for (i = 0; i < ncols; ++i) { l = ncols-1-i; - if (pivs[l] != NULL) { + hm_t *piv = atomic_load_explicit(&pivs[l], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = bs->cf_32[pivs[l][COEFFS]]; - const len_t os = pivs[l][PRELOOP]; - const len_t len = pivs[l][LENGTH]; - const len_t bi = pivs[l][BINDEX]; - const len_t mh = pivs[l][MULT]; - const hm_t * const ds = pivs[l] + OFFSET; + cfs = bs->cf_32[piv[COEFFS]]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -4315,11 +4323,12 @@ static void interreduce_matrix_rows_ff_32( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[l]); - pivs[l] = NULL; - pivs[l] = mat->tr[k--] = + free(piv); + atomic_store_explicit(&pivs[l], NULL, memory_order_relaxed); + mat->tr[k] = reduce_dense_row_by_known_pivots_sparse_ff_32( dr, mat, bs, pivs, sc, l, mh, bi, 0, st); + atomic_store_explicit(&pivs[l], mat->tr[k--], memory_order_relaxed); } } if (free_basis != 0) { diff --git a/src/neogb/la_ff_8.c b/src/neogb/la_ff_8.c index 2f740081..3da7b042 100644 --- a/src/neogb/la_ff_8.c +++ b/src/neogb/la_ff_8.c @@ -98,7 +98,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_8( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -154,7 +154,9 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_8( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -163,7 +165,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_8( } /* found reducer row, get multiplier */ const uint32_t mul= (uint32_t)(fc - dr[i]); - dts = pivs[i]; if (i < ncl) { /* set corresponding bit of reducer in reducer bit array */ if (tr > 0) { @@ -565,7 +566,7 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_8( int64_t *dr, mat_t *mat, const bs_t * const bs, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos, /* position of new coeffs array in tmpcf */ const len_t mh, /* multiplier hash for tracing */ @@ -590,7 +591,9 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_8( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { np = i; } @@ -599,7 +602,6 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_ff_8( } /* found reducer row, get multiplier */ const uint32_t mul= (uint32_t)(fc - dr[i]); - dts = pivs[i]; if (i < ncl) { cfs = bs->cf_8[dts[COEFFS]]; /* set corresponding bit of reducer in reducer bit array */ @@ -652,7 +654,7 @@ static cf8_t *reduce_dense_row_by_all_pivots_ff_8( const bs_t * const bs, len_t *pc, hm_t *const * const pivs, - cf8_t *const * const dpivs, + _Atomic(cf8_t *) * const dpivs, const uint32_t fc ) { @@ -702,15 +704,14 @@ static cf8_t *reduce_dense_row_by_all_pivots_ff_8( if (dr[i] == 0) { continue; } - if (dpivs[i-ncl] == NULL) { + red = atomic_load_explicit(&dpivs[i-ncl], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } k++; continue; } - - red = dpivs[i-ncl]; const uint32_t mul= (uint32_t)(fc - dr[i]); const len_t os = (ncols - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { @@ -810,7 +811,7 @@ static cf8_t *reduce_dense_row_by_old_pivots_ff_8( static cf8_t *reduce_dense_row_by_dense_new_pivots_ff_8( int64_t *dr, len_t *pc, - cf8_t * const * const pivs, + _Atomic(cf8_t *) * const pivs, const len_t ncr, const uint32_t fc ) @@ -826,7 +827,8 @@ static cf8_t *reduce_dense_row_by_dense_new_pivots_ff_8( if (dr[i] == 0) { continue; } - if (pivs[i] == NULL) { + const cf8_t * const red = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (red == NULL) { if (np == -1) { np = i; } @@ -837,13 +839,13 @@ static cf8_t *reduce_dense_row_by_dense_new_pivots_ff_8( const uint32_t mul= (uint32_t)(fc - dr[i]); const len_t os = (ncr - i) % UNROLL; for (l = 0, j = i; l < os; ++l, ++j) { - dr[j] += mul * pivs[i][l]; + dr[j] += mul * red[l]; } for (; j < ncr; l += UNROLL, j += UNROLL) { - dr[j] += mul * pivs[i][l]; - dr[j+1] += mul * pivs[i][l+1]; - dr[j+2] += mul * pivs[i][l+2]; - dr[j+3] += mul * pivs[i][l+3]; + dr[j] += mul * red[l]; + dr[j+1] += mul * red[l+1]; + dr[j+2] += mul * red[l+2]; + dr[j+3] += mul * red[l+3]; } } if (k == 0) { @@ -879,8 +881,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( const len_t ncl = mat->ncl; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); j = nrl; for (i = 0; i < mat->nru; ++i) { mat->cf_8[j] = bs->cf_8[mat->rr[i][COEFFS]]; @@ -982,7 +983,10 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( } cfs = mat->cf_8[npiv[COEFFS]]; sc = npiv[OFFSET]; - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); } while (!k); bctr++; } @@ -997,8 +1001,8 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -1012,15 +1016,16 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( hm_t sc, cfp; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_8[pivs[k][COEFFS]]; - cfp = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_8[piv[COEFFS]]; + cfp = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -1031,12 +1036,13 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_8( dr, mat, bs, pivs, sc, cfp, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(mat->rr); @@ -1069,8 +1075,7 @@ static int exact_application_sparse_reduced_echelon_form_ff_8( const int32_t nthrds = st->in_final_reduction_step == 1 ? 1 : st->nthrds; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); /* unkown pivot rows we have to reduce with the known pivots first */ hm_t **upivs = mat->tr; @@ -1124,7 +1129,10 @@ static int exact_application_sparse_reduced_echelon_form_ff_8( normalize_sparse_matrix_row_ff_8( mat->cf_8[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_8[npiv[COEFFS]]; } while (!k); } @@ -1135,8 +1143,8 @@ static int exact_application_sparse_reduced_echelon_form_ff_8( } /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -1149,15 +1157,16 @@ static int exact_application_sparse_reduced_echelon_form_ff_8( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_8[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_8[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -1168,12 +1177,13 @@ static int exact_application_sparse_reduced_echelon_form_ff_8( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_8( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(pivs); @@ -1205,8 +1215,7 @@ static void exact_trace_sparse_reduced_echelon_form_ff_8( const int32_t nthrds = st->in_final_reduction_step == 1 ? 1 : st->nthrds; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); /* unkown pivot rows we have to reduce with the known pivots first */ hm_t **upivs = mat->tr; @@ -1258,7 +1267,10 @@ static void exact_trace_sparse_reduced_echelon_form_ff_8( normalize_sparse_matrix_row_ff_8( mat->cf_8[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_8[npiv[COEFFS]]; } while (!k); } @@ -1268,8 +1280,8 @@ static void exact_trace_sparse_reduced_echelon_form_ff_8( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -1282,15 +1294,16 @@ static void exact_trace_sparse_reduced_echelon_form_ff_8( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_8[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_8[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -1301,12 +1314,13 @@ static void exact_trace_sparse_reduced_echelon_form_ff_8( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_8( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } free(pivs); @@ -1338,12 +1352,14 @@ static void exact_sparse_reduced_echelon_form_ff_8( len_t bad_prime = 0; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = NULL; if (st->in_final_reduction_step == 0) { - memcpy(pivs, mat->rr, (uint64_t)mat->nru * sizeof(hm_t *)); + pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); } else { + pivs = allocate_atomic_pivots(ncols, NULL, 0); for (i = 0; i < mat->nru; ++i) { - pivs[mat->rr[i][OFFSET]] = mat->rr[i]; + atomic_store_explicit(&pivs[mat->rr[i][OFFSET]], mat->rr[i], + memory_order_relaxed); } } j = nrl; @@ -1416,7 +1432,10 @@ static void exact_sparse_reduced_echelon_form_ff_8( normalize_sparse_matrix_row_ff_8( mat->cf_8[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH], st->fc); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_8[npiv[COEFFS]]; } } while (!k); @@ -1425,8 +1444,8 @@ static void exact_sparse_reduced_echelon_form_ff_8( if (bad_prime == 1) { for (i = 0; i < ncl+ncr; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } mat->np = 0; if (st->info_level > 0) { @@ -1442,8 +1461,8 @@ static void exact_sparse_reduced_echelon_form_ff_8( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -1457,15 +1476,16 @@ static void exact_sparse_reduced_echelon_form_ff_8( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = mat->cf_8[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const len_t bi = pivs[k][BINDEX]; - const len_t mh = pivs[k][MULT]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_8[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -1476,12 +1496,13 @@ static void exact_sparse_reduced_echelon_form_ff_8( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs++] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_ff_8( dr, mat, bs, pivs, sc, cf_array_pos, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[k], mat->tr[npivs++], memory_order_relaxed); } } mat->tr = realloc(mat->tr, (uint64_t)npivs * sizeof(hi_t *)); @@ -1577,8 +1598,8 @@ static cf8_t **sparse_AB_CD_linear_algebra_ff_8( return drs; } -static cf8_t **interreduce_dense_matrix_ff_8( - cf8_t **dm, +static _Atomic(cf8_t *) *interreduce_dense_matrix_ff_8( + _Atomic(cf8_t *) *dm, const len_t ncr, const uint32_t fc ) @@ -1588,32 +1609,34 @@ static cf8_t **interreduce_dense_matrix_ff_8( for (i = 0; i < ncr; ++i) { k = ncr-1-i; - if (dm[k]) { + cf8_t *row = atomic_load_explicit(&dm[k], memory_order_relaxed); + if (row) { memset(dr, 0, (uint64_t)ncr * sizeof(int64_t)); const len_t npc = ncr - k; const len_t os = npc % UNROLL; for (j = k, l = 0; l < os; ++j, ++l) { - dr[j] = (int64_t)dm[k][l]; + dr[j] = (int64_t)row[l]; } for (; l < npc; j += UNROLL, l += UNROLL) { - dr[j] = (int64_t)dm[k][l]; - dr[j+1] = (int64_t)dm[k][l+1]; - dr[j+2] = (int64_t)dm[k][l+2]; - dr[j+3] = (int64_t)dm[k][l+3]; + dr[j] = (int64_t)row[l]; + dr[j+1] = (int64_t)row[l+1]; + dr[j+2] = (int64_t)row[l+2]; + dr[j+3] = (int64_t)row[l+3]; } - free(dm[k]); - dm[k] = NULL; + free(row); + atomic_store_explicit(&dm[k], NULL, memory_order_relaxed); /* start with previous pivot the reduction process, so keep the * pivot element as it is */ - dm[k] = reduce_dense_row_by_dense_new_pivots_ff_8( + row = reduce_dense_row_by_dense_new_pivots_ff_8( dr, &k, dm, ncr, fc); + atomic_store_explicit(&dm[k], row, memory_order_relaxed); } } free(dr); return dm; } -static cf8_t **exact_dense_linear_algebra_ff_8( +static _Atomic(cf8_t *) *exact_dense_linear_algebra_ff_8( cf8_t **dm, mat_t *mat, md_t *st @@ -1625,7 +1648,11 @@ static cf8_t **exact_dense_linear_algebra_ff_8( const len_t ncr = mat->ncr; /* rows already representing new pivots */ - cf8_t **nps = (cf8_t **)calloc((uint64_t)ncr, sizeof(cf8_t *)); + _Atomic(cf8_t *) *nps = (_Atomic(cf8_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf8_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } /* rows to be further reduced */ cf8_t **tbr = (cf8_t **)calloc((uint64_t)nrows, sizeof(cf8_t *)); int64_t *dr = (int64_t *)malloc( @@ -1641,15 +1668,15 @@ static cf8_t **exact_dense_linear_algebra_ff_8( while (dm[i][k] == 0) { ++k; } - if (nps[k] == NULL) { + if (atomic_load_explicit(&nps[k], memory_order_relaxed) == NULL) { /* we have a pivot, cut the dense row down to start * at the first nonzero entry */ memmove(dm[i], dm[i]+k, (uint64_t)(ncr-k) * sizeof(cf8_t)); dm[i] = realloc(dm[i], (uint64_t)(ncr-k) * sizeof(cf8_t)); - nps[k] = dm[i]; - if (nps[k][0] != 1) { - nps[k] = normalize_dense_matrix_row_ff_8(nps[k], ncr-k, st->fc); + if (dm[i][0] != 1) { + dm[i] = normalize_dense_matrix_row_ff_8(dm[i], ncr-k, st->fc); } + atomic_store_explicit(&nps[k], dm[i], memory_order_relaxed); } else { tbr[j++] = dm[i]; } @@ -1691,29 +1718,17 @@ static cf8_t **exact_dense_linear_algebra_ff_8( if (npc == -1) { break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, npiv); + cf8_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, npiv, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); } /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -1725,7 +1740,7 @@ static cf8_t **exact_dense_linear_algebra_ff_8( return nps; } -static cf8_t **probabilistic_dense_linear_algebra_ff_8( +static _Atomic(cf8_t *) *probabilistic_dense_linear_algebra_ff_8( cf8_t **dm, mat_t *mat, md_t *st @@ -1739,7 +1754,11 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( const len_t ncr = mat->ncr; /* rows already representing new pivots */ - cf8_t **nps = (cf8_t **)calloc((uint64_t)ncr, sizeof(cf8_t *)); + _Atomic(cf8_t *) *nps = (_Atomic(cf8_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf8_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } /* rows to be further reduced */ cf8_t **tbr = (cf8_t **)calloc((uint64_t)nrows, sizeof(cf8_t *)); @@ -1753,15 +1772,15 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( while (dm[i][k] == 0) { ++k; } - if (nps[k] == NULL) { + if (atomic_load_explicit(&nps[k], memory_order_relaxed) == NULL) { /* we have a pivot, cut the dense row down to start * at the first nonzero entry */ memmove(dm[i], dm[i]+k, (uint64_t)(ncr-k) * sizeof(cf8_t)); dm[i] = realloc(dm[i], (uint64_t)(ncr-k) * sizeof(cf8_t)); - nps[k] = dm[i]; - if (nps[k][0] != 1) { - nps[k] = normalize_dense_matrix_row_ff_8(nps[k], ncr-k, st->fc); + if (dm[i][0] != 1) { + dm[i] = normalize_dense_matrix_row_ff_8(dm[i], ncr-k, st->fc); } + atomic_store_explicit(&nps[k], dm[i], memory_order_relaxed); } else { tbr[j++] = dm[i]; } @@ -1843,7 +1862,10 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( bctr = nrbl; break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, tmp); + cf8_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, tmp, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); @@ -1856,23 +1878,8 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( } } /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -1885,7 +1892,7 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( return nps; } -static cf8_t **probabilistic_sparse_dense_echelon_form_ff_8( +static _Atomic(cf8_t *) *probabilistic_sparse_dense_echelon_form_ff_8( mat_t *mat, const bs_t * const bs, md_t *st @@ -1906,7 +1913,11 @@ static cf8_t **probabilistic_sparse_dense_echelon_form_ff_8( hm_t **upivs = mat->tr; /* rows already representing new pivots */ - cf8_t **nps = (cf8_t **)calloc((uint64_t)ncr, sizeof(cf8_t *)); + _Atomic(cf8_t *) *nps = (_Atomic(cf8_t *) *)malloc( + (uint64_t)ncr * sizeof(_Atomic(cf8_t *))); + for (i = 0; i < ncr; ++i) { + atomic_init(&nps[i], NULL); + } const uint32_t fc = st->fc; const int64_t mod2 = (int64_t)fc * fc; @@ -1981,7 +1992,10 @@ static cf8_t **probabilistic_sparse_dense_echelon_form_ff_8( bctr = nrbl; break; } - k = __sync_bool_compare_and_swap(&nps[npc], NULL, tmp); + cf8_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &nps[npc], &expected, tmp, + memory_order_release, memory_order_relaxed); /* some other thread has already added a pivot so we have to * recall the dense reduction process */ } while (!k); @@ -1995,23 +2009,8 @@ static cf8_t **probabilistic_sparse_dense_echelon_form_ff_8( } npivs = 0; /* count number of pivots */ - const len_t os = ncr % UNROLL; - for (i = 0; i < os; ++i) { - if (nps[i] != NULL) { - npivs++; - } - } - for (; i < ncr; i += UNROLL) { - if (nps[i] != NULL) { - npivs++; - } - if (nps[i+1] != NULL) { - npivs++; - } - if (nps[i+2] != NULL) { - npivs++; - } - if (nps[i+3] != NULL) { + for (i = 0; i < ncr; ++i) { + if (atomic_load_explicit(&nps[i], memory_order_relaxed) != NULL) { npivs++; } } @@ -2033,7 +2032,7 @@ static cf8_t **probabilistic_sparse_dense_echelon_form_ff_8( static void convert_to_sparse_matrix_rows_ff_8( mat_t *mat, - cf8_t * const * const dm + _Atomic(cf8_t *) * const dm ) { if (mat->np == 0) { @@ -2054,7 +2053,8 @@ static void convert_to_sparse_matrix_rows_ff_8( l = 0; for (i = 0; i < ncr; ++i) { m = ncr-1-i; - if (dm[m] != NULL) { + const cf8_t * const row = atomic_load_explicit(&dm[m], memory_order_relaxed); + if (row != NULL) { cfs = malloc((uint64_t)(ncr-m) * sizeof(cf8_t)); dts = malloc((uint64_t)(ncr-m+OFFSET) * sizeof(hm_t)); const hm_t len = ncr-m; @@ -2063,26 +2063,26 @@ static void convert_to_sparse_matrix_rows_ff_8( dss = dts + OFFSET; for (k = 0, j = 0; j < os; ++j) { - if (dm[m][j] != 0) { - cfs[k] = dm[m][j]; + if (row[j] != 0) { + cfs[k] = row[j]; dss[k++] = j+shift; } } for (; j < len; j += UNROLL) { - if (dm[m][j] != 0) { - cfs[k] = dm[m][j]; + if (row[j] != 0) { + cfs[k] = row[j]; dss[k++] = j+shift; } - if (dm[m][j+1] != 0) { - cfs[k] = dm[m][j+1]; + if (row[j+1] != 0) { + cfs[k] = row[j+1]; dss[k++] = j+1+shift; } - if (dm[m][j+2] != 0) { - cfs[k] = dm[m][j+2]; + if (row[j+2] != 0) { + cfs[k] = row[j+2]; dss[k++] = j+2+shift; } - if (dm[m][j+3] != 0) { - cfs[k] = dm[m][j+3]; + if (row[j+3] != 0) { + cfs[k] = row[j+3]; dss[k++] = j+3+shift; } } @@ -2254,10 +2254,10 @@ static void exact_sparse_dense_linear_algebra_ff_8( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf8_t **dm; - dm = sparse_AB_CD_linear_algebra_ff_8(mat, bs, st); + _Atomic(cf8_t *) *dm = NULL; + cf8_t **drs = sparse_AB_CD_linear_algebra_ff_8(mat, bs, st); if (mat->np > 0) { - dm = exact_dense_linear_algebra_ff_8(dm, mat, st); + dm = exact_dense_linear_algebra_ff_8(drs, mat, st); dm = interreduce_dense_matrix_ff_8(dm, ncr, st->fc); } @@ -2268,7 +2268,7 @@ static void exact_sparse_dense_linear_algebra_ff_8( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -2304,10 +2304,10 @@ static void probabilistic_sparse_dense_linear_algebra_ff_8_2( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf8_t **dm; - dm = sparse_AB_CD_linear_algebra_ff_8(mat, bs, st); + _Atomic(cf8_t *) *dm = NULL; + cf8_t **drs = sparse_AB_CD_linear_algebra_ff_8(mat, bs, st); if (mat->np > 0) { - dm = probabilistic_dense_linear_algebra_ff_8(dm, mat, st); + dm = probabilistic_dense_linear_algebra_ff_8(drs, mat, st); dm = interreduce_dense_matrix_ff_8(dm, mat->ncr, st->fc); } @@ -2318,7 +2318,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_8_2( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -2354,7 +2354,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_8( const len_t ncr = mat->ncr; /* generate updated dense D part via reduction of CD with AB */ - cf8_t **dm = NULL; + _Atomic(cf8_t *) *dm = NULL; mat->np = 0; dm = probabilistic_sparse_dense_echelon_form_ff_8(mat, bs, st); dm = interreduce_dense_matrix_ff_8(dm, mat->ncr, st->fc); @@ -2366,7 +2366,7 @@ static void probabilistic_sparse_dense_linear_algebra_ff_8( /* free dm */ if (dm) { for (i = 0; i < ncr; ++i) { - free(dm[i]); + free(atomic_load_explicit(&dm[i], memory_order_relaxed)); } free(dm); dm = NULL; @@ -2416,12 +2416,13 @@ static void interreduce_matrix_rows_ff_8( mat->cf_8 = realloc(mat->cf_8, (uint64_t)ncols * sizeof(cf32_t *)); memset(mat->cf_8, 0, (uint64_t)ncols * sizeof(cf8_t *)); - hm_t **pivs = (hm_t **)calloc((uint64_t)ncols, sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, NULL, 0); /* copy coefficient arrays from basis in matrix, maybe * several rows need the same coefficient arrays, but we * cannot share them here. */ for (i = 0; i < nrows; ++i) { - pivs[mat->rr[i][OFFSET]] = mat->rr[i]; + atomic_store_explicit(&pivs[mat->rr[i][OFFSET]], mat->rr[i], + memory_order_relaxed); } int64_t *dr = (int64_t *)malloc((uint64_t)ncols * sizeof(int64_t)); @@ -2432,14 +2433,15 @@ static void interreduce_matrix_rows_ff_8( k = nrows - 1; for (i = 0; i < ncols; ++i) { l = ncols-1-i; - if (pivs[l] != NULL) { + hm_t *piv = atomic_load_explicit(&pivs[l], memory_order_relaxed); + if (piv) { memset(dr, 0, (uint64_t)ncols * sizeof(int64_t)); - cfs = bs->cf_8[pivs[l][COEFFS]]; - const len_t bi = pivs[l][BINDEX]; - const len_t mh = pivs[l][MULT]; - const len_t os = pivs[l][PRELOOP]; - const len_t len = pivs[l][LENGTH]; - const hm_t * const ds = pivs[l] + OFFSET; + cfs = bs->cf_8[piv[COEFFS]]; + const len_t bi = piv[BINDEX]; + const len_t mh = piv[MULT]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { dr[ds[j]] = (int64_t)cfs[j]; @@ -2450,11 +2452,12 @@ static void interreduce_matrix_rows_ff_8( dr[ds[j+2]] = (int64_t)cfs[j+2]; dr[ds[j+3]] = (int64_t)cfs[j+3]; } - free(pivs[l]); - pivs[l] = NULL; - pivs[l] = mat->tr[k--] = + free(piv); + atomic_store_explicit(&pivs[l], NULL, memory_order_relaxed); + mat->tr[k] = reduce_dense_row_by_known_pivots_sparse_ff_8( dr, mat, bs, pivs, sc, l, mh, bi, 0, st->fc); + atomic_store_explicit(&pivs[l], mat->tr[k--], memory_order_relaxed); } } for (i = 0; i < ncols; ++i) { diff --git a/src/neogb/la_qq.c b/src/neogb/la_qq.c index f94a8c89..72c7631e 100644 --- a/src/neogb/la_qq.c +++ b/src/neogb/la_qq.c @@ -72,7 +72,7 @@ static inline mpz_t *remove_content_of_sparse_matrix_row_qq( static hm_t *reduce_dense_row_by_known_pivots_sparse_only_ab_qq( mpz_t *dr, mat_t *mat, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos /* position of new coeffs array in tmpcf */ ) @@ -95,7 +95,9 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_only_ab_qq( if (mpz_sgn(dr[i]) == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { row = (hm_t *)malloc( (unsigned long)(ncols-i+OFFSET) * sizeof(hm_t)); @@ -110,7 +112,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_only_ab_qq( continue; } /* found reducer row, get multiplier */ - dts = pivs[i]; cfs = mat->cf_ab_qq[dts[COEFFS]]; const len_t os = dts[PRELOOP]; const len_t len = dts[LENGTH]; @@ -167,7 +168,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_only_ab_qq( static hm_t *reduce_dense_row_by_known_pivots_sparse_ab_first_qq( mpz_t *dr, mat_t *mat, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos /* position of new coeffs array in tmpcf */ ) @@ -191,7 +192,9 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ab_first_qq( if (mpz_sgn(dr[i]) == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { row = (hm_t *)malloc( (unsigned long)(ncols-i+OFFSET) * sizeof(hm_t)); @@ -206,7 +209,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ab_first_qq( continue; } /* found reducer row, get multiplier */ - dts = pivs[i]; if (i < ncl) { cfs = mat->cf_ab_qq[dts[COEFFS]]; } else { @@ -259,7 +261,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_qq( mpz_t *dr, mat_t *mat, const bs_t * const bs, - hm_t * const * const pivs, + _Atomic(hm_t *) * const pivs, const hi_t dpiv, /* pivot of dense row at the beginning */ const hm_t tmp_pos /* position of new coeffs array in tmpcf */ ) @@ -283,7 +285,9 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_qq( if (mpz_sgn(dr[i]) == 0) { continue; } - if (pivs[i] == NULL) { + /* pairs with the release of the thread publishing the pivot */ + dts = atomic_load_explicit(&pivs[i], memory_order_acquire); + if (dts == NULL) { if (np == -1) { row = (hm_t *)malloc( (unsigned long)(ncols-i+OFFSET) * sizeof(hm_t)); @@ -298,7 +302,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_qq( continue; } /* found reducer row, get multiplier */ - dts = pivs[i]; if (i < ncl) { cfs = bs->cf_qq[dts[COEFFS]]; } else { @@ -366,8 +369,7 @@ static void exact_sparse_reduced_echelon_form_ab_first_qq( hm_t cf_array_pos; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((unsigned long)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (unsigned long)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); /* unkown pivot rows we have to reduce with the known pivots first */ hm_t **upivs = mat->tr; @@ -377,21 +379,23 @@ static void exact_sparse_reduced_echelon_form_ab_first_qq( mpz_init(dr[i]); } + hm_t *piv = atomic_load_explicit(&pivs[nru-1], memory_order_relaxed); mat->cf_ab_qq[nru-1] = (mpz_t *)malloc( - (unsigned long)pivs[nru-1][LENGTH] * sizeof(mpz_t)); - for (i = 0 ; i < pivs[nru-1][LENGTH]; ++i) { - mpz_init_set(mat->cf_ab_qq[nru-1][i], bs->cf_qq[pivs[nru-1][COEFFS]][i]); + (unsigned long)piv[LENGTH] * sizeof(mpz_t)); + for (i = 0 ; i < piv[LENGTH]; ++i) { + mpz_init_set(mat->cf_ab_qq[nru-1][i], bs->cf_qq[piv[COEFFS]][i]); } - pivs[nru-1][COEFFS] = nru-1; + piv[COEFFS] = nru-1; for (i = 0; i < nru-1; ++i) { k = nru-2-i; for (j = 0; j < ncols; ++j) { mpz_set_si(dr[j], 0); } - cfs = bs->cf_qq[pivs[k][COEFFS]]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + cfs = bs->cf_qq[piv[COEFFS]]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { mpz_set(dr[ds[j]], cfs[j]); @@ -402,14 +406,15 @@ static void exact_sparse_reduced_echelon_form_ab_first_qq( mpz_set(dr[ds[j+2]], cfs[j+2]); mpz_set(dr[ds[j+3]], cfs[j+3]); } - free(pivs[k]); + free(piv); cfs = NULL;; - pivs[k] = NULL; - pivs[k] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + piv = reduce_dense_row_by_known_pivots_sparse_only_ab_qq( dr, mat, pivs, sc, k); + atomic_store_explicit(&pivs[k], piv, memory_order_relaxed); remove_content_of_sparse_matrix_row_qq( - mat->cf_ab_qq[pivs[k][COEFFS]], pivs[k][PRELOOP], pivs[k][LENGTH]); + mat->cf_ab_qq[piv[COEFFS]], piv[PRELOOP], piv[LENGTH]); } const len_t drlen = st->nthrds * ncols; @@ -487,7 +492,10 @@ static void exact_sparse_reduced_echelon_form_ab_first_qq( remove_content_of_sparse_matrix_row_qq( mat->cf_qq[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH]); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_qq[npiv[COEFFS]]; } while (k == 0); cfs = NULL; @@ -495,13 +503,14 @@ static void exact_sparse_reduced_echelon_form_ab_first_qq( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - for (j = 0; j < pivs[i][LENGTH]; ++j) { - mpz_clear(mat->cf_ab_qq[pivs[i][COEFFS]][j]); + piv = atomic_load_explicit(&pivs[i], memory_order_relaxed); + for (j = 0; j < piv[LENGTH]; ++j) { + mpz_clear(mat->cf_ab_qq[piv[COEFFS]][j]); } - free(mat->cf_ab_qq[pivs[i][COEFFS]]); - mat->cf_ab_qq[pivs[i][COEFFS]] = NULL; - free(pivs[i]); - pivs[i] = NULL; + free(mat->cf_ab_qq[piv[COEFFS]]); + mat->cf_ab_qq[piv[COEFFS]] = NULL; + free(piv); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -515,15 +524,16 @@ static void exact_sparse_reduced_echelon_form_ab_first_qq( /* interreduce new pivots */ for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { for (j = 0; j < ncols; ++j) { mpz_set_si(dr[j], 0); } - cfs = mat->cf_qq[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_qq[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { mpz_swap(dr[ds[j]], cfs[j]); @@ -539,12 +549,13 @@ static void exact_sparse_reduced_echelon_form_ab_first_qq( mpz_swap(dr[ds[j+3]], cfs[j+3]); mpz_clear(cfs[j+3]); } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_qq( dr, mat, bs, pivs, sc, cf_array_pos); + atomic_store_explicit(&pivs[k], mat->tr[npivs], memory_order_relaxed); remove_content_of_sparse_matrix_row_qq( mat->cf_qq[mat->tr[npivs][COEFFS]], mat->tr[npivs][PRELOOP], @@ -579,8 +590,7 @@ static void exact_sparse_reduced_echelon_form_qq( const len_t ncl = mat->ncl; /* we fill in all known lead terms in pivs */ - hm_t **pivs = (hm_t **)calloc((unsigned long)ncols, sizeof(hm_t *)); - memcpy(pivs, mat->rr, (unsigned long)mat->nru * sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, mat->rr, mat->nru); /* unkown pivot rows we have to reduce with the known pivots first */ hm_t **upivs = mat->tr; @@ -661,7 +671,10 @@ static void exact_sparse_reduced_echelon_form_qq( remove_content_of_sparse_matrix_row_qq( mat->cf_qq[npiv[COEFFS]], npiv[PRELOOP], npiv[LENGTH]); } - k = __sync_bool_compare_and_swap(&pivs[npiv[OFFSET]], NULL, npiv); + hm_t *expected = NULL; + k = atomic_compare_exchange_strong_explicit( + &pivs[npiv[OFFSET]], &expected, npiv, + memory_order_release, memory_order_relaxed); cfs = mat->cf_qq[npiv[COEFFS]]; } while (k == 0); cfs = NULL; @@ -669,8 +682,8 @@ static void exact_sparse_reduced_echelon_form_qq( /* we do not need the old pivots anymore */ for (i = 0; i < ncl; ++i) { - free(pivs[i]); - pivs[i] = NULL; + free(atomic_load_explicit(&pivs[i], memory_order_relaxed)); + atomic_store_explicit(&pivs[i], NULL, memory_order_relaxed); } len_t npivs = 0; /* number of new pivots */ @@ -686,15 +699,16 @@ static void exact_sparse_reduced_echelon_form_qq( hm_t cf_array_pos; for (i = 0; i < ncr; ++i) { k = ncols-1-i; - if (pivs[k]) { + hm_t *piv = atomic_load_explicit(&pivs[k], memory_order_relaxed); + if (piv) { for (j = 0; j < ncols; ++j) { mpz_set_si(dr[j], 0); } - cfs = mat->cf_qq[pivs[k][COEFFS]]; - cf_array_pos = pivs[k][COEFFS]; - const len_t os = pivs[k][PRELOOP]; - const len_t len = pivs[k][LENGTH]; - const hm_t * const ds = pivs[k] + OFFSET; + cfs = mat->cf_qq[piv[COEFFS]]; + cf_array_pos = piv[COEFFS]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { mpz_swap(dr[ds[j]], cfs[j]); @@ -710,12 +724,13 @@ static void exact_sparse_reduced_echelon_form_qq( mpz_swap(dr[ds[j+3]], cfs[j+3]); mpz_clear(cfs[j+3]); } - free(pivs[k]); + free(piv); free(cfs); - pivs[k] = NULL; - pivs[k] = mat->tr[npivs] = + atomic_store_explicit(&pivs[k], NULL, memory_order_relaxed); + mat->tr[npivs] = reduce_dense_row_by_known_pivots_sparse_qq( dr, mat, bs, pivs, sc, cf_array_pos); + atomic_store_explicit(&pivs[k], mat->tr[npivs], memory_order_relaxed); remove_content_of_sparse_matrix_row_qq( mat->cf_qq[mat->tr[npivs][COEFFS]], mat->tr[npivs][PRELOOP], @@ -819,12 +834,13 @@ static void interreduce_matrix_rows_qq( mat->cf_qq = realloc(mat->cf_qq, (unsigned long)ncols * sizeof(mpz_t *)); memset(mat->cf_qq, 0, (unsigned long)ncols * sizeof(mpz_t *)); - hm_t **pivs = (hm_t **)calloc((unsigned long)ncols, sizeof(hm_t *)); + _Atomic(hm_t *) *pivs = allocate_atomic_pivots(ncols, NULL, 0); /* copy coefficient arrays from basis in matrix, maybe * several rows need the same coefficient arrays, but we * cannot share them here. */ for (i = 0; i < nrows; ++i) { - pivs[mat->rr[i][OFFSET]] = mat->rr[i]; + atomic_store_explicit(&pivs[mat->rr[i][OFFSET]], mat->rr[i], + memory_order_relaxed); } mpz_t *dr = (mpz_t *)malloc((unsigned long)ncols * sizeof(mpz_t)); @@ -838,14 +854,15 @@ static void interreduce_matrix_rows_qq( k = nrows - 1; for (i = 0; i < ncols; ++i) { l = ncols-1-i; - if (pivs[l] != NULL) { + hm_t *piv = atomic_load_explicit(&pivs[l], memory_order_relaxed); + if (piv) { for (j = 0; j < ncols; ++j) { mpz_set_si(dr[j], 0); } - cfs = bs->cf_qq[pivs[l][COEFFS]]; - const len_t os = pivs[l][PRELOOP]; - const len_t len = pivs[l][LENGTH]; - const hm_t * const ds = pivs[l] + OFFSET; + cfs = bs->cf_qq[piv[COEFFS]]; + const len_t os = piv[PRELOOP]; + const len_t len = piv[LENGTH]; + const hm_t * const ds = piv + OFFSET; sc = ds[0]; for (j = 0; j < os; ++j) { mpz_swap(dr[ds[j]], cfs[j]); @@ -856,11 +873,12 @@ static void interreduce_matrix_rows_qq( mpz_swap(dr[ds[j+2]], cfs[j+2]); mpz_swap(dr[ds[j+3]], cfs[j+3]); } - free(pivs[l]); - pivs[l] = NULL; - pivs[l] = mat->tr[k--] = + free(piv); + atomic_store_explicit(&pivs[l], NULL, memory_order_relaxed); + mat->tr[k] = reduce_dense_row_by_known_pivots_sparse_qq( dr, mat, bs, pivs, sc, l); + atomic_store_explicit(&pivs[l], mat->tr[k--], memory_order_relaxed); } } if (free_basis != 0) {