SIMD 寄存器的数学函数

Ste*_*ven 6 gcc simd avx

根据https://sourceware.org/glibc/wiki/libmvec GCC 具有数学函数的矢量实现。它们可以被编译器用于优化,可以在这个例子中看到: https: //godbolt.org/g/IcxtVi,编译器使用一些损坏的正弦函数并一次对 4 个双精度数进行操作

我知道如果我需要数学函数,可以使用 SIMD 数学库,但我仍然感兴趣是否有一种方法可以使用__m256d某种内在函数或任何其他方式手动调用 GCC 中已经存在的向量化数学函数?

wim*_*wim 9

Libmvec 是一个 x86_64 glibc 库,具有 SSE4、AVX、AVX2 和 AVX-512 向量化函数,用于 单精度和双精度的cos、exp、log、sin、pow和sincos。这些函数的精度为 4-ulp 最大相对误差。

通常,gcc 会插入对 Libmvec 函数的调用,例如_ZGVdN4v_cos,而它是自动向量化标量代码,例如使用-ffast-math或 #pragma omp simd。然而,这些向量函数也适用于手动向量化 C 代码。不幸的是,Libmvec 文档 没有给出很多示例来解释如何使用这些函数来手动矢量化 C 代码。

要使用 Libmvec 函数,必须包含函数声明并将代码与-lm(与数学库链接)选项链接。链接器不需要-lmvec,链接标准数学库-lm对于非静态构建就足够了,请参见此处。编译器选项-ffast-math不是必需的。

以下示例显示如何调用 AVX2 向量化函数。此示例在标准 Ubuntu 18.04 系统(glibc 版本 2.27)上成功运行。

/* gcc -Wall -O3 -m64 -march=skylake libmvec_ex.c -lm */ 
#include <immintrin.h>
#include <stdio.h>
#include <stdint.h>
#define _GNU_SOURCE      /* These two lines are optional. They are only needed to use the scalar */
#include <math.h>        /* functions sin, exp, sincos,...                                       */

__m256d _ZGVdN4v_cos(__m256d x);                    
__m256d _ZGVdN4v_exp(__m256d x);                    
__m256d _ZGVdN4v_log(__m256d x);                    
__m256d _ZGVdN4v_sin(__m256d x);                    
__m256d _ZGVdN4vv_pow(__m256d x, __m256d y); 
void    _ZGVdN4vvv_sincos(__m256d x, __m256i ptrs, __m256i ptrc);

__m256  _ZGVdN8v_cosf(__m256 x);            
__m256  _ZGVdN8v_expf(__m256 x);            
__m256  _ZGVdN8v_logf(__m256 x);            
__m256  _ZGVdN8v_sinf(__m256 x);            
__m256  _ZGVdN8vv_powf(__m256 x, __m256 y); 
void    _ZGVdN8vvv_sincosf(__m256 x, __m256i ptrs_lo, __m256i ptrs_hi, __m256i ptrc_lo, __m256i ptrc_hi);

int mm256_print_pd(__m256d x);
int mm256_print_ps(__m256 x);

