avx512_vectorization.hpp

Go to the documentation of this file.

//=======================================================================
// Copyright (c) 2014-2023 Baptiste Wicht
// Distributed under the terms of the MIT License.
// (See accompanying file LICENSE or copy at
//  http://opensource.org/licenses/MIT)
//=======================================================================

//TODO Implementation of AVX-512 complex multiplication and division

#pragma once

#ifdef __AVX512F__

#include <immintrin.h>

#include "etl/inline.hpp"

#ifdef VECT_DEBUG
#include <iostream>
#endif

#define ETL_INLINE_VEC_VOID ETL_STATIC_INLINE(void)
#define ETL_INLINE_VEC_512 ETL_STATIC_INLINE(__m512)
#define ETL_INLINE_VEC_512D ETL_STATIC_INLINE(__m512d)
#define ETL_OUT_VEC_2512ETL_OUT_INLINE(__m512)
#define ETL_OUT_VEC_512D ETL_OUT_INLINE(__m512d)

namespace etl {

using avx_512_simd_float = simd_pack<vector_mode_t::AVX512, float, __m512>;

using avx_512_simd_double = simd_pack<vector_mode_t::AVX512, double, __m512d>;

template <typename T>
using avx_512_simd_complex_float = simd_pack<vector_mode_t::AVX512, T, __m512>;

template <typename T>
using avx_512_simd_complex_double = simd_pack<vector_mode_t::AVX512, T, __m512d>;

using avx_512_simd_byte = simd_pack<vector_mode_t::AVX512, int8_t, __m512i>;

using avx_512_simd_short = simd_pack<vector_mode_t::AVX512, int16_t, __m512i>;

using avx_512_simd_int = simd_pack<vector_mode_t::AVX512, int32_t, __m512i>;

using avx_512_simd_long = simd_pack<vector_mode_t::AVX512, int64_t, __m512i>;

template <typename T>
struct avx512_intrinsic_traits {
    static constexpr bool vectorizable = false;      
    static constexpr size_t size       = 1;          
    static constexpr size_t alignment  = alignof(T); 

    using intrinsic_type = T; 
};

template <>
struct avx512_intrinsic_traits<float> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 16;   
    static constexpr size_t alignment  = 64;   

    using intrinsic_type = avx_512_simd_float; 
};

template <>
struct avx512_intrinsic_traits<double> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 8;    
    static constexpr size_t alignment  = 64;   

    using intrinsic_type = avx_512_simd_double; 
};

template <>
struct avx512_intrinsic_traits<std::complex<float>> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 8;    
    static constexpr size_t alignment  = 64;   

    using intrinsic_type = avx_512_simd_complex_float<std::complex<float>>; 
};

template <>
struct avx512_intrinsic_traits<std::complex<double>> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 4;    
    static constexpr size_t alignment  = 64;   

    using intrinsic_type = avx_512_simd_complex_double<std::complex<double>>; 
};

template <>
struct avx512_intrinsic_traits<etl::complex<float>> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 8;    
    static constexpr size_t alignment  = 64;   

    using intrinsic_type = avx_512_simd_complex_float<etl::complex<float>>; 
};

template <>
struct avx512_intrinsic_traits<etl::complex<double>> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 4;    
    static constexpr size_t alignment  = 64;   

    using intrinsic_type = avx_512_simd_complex_double<etl::complex<double>>; 
};

template <>
struct avx512_intrinsic_traits<int8_t> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 64;           
    static constexpr size_t alignment  = 64;           

    using intrinsic_type = avx_512_simd_byte; 
};

template <>
struct avx512_intrinsic_traits<int16_t> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 32;           
    static constexpr size_t alignment  = 64;           

    using intrinsic_type = avx_512_simd_short; 
};

template <>
struct avx512_intrinsic_traits<int32_t> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 16;            
    static constexpr size_t alignment  = 64;           

    using intrinsic_type = avx_512_simd_int; 
};

template <>
struct avx512_intrinsic_traits<int64_t> {
    static constexpr bool vectorizable = true; 
    static constexpr size_t size       = 8;            
    static constexpr size_t alignment  = 64;           

    using intrinsic_type = avx_512_simd_long; 
};

