如何在 AVX 中对 16 位压缩整数使用融合乘法和加法
How to use fused multiply and add in AVX for 16 bit packed integers
我知道在 AVX2 中可以使用一条指令进行乘加。我想使用乘加指令,其中每个 256 位 AVX2 变量都包含 16 个 16 位变量。例如,考虑下面的例子,
res=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 ,%%rdx;orq %%rdx,%%rax"
: "=a" (result) :: "%rdx");
return result;
}
我的 gcc -version returns,
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 之外没有 SIMD 整数 FMA / MAC(乘法累加),它将水平添加到更宽的整数中。 (直到 AVX512IFMA _mm_madd52lo_epu64
or AVX512_4VNNIW _mm512_4dpwssd_epi32(__m512i, __m512ix4, __m128i *)
)。
FP-contract 和 -ffast-math
选项与 SIMD-integer 无关;整数数学总是精确的。
我认为您的 "simple" 方法较慢,因为您还修改了输入数组,而这并没有得到优化,例如
a0_avx[i] = _mm256_mullo_epi16(a0_avx[i], b0_avx[i]);
以及更新 res_avx[i]
.
如果编译器没有优化它,那些额外的存储可能正是它比您的 mul_add
函数慢的原因。 rdtsc
没有序列化指令甚至不必等待更早的指令执行,更不用说将存储退回或提交到 L1d 缓存,但前端的额外微指令仍然需要仔细考虑。每个时钟吞吐量只有 1 个存储,这很容易成为新的瓶颈。
仅供参考,您不需要将输入复制到 __m256i
的数组中。通常,您只需对常规数据使用 SIMD 加载。这并不比 __m256i
的索引数组慢。您的数组太大,编译器无法完全展开并将所有内容保存在寄存器中(就像标量 __m256i
变量一样)。
如果你只是在循环中使用 __m256i a0 = _mm256_loadu_si256(...)
那么你可以更新 a0
而不会减慢你的代码,因为它只是一个可以保存在 YMM 中的局部变量注册
但我发现在大多数步骤中使用新命名的 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 中可以使用一条指令进行乘加。我想使用乘加指令,其中每个 256 位 AVX2 变量都包含 16 个 16 位变量。例如,考虑下面的例子,
res=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 ,%%rdx;orq %%rdx,%%rax"
: "=a" (result) :: "%rdx");
return result;
}
我的 gcc -version returns,
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 之外没有 SIMD 整数 FMA / MAC(乘法累加),它将水平添加到更宽的整数中。 (直到 AVX512IFMA _mm_madd52lo_epu64
or AVX512_4VNNIW _mm512_4dpwssd_epi32(__m512i, __m512ix4, __m128i *)
)。
FP-contract 和 -ffast-math
选项与 SIMD-integer 无关;整数数学总是精确的。
我认为您的 "simple" 方法较慢,因为您还修改了输入数组,而这并没有得到优化,例如
a0_avx[i] = _mm256_mullo_epi16(a0_avx[i], b0_avx[i]);
以及更新 res_avx[i]
.
如果编译器没有优化它,那些额外的存储可能正是它比您的 mul_add
函数慢的原因。 rdtsc
没有序列化指令甚至不必等待更早的指令执行,更不用说将存储退回或提交到 L1d 缓存,但前端的额外微指令仍然需要仔细考虑。每个时钟吞吐量只有 1 个存储,这很容易成为新的瓶颈。
仅供参考,您不需要将输入复制到 __m256i
的数组中。通常,您只需对常规数据使用 SIMD 加载。这并不比 __m256i
的索引数组慢。您的数组太大,编译器无法完全展开并将所有内容保存在寄存器中(就像标量 __m256i
变量一样)。
如果你只是在循环中使用 __m256i a0 = _mm256_loadu_si256(...)
那么你可以更新 a0
而不会减慢你的代码,因为它只是一个可以保存在 YMM 中的局部变量注册
但我发现在大多数步骤中使用新命名的 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.