int main()
{
  __m256d x, y, z;
  __m256 xf, yf, zf;
  double z_s[4], z_c[4];          /* sincos output                              */
  float zf_s[8], zf_c[8];         /* sincosf output                             */
  __m256i ptrs, ptrc;             /* Pointers to the elements of z_s and z_c    */ 
  __m256i ptrs_lo, ptrs_hi, ptrc_lo, ptrc_hi;

  x = _mm256_set_pd (0.04, 0.03, 0.02, 0.01);
  y = _mm256_set_pd (2.0, 1.5, 1.0, 0.5);
  xf = _mm256_set_ps (0.08f, 0.07f, 0.06f, 0.05f, 0.04f, 0.03f, 0.02f, 0.01f);
  yf = _mm256_set_ps (4.0f, 3.5f, 3.0f, 2.5f, 2.0f, 1.5f, 1.0f, 0.5f);

  printf("AVX2 Double precision examples\n");
                             printf("x             "); mm256_print_pd(x);
                             printf("y             "); mm256_print_pd(y);
  z =_ZGVdN4v_cos(x);        printf("cos(x)        "); mm256_print_pd(z);
  z =_ZGVdN4v_exp(x);        printf("exp(x)        "); mm256_print_pd(z);
  z =_ZGVdN4v_log(x);        printf("log(x)        "); mm256_print_pd(z);
  z =_ZGVdN4v_sin(x);        printf("sin(x)        "); mm256_print_pd(z);
  z =_ZGVdN4vv_pow(x, y);    printf("pow(x,y)      "); mm256_print_pd(z);

  ptrs = _mm256_set_epi64x((uint64_t)&z_s[3],(uint64_t)&z_s[2],(uint64_t)&z_s[1],(uint64_t)&z_s[0]);
  ptrc = _mm256_set_epi64x((uint64_t)&z_c[3],(uint64_t)&z_c[2],(uint64_t)&z_c[1],(uint64_t)&z_c[0]);
/* Alternative: ptrs = _mm256_add_epi64(_mm256_set1_epi64x((uint64_t)&z_s[0]),_mm256_set_epi64x(24,16,8,0));  */
/* This might be more efficient if the destination addresses are contiguous in memory.                        */
  _ZGVdN4vvv_sincos(x, ptrs, ptrc);                                 /* The results of _ZGVdN4vvv_sincos are scattered into the adresses in ptrs and ptrc */ 
  printf("sincos cos(x) %12.8f %12.8f %12.8f %12.8f  \n", z_c[3], z_c[2], z_c[1], z_c[0]);
  printf("sincos sin(x) %12.8f %12.8f %12.8f %12.8f\n\n", z_s[3], z_s[2], z_s[1], z_s[0]);

  printf("AVX2 Single precision examples\n");
                             printf("x             "); mm256_print_ps(xf);
                             printf("y             "); mm256_print_ps(yf);
  zf =_ZGVdN8v_cosf(xf);     printf("cosf(x)       "); mm256_print_ps(zf);
  zf =_ZGVdN8v_expf(xf);     printf("expf(x)       "); mm256_print_ps(zf);
  zf =_ZGVdN8v_logf(xf);     printf("logf(x)       "); mm256_print_ps(zf);
  zf =_ZGVdN8v_sinf(xf);     printf("sinf(x)       "); mm256_print_ps(zf);
  zf =_ZGVdN8vv_powf(xf, yf);printf("powf(x,y)     "); mm256_print_ps(zf);

  ptrs_lo = _mm256_set_epi64x((uint64_t)&zf_s[3],(uint64_t)&zf_s[2],(uint64_t)&zf_s[1],(uint64_t)&zf_s[0]);
  ptrs_hi = _mm256_set_epi64x((uint64_t)&zf_s[7],(uint64_t)&zf_s[6],(uint64_t)&zf_s[5],(uint64_t)&zf_s[4]);
  ptrc_lo = _mm256_set_epi64x((uint64_t)&zf_c[3],(uint64_t)&zf_c[2],(uint64_t)&zf_c[1],(uint64_t)&zf_c[0]);
  ptrc_hi = _mm256_set_epi64x((uint64_t)&zf_c[7],(uint64_t)&zf_c[6],(uint64_t)&zf_c[5],(uint64_t)&zf_c[4]);
  _ZGVdN8vvv_sincosf(xf, ptrs_lo, ptrs_hi, ptrc_lo, ptrc_hi);       /* The results of _ZGVdN8vvv_sincosf are scattered to the adresses in ptrs_lo, ptrs_hi, ptrc_lo, and ptrc_hi */ 
  printf("sincosf cos(x)%12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f \n", 
           zf_c[7], zf_c[6], zf_c[5], zf_c[4], zf_c[3], zf_c[2], zf_c[1], zf_c[0]);
  printf("sincosf sin(x)%12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f  \n",
           zf_s[7], zf_s[6], zf_s[5], zf_s[4], zf_s[3], zf_s[2], zf_s[1], zf_s[0]);

  return 0;
}


__attribute__ ((noinline)) int mm256_print_pd(__m256d x){
    double vec_x[4];
    _mm256_storeu_pd(vec_x,x);
    printf("%12.8f %12.8f %12.8f %12.8f  \n", vec_x[3], vec_x[2], vec_x[1], vec_x[0]);
    return 0;
}


