From cb0727edd1a34fda7da0821769644ea42e386184 Mon Sep 17 00:00:00 2001 From: Weijia Wang <9713184+wegank@users.noreply.github.com> Date: Thu, 1 Oct 2026 23:36:18 +0200 Subject: [PATCH] Detect AVX2 and AVX-512 support at runtime So far, configure inspected the CPU of the build machine and compiled all of msolve with the SIMD extensions found there. Binaries built on a recent machine thus crashed with an illegal instruction on older x86-64 CPUs, and Linux distributions, which must target the baseline architecture, could not ship the vectorized code at all (see #71). Now, the AVX2 and AVX-512 kernels are always compiled when the compiler supports them, each marked with a function-level target attribute, while the remaining code is compiled for the baseline architecture. The kernel to use is chosen at runtime via __builtin_cpu_supports. The new header src/msolve/cpu_features.h provides the target attributes, the runtime checks and the SIMD intrinsics includes; with e.g. -march=native the checks are resolved at compile time. In neogb, the inline SIMD branches of the row reduction routines are moved into separate scalar/NEON/AVX2/AVX-512 kernels, without changing which instruction sets each routine uses. The fglm kernels already were separate functions and only get the target attributes. configure now only checks whether the compiler can build these kernels and defines HAVE_AVX2_KERNELS and HAVE_AVX512_KERNELS accordingly. AX_EXT and the m4 macros only used by it are removed. Co-Authored-By: Claude Opus 5.5 --- INSTALL | 9 + Makefile.am | 2 +- configure.ac | 48 +- m4/ax_check_compile_flag.m4 | 53 -- m4/ax_cpu_vendor.m4 | 64 -- m4/ax_ext.m4 | 328 --------- m4/ax_gcc_x86_avx_xgetbv.m4 | 79 --- m4/ax_gcc_x86_cpuid.m4 | 89 --- src/fglm/Makefile.am | 2 +- src/fglm/fglm_core.c | 67 +- src/fglm/linalg-fglm.c | 58 +- src/fglm/matrix-mult.c | 14 +- src/msolve/Makefile.am | 3 +- src/msolve/cpu_features.h | 74 +++ src/neogb/Makefile.am | 2 +- src/neogb/la_ff_16.c | 550 ++++++++------- src/neogb/la_ff_32.c | 1253 +++++++++++++++-------------------- src/neogb/la_ff_8.c | 842 ++++++++++++----------- src/neogb/symbol.c | 4 - src/usolve/Makefile.am | 2 +- 20 files changed, 1513 insertions(+), 2030 deletions(-) delete mode 100644 m4/ax_check_compile_flag.m4 delete mode 100644 m4/ax_cpu_vendor.m4 delete mode 100644 m4/ax_ext.m4 delete mode 100644 m4/ax_gcc_x86_avx_xgetbv.m4 delete mode 100644 m4/ax_gcc_x86_cpuid.m4 create mode 100644 src/msolve/cpu_features.h diff --git a/INSTALL b/INSTALL index b8a96dc2..e6eace18 100644 --- a/INSTALL +++ b/INSTALL @@ -30,6 +30,15 @@ via github releases) 4. Run `make install` in order to globally install the library and the binary of msolve. +Note on SIMD instructions +========================= + +On x86-64, msolve contains AVX2 and AVX-512 variants of its linear algebra +kernels and chooses the best one supported by the CPU at runtime. Thus the +binaries run on any x86-64 CPU, independently of the CPU of the build machine. +If msolve is only used on the machine it is built on, adding `-march=native` to +CFLAGS lets the compiler optimize the remaining code for this CPU as well. + Note to macOS users =================== diff --git a/Makefile.am b/Makefile.am index cb3948a8..4f64cdcd 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,7 +1,7 @@ AUTOMAKE_OPTIONS = color-tests SUBDIRS = src/usolve src/fglm src/neogb src/msolve -AM_CFLAGS = $(SIMD_FLAGS) $(CPUEXT_FLAGS) $(OPENMP_CFLAGS) +AM_CFLAGS = $(OPENMP_CFLAGS) LDADD = src/neogb/libneogb.la src/fglm/libfglm.la src/usolve/libusolve.la bin_PROGRAMS = msolve diff --git a/configure.ac b/configure.ac index 2221b5a0..791da899 100644 --- a/configure.ac +++ b/configure.ac @@ -70,7 +70,53 @@ AC_TYPE_UINT8_T AX_COUNT_CPUS -AX_EXT + +# AVX2 and AVX-512 kernels are compiled whenever the compiler supports them +# and are selected at runtime depending on the CPU msolve runs on, see +# src/msolve/cpu_features.h. Thus the binaries do not depend on the CPU of +# the build machine. +case $host_cpu in + x86_64*|amd64*) + AC_CACHE_CHECK([whether $CC can build AVX2 kernels], + [msolve_cv_avx2_kernels], + [AC_LINK_IFELSE([AC_LANG_PROGRAM([[ +#include +__attribute__((target("avx2"))) +static long long f(long long x) +{ + __m256i v = _mm256_set1_epi64x(x); + v = _mm256_add_epi64(_mm256_mul_epu32(v, v), _mm256_cmpgt_epi64(v, v)); + return _mm256_extract_epi64(v, 0); +}]], [[return __builtin_cpu_supports("avx2") ? (int)f(3) : 0;]])], + [msolve_cv_avx2_kernels=yes], + [msolve_cv_avx2_kernels=no])]) + if test "x$msolve_cv_avx2_kernels" = xyes; then + AC_DEFINE([HAVE_AVX2_KERNELS], [1], + [Define to 1 if AVX2 kernels are built (selected at runtime).]) + fi + + AC_CACHE_CHECK([whether $CC can build AVX-512 kernels], + [msolve_cv_avx512_kernels], + [AC_LINK_IFELSE([AC_LANG_PROGRAM([[ +#include +__attribute__((target("avx512f,avx512bw"))) +static long long f(long long x) +{ + __m512i v = _mm512_set1_epi64(x); + __mmask8 m = _mm512_cmpgt_epi64_mask(v, _mm512_setzero_si512()); + v = _mm512_mask_add_epi64(v, m, v, _mm512_mullo_epi16(v, v)); + return _mm512_reduce_add_epi64(v); +}]], [[return __builtin_cpu_supports("avx512f") + && __builtin_cpu_supports("avx512bw") ? (int)f(3) : 0;]])], + [msolve_cv_avx512_kernels=yes], + [msolve_cv_avx512_kernels=no])]) + if test "x$msolve_cv_avx512_kernels" = xyes; then + AC_DEFINE([HAVE_AVX512_KERNELS], [1], + [Define to 1 if AVX-512 kernels are built (selected at runtime).]) + fi + ;; +esac + AX_GCC_BUILTIN([__builtin_constant_p]) AX_GCC_BUILTIN([__builtin_clzll]) AX_GCC_BUILTIN([__builtin_clzl]) diff --git a/m4/ax_check_compile_flag.m4 b/m4/ax_check_compile_flag.m4 deleted file mode 100644 index bd753b34..00000000 --- a/m4/ax_check_compile_flag.m4 +++ /dev/null @@ -1,53 +0,0 @@ -# =========================================================================== -# https://www.gnu.org/software/autoconf-archive/ax_check_compile_flag.html -# =========================================================================== -# -# SYNOPSIS -# -# AX_CHECK_COMPILE_FLAG(FLAG, [ACTION-SUCCESS], [ACTION-FAILURE], [EXTRA-FLAGS], [INPUT]) -# -# DESCRIPTION -# -# Check whether the given FLAG works with the current language's compiler -# or gives an error. (Warnings, however, are ignored) -# -# ACTION-SUCCESS/ACTION-FAILURE are shell commands to execute on -# success/failure. -# -# If EXTRA-FLAGS is defined, it is added to the current language's default -# flags (e.g. CFLAGS) when the check is done. The check is thus made with -# the flags: "CFLAGS EXTRA-FLAGS FLAG". This can for example be used to -# force the compiler to issue an error when a bad flag is given. -# -# INPUT gives an alternative input source to AC_COMPILE_IFELSE. -# -# NOTE: Implementation based on AX_CFLAGS_GCC_OPTION. Please keep this -# macro in sync with AX_CHECK_{PREPROC,LINK}_FLAG. -# -# LICENSE -# -# Copyright (c) 2008 Guido U. Draheim -# Copyright (c) 2011 Maarten Bosmans -# -# Copying and distribution of this file, with or without modification, are -# permitted in any medium without royalty provided the copyright notice -# and this notice are preserved. This file is offered as-is, without any -# warranty. - -#serial 6 - -AC_DEFUN([AX_CHECK_COMPILE_FLAG], -[AC_PREREQ(2.64)dnl for _AC_LANG_PREFIX and AS_VAR_IF -AS_VAR_PUSHDEF([CACHEVAR],[ax_cv_check_[]_AC_LANG_ABBREV[]flags_$4_$1])dnl -AC_CACHE_CHECK([whether _AC_LANG compiler accepts $1], CACHEVAR, [ - ax_check_save_flags=$[]_AC_LANG_PREFIX[]FLAGS - _AC_LANG_PREFIX[]FLAGS="$[]_AC_LANG_PREFIX[]FLAGS $4 $1" - AC_COMPILE_IFELSE([m4_default([$5],[AC_LANG_PROGRAM()])], - [AS_VAR_SET(CACHEVAR,[yes])], - [AS_VAR_SET(CACHEVAR,[no])]) - _AC_LANG_PREFIX[]FLAGS=$ax_check_save_flags]) -AS_VAR_IF(CACHEVAR,yes, - [m4_default([$2], :)], - [m4_default([$3], :)]) -AS_VAR_POPDEF([CACHEVAR])dnl -])dnl AX_CHECK_COMPILE_FLAGS diff --git a/m4/ax_cpu_vendor.m4 b/m4/ax_cpu_vendor.m4 deleted file mode 100644 index 79eeed34..00000000 --- a/m4/ax_cpu_vendor.m4 +++ /dev/null @@ -1,64 +0,0 @@ -# =========================================================================== -# http://www.gnu.org/software/autoconf-archive/ax_cpu_vendor.html -# =========================================================================== -# -# SYNOPSIS -# -# AX_CPU_VENDOR -# -# DESCRIPTION -# -# Find your CPU's vendor by requesting cpuid and define "ax_cpu_vendor" -# accordingly. This macro depends on AX_GCC_X86_CPUID. -# -# LICENSE -# -# Copyright (c) 2008 Christophe Tournayre -# -# Copying and distribution of this file, with or without modification, are -# permitted in any medium without royalty provided the copyright notice -# and this notice are preserved. This file is offered as-is, without any -# warranty. - -#serial 7 - -AC_DEFUN([AX_CPU_VENDOR], -[ - AC_REQUIRE([AX_GCC_X86_CPUID]) - AX_GCC_X86_CPUID(0x0) - - AC_CACHE_CHECK(for the processor vendor, ax_cpu_vendor, - [ - vendor=`echo $ax_cv_gcc_x86_cpuid_0x0 | cut -d ":" -f 2` - - case $vendor in - 756e6547*) - ax_cpu_vendor="Intel" - ;; - 68747541*) - ax_cpu_vendor="AMD" - ;; - 69727943*) - ax_cpu_vendor="Cyrix" - ;; - 746e6543*) - ax_cpu_vendor="IDT" - ;; - 646f6547*) - ax_cpu_vendor="Natsemi Geode" - ;; - 52697365*) - ax_cpu_vendor="Rise" - ;; - 65736952*) - ax_cpu_vendor="Rise" - ;; - 20536953*) - ax_cpu_vendor="SiS" - ;; - *) - ax_cpu_vendor="Unknown" - ;; - esac - ]) -]) diff --git a/m4/ax_ext.m4 b/m4/ax_ext.m4 deleted file mode 100644 index 449daf21..00000000 --- a/m4/ax_ext.m4 +++ /dev/null @@ -1,328 +0,0 @@ -# =========================================================================== -# https://www.gnu.org/software/autoconf-archive/ax_ext.html -# =========================================================================== -# -# SYNOPSIS -# -# AX_EXT -# -# DESCRIPTION -# -# Find supported SIMD extensions by requesting cpuid. When a SIMD -# extension is found, the -m"simdextensionname" is added to SIMD_FLAGS if -# compiler supports it. For example, if "sse2" is available then "-msse2" -# is added to SIMD_FLAGS. -# -# Find other supported CPU extensions by requesting cpuid. When a -# processor extension is found, the -m"extensionname" is added to -# CPUEXT_FLAGS if compiler supports it. For example, if "bmi2" is -# available then "-mbmi2" is added to CPUEXT_FLAGS. -# -# This macro calls: -# -# AC_SUBST(SIMD_FLAGS) -# AC_SUBST(CPUEXT_FLAGS) -# -# And defines: -# -# HAVE_RDRND / HAVE_BMI1 / HAVE_BMI2 / HAVE_ADX / HAVE_MPX -# HAVE_PREFETCHWT1 / HAVE_ABM / HAVE_MMX / HAVE_SSE / HAVE_SSE2 -# HAVE_SSE3 / HAVE_SSSE3 / HAVE_SSE4_1 / HAVE_SSE4_2 / HAVE_SSE4a -# HAVE_SHA / HAVE_AES / HAVE_AVX / HAVE_FMA3 / HAVE_FMA4 / HAVE_XOP -# HAVE_AVX2 / HAVE_AVX512_F / HAVE_AVX512_CD / HAVE_AVX512_PF -# HAVE_AVX512_ER / HAVE_AVX512_VL / HAVE_AVX512_BW / HAVE_AVX512_DQ -# HAVE_AVX512_IFMA / HAVE_AVX512_VBMI / HAVE_ALTIVEC / HAVE_VSX -# -# LICENSE -# -# Copyright (c) 2007 Christophe Tournayre -# Copyright (c) 2013,2015 Michael Petch -# Copyright (c) 2017 Rafael de Lucena Valle -# -# Copying and distribution of this file, with or without modification, are -# permitted in any medium without royalty provided the copyright notice -# and this notice are preserved. This file is offered as-is, without any -# warranty. - -#serial 18 - -AC_DEFUN([AX_EXT], -[ - AC_REQUIRE([AC_CANONICAL_HOST]) - AC_REQUIRE([AC_PROG_CC]) - - CPUEXT_FLAGS="" - SIMD_FLAGS="" - - case $host_cpu in - powerpc*) - AC_CACHE_CHECK([whether altivec is supported for old distros], [ax_cv_have_altivec_old_ext], - [ - if test `/usr/sbin/sysctl -a 2>/dev/null| grep -c hw.optional.altivec` != 0; then - if test `/usr/sbin/sysctl -n hw.optional.altivec` = 1; then - ax_cv_have_altivec_old_ext=yes - fi - fi - ]) - - if test "$ax_cv_have_altivec_old_ext" = yes; then - AC_DEFINE(HAVE_ALTIVEC,,[Support Altivec instructions]) - AX_CHECK_COMPILE_FLAG(-faltivec, SIMD_FLAGS="$SIMD_FLAGS -faltivec", []) - fi - - AC_CACHE_CHECK([whether altivec is supported], [ax_cv_have_altivec_ext], - [ - if test `LD_SHOW_AUXV=1 /bin/true 2>/dev/null|grep -c altivec` != 0; then - ax_cv_have_altivec_ext=yes - fi - ]) - - if test "$ax_cv_have_altivec_ext" = yes; then - AC_DEFINE(HAVE_ALTIVEC,,[Support Altivec instructions]) - AX_CHECK_COMPILE_FLAG(-maltivec, SIMD_FLAGS="$SIMD_FLAGS -maltivec", []) - fi - - AC_CACHE_CHECK([whether vsx is supported], [ax_cv_have_vsx_ext], - [ - if test `LD_SHOW_AUXV=1 /bin/true 2>/dev/null|grep -c vsx` != 0; then - ax_cv_have_vsx_ext=yes - fi - ]) - - if test "$ax_cv_have_vsx_ext" = yes; then - AC_DEFINE(HAVE_VSX,,[Support VSX instructions]) - AX_CHECK_COMPILE_FLAG(-mvsx, SIMD_FLAGS="$SIMD_FLAGS -mvsx", []) - fi - ;; - - i[[3456]]86*|x86_64*|amd64*) - - AC_REQUIRE([AX_GCC_X86_CPUID]) - AC_REQUIRE([AX_GCC_X86_CPUID_COUNT]) - AC_REQUIRE([AX_GCC_X86_AVX_XGETBV]) - - eax_cpuid0=0 - AX_GCC_X86_CPUID(0x00000000) - if test "$ax_cv_gcc_x86_cpuid_0x00000000" != "unknown"; - then - eax_cpuid0=`echo $ax_cv_gcc_x86_cpuid_0x00000000 | cut -d ":" -f 1` - fi - - eax_cpuid80000000=0 - AX_GCC_X86_CPUID(0x80000000) - if test "$ax_cv_gcc_x86_cpuid_0x80000000" != "unknown"; - then - eax_cpuid80000000=`echo $ax_cv_gcc_x86_cpuid_0x80000000 | cut -d ":" -f 1` - fi - - ecx_cpuid1=0 - edx_cpuid1=0 - if test "$((0x$eax_cpuid0))" -ge 1 ; then - AX_GCC_X86_CPUID(0x00000001) - if test "$ax_cv_gcc_x86_cpuid_0x00000001" != "unknown"; - then - ecx_cpuid1=`echo $ax_cv_gcc_x86_cpuid_0x00000001 | cut -d ":" -f 3` - edx_cpuid1=`echo $ax_cv_gcc_x86_cpuid_0x00000001 | cut -d ":" -f 4` - fi - fi - - ebx_cpuid7=0 - ecx_cpuid7=0 - if test "$((0x$eax_cpuid0))" -ge 7 ; then - AX_GCC_X86_CPUID_COUNT(0x00000007, 0x00) - if test "$ax_cv_gcc_x86_cpuid_0x00000007" != "unknown"; - then - ebx_cpuid7=`echo $ax_cv_gcc_x86_cpuid_0x00000007 | cut -d ":" -f 2` - ecx_cpuid7=`echo $ax_cv_gcc_x86_cpuid_0x00000007 | cut -d ":" -f 3` - fi - fi - - ecx_cpuid80000001=0 - edx_cpuid80000001=0 - if test "$((0x$eax_cpuid80000000))" -ge "$((0x80000001))" ; then - AX_GCC_X86_CPUID(0x80000001) - if test "$ax_cv_gcc_x86_cpuid_0x80000001" != "unknown"; - then - ecx_cpuid80000001=`echo $ax_cv_gcc_x86_cpuid_0x80000001 | cut -d ":" -f 3` - edx_cpuid80000001=`echo $ax_cv_gcc_x86_cpuid_0x80000001 | cut -d ":" -f 4` - fi - fi - - AC_CACHE_VAL([ax_cv_have_mmx_os_support_ext], - [ - ax_cv_have_mmx_os_support_ext=yes - ]) - - ax_cv_have_none_os_support_ext=yes - - AC_CACHE_VAL([ax_cv_have_sse_os_support_ext], - [ - ax_cv_have_sse_os_support_ext=no, - if test "$((0x$edx_cpuid1>>25&0x01))" = 1; then - AC_LANG_PUSH([C]) - AC_RUN_IFELSE([AC_LANG_SOURCE([[ -#include -#include - /* No way at ring1 to ring3 in protected mode to check the CR0 and CR4 - control registers directly. Execute an SSE instruction. - If it raises SIGILL then OS doesn't support SSE based instructions */ - void sig_handler(int signum){ exit(1); } - int main(){ - signal(SIGILL, sig_handler); - /* SSE instruction xorps %xmm0,%xmm0 */ - __asm__ __volatile__ (".byte 0x0f, 0x57, 0xc0"); - return 0; - }]])], - [ax_cv_have_sse_os_support_ext=yes], - [ax_cv_have_sse_os_support_ext=no], - [ax_cv_have_sse_os_support_ext=no]) - AC_LANG_POP([C]) - fi - ]) - - xgetbv_eax=0 - if test "$((0x$ecx_cpuid1>>28&0x01))" = 1; then - AX_GCC_X86_AVX_XGETBV(0x00000000) - - if test x"$ax_cv_gcc_x86_avx_xgetbv_0x00000000" != x"unknown"; then - xgetbv_eax=`echo $ax_cv_gcc_x86_avx_xgetbv_0x00000000 | cut -d ":" -f 1` - fi - - AC_CACHE_VAL([ax_cv_have_avx_os_support_ext], - [ - ax_cv_have_avx_os_support_ext=no - if test "$((0x$ecx_cpuid1>>27&0x01))" = 1; then - if test "$((0x$xgetbv_eax&0x6))" = 6; then - ax_cv_have_avx_os_support_ext=yes - fi - fi - ]) - fi - - AC_CACHE_VAL([ax_cv_have_avx512_os_support_ext], - [ - ax_cv_have_avx512_os_support_ext=no - if test "$ax_cv_have_avx_os_support_ext" = yes; then - if test "$((0x$xgetbv_eax&0xe6))" = "$((0xe6))"; then - ax_cv_have_avx512_os_support_ext=yes - fi - fi - ]) - - for ac_instr_info dnl - in "none;rdrnd;RDRND;ecx_cpuid1,30;-mrdrnd;HAVE_RDRND;CPUEXT_FLAGS" dnl - "none;bmi1;BMI1;ebx_cpuid7,3;-mbmi;HAVE_BMI1;CPUEXT_FLAGS" dnl - "none;bmi2;BMI2;ebx_cpuid7,8;-mbmi2;HAVE_BMI2;CPUEXT_FLAGS" dnl - "none;adx;ADX;ebx_cpuid7,19;-madx;HAVE_ADX;CPUEXT_FLAGS" dnl - dnl "none;mpx;MPX;ebx_cpuid7,14;-mmpx;HAVE_MPX;CPUEXT_FLAGS" dnl - "none;prefetchwt1;PREFETCHWT1;ecx_cpuid7,0;-mprefetchwt1;HAVE_PREFETCHWT1;CPUEXT_FLAGS" dnl - "none;abm;ABM;ecx_cpuid80000001,5;-mabm;HAVE_ABM;CPUEXT_FLAGS" dnl - "mmx;mmx;MMX;edx_cpuid1,23;-mmmx;HAVE_MMX;SIMD_FLAGS" dnl - "sse;sse;SSE;edx_cpuid1,25;-msse;HAVE_SSE;SIMD_FLAGS" dnl - "sse;sse2;SSE2;edx_cpuid1,26;-msse2;HAVE_SSE2;SIMD_FLAGS" dnl - "sse;sse3;SSE3;ecx_cpuid1,1;-msse3;HAVE_SSE3;SIMD_FLAGS" dnl - "sse;ssse3;SSSE3;ecx_cpuid1,9;-mssse3;HAVE_SSSE3;SIMD_FLAGS" dnl - "sse;sse41;SSE4.1;ecx_cpuid1,19;-msse4.1;HAVE_SSE4_1;SIMD_FLAGS" dnl - "sse;sse42;SSE4.2;ecx_cpuid1,20;-msse4.2;HAVE_SSE4_2;SIMD_FLAGS" dnl - "sse;sse4a;SSE4a;ecx_cpuid80000001,6;-msse4a;HAVE_SSE4a;SIMD_FLAGS" dnl - "sse;sha;SHA;ebx_cpuid7,29;-msha;HAVE_SHA;SIMD_FLAGS" dnl - "sse;aes;AES;ecx_cpuid1,25;-maes;HAVE_AES;SIMD_FLAGS" dnl - "avx;avx;AVX;ecx_cpuid1,28;-mavx;HAVE_AVX;SIMD_FLAGS" dnl - "avx;fma3;FMA3;ecx_cpuid1,12;-mfma;HAVE_FMA3;SIMD_FLAGS" dnl - "avx;fma4;FMA4;ecx_cpuid80000001,16;-mfma4;HAVE_FMA4;SIMD_FLAGS" dnl - "avx;xop;XOP;ecx_cpuid80000001,11;-mxop;HAVE_XOP;SIMD_FLAGS" dnl - "avx;avx2;AVX2;ebx_cpuid7,5;-mavx2;HAVE_AVX2;SIMD_FLAGS" dnl - "avx512;avx512f;AVX512-F;ebx_cpuid7,16;-mavx512f;HAVE_AVX512_F;SIMD_FLAGS" dnl - "avx512;avx512cd;AVX512-CD;ebx_cpuid7,28;-mavx512cd;HAVE_AVX512_CD;SIMD_FLAGS" dnl - "avx512;avx512pf;AVX512-PF;ebx_cpuid7,26;-mavx512pf;HAVE_AVX512_PF;SIMD_FLAGS" dnl - "avx512;avx512er;AVX512-ER;ebx_cpuid7,27;-mavx512er;HAVE_AVX512_ER;SIMD_FLAGS" dnl - "avx512;avx512vl;AVX512-VL;ebx_cpuid7,31;-mavx512vl;HAVE_AVX512_VL;SIMD_FLAGS" dnl - "avx512;avx512bw;AVX512-BW;ebx_cpuid7,30;-mavx512bw;HAVE_AVX512_BW;SIMD_FLAGS" dnl - "avx512;avx512dq;AVX512-DQ;ebx_cpuid7,17;-mavx512dq;HAVE_AVX512_DQ;SIMD_FLAGS" dnl - "avx512;avx512ifma;AVX512-IFMA;ebx_cpuid7,21;-mavx512ifma;HAVE_AVX512_IFMA;SIMD_FLAGS" dnl - "avx512;avx512vbmi;AVX512-VBMI;ecx_cpuid7,1;-mavx512vbmi;HAVE_AVX512_VBMI;SIMD_FLAGS" dnl - # - do ac_instr_os_support=$(eval echo \$ax_cv_have_$(echo $ac_instr_info | cut -d ";" -f 1)_os_support_ext) - ac_instr_acvar=$(echo $ac_instr_info | cut -d ";" -f 2) - ac_instr_shortname=$(echo $ac_instr_info | cut -d ";" -f 3) - ac_instr_chk_loc=$(echo $ac_instr_info | cut -d ";" -f 4) - ac_instr_chk_reg=0x$(eval echo \$$(echo $ac_instr_chk_loc | cut -d "," -f 1)) - ac_instr_chk_bit=$(echo $ac_instr_chk_loc | cut -d "," -f 2) - ac_instr_compiler_flags=$(echo $ac_instr_info | cut -d ";" -f 5) - ac_instr_have_define=$(echo $ac_instr_info | cut -d ";" -f 6) - ac_instr_flag_type=$(echo $ac_instr_info | cut -d ";" -f 7) - - AC_CACHE_CHECK([whether ${ac_instr_shortname} is supported by the processor], [ax_cv_have_${ac_instr_acvar}_cpu_ext], - [ - eval ax_cv_have_${ac_instr_acvar}_cpu_ext=no - if test "$((${ac_instr_chk_reg}>>${ac_instr_chk_bit}&0x01))" = 1 ; then - eval ax_cv_have_${ac_instr_acvar}_cpu_ext=yes - fi - ]) - - if test x"$(eval echo \$ax_cv_have_${ac_instr_acvar}_cpu_ext)" = x"yes"; then - AC_CACHE_CHECK([whether ${ac_instr_shortname} is supported by the processor and OS], [ax_cv_have_${ac_instr_acvar}_ext], - [ - eval ax_cv_have_${ac_instr_acvar}_ext=no - if test x"${ac_instr_os_support}" = x"yes"; then - eval ax_cv_have_${ac_instr_acvar}_ext=yes - fi - ]) - - if test "$(eval echo \$ax_cv_have_${ac_instr_acvar}_ext)" = yes; then - AX_CHECK_COMPILE_FLAG(${ac_instr_compiler_flags}, eval ax_cv_support_${ac_instr_acvar}_ext=yes, - eval ax_cv_support_${ac_instr_acvar}_ext=no) - if test x"$(eval echo \$ax_cv_support_${ac_instr_acvar}_ext)" = x"yes"; then - eval ${ac_instr_flag_type}=\"\$${ac_instr_flag_type} ${ac_instr_compiler_flags}\" - AC_DEFINE_UNQUOTED([${ac_instr_have_define}]) - else - AC_MSG_WARN([Your processor and OS supports ${ac_instr_shortname} instructions but not your compiler, can you try another compiler?]) - fi - else - if test x"${ac_instr_os_support}" = x"no"; then - AC_CACHE_VAL(ax_cv_support_${ac_instr_acvar}_ext, eval ax_cv_support_${ac_instr_acvar}_ext=no) - AC_MSG_WARN([Your processor supports ${ac_instr_shortname}, but your OS doesn't]) - fi - fi - else - AC_CACHE_VAL(ax_cv_have_${ac_instr_acvar}_ext, eval ax_cv_have_${ac_instr_acvar}_ext=no) - AC_CACHE_VAL(ax_cv_support_${ac_instr_acvar}_ext, eval ax_cv_support_${ac_instr_acvar}_ext=no) - fi - done - ;; - esac - - AH_TEMPLATE([HAVE_RDRND],[Define to 1 to support Digital Random Number Generator]) - AH_TEMPLATE([HAVE_BMI1],[Define to 1 to support Bit Manipulation Instruction Set 1]) - AH_TEMPLATE([HAVE_BMI2],[Define to 1 to support Bit Manipulation Instruction Set 2]) - AH_TEMPLATE([HAVE_ADX],[Define to 1 to support Multi-Precision Add-Carry Instruction Extensions]) - AH_TEMPLATE([HAVE_MPX],[Define to 1 to support Memory Protection Extensions]) - AH_TEMPLATE([HAVE_PREFETCHWT1],[Define to 1 to support Prefetch Vector Data Into Caches WT1]) - AH_TEMPLATE([HAVE_ABM],[Define to 1 to support Advanced Bit Manipulation]) - AH_TEMPLATE([HAVE_MMX],[Define to 1 to support Multimedia Extensions]) - AH_TEMPLATE([HAVE_SSE],[Define to 1 to support Streaming SIMD Extensions]) - AH_TEMPLATE([HAVE_SSE2],[Define to 1 to support Streaming SIMD Extensions]) - AH_TEMPLATE([HAVE_SSE3],[Define to 1 to support Streaming SIMD Extensions 3]) - AH_TEMPLATE([HAVE_SSSE3],[Define to 1 to support Supplemental Streaming SIMD Extensions 3]) - AH_TEMPLATE([HAVE_SSE4_1],[Define to 1 to support Streaming SIMD Extensions 4.1]) - AH_TEMPLATE([HAVE_SSE4_2],[Define to 1 to support Streaming SIMD Extensions 4.2]) - AH_TEMPLATE([HAVE_SSE4a],[Define to 1 to support AMD Streaming SIMD Extensions 4a]) - AH_TEMPLATE([HAVE_SHA],[Define to 1 to support Secure Hash Algorithm Extension]) - AH_TEMPLATE([HAVE_AES],[Define to 1 to support Advanced Encryption Standard New Instruction Set (AES-NI)]) - AH_TEMPLATE([HAVE_AVX],[Define to 1 to support Advanced Vector Extensions]) - AH_TEMPLATE([HAVE_FMA3],[Define to 1 to support Fused Multiply-Add Extensions 3]) - AH_TEMPLATE([HAVE_FMA4],[Define to 1 to support Fused Multiply-Add Extensions 4]) - AH_TEMPLATE([HAVE_XOP],[Define to 1 to support eXtended Operations Extensions]) - AH_TEMPLATE([HAVE_AVX2],[Define to 1 to support Advanced Vector Extensions 2]) - AH_TEMPLATE([HAVE_AVX512_F],[Define to 1 to support AVX-512 Foundation Extensions]) - AH_TEMPLATE([HAVE_AVX512_CD],[Define to 1 to support AVX-512 Conflict Detection Instructions]) - AH_TEMPLATE([HAVE_AVX512_PF],[Define to 1 to support AVX-512 Conflict Prefetch Instructions]) - AH_TEMPLATE([HAVE_AVX512_ER],[Define to 1 to support AVX-512 Exponential & Reciprocal Instructions]) - AH_TEMPLATE([HAVE_AVX512_VL],[Define to 1 to support AVX-512 Vector Length Extensions]) - AH_TEMPLATE([HAVE_AVX512_BW],[Define to 1 to support AVX-512 Byte and Word Instructions]) - AH_TEMPLATE([HAVE_AVX512_DQ],[Define to 1 to support AVX-512 Doubleword and Quadword Instructions]) - AH_TEMPLATE([HAVE_AVX512_IFMA],[Define to 1 to support AVX-512 Integer Fused Multiply Add Instructions]) - AH_TEMPLATE([HAVE_AVX512_VBMI],[Define to 1 to support AVX-512 Vector Byte Manipulation Instructions]) - AC_SUBST(SIMD_FLAGS) - AC_SUBST(CPUEXT_FLAGS) -]) diff --git a/m4/ax_gcc_x86_avx_xgetbv.m4 b/m4/ax_gcc_x86_avx_xgetbv.m4 deleted file mode 100644 index a57fc5ee..00000000 --- a/m4/ax_gcc_x86_avx_xgetbv.m4 +++ /dev/null @@ -1,79 +0,0 @@ -# =========================================================================== -# https://www.gnu.org/software/autoconf-archive/ax_gcc_x86_avx_xgetbv.html -# =========================================================================== -# -# SYNOPSIS -# -# AX_GCC_X86_AVX_XGETBV -# -# DESCRIPTION -# -# On later x86 processors with AVX SIMD support, with gcc or a compiler -# that has a compatible syntax for inline assembly instructions, run a -# small program that executes the xgetbv instruction with input OP. This -# can be used to detect if the OS supports AVX instruction usage. -# -# On output, the values of the eax and edx registers are stored as -# hexadecimal strings as "eax:edx" in the cache variable -# ax_cv_gcc_x86_avx_xgetbv. -# -# If the xgetbv instruction fails (because you are running a -# cross-compiler, or because you are not using gcc, or because you are on -# a processor that doesn't have this instruction), -# ax_cv_gcc_x86_avx_xgetbv_OP is set to the string "unknown". -# -# This macro mainly exists to be used in AX_EXT. -# -# LICENSE -# -# Copyright (c) 2013 Michael Petch -# -# This program is free software: you can redistribute it and/or modify it -# under the terms of the GNU General Public License as published by the -# Free Software Foundation, either version 3 of the License, or (at your -# option) any later version. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General -# Public License for more details. -# -# You should have received a copy of the GNU General Public License along -# with this program. If not, see . -# -# As a special exception, the respective Autoconf Macro's copyright owner -# gives unlimited permission to copy, distribute and modify the configure -# scripts that are the output of Autoconf when processing the Macro. You -# need not follow the terms of the GNU General Public License when using -# or distributing such scripts, even though portions of the text of the -# Macro appear in them. The GNU General Public License (GPL) does govern -# all other use of the material that constitutes the Autoconf Macro. -# -# This special exception to the GPL applies to versions of the Autoconf -# Macro released by the Autoconf Archive. When you make and distribute a -# modified version of the Autoconf Macro, you may extend this special -# exception to the GPL to apply to your modified version as well. - -#serial 3 - -AC_DEFUN([AX_GCC_X86_AVX_XGETBV], -[AC_REQUIRE([AC_PROG_CC]) -AC_LANG_PUSH([C]) -AC_CACHE_CHECK(for x86-AVX xgetbv $1 output, ax_cv_gcc_x86_avx_xgetbv_$1, - [AC_RUN_IFELSE([AC_LANG_PROGRAM([#include ], [ - int op = $1, eax, edx; - FILE *f; - /* Opcodes for xgetbv */ - __asm__ __volatile__ (".byte 0x0f, 0x01, 0xd0" - : "=a" (eax), "=d" (edx) - : "c" (op)); - f = fopen("conftest_xgetbv", "w"); if (!f) return 1; - fprintf(f, "%x:%x\n", eax, edx); - fclose(f); - return 0; -])], - [ax_cv_gcc_x86_avx_xgetbv_$1=`cat conftest_xgetbv`; rm -f conftest_xgetbv], - [ax_cv_gcc_x86_avx_xgetbv_$1=unknown; rm -f conftest_xgetbv], - [ax_cv_gcc_x86_avx_xgetbv_$1=unknown])]) -AC_LANG_POP([C]) -]) diff --git a/m4/ax_gcc_x86_cpuid.m4 b/m4/ax_gcc_x86_cpuid.m4 deleted file mode 100644 index df954658..00000000 --- a/m4/ax_gcc_x86_cpuid.m4 +++ /dev/null @@ -1,89 +0,0 @@ -# =========================================================================== -# https://www.gnu.org/software/autoconf-archive/ax_gcc_x86_cpuid.html -# =========================================================================== -# -# SYNOPSIS -# -# AX_GCC_X86_CPUID(OP) -# AX_GCC_X86_CPUID_COUNT(OP, COUNT) -# -# DESCRIPTION -# -# On Pentium and later x86 processors, with gcc or a compiler that has a -# compatible syntax for inline assembly instructions, run a small program -# that executes the cpuid instruction with input OP. This can be used to -# detect the CPU type. AX_GCC_X86_CPUID_COUNT takes an additional COUNT -# parameter that gets passed into register ECX before calling cpuid. -# -# On output, the values of the eax, ebx, ecx, and edx registers are stored -# as hexadecimal strings as "eax:ebx:ecx:edx" in the cache variable -# ax_cv_gcc_x86_cpuid_OP. -# -# If the cpuid instruction fails (because you are running a -# cross-compiler, or because you are not using gcc, or because you are on -# a processor that doesn't have this instruction), ax_cv_gcc_x86_cpuid_OP -# is set to the string "unknown". -# -# This macro mainly exists to be used in AX_GCC_ARCHFLAG. -# -# LICENSE -# -# Copyright (c) 2008 Steven G. Johnson -# Copyright (c) 2008 Matteo Frigo -# Copyright (c) 2015 Michael Petch -# -# This program is free software: you can redistribute it and/or modify it -# under the terms of the GNU General Public License as published by the -# Free Software Foundation, either version 3 of the License, or (at your -# option) any later version. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General -# Public License for more details. -# -# You should have received a copy of the GNU General Public License along -# with this program. If not, see . -# -# As a special exception, the respective Autoconf Macro's copyright owner -# gives unlimited permission to copy, distribute and modify the configure -# scripts that are the output of Autoconf when processing the Macro. You -# need not follow the terms of the GNU General Public License when using -# or distributing such scripts, even though portions of the text of the -# Macro appear in them. The GNU General Public License (GPL) does govern -# all other use of the material that constitutes the Autoconf Macro. -# -# This special exception to the GPL applies to versions of the Autoconf -# Macro released by the Autoconf Archive. When you make and distribute a -# modified version of the Autoconf Macro, you may extend this special -# exception to the GPL to apply to your modified version as well. - -#serial 10 - -AC_DEFUN([AX_GCC_X86_CPUID], -[AX_GCC_X86_CPUID_COUNT($1, 0) -]) - -AC_DEFUN([AX_GCC_X86_CPUID_COUNT], -[AC_REQUIRE([AC_PROG_CC]) -AC_LANG_PUSH([C]) -AC_CACHE_CHECK(for x86 cpuid $1 output, ax_cv_gcc_x86_cpuid_$1, - [AC_RUN_IFELSE([AC_LANG_PROGRAM([#include ], [ - int op = $1, level = $2, eax, ebx, ecx, edx; - FILE *f; - __asm__ __volatile__ ("xchg %%ebx, %1\n" - "cpuid\n" - "xchg %%ebx, %1\n" - : "=a" (eax), "=r" (ebx), "=c" (ecx), "=d" (edx) - : "a" (op), "2" (level)); - - f = fopen("conftest_cpuid", "w"); if (!f) return 1; - fprintf(f, "%x:%x:%x:%x\n", eax, ebx, ecx, edx); - fclose(f); - return 0; -])], - [ax_cv_gcc_x86_cpuid_$1=`cat conftest_cpuid`; rm -f conftest_cpuid], - [ax_cv_gcc_x86_cpuid_$1=unknown; rm -f conftest_cpuid], - [ax_cv_gcc_x86_cpuid_$1=unknown])]) -AC_LANG_POP([C]) -]) diff --git a/src/fglm/Makefile.am b/src/fglm/Makefile.am index b5e8767a..346e0447 100644 --- a/src/fglm/Makefile.am +++ b/src/fglm/Makefile.am @@ -2,7 +2,7 @@ noinst_LTLIBRARIES = libfglm.la libfglm_la_SOURCES = fglm_core.c libfglm_ladir = $(includedir)/msolve/fglm libfglm_la_HEADERS = fglm.h -libfglm_la_CFLAGS = $(SIMD_FLAGS) $(CPUEXT_FLAGS) $(OPENMP_CFLAGS) -Wall -Wextra +libfglm_la_CFLAGS = $(OPENMP_CFLAGS) -Wall -Wextra EXTRA_DIST = fglm.h \ libfglm.h \ diff --git a/src/fglm/fglm_core.c b/src/fglm/fglm_core.c index 1aac58ed..ad88a76c 100644 --- a/src/fglm/fglm_core.c +++ b/src/fglm/fglm_core.c @@ -423,26 +423,8 @@ static inline void sparse_mat_fglm_mult_vec(CF_t *res, sp_matfglm_t *mat, for(szmat_t i = 0; i < ntriv; i++){ res[mat->triv_idx[i]] = vec[mat->triv_pos[i]]; } -#if defined(HAVE_AVX512_F) - nmod_t mod; - uint64_t pow2_precomp; - nmod_init(&mod, (uint64_t)prime); - NMOD_RED(pow2_precomp, (UINT64_C(1) << __DOT_SPLIT_BITS), mod); - - _avx512_matrix_vector_product(vres, mat->dense_mat, vec, mat->dst, - ncols, nrows, mod, pow2_precomp, st); -#elif defined(HAVE_AVX2) - nmod_t mod; - uint64_t pow2_precomp; - nmod_init(&mod, (uint64_t)prime); - NMOD_RED(pow2_precomp, (UINT64_C(1) << __DOT_SPLIT_BITS), mod); - - _avx2_matrix_vector_product(vres, mat->dense_mat, vec, mat->dst, - ncols, nrows, mod, pow2_precomp, st); -#else - non_avx_matrix_vector_product(vres, mat->dense_mat, vec, - ncols, nrows, prime); -#endif + dense_matrix_vector_product(vres, mat->dense_mat, vec, mat->dst, + ncols, nrows, prime, st); for(szmat_t i = 0; i < nrows; i++){ res[mat->dense_idx[i]] = vres[i]; @@ -478,26 +460,8 @@ static inline void sparse_mat_fglm_colon_mult_vec(CF_t *res, sp_matfglmcol_t *ma } /* printf ("zero\n"); */ /* printf("ncols %u\n", ncols); */ -#if defined(HAVE_AVX512_F) - nmod_t mod; - uint64_t pow2_precomp; - nmod_init(&mod, (uint64_t)prime); - NMOD_RED(pow2_precomp, (UINT64_C(1) << __DOT_SPLIT_BITS), mod); - - _avx512_matrix_vector_product(vres, mat->dense_mat, vec, mat->dst, - ncols, nrows, mod, pow2_precomp, st); -#elif defined(HAVE_AVX2) - nmod_t mod; - uint64_t pow2_precomp; - nmod_init(&mod, (uint64_t)prime); - NMOD_RED(pow2_precomp, (UINT64_C(1) << __DOT_SPLIT_BITS), mod); - - _avx2_matrix_vector_product(vres, mat->dense_mat, vec, mat->dst, - ncols, nrows, mod, pow2_precomp, st); -#else - non_avx_matrix_vector_product(vres, mat->dense_mat, vec, - ncols, nrows, prime); -#endif + dense_matrix_vector_product(vres, mat->dense_mat, vec, mat->dst, + ncols, nrows, prime, st); for(szmat_t i = 0; i < nrows; i++){ res[mat->dense_idx[i]] = vres[i]; } @@ -603,17 +567,20 @@ static void generate_matrix_sequence(sp_matfglm_t *matxn, fglm_data_t *data, szmat_t nb = 2 * matxn->ncols / BL; for(szmat_t i = 0; i < nb; i++){ -#ifdef HAVE_AVX2 - sparse_matfglm_mul(res, matxn, Rmat, - tres, - BL, - prime, preinv, - RED_32, - RED_64); -#else - fprintf(ERRSTREAM, "Not implemented yet\n"); - exit(1); +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + sparse_matfglm_mul(res, matxn, Rmat, + tres, + BL, + prime, preinv, + RED_32, + RED_64); + } else #endif + { + fprintf(ERRSTREAM, "Not implemented yet\n"); + exit(1); + } } posix_memalign_free(Rmat); posix_memalign_free(res); diff --git a/src/fglm/linalg-fglm.c b/src/fglm/linalg-fglm.c index 01463d2d..6ac305be 100644 --- a/src/fglm/linalg-fglm.c +++ b/src/fglm/linalg-fglm.c @@ -26,9 +26,7 @@ #include #endif -#ifdef HAVE_AVX2 -#include -#endif +#include "../msolve/cpu_features.h" // for dot product, simultaneous dot products, and matrix-vector products, // vectorized functions below accumulate 8 terms: @@ -199,10 +197,10 @@ static inline void non_avx_matrix_vector_product(uint32_t* vec_res, const uint32 /* vectorized (AVX2) matrix vector product */ /*-----------------------------------------*/ -#ifdef HAVE_AVX2 +#ifdef HAVE_AVX2_KERNELS // avx2 horizontal sum -static inline +static inline TARGET_AVX2 uint64_t _mm256_hsum(__m256i a) { __m256i a_hi = _mm256_shuffle_epi32(a, 14); // 14 == 0b00001110 @@ -218,6 +216,7 @@ uint64_t _mm256_hsum(__m256i a) #endif } +TARGET_AVX2 uint32_t _nmod32_vec_dot_split_avx2(const uint32_t * vec1, const uint32_t * vec2, int64_t len, nmod_t mod, uint64_t pow2_precomp) { @@ -310,6 +309,7 @@ uint32_t _nmod32_vec_dot_split_avx2(const uint32_t * vec1, const uint32_t * vec2 return (uint32_t)res; } +TARGET_AVX2 void _nmod32_vec_dot2_split_avx2(uint32_t * res0, uint32_t * res1, const uint32_t * vec1, const uint32_t * vec2_0, const uint32_t * vec2_1, int64_t len, nmod_t mod, uint64_t pow2_precomp) @@ -440,6 +440,7 @@ void _nmod32_vec_dot2_split_avx2(uint32_t * res0, uint32_t * res1, NMOD_RED(*res1, pow2_precomp * hsum_hi1 + hsum_lo1, mod); } +TARGET_AVX2 void _nmod32_vec_dot3_split_avx2(uint32_t * res, const uint32_t * vec1, const uint32_t * vec2_0, const uint32_t * vec2_1, const uint32_t * vec2_2, int64_t len, nmod_t mod, uint64_t pow2_precomp) @@ -662,14 +663,15 @@ static inline void _avx2_matrix_vector_product(uint32_t * vec_res, /* vectorized (AVX512) matrix vector product */ /*-------------------------------------------*/ -#ifdef HAVE_AVX512_F +#ifdef HAVE_AVX512_KERNELS // avx512 horizontal sum -FLINT_FORCE_INLINE uint64_t _mm512_hsum(__m512i a) +FLINT_FORCE_INLINE TARGET_AVX512 uint64_t _mm512_hsum(__m512i a) { return _mm512_reduce_add_epi64(a); } +TARGET_AVX512 uint32_t _nmod32_vec_dot_split_avx512(const uint32_t * vec1, const uint32_t * vec2, int64_t len, nmod_t mod, uint64_t pow2_precomp) { const __m512i low_bits = _mm512_set1_epi64(__DOT_SPLIT_MASK); @@ -772,6 +774,7 @@ uint32_t _nmod32_vec_dot_split_avx512(const uint32_t * vec1, const uint32_t * ve return (uint32_t)res; } +TARGET_AVX512 void _nmod32_vec_dot2_split_avx512(uint32_t * res0, uint32_t * res1, const uint32_t * vec1, const uint32_t * vec2_0, const uint32_t * vec2_1, int64_t len, nmod_t mod, uint64_t pow2_precomp) @@ -915,6 +918,7 @@ void _nmod32_vec_dot2_split_avx512(uint32_t * res0, uint32_t * res1, NMOD_RED(*res1, pow2_precomp * hsum_hi1 + hsum_lo1, mod); } +TARGET_AVX512 void _nmod32_vec_dot3_split_avx512(uint32_t * res, const uint32_t * vec1, const uint32_t * vec2_0, const uint32_t * vec2_1, const uint32_t * vec2_2, int64_t len, nmod_t mod, uint64_t pow2_precomp) @@ -1130,3 +1134,43 @@ static inline void _avx512_matrix_vector_product(uint32_t * vec_res, } #endif + +/*-----------------------------------------------------------*/ +/* matrix vector product, choosing the kernel at runtime */ +/*-----------------------------------------------------------*/ + +static inline void dense_matrix_vector_product(uint32_t * vec_res, + const uint32_t * mat, + const uint32_t * vec, + const uint32_t * dst, + const uint32_t ncols, + const uint32_t nrows, + const uint32_t prime, + md_t *st) +{ +#if defined(HAVE_AVX2_KERNELS) || defined(HAVE_AVX512_KERNELS) + nmod_t mod; + uint64_t pow2_precomp; + nmod_init(&mod, (uint64_t)prime); + NMOD_RED(pow2_precomp, (UINT64_C(1) << __DOT_SPLIT_BITS), mod); +#else + (void)dst; + (void)st; +#endif + +#ifdef HAVE_AVX512_KERNELS + if (cpu_has_avx512()) { + _avx512_matrix_vector_product(vec_res, mat, vec, dst, + ncols, nrows, mod, pow2_precomp, st); + return; + } +#endif +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + _avx2_matrix_vector_product(vec_res, mat, vec, dst, + ncols, nrows, mod, pow2_precomp, st); + return; + } +#endif + non_avx_matrix_vector_product(vec_res, mat, vec, ncols, nrows, prime); +} diff --git a/src/fglm/matrix-mult.c b/src/fglm/matrix-mult.c index e5d47d98..5b3abfa8 100644 --- a/src/fglm/matrix-mult.c +++ b/src/fglm/matrix-mult.c @@ -21,8 +21,9 @@ //-1 sur 32 bits #define MONE32 ((uint32_t)0xFFFFFFFF) -#ifdef HAVE_AVX2 -#include +#include "../msolve/cpu_features.h" + +#ifdef HAVE_AVX2_KERNELS #define AVX2LOAD(A) _mm256_load_si256((__m256i*)(A)) #define AVX2LOADU(A) _mm256_loadu_si256((__m256i*)(A)) @@ -112,8 +113,8 @@ static inline uint64_t ADDMODRED32(uint64_t a, uint64_t b, * } */ -#if HAVE_AVX2 -static inline void REDUCE(uint64_t *acc64, uint64_t *acc4x64, +#ifdef HAVE_AVX2_KERNELS +static inline TARGET_AVX2 void REDUCE(uint64_t *acc64, uint64_t *acc4x64, __m256i acc_low, __m256i acc_high, const uint32_t fc, const uint32_t preinv, const uint32_t RED_32, const uint32_t RED_64){ @@ -167,7 +168,7 @@ static inline void make_zero(uint32_t ** D, } -#ifdef HAVE_AVX2 +#ifdef HAVE_AVX2_KERNELS /* m = number of rows in A l = number of cols in A = number of rows in B @@ -176,6 +177,7 @@ static inline void make_zero(uint32_t ** D, D = A*B */ +TARGET_AVX2 void _mod_mat_addmul_transpose_op(uint32_t *D, uint32_t *A, uint32_t *B, const uint32_t m, const uint32_t l, @@ -267,7 +269,7 @@ void _mod_mat_addmul_transpose_op(uint32_t *D, #endif /* not yet implemented for non-AVX machine */ -#ifdef HAVE_AVX2 +#ifdef HAVE_AVX2_KERNELS static inline void sparse_matfglm_mul(CF_t *res, sp_matfglm_t *matxn, CF_t *R, CF_t *tres, const int nc, diff --git a/src/msolve/Makefile.am b/src/msolve/Makefile.am index 1f41ac4e..f2c50856 100644 --- a/src/msolve/Makefile.am +++ b/src/msolve/Makefile.am @@ -2,7 +2,7 @@ lib_LTLIBRARIES = libmsolve.la libmsolve_la_SOURCES = libmsolve.c libmsolve_ladir = $(includedir)/msolve/msolve libmsolve_la_HEADERS = msolve-data.h msolve-version.h msolve.h -libmsolve_la_CFLAGS = $(SIMD_FLAGS) $(CPUEXT_FLAGS) $(OPENMP_CFLAGS) -Wall -Wextra -Wno-unused-parameter -Wno-sign-compare +libmsolve_la_CFLAGS = $(OPENMP_CFLAGS) -Wall -Wextra -Wno-unused-parameter -Wno-sign-compare libmsolve_la_LDFLAGS = -version-info $(LT_VERSION) libmsolve_la_LIBADD = ../usolve/libusolve.la ../fglm/libfglm.la ../neogb/libneogb.la @@ -20,6 +20,7 @@ EXTRA_DIST = msolve-data.h \ msolve.c \ primes.c \ streams.h \ + cpu_features.h \ ../crt/longlong.h \ ../crt/ulong_extras.h \ ../crt/mpq_reconstruct.c \ diff --git a/src/msolve/cpu_features.h b/src/msolve/cpu_features.h new file mode 100644 index 00000000..f6409b46 --- /dev/null +++ b/src/msolve/cpu_features.h @@ -0,0 +1,74 @@ +/* This file is part of msolve. + * + * msolve is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * msolve is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with msolve. If not, see + * + * Authors: + * Jérémy Berthomieu + * Christian Eder + * Mohab Safey El Din */ + +#ifndef CPU_FEATURES_H +#define CPU_FEATURES_H + +/* SIMD intrinsics and CPU feature detection for the vectorized kernels. + * + * x86-64: HAVE_AVX2_KERNELS resp. HAVE_AVX512_KERNELS are defined by + * configure whenever the compiler is able to build AVX2 resp. AVX-512 + * kernels, independently of the CPU of the build machine. Such kernels are + * marked with TARGET_AVX2 resp. TARGET_AVX512, the remaining code is + * compiled for the baseline architecture. Callers choose a kernel at runtime + * via cpu_has_avx2() resp. cpu_has_avx512(), so that the same binary runs on + * every x86-64 CPU. If msolve is compiled with these extensions enabled + * globally (e.g. with -march=native) the checks are resolved at compile time. + * + * AArch64: NEON is part of the baseline architecture, so the NEON kernels + * are selected at compile time via __aarch64__. */ + +#if defined HAVE_AVX2_KERNELS || defined HAVE_AVX512_KERNELS +#include +#endif + +#ifdef __aarch64__ +#include +#endif + +#ifdef HAVE_AVX2_KERNELS +#define TARGET_AVX2 __attribute__((target("avx2"))) + +static inline int cpu_has_avx2(void) +{ +#ifdef __AVX2__ + return 1; +#else + return __builtin_cpu_supports("avx2"); +#endif +} +#endif + +#ifdef HAVE_AVX512_KERNELS +/* the 8 and 16 bit kernels need AVX-512 BW for 16 bit multiplications */ +#define TARGET_AVX512 __attribute__((target("avx512f,avx512bw"))) + +static inline int cpu_has_avx512(void) +{ +#if defined __AVX512F__ && defined __AVX512BW__ + return 1; +#else + return __builtin_cpu_supports("avx512f") + && __builtin_cpu_supports("avx512bw"); +#endif +} +#endif + +#endif diff --git a/src/neogb/Makefile.am b/src/neogb/Makefile.am index f9800338..24b40c4e 100644 --- a/src/neogb/Makefile.am +++ b/src/neogb/Makefile.am @@ -5,7 +5,7 @@ libneogb_la_HEADERS =libneogb.h basis.h data.h engine.h f4.h sba.h hash.h \ io.h modular.h nf.h f4sat.h sort_r.h meta_data.h \ tools.h update.h libneogb_la_LDFLAGS = -version-info $(LT_VERSION) -libneogb_la_CFLAGS = $(SIMD_FLAGS) $(CPUEXT_FLAGS) $(OPENMP_CFLAGS) -Wall -Wextra -Wno-unused-parameter -Wno-sign-compare +libneogb_la_CFLAGS = $(OPENMP_CFLAGS) -Wall -Wextra -Wno-unused-parameter -Wno-sign-compare EXTRA_DIST = basis.h \ data.h \ diff --git a/src/neogb/la_ff_16.c b/src/neogb/la_ff_16.c index de7ba538..52d9cfc1 100644 --- a/src/neogb/la_ff_16.c +++ b/src/neogb/la_ff_16.c @@ -21,12 +21,7 @@ #include "data.h" #include "../msolve/streams.h" -/* That's also enough if AVX512 is avaialable on the system */ -#if defined HAVE_AVX2 -#include -#elif defined __aarch64__ -#include -#endif +#include "../msolve/cpu_features.h" static inline cf16_t *normalize_dense_matrix_row_ff_16( cf16_t *row, @@ -94,6 +89,309 @@ static inline cf16_t *normalize_sparse_matrix_row_ff_16( return row; } +static inline void add_mul_row_ff_16_scalar( + int64_t *dr, + const cf16_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + const len_t os = len % UNROLL; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += UNROLL) { + dr[ds[j]] += mul * cfs[j]; + dr[ds[j+1]] += mul * cfs[j+1]; + dr[ds[j+2]] += mul * cfs[j+2]; + dr[ds[j+3]] += mul * cfs[j+3]; + } +} + +#ifdef __aarch64__ +static inline void add_mul_row_ff_16_neon( + int64_t *dr, + const cf16_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + uint64_t tmp[2] __attribute__((aligned(32))); + uint32x4_t prodv; + uint16x8_t redv; + uint64x2_t drv, resv; + + const len_t os = len % 8; + const cf16_t mul16 = (cf16_t)mul; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 8) { + tmp[0] = (uint64_t)dr[ds[j]]; + tmp[1] = (uint64_t)dr[ds[j+1]]; + drv = vld1q_u64(tmp); + redv = vld1q_u16((cf16_t *)(cfs)+j); + + prodv = vmull_n_u16(vget_low_u16(redv), mul16); + resv = vaddw_u32(drv, vget_low_u32(prodv)); + vst1q_u64(tmp, resv); + dr[ds[j]] = (int64_t)tmp[0]; + dr[ds[j+1]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+2]]; + tmp[1] = (uint64_t)dr[ds[j+3]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_high_u32(prodv)); + vst1q_u64(tmp, resv); + dr[ds[j+2]] = (int64_t)tmp[0]; + dr[ds[j+3]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+4]]; + tmp[1] = (uint64_t)dr[ds[j+5]]; + drv = vld1q_u64(tmp); + + prodv = vmull_n_u16( vget_high_u16(redv), mul16); + resv = vaddw_u32(drv, vget_low_u32(prodv)); + vst1q_u64(tmp, resv); + dr[ds[j+4]] = (int64_t)tmp[0]; + dr[ds[j+5]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+6]]; + tmp[1] = (uint64_t)dr[ds[j+7]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_high_u32(prodv)); + vst1q_u64(tmp, resv); + dr[ds[j+6]] = (int64_t)tmp[0]; + dr[ds[j+7]] = (int64_t)tmp[1]; + } +} +#endif + +#ifdef HAVE_AVX2_KERNELS +static TARGET_AVX2 void add_mul_row_ff_16_avx2( + int64_t *dr, + const cf16_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + uint32_t mone32 = (uint32_t)0xFFFFFFFF; + uint16_t mone16 = (uint16_t)0xFFFF; + uint32_t mone16h = (uint32_t)0xFFFF0000; + __m256i mask32 = _mm256_set1_epi64x(mone32); + __m256i mask16 = _mm256_set1_epi32(mone16); + __m256i mask16h = _mm256_set1_epi32(mone16h); + + int64_t res[4] __attribute__((aligned(32))); + __m256i redv, mulv, prodh, prodl, prod, drv, resv; + + const uint16_t mul16 = (uint16_t)mul; + mulv = _mm256_set1_epi16(mul16); + const len_t os = len % 16; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 16) { + redv = _mm256_loadu_si256((__m256i*)(cfs+j)); + prodh = _mm256_mulhi_epu16(mulv, redv); + prodl = _mm256_mullo_epi16(mulv, redv); + prod = _mm256_xor_si256( + _mm256_and_si256(prodh, mask16h), _mm256_srli_epi32(prodl, 16)); + drv = _mm256_setr_epi64x( + dr[ds[j+1]], + dr[ds[j+5]], + dr[ds[j+9]], + dr[ds[j+13]]); + resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask32)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+1]] = res[0]; + dr[ds[j+5]] = res[1]; + dr[ds[j+9]] = res[2]; + dr[ds[j+13]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+3]], + dr[ds[j+7]], + dr[ds[j+11]], + dr[ds[j+15]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(prod, 32)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+3]] = res[0]; + dr[ds[j+7]] = res[1]; + dr[ds[j+11]] = res[2]; + dr[ds[j+15]] = res[3]; + prod = _mm256_xor_si256( + _mm256_slli_epi32(prodh, 16), _mm256_and_si256(prodl, mask16)); + drv = _mm256_setr_epi64x( + dr[ds[j+0]], + dr[ds[j+4]], + dr[ds[j+8]], + dr[ds[j+12]]); + resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask32)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+0]] = res[0]; + dr[ds[j+4]] = res[1]; + dr[ds[j+8]] = res[2]; + dr[ds[j+12]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+2]], + dr[ds[j+6]], + dr[ds[j+10]], + dr[ds[j+14]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(prod, 32)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+2]] = res[0]; + dr[ds[j+6]] = res[1]; + dr[ds[j+10]] = res[2]; + dr[ds[j+14]] = res[3]; + } +} +#endif + +#ifdef HAVE_AVX512_KERNELS +static TARGET_AVX512 void add_mul_row_ff_16_avx512( + int64_t *dr, + const cf16_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + uint32_t mone32 = (uint32_t)0xFFFFFFFF; + uint16_t mone16 = (uint16_t)0xFFFF; + uint32_t mone16h = (uint32_t)0xFFFF0000; + __m512i mask32 = _mm512_set1_epi64(mone32); + __m512i mask16 = _mm512_set1_epi32(mone16); + __m512i mask16h = _mm512_set1_epi32(mone16h); + + int64_t res[8] __attribute__((aligned(64))); + __m512i redv, mulv, prodh, prodl, prod, drv, resv; + + const uint16_t mul16 = (uint16_t)mul; + mulv = _mm512_set1_epi16(mul16); + const len_t os = len % 32; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 32) { + redv = _mm512_loadu_si512((__m512i*)(cfs+j)); + prodh = _mm512_mulhi_epu16(mulv, redv); + prodl = _mm512_mullo_epi16(mulv, redv); + prod = _mm512_xor_si512( + _mm512_and_si512(prodh, mask16h), _mm512_srli_epi32(prodl, 16)); + drv = _mm512_setr_epi64( + dr[ds[j+1]], + dr[ds[j+5]], + dr[ds[j+9]], + dr[ds[j+13]], + dr[ds[j+17]], + dr[ds[j+21]], + dr[ds[j+25]], + dr[ds[j+29]]); + resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask32)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+1]] = res[0]; + dr[ds[j+5]] = res[1]; + dr[ds[j+9]] = res[2]; + dr[ds[j+13]] = res[3]; + dr[ds[j+17]] = res[4]; + dr[ds[j+21]] = res[5]; + dr[ds[j+25]] = res[6]; + dr[ds[j+29]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+3]], + dr[ds[j+7]], + dr[ds[j+11]], + dr[ds[j+15]], + dr[ds[j+19]], + dr[ds[j+23]], + dr[ds[j+27]], + dr[ds[j+31]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(prod, 32)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+3]] = res[0]; + dr[ds[j+7]] = res[1]; + dr[ds[j+11]] = res[2]; + dr[ds[j+15]] = res[3]; + dr[ds[j+19]] = res[4]; + dr[ds[j+23]] = res[5]; + dr[ds[j+27]] = res[6]; + dr[ds[j+31]] = res[7]; + prod = _mm512_xor_si512( + _mm512_slli_epi32(prodh, 16), _mm512_and_si512(prodl, mask16)); + drv = _mm512_setr_epi64( + dr[ds[j+0]], + dr[ds[j+4]], + dr[ds[j+8]], + dr[ds[j+12]], + dr[ds[j+16]], + dr[ds[j+20]], + dr[ds[j+24]], + dr[ds[j+28]]); + resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask32)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+0]] = res[0]; + dr[ds[j+4]] = res[1]; + dr[ds[j+8]] = res[2]; + dr[ds[j+12]] = res[3]; + dr[ds[j+16]] = res[4]; + dr[ds[j+20]] = res[5]; + dr[ds[j+24]] = res[6]; + dr[ds[j+28]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+2]], + dr[ds[j+6]], + dr[ds[j+10]], + dr[ds[j+14]], + dr[ds[j+18]], + dr[ds[j+22]], + dr[ds[j+26]], + dr[ds[j+30]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(prod, 32)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+2]] = res[0]; + dr[ds[j+6]] = res[1]; + dr[ds[j+10]] = res[2]; + dr[ds[j+14]] = res[3]; + dr[ds[j+18]] = res[4]; + dr[ds[j+22]] = res[5]; + dr[ds[j+26]] = res[6]; + dr[ds[j+30]] = res[7]; + } +} +#endif + +/* dr[ds[j]] += mul * cfs[j] for 0 <= j < len */ +static inline void add_mul_row_ff_16( + int64_t *dr, + const cf16_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ +#ifdef HAVE_AVX512_KERNELS + if (cpu_has_avx512()) { + add_mul_row_ff_16_avx512(dr, cfs, ds, len, mul); + return; + } +#endif +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + add_mul_row_ff_16_avx2(dr, cfs, ds, len, mul); + return; + } +#endif +#ifdef __aarch64__ + add_mul_row_ff_16_neon(dr, cfs, ds, len, mul); +#else + add_mul_row_ff_16_scalar(dr, cfs, ds, len, mul); +#endif +} + static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_16( int64_t *dr, mat_t *mat, @@ -122,32 +420,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_16( } else { rba = NULL; } -#if defined HAVE_AVX512_F - uint32_t mone32 = (uint32_t)0xFFFFFFFF; - uint16_t mone16 = (uint16_t)0xFFFF; - uint32_t mone16h = (uint32_t)0xFFFF0000; - __m512i mask32 = _mm512_set1_epi64(mone32); - __m512i mask16 = _mm512_set1_epi32(mone16); - __m512i mask16h = _mm512_set1_epi32(mone16h); - - int64_t res[8] __attribute__((aligned(64))); - __m512i redv, mulv, prodh, prodl, prod, drv, resv; -#elif defined HAVE_AVX2 - uint32_t mone32 = (uint32_t)0xFFFFFFFF; - uint16_t mone16 = (uint16_t)0xFFFF; - uint32_t mone16h = (uint32_t)0xFFFF0000; - __m256i mask32 = _mm256_set1_epi64x(mone32); - __m256i mask16 = _mm256_set1_epi32(mone16); - __m256i mask16h = _mm256_set1_epi32(mone16h); - - int64_t res[4] __attribute__((aligned(32))); - __m256i redv, mulv, prodh, prodl, prod, drv, resv; -#elif defined __aarch64__ - uint64_t tmp[2] __attribute__((aligned(32))); - uint32x4_t prodv; - uint16x8_t redv; - uint64x2_t drv, resv; -#endif k = 0; for (i = dpiv; i < ncols; ++i) { @@ -174,219 +446,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_16( } } cfs = mcf[dts[COEFFS]]; -#if defined HAVE_AVX512_F - const uint16_t mul16 = (uint16_t)(fc - dr[i]); - mulv = _mm512_set1_epi16(mul16); - const len_t len = dts[LENGTH]; - const len_t os = len % 32; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 32) { - redv = _mm512_loadu_si512((__m512i*)(cfs+j)); - prodh = _mm512_mulhi_epu16(mulv, redv); - prodl = _mm512_mullo_epi16(mulv, redv); - prod = _mm512_xor_si512( - _mm512_and_si512(prodh, mask16h), _mm512_srli_epi32(prodl, 16)); - drv = _mm512_setr_epi64( - dr[ds[j+1]], - dr[ds[j+5]], - dr[ds[j+9]], - dr[ds[j+13]], - dr[ds[j+17]], - dr[ds[j+21]], - dr[ds[j+25]], - dr[ds[j+29]]); - resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask32)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+1]] = res[0]; - dr[ds[j+5]] = res[1]; - dr[ds[j+9]] = res[2]; - dr[ds[j+13]] = res[3]; - dr[ds[j+17]] = res[4]; - dr[ds[j+21]] = res[5]; - dr[ds[j+25]] = res[6]; - dr[ds[j+29]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+3]], - dr[ds[j+7]], - dr[ds[j+11]], - dr[ds[j+15]], - dr[ds[j+19]], - dr[ds[j+23]], - dr[ds[j+27]], - dr[ds[j+31]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(prod, 32)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+3]] = res[0]; - dr[ds[j+7]] = res[1]; - dr[ds[j+11]] = res[2]; - dr[ds[j+15]] = res[3]; - dr[ds[j+19]] = res[4]; - dr[ds[j+23]] = res[5]; - dr[ds[j+27]] = res[6]; - dr[ds[j+31]] = res[7]; - prod = _mm512_xor_si512( - _mm512_slli_epi32(prodh, 16), _mm512_and_si512(prodl, mask16)); - drv = _mm512_setr_epi64( - dr[ds[j+0]], - dr[ds[j+4]], - dr[ds[j+8]], - dr[ds[j+12]], - dr[ds[j+16]], - dr[ds[j+20]], - dr[ds[j+24]], - dr[ds[j+28]]); - resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask32)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+0]] = res[0]; - dr[ds[j+4]] = res[1]; - dr[ds[j+8]] = res[2]; - dr[ds[j+12]] = res[3]; - dr[ds[j+16]] = res[4]; - dr[ds[j+20]] = res[5]; - dr[ds[j+24]] = res[6]; - dr[ds[j+28]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+2]], - dr[ds[j+6]], - dr[ds[j+10]], - dr[ds[j+14]], - dr[ds[j+18]], - dr[ds[j+22]], - dr[ds[j+26]], - dr[ds[j+30]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(prod, 32)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+2]] = res[0]; - dr[ds[j+6]] = res[1]; - dr[ds[j+10]] = res[2]; - dr[ds[j+14]] = res[3]; - dr[ds[j+18]] = res[4]; - dr[ds[j+22]] = res[5]; - dr[ds[j+26]] = res[6]; - dr[ds[j+30]] = res[7]; - } -#elif defined HAVE_AVX2 - const uint16_t mul16 = (uint16_t)(fc - dr[i]); - mulv = _mm256_set1_epi16(mul16); - const len_t len = dts[LENGTH]; - const len_t os = len % 16; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 16) { - redv = _mm256_loadu_si256((__m256i*)(cfs+j)); - prodh = _mm256_mulhi_epu16(mulv, redv); - prodl = _mm256_mullo_epi16(mulv, redv); - prod = _mm256_xor_si256( - _mm256_and_si256(prodh, mask16h), _mm256_srli_epi32(prodl, 16)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+5]], - dr[ds[j+9]], - dr[ds[j+13]]); - resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask32)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+1]] = res[0]; - dr[ds[j+5]] = res[1]; - dr[ds[j+9]] = res[2]; - dr[ds[j+13]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+3]], - dr[ds[j+7]], - dr[ds[j+11]], - dr[ds[j+15]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(prod, 32)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+3]] = res[0]; - dr[ds[j+7]] = res[1]; - dr[ds[j+11]] = res[2]; - dr[ds[j+15]] = res[3]; - prod = _mm256_xor_si256( - _mm256_slli_epi32(prodh, 16), _mm256_and_si256(prodl, mask16)); - drv = _mm256_setr_epi64x( - dr[ds[j+0]], - dr[ds[j+4]], - dr[ds[j+8]], - dr[ds[j+12]]); - resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask32)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+0]] = res[0]; - dr[ds[j+4]] = res[1]; - dr[ds[j+8]] = res[2]; - dr[ds[j+12]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+2]], - dr[ds[j+6]], - dr[ds[j+10]], - dr[ds[j+14]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(prod, 32)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+2]] = res[0]; - dr[ds[j+6]] = res[1]; - dr[ds[j+10]] = res[2]; - dr[ds[j+14]] = res[3]; - } -#elif defined __aarch64__ - const len_t len = dts[LENGTH]; - const len_t os = len % 8; - const hm_t * const ds = dts + OFFSET; - const cf16_t mul16 = (cf16_t)(mod - dr[i]); - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 8) { - tmp[0] = (uint64_t)dr[ds[j]]; - tmp[1] = (uint64_t)dr[ds[j+1]]; - drv = vld1q_u64(tmp); - redv = vld1q_u16((cf16_t *)(cfs)+j); - - prodv = vmull_n_u16(vget_low_u16(redv), mul16); - resv = vaddw_u32(drv, vget_low_u32(prodv)); - vst1q_u64(tmp, resv); - dr[ds[j]] = (int64_t)tmp[0]; - dr[ds[j+1]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+2]]; - tmp[1] = (uint64_t)dr[ds[j+3]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_high_u32(prodv)); - vst1q_u64(tmp, resv); - dr[ds[j+2]] = (int64_t)tmp[0]; - dr[ds[j+3]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+4]]; - tmp[1] = (uint64_t)dr[ds[j+5]]; - drv = vld1q_u64(tmp); - - prodv = vmull_n_u16( vget_high_u16(redv), mul16); - resv = vaddw_u32(drv, vget_low_u32(prodv)); - vst1q_u64(tmp, resv); - dr[ds[j+4]] = (int64_t)tmp[0]; - dr[ds[j+5]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+6]]; - tmp[1] = (uint64_t)dr[ds[j+7]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_high_u32(prodv)); - vst1q_u64(tmp, resv); - dr[ds[j+6]] = (int64_t)tmp[0]; - dr[ds[j+7]] = (int64_t)tmp[1]; - } -#else - const len_t os = dts[PRELOOP]; - const len_t len = dts[LENGTH]; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] += mul * cfs[j]; - dr[ds[j+1]] += mul * cfs[j+1]; - dr[ds[j+2]] += mul * cfs[j+2]; - dr[ds[j+3]] += mul * cfs[j+3]; - } -#endif + add_mul_row_ff_16(dr, cfs, dts + OFFSET, dts[LENGTH], mul); dr[i] = 0; } if (k == 0) { diff --git a/src/neogb/la_ff_32.c b/src/neogb/la_ff_32.c index c5794cf8..819bb396 100644 --- a/src/neogb/la_ff_32.c +++ b/src/neogb/la_ff_32.c @@ -20,13 +20,7 @@ #include "data.h" #include "../msolve/streams.h" - -/* That's also enough if AVX512 is available on the system */ -#if defined HAVE_AVX2 -#include -#elif defined __aarch64__ -#include -#endif +#include "../msolve/cpu_features.h" static inline cf32_t *normalize_dense_matrix_row_ff_32( cf32_t *row, @@ -305,6 +299,259 @@ static inline cf32_t *multiply_sparse_matrix_row_ff_32( return row; } +static inline void add_mul_row_17_bit_scalar( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul + ) +{ + len_t j; + const len_t os = len % UNROLL; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += UNROLL) { + dr[ds[j]] += mul * cfs[j]; + dr[ds[j+1]] += mul * cfs[j+1]; + dr[ds[j+2]] += mul * cfs[j+2]; + dr[ds[j+3]] += mul * cfs[j+3]; + } +} + +#ifdef __aarch64__ +static inline void add_mul_row_17_bit_neon( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul + ) +{ + len_t j; + uint64_t tmp[2] __attribute__((aligned(32))); + uint32x4_t redv; + uint64x2_t drv, resv; + + const len_t os = len % 16; + const cf32_t mul32 = (cf32_t)mul; + const uint32x2_t mulv = vmov_n_u32(mul32); + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 16) { + tmp[0] = (uint64_t)dr[ds[j]]; + tmp[1] = (uint64_t)dr[ds[j+1]]; + drv = vld1q_u64(tmp); + redv = vld1q_u32((cf32_t *)(cfs)+j); + resv = vmlal_u32(drv, vget_low_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j]] = (int64_t)tmp[0]; + dr[ds[j+1]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+2]]; + tmp[1] = (uint64_t)dr[ds[j+3]]; + drv = vld1q_u64(tmp); + resv = vmlal_u32(drv, vget_high_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j+2]] = (int64_t)tmp[0]; + dr[ds[j+3]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+4]]; + tmp[1] = (uint64_t)dr[ds[j+5]]; + drv = vld1q_u64(tmp); + redv = vld1q_u32((cf32_t *)(cfs)+j+4); + resv = vmlal_u32(drv, vget_low_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j+4]] = (int64_t)tmp[0]; + dr[ds[j+5]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+6]]; + tmp[1] = (uint64_t)dr[ds[j+7]]; + drv = vld1q_u64(tmp); + resv = vmlal_u32(drv, vget_high_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j+6]] = (int64_t)tmp[0]; + dr[ds[j+7]] = (int64_t)tmp[1]; + + tmp[0] = (uint64_t)dr[ds[j+8]]; + tmp[1] = (uint64_t)dr[ds[j+9]]; + drv = vld1q_u64(tmp); + redv = vld1q_u32((cf32_t *)(cfs)+j+8); + resv = vmlal_u32(drv, vget_low_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j+8]] = (int64_t)tmp[0]; + dr[ds[j+9]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+10]]; + tmp[1] = (uint64_t)dr[ds[j+11]]; + drv = vld1q_u64(tmp); + resv = vmlal_u32(drv, vget_high_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j+10]] = (int64_t)tmp[0]; + dr[ds[j+11]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+12]]; + tmp[1] = (uint64_t)dr[ds[j+13]]; + drv = vld1q_u64(tmp); + redv = vld1q_u32((cf32_t *)(cfs)+j+12); + resv = vmlal_u32(drv, vget_low_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j+12]] = (int64_t)tmp[0]; + dr[ds[j+13]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+14]]; + tmp[1] = (uint64_t)dr[ds[j+15]]; + drv = vld1q_u64(tmp); + resv = vmlal_u32(drv, vget_high_u32(redv), mulv); + vst1q_u64(tmp, resv); + dr[ds[j+14]] = (int64_t)tmp[0]; + dr[ds[j+15]] = (int64_t)tmp[1]; + } +} +#endif + +#ifdef HAVE_AVX2_KERNELS +static TARGET_AVX2 void add_mul_row_17_bit_avx2( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul + ) +{ + len_t j; + int64_t res[4] __attribute__((aligned(32))); + __m256i redv, mulv, prodv, drv, resv; + + const len_t os = len % 8; + const uint32_t mul32 = (int32_t)mul; + mulv = _mm256_set1_epi32(mul32); + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 8) { + redv = _mm256_lddqu_si256((__m256i*)(cfs+j)); + drv = _mm256_setr_epi64x( + dr[ds[j+1]], + dr[ds[j+3]], + dr[ds[j+5]], + dr[ds[j+7]]); + /* first four mult-adds -- lower */ + prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); + resv = _mm256_add_epi64(drv, prodv); + _mm256_store_si256((__m256i*)(res), resv); + dr[ds[j+1]] = res[0]; + dr[ds[j+3]] = res[1]; + dr[ds[j+5]] = res[2]; + dr[ds[j+7]] = res[3]; + /* second four mult-adds -- higher */ + prodv = _mm256_mul_epu32(mulv, redv); + drv = _mm256_setr_epi64x( + dr[ds[j]], + dr[ds[j+2]], + dr[ds[j+4]], + dr[ds[j+6]]); + resv = _mm256_add_epi64(drv, prodv); + _mm256_store_si256((__m256i*)(res), resv); + dr[ds[j]] = res[0]; + dr[ds[j+2]] = res[1]; + dr[ds[j+4]] = res[2]; + dr[ds[j+6]] = res[3]; + } +} +#endif + +#ifdef HAVE_AVX512_KERNELS +static TARGET_AVX512 void add_mul_row_17_bit_avx512( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul + ) +{ + len_t j; + int64_t res[8] __attribute__((aligned(64))); + __m512i redv, mulv, prodv, drv, resv; + + const len_t os = len % 16; + const uint32_t mul32 = (int32_t)mul; + mulv = _mm512_set1_epi32(mul32); + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 16) { + redv = _mm512_loadu_si512((__m512i*)(cfs+j)); + drv = _mm512_setr_epi64( + dr[ds[j+1]], + dr[ds[j+3]], + dr[ds[j+5]], + dr[ds[j+7]], + dr[ds[j+9]], + dr[ds[j+11]], + dr[ds[j+13]], + dr[ds[j+15]]); + /* first four mult-adds -- lower */ + prodv = _mm512_mul_epu32(mulv, _mm512_srli_epi64(redv, 32)); + resv = _mm512_add_epi64(drv, prodv); + _mm512_store_si512((__m512*)(res), resv); + dr[ds[j+1]] = res[0]; + dr[ds[j+3]] = res[1]; + dr[ds[j+5]] = res[2]; + dr[ds[j+7]] = res[3]; + dr[ds[j+9]] = res[4]; + dr[ds[j+11]] = res[5]; + dr[ds[j+13]] = res[6]; + dr[ds[j+15]] = res[7]; + /* second four mult-adds -- higher */ + prodv = _mm512_mul_epu32(mulv, redv); + drv = _mm512_setr_epi64( + dr[ds[j]], + dr[ds[j+2]], + dr[ds[j+4]], + dr[ds[j+6]], + dr[ds[j+8]], + dr[ds[j+10]], + dr[ds[j+12]], + dr[ds[j+14]]); + resv = _mm512_add_epi64(drv, prodv); + _mm512_store_si512((__m512i*)(res), resv); + dr[ds[j]] = res[0]; + dr[ds[j+2]] = res[1]; + dr[ds[j+4]] = res[2]; + dr[ds[j+6]] = res[3]; + dr[ds[j+8]] = res[4]; + dr[ds[j+10]] = res[5]; + dr[ds[j+12]] = res[6]; + dr[ds[j+14]] = res[7]; + } +} +#endif + +/* dr[ds[j]] += mul * cfs[j] for 0 <= j < len */ +static inline void add_mul_row_17_bit( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul + ) +{ +#ifdef HAVE_AVX512_KERNELS + if (cpu_has_avx512()) { + add_mul_row_17_bit_avx512(dr, cfs, ds, len, mul); + return; + } +#endif +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + add_mul_row_17_bit_avx2(dr, cfs, ds, len, mul); + return; + } +#endif +#ifdef __aarch64__ + add_mul_row_17_bit_neon(dr, cfs, ds, len, mul); +#else + add_mul_row_17_bit_scalar(dr, cfs, ds, len, mul); +#endif +} + static hm_t *reduce_dense_row_by_known_pivots_sparse_17_bit( int64_t *dr, mat_t *mat, @@ -333,17 +580,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_17_bit( } else { rba = NULL; } -#if defined HAVE_AVX512_F - int64_t res[8] __attribute__((aligned(64))); - __m512i redv, mulv, prodv, drv, resv; -#elif defined HAVE_AVX2 - int64_t res[4] __attribute__((aligned(32))); - __m256i redv, mulv, prodv, drv, resv; -#elif defined __aarch64__ - uint64_t tmp[2] __attribute__((aligned(32))); - uint32x4_t redv; - uint64x2_t drv, resv; -#endif k = 0; for (i = dpiv; i < ncols; ++i) { @@ -371,184 +607,8 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_17_bit( } } cfs = mcf[dts[COEFFS]]; -#if defined HAVE_AVX512_F const len_t len = dts[LENGTH]; - const len_t os = len % 16; - const hm_t * const ds = dts + OFFSET; - const uint32_t mul32 = (int32_t)(mod - dr[i]); - mulv = _mm512_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 16) { - redv = _mm512_loadu_si512((__m512i*)(cfs+j)); - drv = _mm512_setr_epi64( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]], - dr[ds[j+9]], - dr[ds[j+11]], - dr[ds[j+13]], - dr[ds[j+15]]); - /* first four mult-adds -- lower */ - prodv = _mm512_mul_epu32(mulv, _mm512_srli_epi64(redv, 32)); - resv = _mm512_add_epi64(drv, prodv); - _mm512_store_si512((__m512*)(res), resv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - dr[ds[j+9]] = res[4]; - dr[ds[j+11]] = res[5]; - dr[ds[j+13]] = res[6]; - dr[ds[j+15]] = res[7]; - /* second four mult-adds -- higher */ - prodv = _mm512_mul_epu32(mulv, redv); - drv = _mm512_setr_epi64( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]], - dr[ds[j+8]], - dr[ds[j+10]], - dr[ds[j+12]], - dr[ds[j+14]]); - resv = _mm512_add_epi64(drv, prodv); - _mm512_store_si512((__m512i*)(res), resv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - dr[ds[j+8]] = res[4]; - dr[ds[j+10]] = res[5]; - dr[ds[j+12]] = res[6]; - dr[ds[j+14]] = res[7]; - } -#elif defined HAVE_AVX2 - const len_t len = dts[LENGTH]; - const len_t os = len % 8; - const hm_t * const ds = dts + OFFSET; - const uint32_t mul32 = (int32_t)(mod - dr[i]); - mulv = _mm256_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 8) { - redv = _mm256_lddqu_si256((__m256i*)(cfs+j)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]]); - /* first four mult-adds -- lower */ - prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); - resv = _mm256_add_epi64(drv, prodv); - _mm256_store_si256((__m256i*)(res), resv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - /* second four mult-adds -- higher */ - prodv = _mm256_mul_epu32(mulv, redv); - drv = _mm256_setr_epi64x( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]]); - resv = _mm256_add_epi64(drv, prodv); - _mm256_store_si256((__m256i*)(res), resv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - } -#elif defined __aarch64__ - const len_t len = dts[LENGTH]; - const len_t os = len % 16; - const hm_t * const ds = dts + OFFSET; - const cf32_t mul32 = (cf32_t)(mod - dr[i]); - const uint32x2_t mulv = vmov_n_u32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 16) { - tmp[0] = (uint64_t)dr[ds[j]]; - tmp[1] = (uint64_t)dr[ds[j+1]]; - drv = vld1q_u64(tmp); - redv = vld1q_u32((cf32_t *)(cfs)+j); - resv = vmlal_u32(drv, vget_low_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j]] = (int64_t)tmp[0]; - dr[ds[j+1]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+2]]; - tmp[1] = (uint64_t)dr[ds[j+3]]; - drv = vld1q_u64(tmp); - resv = vmlal_u32(drv, vget_high_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j+2]] = (int64_t)tmp[0]; - dr[ds[j+3]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+4]]; - tmp[1] = (uint64_t)dr[ds[j+5]]; - drv = vld1q_u64(tmp); - redv = vld1q_u32((cf32_t *)(cfs)+j+4); - resv = vmlal_u32(drv, vget_low_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j+4]] = (int64_t)tmp[0]; - dr[ds[j+5]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+6]]; - tmp[1] = (uint64_t)dr[ds[j+7]]; - drv = vld1q_u64(tmp); - resv = vmlal_u32(drv, vget_high_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j+6]] = (int64_t)tmp[0]; - dr[ds[j+7]] = (int64_t)tmp[1]; - - tmp[0] = (uint64_t)dr[ds[j+8]]; - tmp[1] = (uint64_t)dr[ds[j+9]]; - drv = vld1q_u64(tmp); - redv = vld1q_u32((cf32_t *)(cfs)+j+8); - resv = vmlal_u32(drv, vget_low_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j+8]] = (int64_t)tmp[0]; - dr[ds[j+9]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+10]]; - tmp[1] = (uint64_t)dr[ds[j+11]]; - drv = vld1q_u64(tmp); - resv = vmlal_u32(drv, vget_high_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j+10]] = (int64_t)tmp[0]; - dr[ds[j+11]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+12]]; - tmp[1] = (uint64_t)dr[ds[j+13]]; - drv = vld1q_u64(tmp); - redv = vld1q_u32((cf32_t *)(cfs)+j+12); - resv = vmlal_u32(drv, vget_low_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j+12]] = (int64_t)tmp[0]; - dr[ds[j+13]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+14]]; - tmp[1] = (uint64_t)dr[ds[j+15]]; - drv = vld1q_u64(tmp); - resv = vmlal_u32(drv, vget_high_u32(redv), mulv); - vst1q_u64(tmp, resv); - dr[ds[j+14]] = (int64_t)tmp[0]; - dr[ds[j+15]] = (int64_t)tmp[1]; - } -#else - const len_t os = dts[PRELOOP]; - const len_t len = dts[LENGTH]; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] += mul * cfs[j]; - dr[ds[j+1]] += mul * cfs[j+1]; - dr[ds[j+2]] += mul * cfs[j+2]; - dr[ds[j+3]] += mul * cfs[j+3]; - } -#endif + add_mul_row_17_bit(dr, cfs, dts + OFFSET, len, mul); dr[i] = 0; st->application_nr_mult += len / 1000.0; st->application_nr_add += len / 1000.0; @@ -665,6 +725,265 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_17_bit( return row; } +static inline void sub_mul_row_31_bit_scalar( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul, + const int64_t mod2 + ) +{ + len_t j; + const len_t os = len % UNROLL; + for (j = 0; j < os; ++j) { + dr[ds[j]] -= mul * cfs[j]; + dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; + } + for (; j < len; j += UNROLL) { + dr[ds[j]] -= mul * cfs[j]; + dr[ds[j+1]] -= mul * cfs[j+1]; + dr[ds[j+2]] -= mul * cfs[j+2]; + dr[ds[j+3]] -= mul * cfs[j+3]; + dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; + dr[ds[j+1]] += (dr[ds[j+1]] >> 63) & mod2; + dr[ds[j+2]] += (dr[ds[j+2]] >> 63) & mod2; + dr[ds[j+3]] += (dr[ds[j+3]] >> 63) & mod2; + } +} + +#ifdef __aarch64__ +static inline void sub_mul_row_31_bit_neon( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul, + const int64_t mod2 + ) +{ + len_t j; + const int64x2_t mod2v = vmovq_n_s64(mod2); + int64_t tmp[2] __attribute__((aligned(32))); + int32x4_t redv; + int64x2_t drv, mask, resv; + + const len_t os = len % 8; + const int32_t mul32 = (int32_t)mul; + const int32x2_t mulv = vmov_n_s32(mul32); + for (j = 0; j < os; ++j) { + dr[ds[j]] -= mul * cfs[j]; + dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; + } + for (; j < len; j += 8) { + tmp[0] = dr[ds[j]]; + tmp[1] = dr[ds[j+1]]; + drv = vld1q_s64(tmp); + redv = vld1q_s32((int32_t *)(cfs)+j); + /* multiply and subtract */ + resv = vmlsl_s32(drv, vget_low_s32(redv), mulv); + mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); + resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); + vst1q_s64(tmp, resv); + dr[ds[j]] = tmp[0]; + dr[ds[j+1]] = tmp[1]; + tmp[0] = dr[ds[j+2]]; + tmp[1] = dr[ds[j+3]]; + drv = vld1q_s64(tmp); + resv = vmlsl_s32(drv, vget_high_s32(redv), mulv); + mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); + resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); + vst1q_s64(tmp, resv); + dr[ds[j+2]] = tmp[0]; + dr[ds[j+3]] = tmp[1]; + tmp[0] = dr[ds[j+4]]; + tmp[1] = dr[ds[j+5]]; + drv = vld1q_s64(tmp); + redv = vld1q_s32((int32_t *)(cfs)+j+4); + /* multiply and subtract */ + resv = vmlsl_s32(drv, vget_low_s32(redv), mulv); + mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); + resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); + vst1q_s64(tmp, resv); + dr[ds[j+4]] = tmp[0]; + dr[ds[j+5]] = tmp[1]; + tmp[0] = dr[ds[j+6]]; + tmp[1] = dr[ds[j+7]]; + drv = vld1q_s64(tmp); + resv = vmlsl_s32(drv, vget_high_s32(redv), mulv); + mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); + resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); + vst1q_s64(tmp, resv); + dr[ds[j+6]] = tmp[0]; + dr[ds[j+7]] = tmp[1]; + } +} +#endif + +#ifdef HAVE_AVX2_KERNELS +static TARGET_AVX2 void sub_mul_row_31_bit_avx2( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul, + const int64_t mod2 + ) +{ + len_t j; + int64_t res[4] __attribute__((aligned(32))); + __m256i cmpv, redv, drv, mulv, prodv, resv, rresv; + __m256i zerov= _mm256_set1_epi64x(0); + __m256i mod2v = _mm256_set1_epi64x(mod2); + + const len_t os = len % 8; + const uint32_t mul32 = (uint32_t)mul; + mulv = _mm256_set1_epi32(mul32); + for (j = 0; j < os; ++j) { + dr[ds[j]] -= mul * cfs[j]; + dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; + } + for (; j < len; j += 8) { + redv = _mm256_loadu_si256((__m256i*)(cfs+j)); + drv = _mm256_setr_epi64x( + dr[ds[j+1]], + dr[ds[j+3]], + dr[ds[j+5]], + dr[ds[j+7]]); + /* first four mult-adds -- lower */ + prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); + resv = _mm256_sub_epi64(drv, prodv); + cmpv = _mm256_cmpgt_epi64(zerov, resv); + rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); + _mm256_store_si256((__m256i*)(res), rresv); + dr[ds[j+1]] = res[0]; + dr[ds[j+3]] = res[1]; + dr[ds[j+5]] = res[2]; + dr[ds[j+7]] = res[3]; + /* second four mult-adds -- higher */ + prodv = _mm256_mul_epu32(mulv, redv); + drv = _mm256_setr_epi64x( + dr[ds[j]], + dr[ds[j+2]], + dr[ds[j+4]], + dr[ds[j+6]]); + resv = _mm256_sub_epi64(drv, prodv); + cmpv = _mm256_cmpgt_epi64(zerov, resv); + rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); + _mm256_store_si256((__m256i*)(res), rresv); + dr[ds[j]] = res[0]; + dr[ds[j+2]] = res[1]; + dr[ds[j+4]] = res[2]; + dr[ds[j+6]] = res[3]; + } +} +#endif + +#ifdef HAVE_AVX512_KERNELS +static TARGET_AVX512 void sub_mul_row_31_bit_avx512( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul, + const int64_t mod2 + ) +{ + len_t j; + int64_t res[8] __attribute__((aligned(64))); + __m512i redv, drv, mulv, prodv, resv, rresv; + __mmask8 cmpv; + __m512i zerov = _mm512_set1_epi64(0); + __m512i mod2v = _mm512_set1_epi64(mod2); + + const len_t os = len % 16; + const uint32_t mul32 = (int32_t)mul; + mulv = _mm512_set1_epi32(mul32); + for (j = 0; j < os; ++j) { + dr[ds[j]] -= mul * cfs[j]; + dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; + } + for (; j < len; j += 16) { + redv = _mm512_loadu_si512((__m512i*)(cfs+j)); + drv = _mm512_setr_epi64( + dr[ds[j+1]], + dr[ds[j+3]], + dr[ds[j+5]], + dr[ds[j+7]], + dr[ds[j+9]], + dr[ds[j+11]], + dr[ds[j+13]], + dr[ds[j+15]]); + /* first four mult-adds -- lower */ + prodv = _mm512_mul_epu32(mulv, _mm512_srli_epi64(redv, 32)); + resv = _mm512_sub_epi64(drv, prodv); + cmpv = _mm512_cmpgt_epi64_mask(zerov, resv); + rresv = _mm512_mask_add_epi64(resv, cmpv, resv, mod2v); + _mm512_store_si512((__m512*)(res), rresv); + dr[ds[j+1]] = res[0]; + dr[ds[j+3]] = res[1]; + dr[ds[j+5]] = res[2]; + dr[ds[j+7]] = res[3]; + dr[ds[j+9]] = res[4]; + dr[ds[j+11]] = res[5]; + dr[ds[j+13]] = res[6]; + dr[ds[j+15]] = res[7]; + /* second four mult-adds -- higher */ + prodv = _mm512_mul_epu32(mulv, redv); + drv = _mm512_setr_epi64( + dr[ds[j]], + dr[ds[j+2]], + dr[ds[j+4]], + dr[ds[j+6]], + dr[ds[j+8]], + dr[ds[j+10]], + dr[ds[j+12]], + dr[ds[j+14]]); + resv = _mm512_sub_epi64(drv, prodv); + cmpv = _mm512_cmpgt_epi64_mask(zerov, resv); + rresv = _mm512_mask_add_epi64(resv, cmpv, resv, mod2v); + _mm512_store_si512((__m512i*)(res), rresv); + dr[ds[j]] = res[0]; + dr[ds[j+2]] = res[1]; + dr[ds[j+4]] = res[2]; + dr[ds[j+6]] = res[3]; + dr[ds[j+8]] = res[4]; + dr[ds[j+10]] = res[5]; + dr[ds[j+12]] = res[6]; + dr[ds[j+14]] = res[7]; + } +} +#endif + +/* dr[ds[j]] -= mul * cfs[j] for 0 <= j < len, keeping dr[ds[j]] in [0, mod2) */ +static inline void sub_mul_row_31_bit( + int64_t *dr, + const cf32_t *cfs, + const hm_t *ds, + const len_t len, + const int64_t mul, + const int64_t mod2 + ) +{ +#ifdef HAVE_AVX512_KERNELS + if (cpu_has_avx512()) { + sub_mul_row_31_bit_avx512(dr, cfs, ds, len, mul, mod2); + return; + } +#endif +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + sub_mul_row_31_bit_avx2(dr, cfs, ds, len, mul, mod2); + return; + } +#endif +#ifdef __aarch64__ + sub_mul_row_31_bit_neon(dr, cfs, ds, len, mul, mod2); +#else + sub_mul_row_31_bit_scalar(dr, cfs, ds, len, mul, mod2); +#endif +} + static hm_t *reduce_dense_row_by_known_pivots_sparse_up_to_ff_31_bit( int64_t *dr, bs_t *sat, @@ -684,12 +1003,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_up_to_ff_31_bit( int64_t np = -1; const int64_t mod = (int64_t)st->fc; const int64_t mod2 = (int64_t)st->fc * st->fc; -#ifdef HAVE_AVX2 - int64_t res[4] __attribute__((aligned(32))); - __m256i cmpv, redv, drv, mulv, prodv, resv, rresv; - __m256i zerov= _mm256_set1_epi64x(0); - __m256i mod2v = _mm256_set1_epi64x(mod2); -#endif for (i = dpiv; i < end; ++i) { if (dr[i] != 0) { @@ -709,68 +1022,15 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_up_to_ff_31_bit( const int64_t mul = (int64_t)dr[i]; dts = pivs[i]; cfs = bs->cf_32[dts[COEFFS]]; -#ifdef HAVE_AVX2 - const len_t len = dts[LENGTH]; - const len_t os = len % 8; - const hm_t * const ds = dts + OFFSET; - const uint32_t mul32 = (uint32_t)(dr[i]); - mulv = _mm256_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += 8) { - redv = _mm256_loadu_si256((__m256i*)(cfs+j)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]]); - /* first four mult-adds -- lower */ - prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - /* second four mult-adds -- higher */ - prodv = _mm256_mul_epu32(mulv, redv); - drv = _mm256_setr_epi64x( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]]); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - } -#else - const len_t os = dts[PRELOOP]; const len_t len = dts[LENGTH]; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j+1]] -= mul * cfs[j+1]; - dr[ds[j+2]] -= mul * cfs[j+2]; - dr[ds[j+3]] -= mul * cfs[j+3]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - dr[ds[j+1]] += (dr[ds[j+1]] >> 63) & mod2; - dr[ds[j+2]] += (dr[ds[j+2]] >> 63) & mod2; - dr[ds[j+3]] += (dr[ds[j+3]] >> 63) & mod2; - } +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + sub_mul_row_31_bit_avx2(dr, cfs, dts + OFFSET, len, mul, mod2); + } else #endif + { + sub_mul_row_31_bit_scalar(dr, cfs, dts + OFFSET, len, mul, mod2); + } dr[i] = 0; st->application_nr_mult += len / 1000.0; st->application_nr_add += len / 1000.0; @@ -826,12 +1086,6 @@ static hm_t *sba_reduce_dense_row_by_known_pivots_sparse_31_bit( const int64_t mod = (int64_t)st->fc; const int64_t mod2 = (int64_t)st->fc * st->fc; const len_t nc = smat->nc; -#ifdef HAVE_AVX2 - int64_t res[4] __attribute__((aligned(32))); - __m256i cmpv, redv, drv, mulv, prodv, resv, rresv; - __m256i zerov= _mm256_set1_epi64x(0); - __m256i mod2v = _mm256_set1_epi64x(mod2); -#endif k = 0; for (i = dpiv; i < nc; ++i) { @@ -863,68 +1117,15 @@ static hm_t *sba_reduce_dense_row_by_known_pivots_sparse_31_bit( dts = pivs[i]; cfs = smat->cc32[dts[SM_CFS]]; -#ifdef HAVE_AVX2 - const len_t len = dts[SM_LEN]; - const len_t os = len % 8; - const hm_t * const ds = dts + SM_OFFSET; - const uint32_t mul32 = (uint32_t)(dr[i]); - mulv = _mm256_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += 8) { - redv = _mm256_loadu_si256((__m256i*)(cfs+j)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]]); - /* first four mult-adds -- lower */ - prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - /* second four mult-adds -- higher */ - prodv = _mm256_mul_epu32(mulv, redv); - drv = _mm256_setr_epi64x( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]]); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - } -#else - const len_t os = dts[SM_PRE]; const len_t len = dts[SM_LEN]; - const hm_t * const ds = dts + SM_OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j+1]] -= mul * cfs[j+1]; - dr[ds[j+2]] -= mul * cfs[j+2]; - dr[ds[j+3]] -= mul * cfs[j+3]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - dr[ds[j+1]] += (dr[ds[j+1]] >> 63) & mod2; - dr[ds[j+2]] += (dr[ds[j+2]] >> 63) & mod2; - dr[ds[j+3]] += (dr[ds[j+3]] >> 63) & mod2; - } +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + sub_mul_row_31_bit_avx2(dr, cfs, dts + SM_OFFSET, len, mul, mod2); + } else #endif + { + sub_mul_row_31_bit_scalar(dr, cfs, dts + SM_OFFSET, len, mul, mod2); + } dr[i] = 0; st->application_nr_mult += len / 1000.0; st->application_nr_add += len / 1000.0; @@ -993,23 +1194,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_31_bit( } else { rba = NULL; } -#if defined HAVE_AVX512_F - int64_t res[8] __attribute__((aligned(64))); - __m512i redv, drv, mulv, prodv, resv, rresv; - __mmask8 cmpv; - __m512i zerov = _mm512_set1_epi64(0); - __m512i mod2v = _mm512_set1_epi64(mod2); -#elif defined HAVE_AVX2 - int64_t res[4] __attribute__((aligned(32))); - __m256i cmpv, redv, drv, mulv, prodv, resv, rresv; - __m256i zerov= _mm256_set1_epi64x(0); - __m256i mod2v = _mm256_set1_epi64x(mod2); -#elif defined __aarch64__ - const int64x2_t mod2v = vmovq_n_s64(mod2); - int64_t tmp[2] __attribute__((aligned(32))); - int32x4_t redv; - int64x2_t drv, mask, resv; -#endif k = 0; for (i = dpiv; i < ncols; ++i) { @@ -1038,180 +1222,8 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_31_bit( } } cfs = mcf[dts[COEFFS]]; -#if defined HAVE_AVX512_F - const len_t len = dts[LENGTH]; - const len_t os = len % 16; - const hm_t * const ds = dts + OFFSET; - const uint32_t mul32 = (int32_t)(dr[i]); - mulv = _mm512_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += 16) { - redv = _mm512_loadu_si512((__m512i*)(cfs+j)); - drv = _mm512_setr_epi64( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]], - dr[ds[j+9]], - dr[ds[j+11]], - dr[ds[j+13]], - dr[ds[j+15]]); - /* first four mult-adds -- lower */ - prodv = _mm512_mul_epu32(mulv, _mm512_srli_epi64(redv, 32)); - resv = _mm512_sub_epi64(drv, prodv); - cmpv = _mm512_cmpgt_epi64_mask(zerov, resv); - rresv = _mm512_mask_add_epi64(resv, cmpv, resv, mod2v); - _mm512_store_si512((__m512*)(res), rresv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - dr[ds[j+9]] = res[4]; - dr[ds[j+11]] = res[5]; - dr[ds[j+13]] = res[6]; - dr[ds[j+15]] = res[7]; - /* second four mult-adds -- higher */ - prodv = _mm512_mul_epu32(mulv, redv); - drv = _mm512_setr_epi64( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]], - dr[ds[j+8]], - dr[ds[j+10]], - dr[ds[j+12]], - dr[ds[j+14]]); - resv = _mm512_sub_epi64(drv, prodv); - cmpv = _mm512_cmpgt_epi64_mask(zerov, resv); - rresv = _mm512_mask_add_epi64(resv, cmpv, resv, mod2v); - _mm512_store_si512((__m512i*)(res), rresv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - dr[ds[j+8]] = res[4]; - dr[ds[j+10]] = res[5]; - dr[ds[j+12]] = res[6]; - dr[ds[j+14]] = res[7]; - } -#elif defined HAVE_AVX2 const len_t len = dts[LENGTH]; - const len_t os = len % 8; - const hm_t * const ds = dts + OFFSET; - const uint32_t mul32 = (uint32_t)(dr[i]); - mulv = _mm256_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += 8) { - redv = _mm256_loadu_si256((__m256i*)(cfs+j)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]]); - /* first four mult-adds -- lower */ - prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - /* second four mult-adds -- higher */ - prodv = _mm256_mul_epu32(mulv, redv); - drv = _mm256_setr_epi64x( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]]); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - } -#elif defined __aarch64__ - const len_t len = dts[LENGTH]; - const len_t os = len % 8; - const hm_t * const ds = dts + OFFSET; - const int32_t mul32 = (int32_t)(dr[i]); - const int32x2_t mulv = vmov_n_s32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += 8) { - tmp[0] = dr[ds[j]]; - tmp[1] = dr[ds[j+1]]; - drv = vld1q_s64(tmp); - redv = vld1q_s32((int32_t *)(cfs)+j); - /* multiply and subtract */ - resv = vmlsl_s32(drv, vget_low_s32(redv), mulv); - mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); - resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); - vst1q_s64(tmp, resv); - dr[ds[j]] = tmp[0]; - dr[ds[j+1]] = tmp[1]; - tmp[0] = dr[ds[j+2]]; - tmp[1] = dr[ds[j+3]]; - drv = vld1q_s64(tmp); - resv = vmlsl_s32(drv, vget_high_s32(redv), mulv); - mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); - resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); - vst1q_s64(tmp, resv); - dr[ds[j+2]] = tmp[0]; - dr[ds[j+3]] = tmp[1]; - tmp[0] = dr[ds[j+4]]; - tmp[1] = dr[ds[j+5]]; - drv = vld1q_s64(tmp); - redv = vld1q_s32((int32_t *)(cfs)+j+4); - /* multiply and subtract */ - resv = vmlsl_s32(drv, vget_low_s32(redv), mulv); - mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); - resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); - vst1q_s64(tmp, resv); - dr[ds[j+4]] = tmp[0]; - dr[ds[j+5]] = tmp[1]; - tmp[0] = dr[ds[j+6]]; - tmp[1] = dr[ds[j+7]]; - drv = vld1q_s64(tmp); - resv = vmlsl_s32(drv, vget_high_s32(redv), mulv); - mask = vreinterpretq_s64_u64(vcltzq_s64(resv)); - resv = vaddq_s64(resv, vandq_s64(mask, mod2v)); - vst1q_s64(tmp, resv); - dr[ds[j+6]] = tmp[0]; - dr[ds[j+7]] = tmp[1]; - } - -#else - const len_t os = dts[PRELOOP]; - const len_t len = dts[LENGTH]; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j+1]] -= mul * cfs[j+1]; - dr[ds[j+2]] -= mul * cfs[j+2]; - dr[ds[j+3]] -= mul * cfs[j+3]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - dr[ds[j+1]] += (dr[ds[j+1]] >> 63) & mod2; - dr[ds[j+2]] += (dr[ds[j+2]] >> 63) & mod2; - dr[ds[j+3]] += (dr[ds[j+3]] >> 63) & mod2; - } -#endif + sub_mul_row_31_bit(dr, cfs, dts + OFFSET, len, mul, mod2); dr[i] = 0; st->application_nr_mult += len / 1000.0; st->application_nr_add += len / 1000.0; @@ -1265,12 +1277,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_sat_ff_31_bit( int64_t np = -1; const int64_t mod = (int64_t)st->fc; const int64_t mod2 = (int64_t)st->fc * st->fc; -#ifdef HAVE_AVX2 - int64_t res[4] __attribute__((aligned(32))); - __m256i cmpv, redv, drv, mulv, prodv, resv, rresv; - __m256i zerov= _mm256_set1_epi64x(0); - __m256i mod2v = _mm256_set1_epi64x(mod2); -#endif k = 0; for (i = dpiv; i < ncols; ++i) { @@ -1294,127 +1300,17 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_sat_ff_31_bit( dtsm = mulh[dts[MULT]]; cfsm = mulcf[dtsm[COEFFS]]; cfs = pivcf[dts[COEFFS]]; -#ifdef HAVE_AVX2 const len_t len = dts[LENGTH]; - const len_t os = len % 8; - const hm_t * const ds = dts + OFFSET; - const uint32_t mul32 = (uint32_t)(dr[i]); - mulv = _mm256_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += 8) { - redv = _mm256_loadu_si256((__m256i*)(cfs+j)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]]); - /* first four mult-adds -- lower */ - prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - /* second four mult-adds -- higher */ - prodv = _mm256_mul_epu32(mulv, redv); - drv = _mm256_setr_epi64x( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]]); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - } - const len_t lenm = dtsm[LENGTH]; - const len_t osm = lenm % 8; - const hm_t * const dsm = dtsm + OFFSET; - /* const uint32_t mulm32 = (uint32_t)(drm[i]); */ - mulv = _mm256_set1_epi32(mul32); - for (j = 0; j < osm; ++j) { - drm[dsm[j]] -= mul * cfsm[j]; - drm[dsm[j]] += (drm[dsm[j]] >> 63) & mod2; - } - for (; j < lenm; j += 8) { - redv = _mm256_loadu_si256((__m256i*)(cfsm+j)); - drv = _mm256_setr_epi64x( - drm[dsm[j+1]], - drm[dsm[j+3]], - drm[dsm[j+5]], - drm[dsm[j+7]]); - /* first four mult-adds -- lower */ - prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - drm[dsm[j+1]] = res[0]; - drm[dsm[j+3]] = res[1]; - drm[dsm[j+5]] = res[2]; - drm[dsm[j+7]] = res[3]; - /* second four mult-adds -- higher */ - prodv = _mm256_mul_epu32(mulv, redv); - drv = _mm256_setr_epi64x( - drm[dsm[j]], - drm[dsm[j+2]], - drm[dsm[j+4]], - drm[dsm[j+6]]); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - drm[dsm[j]] = res[0]; - drm[dsm[j+2]] = res[1]; - drm[dsm[j+4]] = res[2]; - drm[dsm[j+6]] = res[3]; - } -#else - const len_t os = dts[PRELOOP]; - const len_t len = dts[LENGTH]; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j+1]] -= mul * cfs[j+1]; - dr[ds[j+2]] -= mul * cfs[j+2]; - dr[ds[j+3]] -= mul * cfs[j+3]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - dr[ds[j+1]] += (dr[ds[j+1]] >> 63) & mod2; - dr[ds[j+2]] += (dr[ds[j+2]] >> 63) & mod2; - dr[ds[j+3]] += (dr[ds[j+3]] >> 63) & mod2; - } - const len_t osm = dtsm[PRELOOP]; - const len_t lenm = dtsm[LENGTH]; - const hm_t * const dsm = dtsm + OFFSET; - for (j = 0; j < osm; ++j) { - drm[dsm[j]] -= mul * cfsm[j]; - drm[dsm[j]] += (drm[dsm[j]] >> 63) & mod2; - } - for (; j < lenm; j += UNROLL) { - drm[dsm[j]] -= mul * cfsm[j]; - drm[dsm[j+1]] -= mul * cfsm[j+1]; - drm[dsm[j+2]] -= mul * cfsm[j+2]; - drm[dsm[j+3]] -= mul * cfsm[j+3]; - drm[dsm[j]] += (drm[dsm[j]] >> 63) & mod2; - drm[dsm[j+1]] += (drm[dsm[j+1]] >> 63) & mod2; - drm[dsm[j+2]] += (drm[dsm[j+2]] >> 63) & mod2; - drm[dsm[j+3]] += (drm[dsm[j+3]] >> 63) & mod2; - } +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + sub_mul_row_31_bit_avx2(dr, cfs, dts + OFFSET, len, mul, mod2); + sub_mul_row_31_bit_avx2(drm, cfsm, dtsm + OFFSET, dtsm[LENGTH], mul, mod2); + } else #endif + { + sub_mul_row_31_bit_scalar(dr, cfs, dts + OFFSET, len, mul, mod2); + sub_mul_row_31_bit_scalar(drm, cfsm, dtsm + OFFSET, dtsm[LENGTH], mul, mod2); + } dr[i] = 0; st->application_nr_mult += len / 1000.0; st->application_nr_add += len / 1000.0; @@ -1500,12 +1396,6 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_31_bit( const len_t ncols = mat->nc; const len_t ncl = mat->ncl; cf32_t * const * const mcf = mat->cf_32; -#ifdef HAVE_AVX2 - int64_t res[4] __attribute__((aligned(32))); - __m256i cmpv, redv, drv, mulv, prodv, resv, rresv; - __m256i zerov = _mm256_set1_epi64x(0); - __m256i mod2v = _mm256_set1_epi64x(mod2); -#endif k = 0; for (i = dpiv; i < ncols; ++i) { @@ -1534,68 +1424,15 @@ static hm_t *trace_reduce_dense_row_by_known_pivots_sparse_31_bit( cfs = mcf[dts[COEFFS]]; } -#ifdef HAVE_AVX2 const len_t len = dts[LENGTH]; - const len_t os = len % 8; - const hm_t * const ds = dts + OFFSET; - const uint32_t mul32 = (uint32_t)(dr[i]); - mulv = _mm256_set1_epi32(mul32); - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += 8) { - redv = _mm256_loadu_si256((__m256i*)(cfs+j)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+3]], - dr[ds[j+5]], - dr[ds[j+7]]); - /* first four mult-adds -- lower */ - prodv = _mm256_mul_epu32(mulv, _mm256_srli_epi64(redv, 32)); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j+1]] = res[0]; - dr[ds[j+3]] = res[1]; - dr[ds[j+5]] = res[2]; - dr[ds[j+7]] = res[3]; - /* second four mult-adds -- higher */ - prodv = _mm256_mul_epu32(mulv, redv); - drv = _mm256_setr_epi64x( - dr[ds[j]], - dr[ds[j+2]], - dr[ds[j+4]], - dr[ds[j+6]]); - resv = _mm256_sub_epi64(drv, prodv); - cmpv = _mm256_cmpgt_epi64(zerov, resv); - rresv = _mm256_add_epi64(resv, _mm256_and_si256(cmpv, mod2v)); - _mm256_store_si256((__m256i*)(res), rresv); - dr[ds[j]] = res[0]; - dr[ds[j+2]] = res[1]; - dr[ds[j+4]] = res[2]; - dr[ds[j+6]] = res[3]; - } -#else - const len_t os = dts[PRELOOP]; - const len_t len = dts[LENGTH]; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] -= mul * cfs[j]; - dr[ds[j+1]] -= mul * cfs[j+1]; - dr[ds[j+2]] -= mul * cfs[j+2]; - dr[ds[j+3]] -= mul * cfs[j+3]; - dr[ds[j]] += (dr[ds[j]] >> 63) & mod2; - dr[ds[j+1]] += (dr[ds[j+1]] >> 63) & mod2; - dr[ds[j+2]] += (dr[ds[j+2]] >> 63) & mod2; - dr[ds[j+3]] += (dr[ds[j+3]] >> 63) & mod2; - } +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + sub_mul_row_31_bit_avx2(dr, cfs, dts + OFFSET, len, mul, mod2); + } else #endif + { + sub_mul_row_31_bit_scalar(dr, cfs, dts + OFFSET, len, mul, mod2); + } dr[i] = 0; st->trace_nr_mult += len / 1000.0; st->trace_nr_add += len / 1000.0; diff --git a/src/neogb/la_ff_8.c b/src/neogb/la_ff_8.c index 2f740081..871695da 100644 --- a/src/neogb/la_ff_8.c +++ b/src/neogb/la_ff_8.c @@ -21,12 +21,7 @@ #include "data.h" #include "../msolve/streams.h" -/* That's also enough if AVX512 is avaialable on the system */ -#if defined HAVE_AVX2 -#include -#elif defined __aarch64__ -#include -#endif +#include "../msolve/cpu_features.h" static inline cf8_t *normalize_dense_matrix_row_ff_8( cf8_t *row, @@ -94,6 +89,455 @@ static inline cf8_t *normalize_sparse_matrix_row_ff_8( return row; } +static inline void add_mul_row_ff_8_scalar( + int64_t *dr, + const cf8_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + const len_t os = len % UNROLL; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += UNROLL) { + dr[ds[j]] += mul * cfs[j]; + dr[ds[j+1]] += mul * cfs[j+1]; + dr[ds[j+2]] += mul * cfs[j+2]; + dr[ds[j+3]] += mul * cfs[j+3]; + } +} + +#ifdef __aarch64__ +static inline void add_mul_row_ff_8_neon( + int64_t *dr, + const cf8_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + uint64_t tmp[2] __attribute__((aligned(32))); + uint16x8_t prodv; + uint32x4_t prodvl, prodvh; + uint8x16_t redv; + uint64x2_t drv, resv; + + const len_t os = len % 16; + const cf8_t mul8 = (cf8_t)mul; + const uint8x8_t mulv = vmov_n_u8(mul8); + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 16) { + redv = vld1q_u8(cfs+j); + + prodv = vmull_u8(vget_low_u8(redv), mulv); + prodvl = vmovl_u16(vget_low_u16(prodv)); + tmp[0] = (uint64_t)dr[ds[j]]; + tmp[1] = (uint64_t)dr[ds[j+1]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_low_u32(prodvl)); + vst1q_u64(tmp, resv); + dr[ds[j]] = (int64_t)tmp[0]; + dr[ds[j+1]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+2]]; + tmp[1] = (uint64_t)dr[ds[j+3]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_high_u32(prodvl)); + vst1q_u64(tmp, resv); + dr[ds[j+2]] = (int64_t)tmp[0]; + dr[ds[j+3]] = (int64_t)tmp[1]; + prodvh = vmovl_u16(vget_high_u16(prodv)); + tmp[0] = (uint64_t)dr[ds[j+4]]; + tmp[1] = (uint64_t)dr[ds[j+5]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_low_u32(prodvh)); + vst1q_u64(tmp, resv); + dr[ds[j+4]] = (int64_t)tmp[0]; + dr[ds[j+5]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+6]]; + tmp[1] = (uint64_t)dr[ds[j+7]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_high_u32(prodvh)); + vst1q_u64(tmp, resv); + dr[ds[j+6]] = (int64_t)tmp[0]; + dr[ds[j+7]] = (int64_t)tmp[1]; + + prodv = vmull_u8(vget_high_u8(redv), mulv); + prodvl = vmovl_u16(vget_low_u16(prodv)); + tmp[0] = (uint64_t)dr[ds[j+8]]; + tmp[1] = (uint64_t)dr[ds[j+9]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_low_u32(prodvl)); + vst1q_u64(tmp, resv); + dr[ds[j+8]] = (int64_t)tmp[0]; + dr[ds[j+9]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+10]]; + tmp[1] = (uint64_t)dr[ds[j+11]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_high_u32(prodvl)); + vst1q_u64(tmp, resv); + dr[ds[j+10]] = (int64_t)tmp[0]; + dr[ds[j+11]] = (int64_t)tmp[1]; + prodvh = vmovl_u16(vget_high_u16(prodv)); + tmp[0] = (uint64_t)dr[ds[j+12]]; + tmp[1] = (uint64_t)dr[ds[j+13]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_low_u32(prodvh)); + vst1q_u64(tmp, resv); + dr[ds[j+12]] = (int64_t)tmp[0]; + dr[ds[j+13]] = (int64_t)tmp[1]; + tmp[0] = (uint64_t)dr[ds[j+14]]; + tmp[1] = (uint64_t)dr[ds[j+15]]; + drv = vld1q_u64(tmp); + resv = vaddw_u32(drv, vget_high_u32(prodvh)); + vst1q_u64(tmp, resv); + dr[ds[j+14]] = (int64_t)tmp[0]; + dr[ds[j+15]] = (int64_t)tmp[1]; + } +} +#endif + +#ifdef HAVE_AVX2_KERNELS +static TARGET_AVX2 void add_mul_row_ff_8_avx2( + int64_t *dr, + const cf8_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + __m256i mask1 = _mm256_set1_epi64x(0x000000000000FFFF); + __m256i mask2 = _mm256_set1_epi64x(0x00000000FFFF0000); + __m256i mask3 = _mm256_set1_epi64x(0x0000FFFF00000000); + __m256i mask4 = _mm256_set1_epi64x(0xFFFF000000000000); + __m256i mask8 = _mm256_set1_epi16(0x00FF); + int64_t res[4] __attribute__((aligned(32))); + __m256i redv, mulv, prod, drv, resv; + + const uint16_t mul8 = (uint16_t)mul; + mulv = _mm256_set1_epi16(mul8); + const len_t os = len % 32; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 32) { + redv = _mm256_loadu_si256((__m256i*)(cfs+j)); + /* first 16 mult-adds -- lower */ + prod = _mm256_mullo_epi16(mulv,_mm256_and_si256(redv, mask8)); + drv = _mm256_setr_epi64x( + dr[ds[j]], + dr[ds[j+8]], + dr[ds[j+16]], + dr[ds[j+24]]); + resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask1)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j]] = res[0]; + dr[ds[j+8]] = res[1]; + dr[ds[j+16]] = res[2]; + dr[ds[j+24]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+2]], + dr[ds[j+10]], + dr[ds[j+18]], + dr[ds[j+26]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask2), 16)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+2]] = res[0]; + dr[ds[j+10]] = res[1]; + dr[ds[j+18]] = res[2]; + dr[ds[j+26]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+4]], + dr[ds[j+12]], + dr[ds[j+20]], + dr[ds[j+28]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask3), 32)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+4]] = res[0]; + dr[ds[j+12]] = res[1]; + dr[ds[j+20]] = res[2]; + dr[ds[j+28]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+6]], + dr[ds[j+14]], + dr[ds[j+22]], + dr[ds[j+30]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask4), 48)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+6]] = res[0]; + dr[ds[j+14]] = res[1]; + dr[ds[j+22]] = res[2]; + dr[ds[j+30]] = res[3]; + /* second 16 mult-adds -- higher */ + prod = _mm256_mullo_epi16(mulv, _mm256_srli_epi16(redv, 8)); + drv = _mm256_setr_epi64x( + dr[ds[j+1]], + dr[ds[j+9]], + dr[ds[j+17]], + dr[ds[j+25]]); + resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask1)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+1]] = res[0]; + dr[ds[j+9]] = res[1]; + dr[ds[j+17]] = res[2]; + dr[ds[j+25]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+3]], + dr[ds[j+11]], + dr[ds[j+19]], + dr[ds[j+27]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask2), 16)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+3]] = res[0]; + dr[ds[j+11]] = res[1]; + dr[ds[j+19]] = res[2]; + dr[ds[j+27]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+5]], + dr[ds[j+13]], + dr[ds[j+21]], + dr[ds[j+29]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask3), 32)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+5]] = res[0]; + dr[ds[j+13]] = res[1]; + dr[ds[j+21]] = res[2]; + dr[ds[j+29]] = res[3]; + drv = _mm256_setr_epi64x( + dr[ds[j+7]], + dr[ds[j+15]], + dr[ds[j+23]], + dr[ds[j+31]]); + resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask4), 48)); + _mm256_store_si256((__m256i*)(res),resv); + dr[ds[j+7]] = res[0]; + dr[ds[j+15]] = res[1]; + dr[ds[j+23]] = res[2]; + dr[ds[j+31]] = res[3]; + } +} +#endif + +#ifdef HAVE_AVX512_KERNELS +static TARGET_AVX512 void add_mul_row_ff_8_avx512( + int64_t *dr, + const cf8_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ + len_t j; + __m512i mask1 = _mm512_set1_epi64(0x000000000000FFFF); + __m512i mask2 = _mm512_set1_epi64(0x00000000FFFF0000); + __m512i mask3 = _mm512_set1_epi64(0x0000FFFF00000000); + __m512i mask4 = _mm512_set1_epi64(0xFFFF000000000000); + __m512i mask8 = _mm512_set1_epi16(0x00FF); + int64_t res[8] __attribute__((aligned(64))); + __m512i redv, mulv, prod, drv, resv; + + const uint16_t mul8 = (uint16_t)mul; + mulv = _mm512_set1_epi16(mul8); + const len_t os = len % 64; + for (j = 0; j < os; ++j) { + dr[ds[j]] += mul * cfs[j]; + } + for (; j < len; j += 64) { + redv = _mm512_loadu_si512((__m512i*)(cfs+j)); + /* first 16 mult-adds -- lower */ + prod = _mm512_mullo_epi16(mulv,_mm512_and_si512(redv, mask8)); + drv = _mm512_setr_epi64( + dr[ds[j]], + dr[ds[j+8]], + dr[ds[j+16]], + dr[ds[j+24]], + dr[ds[j+32]], + dr[ds[j+40]], + dr[ds[j+48]], + dr[ds[j+56]]); + resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask1)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j]] = res[0]; + dr[ds[j+8]] = res[1]; + dr[ds[j+16]] = res[2]; + dr[ds[j+24]] = res[3]; + dr[ds[j+32]] = res[4]; + dr[ds[j+40]] = res[5]; + dr[ds[j+48]] = res[6]; + dr[ds[j+56]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+2]], + dr[ds[j+10]], + dr[ds[j+18]], + dr[ds[j+26]], + dr[ds[j+34]], + dr[ds[j+42]], + dr[ds[j+50]], + dr[ds[j+58]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask2), 16)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+2]] = res[0]; + dr[ds[j+10]] = res[1]; + dr[ds[j+18]] = res[2]; + dr[ds[j+26]] = res[3]; + dr[ds[j+34]] = res[4]; + dr[ds[j+42]] = res[5]; + dr[ds[j+50]] = res[6]; + dr[ds[j+58]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+4]], + dr[ds[j+12]], + dr[ds[j+20]], + dr[ds[j+28]], + dr[ds[j+36]], + dr[ds[j+44]], + dr[ds[j+52]], + dr[ds[j+60]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask3), 32)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+4]] = res[0]; + dr[ds[j+12]] = res[1]; + dr[ds[j+20]] = res[2]; + dr[ds[j+28]] = res[3]; + dr[ds[j+36]] = res[4]; + dr[ds[j+44]] = res[5]; + dr[ds[j+52]] = res[6]; + dr[ds[j+60]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+6]], + dr[ds[j+14]], + dr[ds[j+22]], + dr[ds[j+30]], + dr[ds[j+38]], + dr[ds[j+46]], + dr[ds[j+54]], + dr[ds[j+62]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask4), 48)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+6]] = res[0]; + dr[ds[j+14]] = res[1]; + dr[ds[j+22]] = res[2]; + dr[ds[j+30]] = res[3]; + dr[ds[j+38]] = res[4]; + dr[ds[j+46]] = res[5]; + dr[ds[j+54]] = res[6]; + dr[ds[j+62]] = res[7]; + /* second 16 mult-adds -- higher */ + prod = _mm512_mullo_epi16(mulv, _mm512_srli_epi16(redv, 8)); + drv = _mm512_setr_epi64( + dr[ds[j+1]], + dr[ds[j+9]], + dr[ds[j+17]], + dr[ds[j+25]], + dr[ds[j+33]], + dr[ds[j+41]], + dr[ds[j+49]], + dr[ds[j+57]]); + resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask1)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+1]] = res[0]; + dr[ds[j+9]] = res[1]; + dr[ds[j+17]] = res[2]; + dr[ds[j+25]] = res[3]; + dr[ds[j+33]] = res[4]; + dr[ds[j+41]] = res[5]; + dr[ds[j+49]] = res[6]; + dr[ds[j+57]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+3]], + dr[ds[j+11]], + dr[ds[j+19]], + dr[ds[j+27]], + dr[ds[j+35]], + dr[ds[j+43]], + dr[ds[j+51]], + dr[ds[j+59]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask2), 16)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+3]] = res[0]; + dr[ds[j+11]] = res[1]; + dr[ds[j+19]] = res[2]; + dr[ds[j+27]] = res[3]; + dr[ds[j+35]] = res[4]; + dr[ds[j+43]] = res[5]; + dr[ds[j+51]] = res[6]; + dr[ds[j+59]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+5]], + dr[ds[j+13]], + dr[ds[j+21]], + dr[ds[j+29]], + dr[ds[j+37]], + dr[ds[j+45]], + dr[ds[j+53]], + dr[ds[j+61]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask3), 32)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+5]] = res[0]; + dr[ds[j+13]] = res[1]; + dr[ds[j+21]] = res[2]; + dr[ds[j+29]] = res[3]; + dr[ds[j+37]] = res[4]; + dr[ds[j+45]] = res[5]; + dr[ds[j+53]] = res[6]; + dr[ds[j+61]] = res[7]; + drv = _mm512_setr_epi64( + dr[ds[j+7]], + dr[ds[j+15]], + dr[ds[j+23]], + dr[ds[j+31]], + dr[ds[j+39]], + dr[ds[j+47]], + dr[ds[j+55]], + dr[ds[j+63]]); + resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask4), 48)); + _mm512_store_si512((__m512i*)(res),resv); + dr[ds[j+7]] = res[0]; + dr[ds[j+15]] = res[1]; + dr[ds[j+23]] = res[2]; + dr[ds[j+31]] = res[3]; + dr[ds[j+39]] = res[4]; + dr[ds[j+47]] = res[5]; + dr[ds[j+55]] = res[6]; + dr[ds[j+63]] = res[7]; + } +} +#endif + +/* dr[ds[j]] += mul * cfs[j] for 0 <= j < len */ +static inline void add_mul_row_ff_8( + int64_t *dr, + const cf8_t *cfs, + const hm_t *ds, + const len_t len, + const uint32_t mul + ) +{ +#ifdef HAVE_AVX512_KERNELS + if (cpu_has_avx512()) { + add_mul_row_ff_8_avx512(dr, cfs, ds, len, mul); + return; + } +#endif +#ifdef HAVE_AVX2_KERNELS + if (cpu_has_avx2()) { + add_mul_row_ff_8_avx2(dr, cfs, ds, len, mul); + return; + } +#endif +#ifdef __aarch64__ + add_mul_row_ff_8_neon(dr, cfs, ds, len, mul); +#else + add_mul_row_ff_8_scalar(dr, cfs, ds, len, mul); +#endif +} + static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_8( int64_t *dr, mat_t *mat, @@ -122,29 +566,6 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_8( } else { rba = NULL; } -#if defined HAVE_AVX512_F - __m512i mask1 = _mm512_set1_epi64(0x000000000000FFFF); - __m512i mask2 = _mm512_set1_epi64(0x00000000FFFF0000); - __m512i mask3 = _mm512_set1_epi64(0x0000FFFF00000000); - __m512i mask4 = _mm512_set1_epi64(0xFFFF000000000000); - __m512i mask8 = _mm512_set1_epi16(0x00FF); - int64_t res[8] __attribute__((aligned(64))); - __m512i redv, mulv, prod, drv, resv; -#elif defined HAVE_AVX2 - __m256i mask1 = _mm256_set1_epi64x(0x000000000000FFFF); - __m256i mask2 = _mm256_set1_epi64x(0x00000000FFFF0000); - __m256i mask3 = _mm256_set1_epi64x(0x0000FFFF00000000); - __m256i mask4 = _mm256_set1_epi64x(0xFFFF000000000000); - __m256i mask8 = _mm256_set1_epi16(0x00FF); - int64_t res[4] __attribute__((aligned(32))); - __m256i redv, mulv, prod, drv, resv; -#elif defined __aarch64__ - uint64_t tmp[2] __attribute__((aligned(32))); - uint16x8_t prodv; - uint32x4_t prodvl, prodvh; - uint8x16_t redv; - uint64x2_t drv, resv; -#endif k = 0; for (i = dpiv; i < ncols; ++i) { @@ -171,368 +592,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_8( } } cfs = mcf[dts[COEFFS]]; -#if defined HAVE_AVX512_F - const uint16_t mul8 = (uint16_t)(fc - dr[i]); - mulv = _mm512_set1_epi16(mul8); - const len_t len = dts[LENGTH]; - const len_t os = len % 64; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 64) { - redv = _mm512_loadu_si512((__m512i*)(cfs+j)); - /* first 16 mult-adds -- lower */ - prod = _mm512_mullo_epi16(mulv,_mm512_and_si512(redv, mask8)); - drv = _mm512_setr_epi64( - dr[ds[j]], - dr[ds[j+8]], - dr[ds[j+16]], - dr[ds[j+24]], - dr[ds[j+32]], - dr[ds[j+40]], - dr[ds[j+48]], - dr[ds[j+56]]); - resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask1)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j]] = res[0]; - dr[ds[j+8]] = res[1]; - dr[ds[j+16]] = res[2]; - dr[ds[j+24]] = res[3]; - dr[ds[j+32]] = res[4]; - dr[ds[j+40]] = res[5]; - dr[ds[j+48]] = res[6]; - dr[ds[j+56]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+2]], - dr[ds[j+10]], - dr[ds[j+18]], - dr[ds[j+26]], - dr[ds[j+34]], - dr[ds[j+42]], - dr[ds[j+50]], - dr[ds[j+58]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask2), 16)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+2]] = res[0]; - dr[ds[j+10]] = res[1]; - dr[ds[j+18]] = res[2]; - dr[ds[j+26]] = res[3]; - dr[ds[j+34]] = res[4]; - dr[ds[j+42]] = res[5]; - dr[ds[j+50]] = res[6]; - dr[ds[j+58]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+4]], - dr[ds[j+12]], - dr[ds[j+20]], - dr[ds[j+28]], - dr[ds[j+36]], - dr[ds[j+44]], - dr[ds[j+52]], - dr[ds[j+60]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask3), 32)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+4]] = res[0]; - dr[ds[j+12]] = res[1]; - dr[ds[j+20]] = res[2]; - dr[ds[j+28]] = res[3]; - dr[ds[j+36]] = res[4]; - dr[ds[j+44]] = res[5]; - dr[ds[j+52]] = res[6]; - dr[ds[j+60]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+6]], - dr[ds[j+14]], - dr[ds[j+22]], - dr[ds[j+30]], - dr[ds[j+38]], - dr[ds[j+46]], - dr[ds[j+54]], - dr[ds[j+62]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask4), 48)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+6]] = res[0]; - dr[ds[j+14]] = res[1]; - dr[ds[j+22]] = res[2]; - dr[ds[j+30]] = res[3]; - dr[ds[j+38]] = res[4]; - dr[ds[j+46]] = res[5]; - dr[ds[j+54]] = res[6]; - dr[ds[j+62]] = res[7]; - /* second 16 mult-adds -- higher */ - prod = _mm512_mullo_epi16(mulv, _mm512_srli_epi16(redv, 8)); - drv = _mm512_setr_epi64( - dr[ds[j+1]], - dr[ds[j+9]], - dr[ds[j+17]], - dr[ds[j+25]], - dr[ds[j+33]], - dr[ds[j+41]], - dr[ds[j+49]], - dr[ds[j+57]]); - resv = _mm512_add_epi64(drv, _mm512_and_si512(prod, mask1)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+1]] = res[0]; - dr[ds[j+9]] = res[1]; - dr[ds[j+17]] = res[2]; - dr[ds[j+25]] = res[3]; - dr[ds[j+33]] = res[4]; - dr[ds[j+41]] = res[5]; - dr[ds[j+49]] = res[6]; - dr[ds[j+57]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+3]], - dr[ds[j+11]], - dr[ds[j+19]], - dr[ds[j+27]], - dr[ds[j+35]], - dr[ds[j+43]], - dr[ds[j+51]], - dr[ds[j+59]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask2), 16)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+3]] = res[0]; - dr[ds[j+11]] = res[1]; - dr[ds[j+19]] = res[2]; - dr[ds[j+27]] = res[3]; - dr[ds[j+35]] = res[4]; - dr[ds[j+43]] = res[5]; - dr[ds[j+51]] = res[6]; - dr[ds[j+59]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+5]], - dr[ds[j+13]], - dr[ds[j+21]], - dr[ds[j+29]], - dr[ds[j+37]], - dr[ds[j+45]], - dr[ds[j+53]], - dr[ds[j+61]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask3), 32)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+5]] = res[0]; - dr[ds[j+13]] = res[1]; - dr[ds[j+21]] = res[2]; - dr[ds[j+29]] = res[3]; - dr[ds[j+37]] = res[4]; - dr[ds[j+45]] = res[5]; - dr[ds[j+53]] = res[6]; - dr[ds[j+61]] = res[7]; - drv = _mm512_setr_epi64( - dr[ds[j+7]], - dr[ds[j+15]], - dr[ds[j+23]], - dr[ds[j+31]], - dr[ds[j+39]], - dr[ds[j+47]], - dr[ds[j+55]], - dr[ds[j+63]]); - resv = _mm512_add_epi64(drv, _mm512_srli_epi64(_mm512_and_si512(prod, mask4), 48)); - _mm512_store_si512((__m512i*)(res),resv); - dr[ds[j+7]] = res[0]; - dr[ds[j+15]] = res[1]; - dr[ds[j+23]] = res[2]; - dr[ds[j+31]] = res[3]; - dr[ds[j+39]] = res[4]; - dr[ds[j+47]] = res[5]; - dr[ds[j+55]] = res[6]; - dr[ds[j+63]] = res[7]; - } -#elif defined HAVE_AVX2 - const uint16_t mul8 = (uint16_t)(fc - dr[i]); - mulv = _mm256_set1_epi16(mul8); - const len_t len = dts[LENGTH]; - const len_t os = len % 32; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 32) { - redv = _mm256_loadu_si256((__m256i*)(cfs+j)); - /* first 16 mult-adds -- lower */ - prod = _mm256_mullo_epi16(mulv,_mm256_and_si256(redv, mask8)); - drv = _mm256_setr_epi64x( - dr[ds[j]], - dr[ds[j+8]], - dr[ds[j+16]], - dr[ds[j+24]]); - resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask1)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j]] = res[0]; - dr[ds[j+8]] = res[1]; - dr[ds[j+16]] = res[2]; - dr[ds[j+24]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+2]], - dr[ds[j+10]], - dr[ds[j+18]], - dr[ds[j+26]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask2), 16)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+2]] = res[0]; - dr[ds[j+10]] = res[1]; - dr[ds[j+18]] = res[2]; - dr[ds[j+26]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+4]], - dr[ds[j+12]], - dr[ds[j+20]], - dr[ds[j+28]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask3), 32)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+4]] = res[0]; - dr[ds[j+12]] = res[1]; - dr[ds[j+20]] = res[2]; - dr[ds[j+28]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+6]], - dr[ds[j+14]], - dr[ds[j+22]], - dr[ds[j+30]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask4), 48)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+6]] = res[0]; - dr[ds[j+14]] = res[1]; - dr[ds[j+22]] = res[2]; - dr[ds[j+30]] = res[3]; - /* second 16 mult-adds -- higher */ - prod = _mm256_mullo_epi16(mulv, _mm256_srli_epi16(redv, 8)); - drv = _mm256_setr_epi64x( - dr[ds[j+1]], - dr[ds[j+9]], - dr[ds[j+17]], - dr[ds[j+25]]); - resv = _mm256_add_epi64(drv, _mm256_and_si256(prod, mask1)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+1]] = res[0]; - dr[ds[j+9]] = res[1]; - dr[ds[j+17]] = res[2]; - dr[ds[j+25]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+3]], - dr[ds[j+11]], - dr[ds[j+19]], - dr[ds[j+27]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask2), 16)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+3]] = res[0]; - dr[ds[j+11]] = res[1]; - dr[ds[j+19]] = res[2]; - dr[ds[j+27]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+5]], - dr[ds[j+13]], - dr[ds[j+21]], - dr[ds[j+29]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask3), 32)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+5]] = res[0]; - dr[ds[j+13]] = res[1]; - dr[ds[j+21]] = res[2]; - dr[ds[j+29]] = res[3]; - drv = _mm256_setr_epi64x( - dr[ds[j+7]], - dr[ds[j+15]], - dr[ds[j+23]], - dr[ds[j+31]]); - resv = _mm256_add_epi64(drv, _mm256_srli_epi64(_mm256_and_si256(prod, mask4), 48)); - _mm256_store_si256((__m256i*)(res),resv); - dr[ds[j+7]] = res[0]; - dr[ds[j+15]] = res[1]; - dr[ds[j+23]] = res[2]; - dr[ds[j+31]] = res[3]; - } -#elif defined __aarch64__ - const len_t len = dts[LENGTH]; - const len_t os = len % 16; - const hm_t * const ds = dts + OFFSET; - const cf8_t mul8 = (cf8_t)(mod - dr[i]); - const uint8x8_t mulv = vmov_n_u8(mul8); - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += 16) { - redv = vld1q_u8(cfs+j); - - prodv = vmull_u8(vget_low_u8(redv), mulv); - prodvl = vmovl_u16(vget_low_u16(prodv)); - tmp[0] = (uint64_t)dr[ds[j]]; - tmp[1] = (uint64_t)dr[ds[j+1]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_low_u32(prodvl)); - vst1q_u64(tmp, resv); - dr[ds[j]] = (int64_t)tmp[0]; - dr[ds[j+1]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+2]]; - tmp[1] = (uint64_t)dr[ds[j+3]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_high_u32(prodvl)); - vst1q_u64(tmp, resv); - dr[ds[j+2]] = (int64_t)tmp[0]; - dr[ds[j+3]] = (int64_t)tmp[1]; - prodvh = vmovl_u16(vget_high_u16(prodv)); - tmp[0] = (uint64_t)dr[ds[j+4]]; - tmp[1] = (uint64_t)dr[ds[j+5]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_low_u32(prodvh)); - vst1q_u64(tmp, resv); - dr[ds[j+4]] = (int64_t)tmp[0]; - dr[ds[j+5]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+6]]; - tmp[1] = (uint64_t)dr[ds[j+7]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_high_u32(prodvh)); - vst1q_u64(tmp, resv); - dr[ds[j+6]] = (int64_t)tmp[0]; - dr[ds[j+7]] = (int64_t)tmp[1]; - - prodv = vmull_u8(vget_high_u8(redv), mulv); - prodvl = vmovl_u16(vget_low_u16(prodv)); - tmp[0] = (uint64_t)dr[ds[j+8]]; - tmp[1] = (uint64_t)dr[ds[j+9]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_low_u32(prodvl)); - vst1q_u64(tmp, resv); - dr[ds[j+8]] = (int64_t)tmp[0]; - dr[ds[j+9]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+10]]; - tmp[1] = (uint64_t)dr[ds[j+11]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_high_u32(prodvl)); - vst1q_u64(tmp, resv); - dr[ds[j+10]] = (int64_t)tmp[0]; - dr[ds[j+11]] = (int64_t)tmp[1]; - prodvh = vmovl_u16(vget_high_u16(prodv)); - tmp[0] = (uint64_t)dr[ds[j+12]]; - tmp[1] = (uint64_t)dr[ds[j+13]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_low_u32(prodvh)); - vst1q_u64(tmp, resv); - dr[ds[j+12]] = (int64_t)tmp[0]; - dr[ds[j+13]] = (int64_t)tmp[1]; - tmp[0] = (uint64_t)dr[ds[j+14]]; - tmp[1] = (uint64_t)dr[ds[j+15]]; - drv = vld1q_u64(tmp); - resv = vaddw_u32(drv, vget_high_u32(prodvh)); - vst1q_u64(tmp, resv); - dr[ds[j+14]] = (int64_t)tmp[0]; - dr[ds[j+15]] = (int64_t)tmp[1]; - } -#else - const len_t os = dts[PRELOOP]; - const len_t len = dts[LENGTH]; - const hm_t * const ds = dts + OFFSET; - for (j = 0; j < os; ++j) { - dr[ds[j]] += mul * cfs[j]; - } - for (; j < len; j += UNROLL) { - dr[ds[j]] += mul * cfs[j]; - dr[ds[j+1]] += mul * cfs[j+1]; - dr[ds[j+2]] += mul * cfs[j+2]; - dr[ds[j+3]] += mul * cfs[j+3]; - } -#endif + add_mul_row_ff_8(dr, cfs, dts + OFFSET, dts[LENGTH], mul); dr[i] = 0; } if (k == 0) { diff --git a/src/neogb/symbol.c b/src/neogb/symbol.c index 6f8ab251..b21c147d 100644 --- a/src/neogb/symbol.c +++ b/src/neogb/symbol.c @@ -23,10 +23,6 @@ #include "../msolve/streams.h" -#ifdef HAVE_AVX2 -#include -#endif - /* select_all_pairs() is unused at the moment */ #if 0 static void select_all_spairs( diff --git a/src/usolve/Makefile.am b/src/usolve/Makefile.am index 43216bdc..70f63283 100644 --- a/src/usolve/Makefile.am +++ b/src/usolve/Makefile.am @@ -2,7 +2,7 @@ noinst_LTLIBRARIES = libusolve.la libusolve_la_SOURCES = usolve.c libusolve_ladir = $(includedir)/msolve/usolve libusolve_la_HEADERS = libusolve.h -libusolve_la_CFLAGS = $(SIMD_FLAGS) $(CPUEXT_FLAGS) $(OPENMP_CFLAGS) -Wall -Wextra +libusolve_la_CFLAGS = $(OPENMP_CFLAGS) -Wall -Wextra EXTRA_DIST = libusolve.h \ bisection.h \