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 \