Replace __sync builtins with C11 atomics in linear algebra

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 <stdatomic.h>
+#include <stddef.h>
+]], [[
+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 (<stdatomic.h>) 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

src/neogb/data.c

@@ -19,8 +19,38 @@
  * Mohab Safey El Din */
 
 
+#include <stdatomic.h>
+
 #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,15 +184,15 @@ 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
         );
 
 
 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
         );

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

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
         )