Fallbacks for BLAS sgemm/dgemm and (parts of) machine_vectors

CMakeLists.txt

@@ -206,6 +206,7 @@ endif()
 set(_BUILD_DIRS
     generic_files
     thread_pool                     thread_support
+    machine_vectors
 
     ulong_extras
     long_extras

Makefile.in

@@ -168,6 +168,7 @@ HEADLESS_DIRS := generic_files
 
 HEADER_DIRS	:=	\
 		thread_pool			thread_support					\
+		machine_vectors									\
 															\
 		ulong_extras		long_extras			perm		\
 		double_extras		d_vec				d_mat		\

doc/source/fmpz_mat.rst

@@ -611,9 +611,12 @@ Matrix multiplication
 
 .. function:: int fmpz_mat_mul_blas(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B)
 
-    Tries to set `C = AB` using BLAS and returns `1` for success and `0` for failure.
+    Tries to set `C = AB` by multimodular reduction to floating-point
+    matrix multiplication (:func:`flint_dgemm`, which uses BLAS if FLINT
+    was built with BLAS support and FLINT's own kernels otherwise), and
+    returns `1` for success and `0` for failure.
     Dimensions must be compatible for matrix multiplication. No aliasing is allowed.
-    This function currently will fail if the matrices are empty, their dimensions are too large, or their max bits size is over one million bits.
+    This function currently will fail if the matrices are empty, their dimensions are too large, their max bits size is over one million bits, or FLINT is built with a 32-bit word size.
 
 .. function:: void fmpz_mat_mul_fft(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B)
 

doc/source/machine_vectors.rst

@@ -3,8 +3,46 @@
 **machine_vectors.h** -- SIMD-accelerated operations on fixed-length vectors
 ===============================================================================
 
-This module currently requires building FLINT with support for
-AVX2 or NEON instructions.
+Vector types and operations mapping onto the target's SIMD instructions,
+together with the ``flint_sgemm`` and ``flint_dgemm`` matrix
+multiplication kernels built on top of them.
+
+Backends are selected automatically: AVX2 (and AVX512 for the ``vec8dz``
+family) on x86, NEON on ARM, and otherwise a generic backend, using GNU
+vector extensions where the compiler supports them and plain ISO C
+structs elsewhere. The AVX2 and NEON backends require a 64-bit word
+size, since their integer vectors have :type:`ulong` lanes; a 32-bit
+build uses the generic backends, which provide only the floating-point
+types. The generic backends implement only the subset of the
+interface required by ``flint_sgemm``/``flint_dgemm``; the full
+interface, as used by ``fft_small``, still requires AVX2 or NEON.
+
+For the vector operations to use the target's instructions, FLINT must
+be built with appropriate compiler flags. ``configure`` chooses these
+from the detected CPU; ``--enable-avx2`` and ``--enable-avx512`` force
+them on.
+
+Defining ``FLINT_MACHINE_VECTORS_FORCE_GENERIC`` before including this
+header selects the generic backend even on a target with AVX2 or NEON,
+and ``FLINT_MACHINE_VECTORS_STRICT_C`` additionally selects the ISO C
+tier over GNU vector extensions. These are intended for testing and
+profiling the portable code paths.
+
+The strict ISO C tier has no way to express a vector operation, so
+whether its operations become SIMD instructions is entirely up to the
+compiler's SLP vectorizer; its performance therefore varies between
+compilers and compiler versions, unlike the other backends. It exists
+so that the header works on compilers without GNU vector extensions,
+and is not the fallback used on GCC or clang.
+
+The generic backends express a fused multiply-add as ``a * b + c``,
+which a compiler fuses into an FMA instruction only when floating-point
+contraction is enabled. GCC in a strict ISO mode, which is how FLINT is
+built, does not contract by default; code using these operations in a
+performance-critical loop should request contraction, as
+``machine_vectors/gemm.c`` does with ``#pragma GCC optimize
+("fp-contract=fast")``. The AVX2, AVX512 and NEON backends use fused
+intrinsics and are unaffected.
 
 Some functions may require that vectors are aligned in memory.
 
@@ -25,6 +63,27 @@ Types
 
     Vector with 1, 2, 4, or 8 ``double`` entries.
 
+.. type:: vec1f
+          vec4f
+          vec8f
+          vec16f
+
+    Vector with 1, 4, 8, or 16 ``float`` entries.
+
+.. type:: vec8dz
+          vec16dz
+          vec16fz
+          vec32fz
+          vec8nz
+
+    Vectors backed by AVX512 registers, available only when building
+    with AVX512F support: 8 or 16 ``double`` entries, 16 or 32 ``float``
+    entries, and 8 :type:`ulong` entries respectively. The ``z`` suffix
+    distinguishes these from the equally-named types built from pairs of
+    narrower registers; ``vec8d``, for instance, remains a pair of AVX2
+    registers, which gives more instruction level parallelism in
+    existing code.
+
 Printing
 -------------------------------------------------------------------------------
 
@@ -282,3 +341,55 @@ Other assumptions are not yet documented.
               vec8n vec8n_addmod_limited(vec8n a, vec8n b, vec8n n)
 
     Return `a + b \bmod n` in `[0,n)`, assuming that `n < 2^{63}`.
+
+Matrix multiplication
+-------------------------------------------------------------------------------
+
+These functions compute a matrix product in single or double precision.
+They are always available: when FLINT is built with BLAS, the default is
+to call it, and otherwise FLINT's own kernels are used. The intended use
+is as a building block for exact linear algebra over `\mathbb{Z}` and
+`\mathbb{Z}/n\mathbb{Z}`, for example in :func:`nmod_mat_mul_blas` and
+:func:`fmpz_mat_mul_blas`.
+
+All of these functions compute `C = AB` for row-major matrices, where
+*C* is *m* by *n*, *A* is *m* by *k* and *B* is *k* by *n*, with
+*ldc*, *lda* and *ldb* the respective leading dimensions (the number of
+entries between the start of consecutive rows, which must be at least
+the number of columns). No transposition or accumulation is performed:
+the previous contents of *C* are overwritten, and `k = 0` sets *C* to
+zero. This is equivalent to ``cblas_sgemm`` or ``cblas_dgemm`` called
+with ``CblasRowMajor``, ``CblasNoTrans``, ``CblasNoTrans``, ``alpha``
+equal to 1 and ``beta`` equal to 0. Aliasing of *C* with *A* or *B* is
+not allowed.
+
+The FLINT kernels are multithreaded internally according to
+:func:`flint_get_num_threads`, using FLINT's thread pool. They handle
+arbitrary dimensions, including thin and unbalanced shapes, without
+requiring any padding or alignment of the input.
+
+.. function:: void flint_sgemm(slong m, slong k, slong n, const float * A, slong lda, const float * B, slong ldb, float * C, slong ldc)
+              void flint_dgemm(slong m, slong k, slong n, const double * A, slong lda, const double * B, slong ldb, double * C, slong ldc)
+
+    Sets `C = AB`, calling either the BLAS or the FLINT implementation
+    according to :var:`flint_gemm_use_blas`.
+
+.. function:: void flint_sgemm_blas(slong m, slong k, slong n, const float * A, slong lda, const float * B, slong ldb, float * C, slong ldc)
+              void flint_dgemm_blas(slong m, slong k, slong n, const double * A, slong lda, const double * B, slong ldb, double * C, slong ldc)
+
+    Sets `C = AB` using ``cblas_sgemm`` or ``cblas_dgemm``. These raise
+    an exception if FLINT was built without BLAS support.
+
+.. function:: void flint_sgemm_fallback(slong m, slong k, slong n, const float * A, slong lda, const float * B, slong ldb, float * C, slong ldc)
+              void flint_dgemm_fallback(slong m, slong k, slong n, const double * A, slong lda, const double * B, slong ldb, double * C, slong ldc)
+
+    Sets `C = AB` using FLINT's own kernels, which are always available.
+
+.. var:: int flint_gemm_use_blas
+
+    Selects the implementation used by :func:`flint_sgemm` and
+    :func:`flint_dgemm`. It is initialized to 1 if FLINT was built with
+    BLAS support and 0 otherwise, and may be set to either value at
+    runtime, for example to compare the two implementations. Setting it
+    to 1 in a build without BLAS support will result in an exception
+    when a gemm is attempted.

doc/source/nmod_mat.rst

@@ -382,7 +382,14 @@ Matrix multiplication
 
 .. function:: int nmod_mat_mul_blas(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B)
 
-    Tries to set `C = AB` using BLAS and returns `1` for success and `0` for failure. Dimensions must be compatible for matrix multiplication.
+    Tries to set `C = AB` by lifting to floating-point matrix
+    multiplication (:func:`flint_sgemm` or :func:`flint_dgemm`, which use
+    BLAS if FLINT was built with BLAS support and FLINT's own kernels
+    otherwise), with multimodular reduction and CRT when the entries do
+    not fit directly. Returns `1` for success and `0` for failure;
+    failure occurs when the dimensions or the modulus are too large, or
+    when FLINT is built with a 32-bit word size.
+    Dimensions must be compatible for matrix multiplication.
 
 .. function:: void nmod_mat_addmul(nmod_mat_t D, const nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B)
 

src/fmpz_mat/mul.c

@@ -192,7 +192,7 @@ fmpz_mat_mul(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B)
             return;
     }
 
-#if FLINT_USES_BLAS && FLINT_BITS == 64
+#if FLINT_BITS == 64
     if (dim > 50)
     {
         if (cbits <= 53)

src/fmpz_mat/mul_blas.c

@@ -14,11 +14,11 @@
 
 /* todo: squaring optimizations */
 
-#if FLINT_USES_BLAS && FLINT_BITS == 64
+#if FLINT_BITS == 64
 
 #include <stdint.h>
 #include <limits.h>
-#include <cblas.h>
+#include "machine_vectors.h"
 #include "nmod.h"
 #include "fmpz.h"
 #include "thread_pool.h"
@@ -141,8 +141,7 @@ static int _fmpz_mat_mul_blas_direct(
         flint_free(args);
     }
 
-    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
-                                       m, n, k, 1.0, dA, k, dB, n, 0.0, dC, n);
+    flint_dgemm(m, k, n, dA, k, dB, n, dC, n);
 
     for (i = 0; i < m; i++)
         for (j = 0; j < n; j++)
@@ -355,7 +354,7 @@ static void _fromd_worker(void * arg_ptr)
         {
             ulong r;
             slong a = (slong) dC[i*n + j];
-            ulong b = (a < 0) ? a + shift : a;
+            ulong b = (a < 0) ? (ulong) a + shift : (ulong) a;
             NMOD_RED(r, b, mod);
             bigC[n*(num_primes*i + l) + j] = r;
         }
@@ -475,7 +474,7 @@ int _fmpz_mat_mul_blas(
     slong num_primes;
     fmpz_comb_t comb;
     thread_pool_handle * handles;
-    slong num_workers;
+    slong num_workers, max_workers;
     _worker_arg * args;
 
     FLINT_ASSERT(sign == 0 || sign == 1);
@@ -509,6 +508,7 @@ int _fmpz_mat_mul_blas(
     dC = (double *) flint_calloc(m*n, sizeof(double));
 
     num_workers = flint_request_threads(&handles, INT_MAX);
+    max_workers = num_workers;
 
     args = FLINT_ARRAY_ALLOC(num_workers + 1, _worker_arg);
     for (start = 0, i = 0; i <= num_workers; start = stop, i++)
@@ -567,8 +567,40 @@ int _fmpz_mat_mul_blas(
         for (i = 0; i < num_workers; i++)
             thread_pool_wait(global_thread_pool, handles[i]);
 
-        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
-                                       m, n, k, 1.0, dA, k, dB, n, 0.0, dC, n);
+        /*
+            The fallback gemm threads through FLINT's pool itself, so
+            the workers held for the conversions must be returned
+            before the call or it finds the pool empty and runs on one
+            thread. (An external BLAS schedules its own threads and
+            does not care.) The re-request is capped by the first
+            grant, so the args array stays large enough; if fewer
+            workers come back, the work is redistributed below.
+        */
+        {
+            slong prev_workers = num_workers;
+
+            flint_give_back_threads(handles, num_workers);
+
+            flint_dgemm(m, k, n, dA, k, dB, n, dC, n);
+
+            num_workers = flint_request_threads(&handles, max_workers + 1);
+
+            if (num_workers != prev_workers)
+            {
+                for (start = 0, i = 0; i <= num_workers; start = stop, i++)
+                {
+                    args[i].l = l;
+                    args[i].prime = primes[l];
+                    args[i].Cstartrow = ((i + 0)*m)/(num_workers + 1);
+                    args[i].Cstoprow  = ((i + 1)*m)/(num_workers + 1);
+                    stop = _thread_pool_find_work_2(m, k, k, n,
+                                                    i + 1, num_workers + 1);
+                    _thread_pool_distribute_work_2(start, stop,
+                                  &args[i].Astartrow, &args[i].Astoprow, m,
+                                  &args[i].Bstartrow, &args[i].Bstoprow, k);
+                }
+            }
+        }
 
         for (i = 0; i < num_workers; i++)
             thread_pool_wake(global_thread_pool, handles[i], 0, _fromd_worker, &args[i]);

src/fmpz_mat/profile/p-mul_blas_v_mul.c

@@ -11,8 +11,6 @@
 
 #include "flint.h"
 
-#if FLINT_USES_BLAS
-#include <cblas.h>
 #include "longlong.h"  // for FLINT_BIT_COUNT
 #include "fmpz_mat.h"
 #include "profiler.h"
@@ -96,10 +94,3 @@ int main(void)
     FLINT_TEST_CLEAR(state);
     return 0;
 }
-
-#else
-int main(void)
-{
-    return 0;
-}
-#endif

src/fmpz_mat/test/t-mul_blas.c

@@ -54,10 +54,10 @@ TEST_FUNCTION_START(fmpz_mat_mul_blas, state)
                 flint_abort();
             }
         }
-#if FLINT_USES_BLAS && FLINT_BITS == 64
+#if FLINT_BITS == 64
         else
         {
-            flint_printf("FAIL: blas should have worked\n");
+            flint_printf("FAIL: mul_blas should have worked\n");
             fflush(stdout);
             flint_abort();
         }
@@ -102,10 +102,10 @@ TEST_FUNCTION_START(fmpz_mat_mul_blas, state)
                 flint_abort();
             }
         }
-#if FLINT_USES_BLAS && FLINT_BITS == 64
+#if FLINT_BITS == 64
         else
         {
-            flint_printf("FAIL: blas should have worked\n");
+            flint_printf("FAIL: mul_blas should have worked\n");
             fflush(stdout);
             flint_abort();
         }