如何在 AVX 中使用融合乘法和加法来处理 16 位压缩整数

Posted 2023-02-16

【中文标题】如何在 AVX 中使用融合乘法和加法来处理 16 位压缩整数

【英文标题】：How to use fused multiply and add in AVX for 16 bit packed integers

【发布时间】：2019-07-31 09:11:57

【问题描述】：

我知道在 AVX2 中可以使用一条指令进行乘加。我想使用乘加指令，其中每个 256 位 AVX2 变量都包含 16 个 16 位变量。例如，考虑下面的例子，

分辨率=a0*b0+a1*b1+a2*b2+a3*b3

这里的 res、a0、a1、a2、a3、b0、b1、b2、b3 都是 16 位变量。 我密切关注discussion。请在下面找到我的代码来计算上面显示的示例，

#include<stdio.h>
#include<stdint.h>
#include <immintrin.h>
#include<time.h>
#include "cpucycles.c"

#pragma STDC FP_CONTRACT ON

#define AVX_LEN 16

inline __m256i mul_add(__m256i a, __m256i b, __m256i c)  
    return _mm256_add_epi16(_mm256_mullo_epi16(a, b), c);


void fill_random(int16_t *a, int32_t len)  //to fill up the random array

    int32_t i;

    for(i=0;i<len;i++)     
        a[i]=(int16_t)rand()&0xffff;
    


void main()


    int16_t a0[16*AVX_LEN], b0[16*AVX_LEN];
    int16_t a1[16*AVX_LEN], b1[16*AVX_LEN];
    int16_t a2[16*AVX_LEN], b2[16*AVX_LEN];
    int16_t a3[16*AVX_LEN], b3[16*AVX_LEN];
    int16_t res[16*AVX_LEN];


    __m256i a0_avx[AVX_LEN], b0_avx[AVX_LEN];
    __m256i a1_avx[AVX_LEN], b1_avx[AVX_LEN];
    __m256i a2_avx[AVX_LEN], b2_avx[AVX_LEN];
    __m256i a3_avx[AVX_LEN], b3_avx[AVX_LEN];

    __m256i res_avx[AVX_LEN];

    int16_t res_avx_check[16*AVX_LEN];
    int32_t i,j;

    uint64_t mask_ar[4]; //for unloading AVX variables
    mask_ar[0]=~(0UL);mask_ar[1]=~(0UL);mask_ar[2]=~(0UL);mask_ar[3]=~(0UL);
    __m256i mask;
    mask = _mm256_loadu_si256 ((__m256i const *)mask_ar);

    time_t t;
    srand((unsigned) time(&t));


    int32_t repeat=100000;

    uint64_t clock1, clock2, fma_clock;

    clock1=clock2=fma_clock=0;

    for(j=0;j<repeat;j++)
        printf("j : %d\n",j);

        fill_random(a0,16*AVX_LEN);// Genrate random data
        fill_random(a1,16*AVX_LEN);
        fill_random(a2,16*AVX_LEN);
        fill_random(a3,16*AVX_LEN);

        fill_random(b0,16*AVX_LEN);
        fill_random(b1,16*AVX_LEN);
        fill_random(b2,16*AVX_LEN);
        fill_random(b3,16*AVX_LEN);


        for(i=0;i<AVX_LEN;i++) //Load values in AVX variables

            a0_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&a0[i*16]));
            a1_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&a1[i*16]));
            a2_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&a2[i*16]));
            a3_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&a3[i*16]));

            b0_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&b0[i*16]));
            b1_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&b1[i*16]));
            b2_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&b2[i*16]));
            b3_avx[i] = _mm256_loadu_si256 ((__m256i const *) (&b3[i*16]));
        

        for(i=0;i<AVX_LEN;i++)
            res_avx[i]= _mm256_set_epi64x(0, 0, 0, 0);
        

        //to calculate a0*b0 + a1*b1 + a2*b2 + a3*b3

        //----standard calculation----
        for(i=0;i<16*AVX_LEN;i++)
            res[i]=a0[i]*b0[i] + a1[i]*b1[i] + a2[i]*b2[i] + a3[i]*b3[i];
        


        //-----AVX-----

        clock1=cpucycles();


        for(i=0;i<AVX_LEN;i++) //simple approach

            a0_avx[i]=_mm256_mullo_epi16(a0_avx[i], b0_avx[i]);
            res_avx[i]=_mm256_add_epi16(a0_avx[i], res_avx[i]);

            a1_avx[i]=_mm256_mullo_epi16(a1_avx[i], b1_avx[i]);
            res_avx[i]=_mm256_add_epi16(a1_avx[i], res_avx[i]);

            a2_avx[i]=_mm256_mullo_epi16(a2_avx[i], b2_avx[i]);
            res_avx[i]=_mm256_add_epi16(a2_avx[i], res_avx[i]);

            a3_avx[i]=_mm256_mullo_epi16(a3_avx[i], b3_avx[i]);
            res_avx[i]=_mm256_add_epi16(a3_avx[i], res_avx[i]);

        

        /*
        for(i=0;i<AVX_LEN;i++) //FMA approach

            res_avx[i]=mul_add(a0_avx[i], b0_avx[i], res_avx[i]);

            res_avx[i]=mul_add(a1_avx[i], b1_avx[i], res_avx[i]);
            res_avx[i]=mul_add(a2_avx[i], b2_avx[i], res_avx[i]);

            res_avx[i]=mul_add(a3_avx[i], b3_avx[i], res_avx[i]);

        
        */

        clock2=cpucycles();
        fma_clock = fma_clock + (clock2-clock1);

        //-----Check----

        for(i=0;i<AVX_LEN;i++) //store avx results for comparison
            _mm256_maskstore_epi64 (res_avx_check + i*16, mask, res_avx[i]);
        

        for(i=0;i<16*AVX_LEN;i++)
            if(res[i]!=res_avx_check[i])

                printf("\n--ERROR--\n");
                return;
               

        
    


    printf("Total time taken is :%llu\n", fma_clock/repeat);