struct avx512_vec {
    template <typename T>
    using traits = avx512_intrinsic_traits<T>;

    template <typename T>
    using vec_type = typename traits<T>::intrinsic_type;

#ifdef VEC_DEBUG

    template <typename T>
    static std::string debug_d(T value) {
        union test {
            __m512d vec;
            double array[8];
            test(__m512d vec) : vec(vec) {}
        };

        test u_value = value;
        std::cout << "[" << u_value.array[0] << "," << u_value.array[1] << "," << u_value.array[2] << "," << u_value.array[3] << "," << u_value.array[4] << ","
                  << u_value.array[5] << "," << u_value.array[6] << "," << u_value.array[7] << "]" << std::endl;
    }

    template <typename T>
    static std::string debug_s(T value) {
        union test {
            __m512 vec;
            float array[16];
            test(__m512 vec) : vec(vec) {}
        };

        test u_value = value;
        std::cout << "[" << u_value.array[0] << "," << u_value.array[1] << "," << u_value.array[2] << "," << u_value.array[3] << "," << u_value.array[4] << ","
                  << u_value.array[5] << "," << u_value.array[6] << "," << u_value.array[7] << "," << u_value.array[8] << "," << u_value.array[9] << ","
                  << u_value.array[10] << "," << u_value.array[11] << "," << u_value.array[12] << "," << u_value.array[13] << "," << u_value.array[14] << ","
                  << u_value.array[15] << "]" << std::endl;
    }

#else

    template <typename T>
    static std::string debug_d(T) {
        return "";
    }

    template <typename T>
    static std::string debug_s(T) {
        return "";
    }

#endif

    ETL_INLINE_VEC_VOID storeu(float* memory, avx_512_simd_float value) {
        _mm512_storeu_ps(memory, value.value);
    }

    ETL_INLINE_VEC_VOID storeu(double* memory, avx_512_simd_double value) {
        _mm512_storeu_pd(memory, value.value);
    }

