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) {