cpucycles 代码来自ECRYPT，如下所示，

#include "cpucycles.h"

long long cpucycles(void)

  unsigned long long result;
  asm volatile(".byte 15;.byte 49;shlq $32,%%rdx;orq %%rdx,%%rax"
    : "=a" (result) ::  "%rdx");
  return result;

我的 gcc -version 返回，

gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-36)

我正在使用

 Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz

当我在我的计算机上运行它时，我分别得到以下 fma 方法和简单方法的循环

FMA approach : Total time taken is :109
Simple approach : Total time taken is :141

如您所见，FMA 方法稍微快一些，但我预计会更快。我知道在我的示例代码中有许多内存访问，这可能是性能下降的原因。但是，

当我转储程序集时，我看到两种方法的说明几乎相似。我在 FMA 版本中没有看到任何 fma 指令。我不明白原因。是因为_mm256_mullo_epi16指令吗？

我的方法正确吗？

你能帮我解决这个问题吗？

我是 AVX2 编程的新手，所以我很有可能做了一些不是很标准的事情，但我很乐意回答一些不清楚的事情。 提前感谢大家的帮助。

【问题讨论】：

x86 除了水平 pmaddubsw / pmaddwd 之外没有整数 FMA / MAC（乘法累加）。 FP 合约期权与整数无关；整数数学总是精确的。

@PeterCordes 我不明白。我说的是向量指令。

[v]pmaddubsw 和 [v]pmaddwd 是 SIMD SSE / AVX2 指令。我的意思是 SIMD 整数，当然，不是 GP 寄存器中的标量。

@PeterCordes 谢谢。我现在知道了。我不明白当我使用 FMA 方法时为什么时钟周期会减少

您使用了哪些编译器选项？您是否忘记启用优化？ gcc -O3 -march=native。实际上你的编译器太老了，不知道-march=skylake 所以-march=native 不会正确设置-mtune 设置。但无论如何，优化后应该没有区别。

【参考方案1】：

x86 除了水平 pmaddubsw / pmaddwd 之外没有 SIMD 整数 FMA / MAC（乘法累加），它将水平添加到更宽的整数中。 （直到 AVX512IFMA _mm_madd52lo_epu64 或 AVX512_4VNNIW _mm512_4dpwssd_epi32(__m512i, __m512ix4, __m128i *)）。

FP-contract 和-ffast-math 选项与SIMD-integer 无关；整数数学总是精确的。

---

我认为您的“简单”方法较慢，因为您还在修改输入数组，而这并没有得到优化，例如

a0_avx[i] = _mm256_mullo_epi16(a0_avx[i], b0_avx[i]);

以及更新res_avx[i]。

如果编译器不对其进行优化，那么这些额外的存储可能正是它比您的 mul_add 函数慢的原因。 rdtsc 没有序列化指令甚至不必等待更早的指令执行，更不用说退休或提交存储到 L1d 缓存，但是前端的额外 uops 仍然需要更多的咀嚼。每个时钟吞吐量只有 1 个存储，这很容易成为新的瓶颈。

---

仅供参考，您无需将输入复制到 __m256i 的数组中。通常，您只需对常规数据使用 SIMD 加载。这并不比索引__m256i 的数组慢。您的数组太大，编译器无法完全展开并将所有内容保存在寄存器中（就像标量 __m256i 变量一样）。

如果您只是在循环中使用 __m256i a0 = _mm256_loadu_si256(...)，那么您可以更新 a0 而不会减慢您的代码速度，因为它只是一个可以保存在 YMM reg 中的局部变量。

但我发现在大多数步骤中使用新命名的 tmp 变量是一种很好的风格，可以使代码更具自我记录性。喜欢__m256i ab = ... 或sum = ...。您可以为每个 a0+b0 和 a1+b1 重复使用相同的 sum 临时。

您也可以对结果向量使用临时变量，而不是让编译器优化 res_avx[i] 处的内存更新，直到最后一个。

您可以使用alignas(32) int16_t a0[...]; 使普通数组与_mm256_load 对齐，而不是loadu。

---

您的cpucycles() RDTSC 函数不需要使用内联汇编。 Use __rdtsc() instead.

【讨论】：

感谢您的图解回答。你已经写了“你也在修改输入数组，这并没有得到优化”，你能告诉我像这样更新 AVX2 数组的最佳方法是什么。感谢您的帮助。

@Rick：首先一般不要使用__m256d 的数组，只需使用一些 __m256d 局部变量，以便将它们保存在寄存器中。如果您确实有数组（无论它们是 double 还是 __m256d），请使用 tmp var 而不是修改它们，除非您真的想要那种副作用。