    ETL_INLINE_VEC_VOID storeu(std::complex<float>* memory, avx_512_simd_complex_float<std::complex<float>> value) {
        _mm512_storeu_ps(reinterpret_cast<float*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID storeu(std::complex<double>* memory, avx_512_simd_complex_double<std::complex<double>> value) {
        _mm512_storeu_pd(reinterpret_cast<double*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID storeu(etl::complex<float>* memory, avx_512_simd_complex_float<etl::complex<float>> value) {
        _mm512_storeu_ps(reinterpret_cast<float*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID storeu(etl::complex<double>* memory, avx_512_simd_complex_double<etl::complex<double>> value) {
        _mm512_storeu_pd(reinterpret_cast<double*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) store(int8_t* memory, avx_512_simd_byte value) {
        _mm512_store_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) store(int16_t* memory, avx_512_simd_short value) {
        _mm512_store_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) store(int32_t* memory, avx_512_simd_int value) {
        _mm512_store_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) store(int64_t* memory, avx_512_simd_long value) {
        _mm512_store_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID store(float* memory, avx_512_simd_float value) {
        _mm512_store_ps(memory, value.value);
    }

    ETL_INLINE_VEC_VOID store(double* memory, avx_512_simd_double value) {
        _mm512_store_pd(memory, value.value);
    }

    ETL_INLINE_VEC_VOID store(std::complex<float>* memory, avx_512_simd_complex_float<std::complex<float>> value) {
        _mm512_store_ps(reinterpret_cast<float*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID store(std::complex<double>* memory, avx_512_simd_complex_double<std::complex<double>> value) {
        _mm512_store_pd(reinterpret_cast<double*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID store(etl::complex<float>* memory, avx_512_simd_complex_float<etl::complex<float>> value) {
        _mm512_store_ps(reinterpret_cast<float*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID store(etl::complex<double>* memory, avx_512_simd_complex_double<etl::complex<double>> value) {
        _mm512_store_pd(reinterpret_cast<double*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) stream(int8_t* memory, avx_512_simd_byte value) {
        _mm512_stream_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) stream(int16_t* memory, avx_512_simd_short value) {
        _mm512_stream_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) stream(int32_t* memory, avx_512_simd_int value) {
        _mm512_stream_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_STATIC_INLINE(void) stream(int64_t* memory, avx_512_simd_long value) {
        _mm512_stream_si512(reinterpret_cast<__m512i*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID stream(float* memory, avx_512_simd_float value) {
        _mm512_stream_ps(memory, value.value);
    }

    ETL_INLINE_VEC_VOID stream(double* memory, avx_512_simd_double value) {
        _mm512_stream_pd(memory, value.value);
    }

    ETL_INLINE_VEC_VOID stream(std::complex<float>* memory, avx_512_simd_complex_float<std::complex<float>> value) {
        _mm512_stream_ps(reinterpret_cast<float*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID stream(std::complex<double>* memory, avx_512_simd_complex_double<std::complex<double>> value) {
        _mm512_stream_pd(reinterpret_cast<double*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID stream(etl::complex<float>* memory, avx_512_simd_complex_float<etl::complex<float>> value) {
        _mm512_stream_ps(reinterpret_cast<float*>(memory), value.value);
    }

    ETL_INLINE_VEC_VOID stream(etl::complex<double>* memory, avx_512_simd_complex_double<etl::complex<double>> value) {
        _mm512_stream_pd(reinterpret_cast<double*>(memory), value.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_byte) load(const int8_t* memory) {
        return _mm512_load_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_short) load(const int16_t* memory) {
        return _mm512_load_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_int) load(const int32_t* memory) {
        return _mm512_load_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_long) load(const int64_t* memory) {
        return _mm512_load_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_float) load(const float* memory) {
        return _mm512_load_ps(memory);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) load(const double* memory) {
        return _mm512_load_pd(memory);
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_float<std::complex<float>>) load(const std::complex<float>* memory) {
        return _mm512_load_ps(reinterpret_cast<const float*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_double<std::complex<double>>) load(const std::complex<double>* memory) {
        return _mm512_load_pd(reinterpret_cast<const double*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_float<etl::complex<float>>) load(const etl::complex<float>* memory) {
        return _mm512_load_ps(reinterpret_cast<const float*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_double<etl::complex<double>>) load(const etl::complex<double>* memory) {
        return _mm512_load_pd(reinterpret_cast<const double*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_byte) loadu(const int8_t* memory) {
        return _mm512_loadu_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_short) loadu(const int16_t* memory) {
        return _mm512_loadu_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_int) loadu(const int32_t* memory) {
        return _mm512_loadu_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_long) loadu(const int64_t* memory) {
        return _mm512_loadu_si512(reinterpret_cast<const __m512i*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_float) loadu(const float* memory) {
        return _mm512_loadu_ps(memory);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) loadu(const double* memory) {
        return _mm512_loadu_pd(memory);
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_float<std::complex<float>>) loadu(const std::complex<float>* memory) {
        return _mm512_loadu_ps(reinterpret_cast<const float*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_double<std::complex<double>>) loadu(const std::complex<double>* memory) {
        return _mm512_loadu_pd(reinterpret_cast<const double*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_float<etl::complex<float>>) loadu(const etl::complex<float>* memory) {
        return _mm512_loadu_ps(reinterpret_cast<const float*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_double<etl::complex<double>>) loadu(const etl::complex<double>* memory) {
        return _mm512_loadu_pd(reinterpret_cast<const double*>(memory));
    }

    ETL_STATIC_INLINE(avx_512_simd_byte) set(int8_t value) {
        return _mm512_set1_epi8(value);
    }

    ETL_STATIC_INLINE(avx_512_simd_short) set(int16_t value) {
        return _mm512_set1_epi16(value);
    }

    ETL_STATIC_INLINE(avx_512_simd_int) set(int32_t value) {
        return _mm512_set1_epi32(value);
    }

    ETL_STATIC_INLINE(avx_512_simd_long) set(int64_t value) {
        return _mm512_set1_epi64(value);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) set(double value) {
        return _mm512_set1_pd(value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) set(float value) {
        return _mm512_set1_ps(value);
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_float<std::complex<float>>) set(std::complex<float> value) {
        std::complex<float> tmp[]{value, value, value, value, value, value, value, value};
        return loadu(tmp);
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_double<std::complex<double>>) set(std::complex<double> value) {
        std::complex<double> tmp[]{value, value, value, value};
        return loadu(tmp);
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_float<etl::complex<float>>) set(etl::complex<float> value) {
        etl::complex<float> tmp[]{value, value, value, value, value, value, value, value};
        return loadu(tmp);
    }

    ETL_STATIC_INLINE(avx_512_simd_complex_double<etl::complex<double>>) set(etl::complex<double> value) {
        etl::complex<double> tmp[]{value, value, value, value};
        return loadu(tmp);
    }

    template <typename T>
    ETL_TMP_INLINE(typename avx512_intrinsic_traits<T>::intrinsic_type)
    zero();

    ETL_STATIC_INLINE(avx_512_simd_float) round_up(avx_512_simd_float x) {
        return _mm512_roundscale_round_ps(x.value, _MM_FROUND_TO_POS_INF, _MM_FROUND_NO_EXC);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) round_up(avx_512_simd_double x) {
        return _mm512_roundscale_round_pd(x.value, _MM_FROUND_TO_POS_INF,  _MM_FROUND_NO_EXC);
    }

    ETL_STATIC_INLINE(avx_512_simd_byte) add(avx_512_simd_byte lhs, avx_512_simd_byte rhs) {
        return _mm512_add_epi8(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_short) add(avx_512_simd_short lhs, avx_512_simd_short rhs) {
        return _mm512_add_epi16(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_int) add(avx_512_simd_int lhs, avx_512_simd_int rhs) {
        return _mm512_add_epi32(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_long) add(avx_512_simd_long lhs, avx_512_simd_long rhs) {
        return _mm512_add_epi64(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) add(avx_512_simd_float lhs, avx_512_simd_float rhs) {
        return _mm512_add_ps(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) add(avx_512_simd_double lhs, avx_512_simd_double rhs) {
        return _mm512_add_pd(lhs.value, rhs.value);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_float<T>)
    add(avx_512_simd_complex_float<T> lhs, avx_512_simd_complex_float<T> rhs) {
        return _mm512_add_ps(lhs.value, rhs.value);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_double<T>)
    add(avx_512_simd_complex_double<T> lhs, avx_512_simd_complex_double<T> rhs) {
        return _mm512_add_pd(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_byte) sub(avx_512_simd_byte lhs, avx_512_simd_byte rhs) {
        return _mm512_sub_epi8(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_short) sub(avx_512_simd_short lhs, avx_512_simd_short rhs) {
        return _mm512_sub_epi16(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_int) sub(avx_512_simd_int lhs, avx_512_simd_int rhs) {
        return _mm512_sub_epi32(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_long) sub(avx_512_simd_long lhs, avx_512_simd_long rhs) {
        return _mm512_sub_epi64(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) sub(avx_512_simd_float lhs, avx_512_simd_float rhs) {
        return _mm512_sub_ps(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) sub(avx_512_simd_double lhs, avx_512_simd_double rhs) {
        return _mm512_sub_pd(lhs.value, rhs.value);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_float<T>)
    sub(avx_512_simd_complex_float<T> lhs, avx_512_simd_complex_float<T> rhs) {
        return _mm512_sub_ps(lhs.value, rhs.value);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_double<T>)
    sub(avx_512_simd_complex_double<T> lhs, avx_512_simd_complex_double<T> rhs) {
        return _mm512_sub_pd(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) sqrt(avx_512_simd_float x) {
        return _mm512_sqrt_ps(x.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) sqrt(avx_512_simd_double x) {
        return _mm512_sqrt_pd(x.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) minus(avx_512_simd_float x) {
        return _mm512_xor_ps(x.value, _mm512_set1_ps(-0.f));
    }

    ETL_STATIC_INLINE(avx_512_simd_double) minus(avx_512_simd_double x) {
        return _mm512_xor_pd(x.value, _mm512_set1_pd(-0.));
    }

    ETL_STATIC_INLINE(avx_512_simd_byte) mul(avx_512_simd_byte lhs, avx_512_simd_byte rhs) {
        // Split in multiple vectors (odd and even)
        __m512i lhs_odd = _mm512_srli_epi16(lhs.value, 8);
        __m512i rhs_odd = _mm512_srli_epi16(rhs.value, 8);
        // Do the multiplication on each side
        __m512i mul_even = _mm512_mullo_epi16(lhs.value, rhs.value);
        __m512i mul_odd  = _mm512_mullo_epi16(lhs_odd, rhs_odd);
        // Combine again
        __m512i temp = _mm512_slli_epi16(mul_odd, 8);
        return _mm512_mask_mov_epi8(mul_even, 0xAAAAAAAAAAAAAAAA, temp);
    }

    ETL_STATIC_INLINE(avx_512_simd_short) mul(avx_512_simd_short lhs, avx_512_simd_short rhs) {
        return _mm512_mullo_epi16(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_int) mul(avx_512_simd_int lhs, avx_512_simd_int rhs) {
        return _mm512_mullo_epi32(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_long) mul(avx_512_simd_long lhs, avx_512_simd_long rhs) {
        return _mm512_mullo_epi64(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) mul(avx_512_simd_float lhs, avx_512_simd_float rhs) {
        return _mm512_mul_ps(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) mul(avx_512_simd_double lhs, avx_512_simd_double rhs) {
        return _mm512_mul_pd(lhs.value, rhs.value);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_float<T>)
    mul(avx_512_simd_complex_float<T> lhs, avx_512_simd_complex_float<T> rhs) {
        //lhs = [x1.real, x1.img, x2.real, x2.img, ...]
        //rhs = [y1.real, y1.img, y2.real, y2.img, ...]

        //ymm1 = [y1.real, y1.real, y2.real, y2.real, ...]
        __m512 ymm1 = _mm512_moveldup_ps(rhs.value);

        //ymm2 = [x1.img, x1.real, x2.img, x2.real]
        __m512 ymm2 = _mm512_permute_ps(lhs.value, 0b10110001);

        //ymm3 = [y1.imag, y1.imag, y2.imag, y2.imag]
        __m512 ymm3 = _mm512_movehdup_ps(rhs.value);

        //ymm4 = ymm2 * ymm3
        __m512 ymm4 = _mm512_mul_ps(ymm2, ymm3);

        //result = [(lhs * ymm1) -+ ymm4];

        return _mm512_fmaddsub_ps(lhs.value, ymm1, ymm4);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_double<T>)
    mul(avx_512_simd_complex_double<T> lhs, avx_512_simd_complex_double<T> rhs) {
        __m512d ymm1 = _mm512_shuffle_pd(rhs.value, rhs.value, 0x55);
        __m512d ymm2 = _mm512_shuffle_pd(lhs.value, lhs.value, 0xFF);
        __m512d ymm3 = _mm512_shuffle_pd(lhs.value, lhs.value, 0);
        __m512d ymm4 = _mm512_mul_pd(ymm2, ymm1);
        return _mm512_fmaddsub_pd(ymm3, rhs.value, ymm4);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) fmadd(avx_512_simd_float a, avx_512_simd_float b, avx_512_simd_float c) {
        return _mm512_fmadd_ps(a.value, b.value, c.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) fmadd(avx_512_simd_double a, avx_512_simd_double b, avx_512_simd_double c) {
        return _mm512_fmadd_pd(a.value, b.value, c.value);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_float<T>)
    fmadd(avx_512_simd_complex_float<T> a, avx_512_simd_complex_float<T> b, avx_512_simd_complex_float<T> c) {
        return add(mul(a, b), c);
    }

    template <typename T>
    ETL_STATIC_INLINE(avx_512_simd_complex_double<T>)
    fmadd(avx_512_simd_complex_double<T> a, avx_512_simd_complex_double<T> b, avx_512_simd_complex_double<T> c) {
        return add(mul(a, b), c);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) div(avx_512_simd_float lhs, avx_512_simd_float rhs) {
        return _mm512_div_ps(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_double) div(avx_512_simd_double lhs, avx_512_simd_double rhs) {
        return _mm512_div_pd(lhs.value, rhs.value);
    }

    //Min

    ETL_STATIC_INLINE(avx_512_simd_double) min(avx_512_simd_double lhs, avx_512_simd_double rhs) {
        return _mm512_min_pd(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) min(avx_512_simd_float lhs, avx_512_simd_float rhs) {
        return _mm512_min_ps(lhs.value, rhs.value);
    }

    //Max

    ETL_STATIC_INLINE(avx_512_simd_double) max(avx_512_simd_double lhs, avx_512_simd_double rhs) {
        return _mm512_max_pd(lhs.value, rhs.value);
    }

    ETL_STATIC_INLINE(avx_512_simd_float) max(avx_512_simd_float lhs, avx_512_simd_float rhs) {
        return _mm512_max_ps(lhs.value, rhs.value);
    }

    // Horizontal sum reductions
    // TODO "Vectorize" these reductions

    ETL_STATIC_INLINE(int8_t) hadd(avx_512_simd_byte in) {
        int8_t acc = 0;
        for (size_t i = 0; i < 64; ++i) {
            acc += in[i];
        }
        return acc;
    }

    ETL_STATIC_INLINE(int16_t) hadd(avx_512_simd_short in) {
        int16_t acc = 0;
        for (size_t i = 0; i < 32; ++i) {
            acc += in[i];
        }
        return acc;
    }

    ETL_STATIC_INLINE(int32_t) hadd(avx_512_simd_int in) {
        int32_t acc = 0;
        for (size_t i = 0; i < 16; ++i) {
            acc += in[i];
        }
        return acc;
    }

    ETL_STATIC_INLINE(int64_t) hadd(avx_512_simd_long in) {
        int64_t acc = 0;
        for (size_t i = 0; i < 8; ++i) {
            acc += in[i];
        }
        return acc;
    }

    ETL_STATIC_INLINE(float) hadd(avx_512_simd_float in) {
        return in[0] + in[1] + in[2] + in[3] + in[4] + in[5] + in[6] + in[7] + in[8] + in[9] + in[10] + in[11] + in[12] + in[13] + in[14] + in[15];
    }

    ETL_STATIC_INLINE(double) hadd(avx_512_simd_double in) {
        return in[0] + in[1] + in[2] + in[3] + in[4] + in[5] + in[6] + in[7];
    }

    template <typename T>
    ETL_STATIC_INLINE(T)
    hadd(avx_512_simd_complex_float<T> in) {
        return in[0] + in[1] + in[2] + in[3] + in[4] + in[5] + in[6] + in[7];
    }

    template <typename T>
    ETL_STATIC_INLINE(T)
    hadd(avx_512_simd_complex_double<T> in) {
        return in[0] + in[1] + in[2] + in[3];
    }
};

template <>
ETL_OUT_INLINE(avx_512_simd_byte)
avx512_vec::zero<int8_t>() {
    return _mm512_setzero_si512();
}

template <>
ETL_OUT_INLINE(avx_512_simd_short)
avx512_vec::zero<int16_t>() {
    return _mm512_setzero_si512();
}

template <>
ETL_OUT_INLINE(avx_512_simd_int)
avx512_vec::zero<int32_t>() {
    return _mm512_setzero_si512();
}

template <>
ETL_OUT_INLINE(avx_512_simd_long)
avx512_vec::zero<int64_t>() {
    return _mm512_setzero_si512();
}

template <>
ETL_OUT_INLINE(avx_512_simd_float)
avx512_vec::zero<float>() {
    return _mm512_setzero_ps();
}

template <>
ETL_OUT_INLINE(avx_512_simd_double)
avx512_vec::zero<double>() {
    return _mm512_setzero_pd();
}

template <>
ETL_OUT_INLINE(avx_512_simd_complex_float<etl::complex<float>>)
avx512_vec::zero<etl::complex<float>>() {
    return _mm512_setzero_ps();
}

template <>
ETL_OUT_INLINE(avx_512_simd_complex_double<etl::complex<double>>)
avx512_vec::zero<etl::complex<double>>() {
    return _mm512_setzero_pd();
}

template <>
ETL_OUT_INLINE(avx_512_simd_complex_float<std::complex<float>>)
avx512_vec::zero<std::complex<float>>() {
    return _mm512_setzero_ps();
}

template <>
ETL_OUT_INLINE(avx_512_simd_complex_double<std::complex<double>>)
avx512_vec::zero<std::complex<double>>() {
    return _mm512_setzero_pd();
}

} //end of namespace etl

#endif //__AVX512F__