Как использовать плавное умножение и добавление в AVX для 16-битных упакованных целых чисел

Я знаю, что в AVX2 можно выполнить умножение и сложение с помощью одной инструкции. Я хочу использовать инструкцию умножения и добавления, в которой каждая 256-битная переменная AVX2 упакована 16 16-битными переменными. Например, рассмотрим пример ниже,

разрешение=a0*b0+a1*b1+a2*b2+a3*b3

здесь каждая из res, a0, a1, a2, a3, b0, b1, b2, b3 — 16-битные переменные. Я внимательно следил за обсуждением . Пожалуйста, найдите мой код ниже, чтобы рассчитать пример, показанный выше,

#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 возвращается,

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 немного быстрее, но я ожидал, что он будет еще быстрее. Я понимаю, что в моем примере кода есть много обращений к памяти, которые могут быть причиной ухудшения производительности. Но,

1. Когда я выгружаю сборку, я вижу почти одинаковые инструкции для обоих подходов. Я не вижу никаких инструкций по fma в версии FMA. Я не понимаю причины. Это из-за инструкций _mm256_mullo_epi16?

2. Верен ли мой подход?

3. Не могли бы вы помочь мне исправить это?

Я новичок в программировании AVX2, поэтому вполне возможно, что я сделал что-то не очень стандартное, но я буду рад ответить на что-то непонятное. Я благодарю вас всех за вашу помощь заранее.