__attribute__ ((noinline)) int mm256_print_ps(__m256 x){
    float vec_x[8];
    _mm256_storeu_ps(vec_x,x);
    printf("%12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f \n", vec_x[7], vec_x[6], vec_x[5], vec_x[4],
                     vec_x[3], vec_x[2], vec_x[1], vec_x[0]);
    return 0;
}
/* 
Run Code Online (Sandbox Code Playgroud)

输出如预期:

$ ./a.out
AVX2 Double precision examples
x               0.04000000   0.03000000   0.02000000   0.01000000  
y               2.00000000   1.50000000   1.00000000   0.50000000  
cos(x)          0.99920011   0.99955003   0.99980001   0.99995000  
exp(x)          1.04081077   1.03045453   1.02020134   1.01005017  
log(x)         -3.21887582  -3.50655790  -3.91202301  -4.60517019  
sin(x)          0.03998933   0.02999550   0.01999867   0.00999983  
pow(x,y)        0.00160000   0.00519615   0.02000000   0.10000000  
sincos cos(x)   0.99920011   0.99955003   0.99980001   0.99995000  
sincos sin(x)   0.03998933   0.02999550   0.01999867   0.00999983

AVX2 Single precision examples
x               0.08000000   0.07000000   0.06000000   0.05000000   0.04000000   0.03000000   0.02000000   0.01000000 
y               4.00000000   3.50000000   3.00000000   2.50000000   2.00000000   1.50000000   1.00000000   0.50000000 
cosf(x)         0.99680173   0.99755102   0.99820060   0.99875027   0.99920011   0.99955004   0.99979997   0.99995005 
expf(x)         1.08328700   1.07250810   1.06183648   1.05127108   1.04081070   1.03045452   1.02020133   1.01005018 
logf(x)        -2.52572870  -2.65926003  -2.81341076  -2.99573231  -3.21887589  -3.50655794  -3.91202307  -4.60517025 
sinf(x)         0.07991469   0.06994285   0.05996400   0.04997917   0.03998933   0.02999550   0.01999867   0.00999983 
powf(x,y)       0.00004096   0.00009075   0.00021600   0.00055902   0.00160000   0.00519615   0.02000000   0.10000000 
sincosf cos(x)  0.99680173   0.99755102   0.99820060   0.99875027   0.99920011   0.99955004   0.99979997   0.99995005 
sincosf sin(x)  0.07991469   0.06994285   0.05996400   0.04997917   0.03998933   0.02999550   0.01999867   0.00999983  
Run Code Online (Sandbox Code Playgroud)

请注意,矢量化sincos函数将结果分散到 4 个(双精度型)或 8 个(浮点型)(不一定是连续的)内存位置。


SSE4、AVX 和 AVX2 的函数声明,以及 Intel SVML 替代方案。

下面的代码块显示了许多不同情况的函数声明,以及来自英特尔 SVML 库的相应 SVML 函数。不知何故,我在互联网上找不到这些带有正确参数列表的声明。

$ ./a.out
AVX2 Double precision examples
x               0.04000000   0.03000000   0.02000000   0.01000000  
y               2.00000000   1.50000000   1.00000000   0.50000000  
cos(x)          0.99920011   0.99955003   0.99980001   0.99995000  
exp(x)          1.04081077   1.03045453   1.02020134   1.01005017  
log(x)         -3.21887582  -3.50655790  -3.91202301  -4.60517019  
sin(x)          0.03998933   0.02999550   0.01999867   0.00999983  
pow(x,y)        0.00160000   0.00519615   0.02000000   0.10000000  
sincos cos(x)   0.99920011   0.99955003   0.99980001   0.99995000  
sincos sin(x)   0.03998933   0.02999550   0.01999867   0.00999983

AVX2 Single precision examples
x               0.08000000   0.07000000   0.06000000   0.05000000   0.04000000   0.03000000   0.02000000   0.01000000 
y               4.00000000   3.50000000   3.00000000   2.50000000   2.00000000   1.50000000   1.00000000   0.50000000 
cosf(x)         0.99680173   0.99755102   0.99820060   0.99875027   0.99920011   0.99955004   0.99979997   0.99995005 
expf(x)         1.08328700   1.07250810   1.06183648   1.05127108   1.04081070   1.03045452   1.02020133   1.01005018 
logf(x)        -2.52572870  -2.65926003  -2.81341076  -2.99573231  -3.21887589  -3.50655794  -3.91202307  -4.60517025 
sinf(x)         0.07991469   0.06994285   0.05996400   0.04997917   0.03998933   0.02999550   0.01999867   0.00999983 
powf(x,y)       0.00004096   0.00009075   0.00021600   0.00055902   0.00160000   0.00519615   0.02000000   0.10000000 
sincosf cos(x)  0.99680173   0.99755102   0.99820060   0.99875027   0.99920011   0.99955004   0.99979997   0.99995005 
sincosf sin(x)  0.07991469   0.06994285   0.05996400   0.04997917   0.03998933   0.02999550   0.01999867   0.00999983  
Run Code Online (Sandbox Code Playgroud)


有用的链接:

https://sourceware.org/glibc/wiki/libmvec

https://sourceware.org/glibc/wiki/libmvec?action=AttachFile&do=view&target=VectorABI.txt

https://github.com/sgallagher/glibc/blob/master/sysdeps/unix/sysv/linux/x86_64/libmvec.abilist

https://software.intel.com/sites/default/files/management/b4/c8/Intel-Vector-Function-ABI.pdf