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

[fredrik-johansson]

We (using Claude Fable 5 and Claude Opus 5) implement flint_sgemm and flint_dgemm and make algorithms like nmod_mat_mul_blas available unconditionally, obsoleting non-BLAS tuning values.

This should make --with-blas unnecessary (actually, undesirable) for most users. On my Zen3 machine with AVX2, flint_sgemm and flint_dgemm match the speed of their OpenBLAS counterparts to within a few % both single-threaded and multi-threaded. Claude also claims that the speed matches OpenBLAS in its AVX512-enabled virtual machine, but someone should verify that, and also on NEON.

For the AVX512 variant, this uses the z extensions of machine_vectors.h proposed by @user202729 in #2715. This PR also adds a plain-C fallback for machine_vectors.h. The fallbacks don't cover enough ground to build fft_small just yet, but they allow flint_sgemm and flint_dgemm to work in any environment. On my machine, forcing the plain-C machine_vectors.h fallback and letting GCC auto-vectorize results in just ~2x worse performance than the specialized version using AVX2 intrinsics.

As noted above, linking to OpenBLAS is undesirable as it uses an OpenMP thread pool rather than FLINT's thread pool so that flint_set_num_threads does not have the intended effect (and parallel FLINT / BLAS work would compete for resources). At least on my machine, the OpenMP thread pool also adds a ridiculous ~1 second of startup CPU time to FLINT programs that link in OpenBLAS whether BLAS is actually used or not.

If FLINT is compiled with --with-blas, use of BLAS in flint_sgemm/flint_dgemm can be toggled with a switch flint_gemm_use_blas.

Some timings comparisons between flint_sgemm/dgemm ("sgemm/dgemm") and OpenBLAS ("blas_s/blas_d"), excerpted from build/machine_vectors/profile/p-gemm:

  128 x   128 x   128  [pow2]
  threads        sgemm   GFLOP/s      blas_s   GFLOP/s   ratio       dgemm   GFLOP/s      blas_d   GFLOP/s   ratio
  1             0.0426     98.48      0.0472     88.88    1.11      0.0814     51.53      0.0840     49.92    1.03
  2             0.0339    123.55      0.0295    142.22    0.87      0.0582     72.02      0.0542     77.32    0.93
  4             0.0212    197.61      0.0191    219.84    0.90      0.0344    121.85      0.0338    123.96    0.98
  8             0.0231    181.69      0.0159    264.55    0.69      0.0366    114.72      0.0232    180.66    0.63

  512 x   512 x   512  [pow2]
  threads        sgemm   GFLOP/s      blas_s   GFLOP/s   ratio       dgemm   GFLOP/s      blas_d   GFLOP/s   ratio
  1             2.4647    108.91      2.4674    108.79    1.00      4.9720     53.99      5.1853     51.77    1.04
  2             1.3332    201.35      1.3455    199.51    1.01      2.6928     99.69      2.8227     95.10    1.05
  4             0.8283    324.09      0.7242    370.68    0.87      1.6634    161.38      1.5701    170.96    0.94
  8             0.5051    531.43      0.5085    527.85    1.01      1.0017    267.97      1.0639    252.32    1.06

 1024 x  1024 x  1024  [pow2]
  threads        sgemm   GFLOP/s      blas_s   GFLOP/s   ratio       dgemm   GFLOP/s      blas_d   GFLOP/s   ratio
  1            20.4587    104.97     19.3027    111.25    0.94     41.1767     52.15     41.3800     51.90    1.00
  2            11.2055    191.65     10.4636    205.23    0.93     23.0270     93.26     22.3943     95.89    0.97
  4             6.2173    345.40      5.4953    390.78    0.88     12.8162    167.56     12.0377    178.40    0.94
  8             4.7759    449.65      3.5060    612.52    0.73      8.3332    257.70      7.6837    279.48    0.92

 2048 x  2048 x  2048  [pow2]
  threads        sgemm   GFLOP/s      blas_s   GFLOP/s   ratio       dgemm   GFLOP/s      blas_d   GFLOP/s   ratio
  1           166.2027    103.37    151.1250    113.68    0.91    330.6360     51.96    326.5663     52.61    0.99
  2            88.6977    193.69     81.8490    209.90    0.92    180.5517     95.15    176.6093     97.28    0.98
  4            52.4010    327.85     43.9717    390.70    0.84    108.1677    158.83     96.8657    177.36    0.90
  8            31.6117    543.47     28.2883    607.31    0.89     63.8087    269.24     62.1223    276.55    0.97

 4096 x  4096 x  4096  [pow2]
  threads        sgemm   GFLOP/s      blas_s   GFLOP/s   ratio       dgemm   GFLOP/s      blas_d   GFLOP/s   ratio
  1          1306.9447    105.16   1184.6850    116.01    0.91   2568.5930     53.51   2583.5690     53.20    1.01
  2           698.3930    196.79    668.1973    205.69    0.96   1424.2030     96.50   1455.6393     94.42    1.02
  4           423.1823    324.77    379.4470    362.21    0.90    863.0640    159.25    855.1230    160.72    0.99
  8           277.0863    496.01    237.0113    579.88    0.86    577.4873    237.99    516.7377    265.97    0.89

Fortunately, it isn't that hard to match BLAS these days: this is a very standard design with a tiling microkernel.

The most complex part of gemm code is the multithreading. There is a complex version using synchronization, invented by Claude, and a simpler version that just splits the work into independent submatrices. The complex version seems to work and is slightly more efficient, but it's only enabled by default on x86-64 for now pending stress-testing in other environments.

[vneiger]

Nice!!

If this can help, here is the output of build/machine_vectors/profile/p-gemm on a zen4 laptop (supports avx512, although with the same throughput as avx2 for many operations), and on an icelake server (supports avx512). Comparison is with BLIS (AOCL-BLAS for zen4, standard BLIS for icelake), and also with openblas for zen4.

The comparison looks very good overall. There is still some gap for the thin cases, e.g. matrix-vector products, or n x n x 8 products (if I remember well, BLIS has specific optimized code for tall/skinny matrices).

gemm_zen4_icelake.zip

It seems that some of the shortest measurements (small sizes with 1 thread, or smallish sizes with a few threads) should be taken with a grain of salt.

[fredrik-johansson]

Agreed, thin cases and matrix-vector products would need dedicated code.

Your Icelake timings look fine for a first iteration of the code. Interestingly sgemm stays near 1.0 up to 2048x2048, only degrading slightly at 4096x4096, while dgemm is closer to 2/3 throughout. Hopefully this gap can be closed by just changing some of the parameters.

[vneiger]

Timings on Apple M4 (with code before the last commit). Nothing new for dgemm: this PR's code and OpenBLAS are mostly on par, except for some "thin" cases where the gap remains quite reasonable. Something new for sgemm: OpenBLAS is often much faster (more than 5x or even 10x), I have not checked carefully but I guess this is due ARMv9 scalable matrix extension, supported by M4 and recently added in OpenBLAS.

mac-m4-openblas.txt

[fredrik-johansson]

Thanks, that looks very good actually. We could open a followup issue about using the ARMv9 scalable matrix extensions.