两个gcc_GCC神坑：-march=native

曾几何时，当我的编译原理还在60分上下徘徊的时候，我从未想过，会有一天，我会用avx指令写程序，然后一脚踏入-march=native这个神坑。

毕竟，当初的编译原理课，我一共就学会了两个命令

一个叫-O3，一个叫-march=native

凭借这两个命令，我把手里的程序加速加速再加速……从未翻车——直到某日我发现电脑CPU不够多——于是，我找老板借了服务器……一台Centos服务器

灾难……降临了

---

在揭晓答案之前，我想先附一段会随机崩溃的程序：这是一个用fma（每一个函数都是查Intel Intrinsics Guide找到的）写的矩阵乘法，大概就是R语言用.Call调用MmV，MmV调用matmul，然后matmul调用真正干活的mmv，matmul到mmv这里用到了循环展开的技巧，把mat先分段，然后每段执行mmv（比如说{{a1,a6,a11},{a2,a7,a12},{a3,a8,a13},{a4,a9,a14},{a5,a10,a15}}，此时vec是一个3维的向量{v1,v2,v3}，我们希望得到结果{a1v1+a6v2+a11v3,a2v1+a7v2+a12v3,...}——这时候我们把mat拆成两部分，一部分长度是4的倍数，另一部分是剩下的，实验&Guide证明，mmv可以处理4的倍数，于是第一种情况解决，对于第二种情况，我们可能会多读几个double数据，不过没关系，为了这玩意我贴心地算了一个mask出来，不管mmv算出来的是什么鬼，反正mask之外的都会被扔掉）

事实上，如果用"/*"注释掉R语言的部分打开main的注释（同时需要一台支持AVX2/FMA的电脑），就会发现，main部分的计算是可以通过测试的，但如果重复若干次（R，重复1000次某一个大矩阵*一个向量）之后，程序会直接爆掉——而神奇的是，在爆掉之前，我们总能得到完全正确的答案。

#include<immintrin.h>
//void _mm256_maskstore_pd (double * mem_addr, __m256i mask, __m256d a)
//__m256d _mm256_loadu_pd (double const * mem_addr)
//__m256d _mm256_setzero_pd (void)
//void _mm256_maskstore_pd (double * mem_addr, __m256i mask, __m256d a)

inline __m256d mmv(register double*mat,register double*vec,register unsigned int matlen,register int veclen){
  register __m256d a,b,c;
  c=_mm256_setzero_pd();
  for(register int i=0;i<veclen;i++){
    a= _mm256_loadu_pd (mat);
    b= _mm256_set1_pd(*vec);
    c= _mm256_fmadd_pd (a,b,c);
    vec++;
    mat+=matlen;
  }
  return c;
}

void matmul(double*res,double*mat,double*vec,unsigned int matlen,unsigned int veclen){
  int c;
  for(c=0;c<(matlen&~3);c+=4){
    _mm256_storeu_pd(res+c,mmv(mat+c,vec,matlen,veclen));
  }
  __int64 d=c-matlen;
  if(d){
    register __m256i mask=_mm256_set_epi64x(2-d,1-d,-d,-1LL);
    _mm256_maskstore_pd(res+c,mask,mmv(mat+c,vec,matlen,veclen));
  }
  _mm256_zeroupper();
}
//#define DEBUG
#ifdef DEBUG
  #include<R.h>
  #define DEBUG_printf Rprintf
#else
  #define DEBUG_printf(x...)
#endif
#include <Rinternals.h>
#include<stdlib.h>
SEXP MmV(SEXP M,SEXP V){
  double *mat=REAL(PROTECT(coerceVector(M, REALSXP)));
  double *vec=REAL(PROTECT(coerceVector(V, REALSXP)));
  int n=length(V);
  int m=length(M)/n;
  SEXP RES=PROTECT(Rf_allocVector(REALSXP, m));
  double *res=REAL(RES);
DEBUG_printf("n=%d,m=%d,mat[1,1]=%f,vec[1,1]=%lfn",n,m,*mat,*vec);
  matmul(res,mat,vec,m,n);
  UNPROTECT(3);
  return RES;
}
/*/
#include<stdio.h>
int main(){
int zero0[100];
  double mat[100];//20*5
int zero1[100];
  double vec[5];
int zero2[100];
  double res[20];
int zero3[100];
  int c,m,n;
  for(c=0;c<100;c++)zero0[c]=zero1[c]=zero2[c]=zero3[c]=0;
printf("m? n?n");
  scanf("%d %d",&m,&n);//20 5
printf("vec?n");
  for(c=0;c<5;c++)scanf("%lf",&vec[c]);
printf("mat?n");
  for(c=0;c<100;c++)scanf("%lf",&mat[c]);
printf("start!n");
  matmul(res,mat,vec,m,n);
printf("end.n");
  for(c=0;c<100;c++)if(!(zero0[c]==0&&zero1[c]==0&&zero2[c]==0&&zero3[c]==0))printf("leak!n");
printf("Checked.n");
  for(c=0;c<20;c++)printf("%.18f ",res[c]);
}*/

---

事情是这样的……自从我成为了我们组里编程最6的存在之后……一路几乎从未翻车，而同学也经常有向我请教程序的时候。

某日，有同学写了一个程序，送去服务器里面的R Open，CPU占用率飙升至1600%，心乡往之，毕竟R的并行从来都是炒鸡难写的，而那个同学又不像是会手写并行的样子。

细细检查，发现……R Open自带的矩阵乘向量的算法就是一个坑

一怒之下，我写了一个AVX

-调试成功-测试失败-调整程序-测试成功-送服务器-服务器炸了……

emmm，服务器是分好几个阶段炸的

比如，服务器的gcc是6.1.0，gfortran是4.4.7……神奇的是，cat /proc/cpuinfo里面显示，CPU是Xeon，支持fma，更仔细的检查说明，这玩意的架构是knights landing，22nm，25MB缓存（应该是L3）……忽然感觉到被嘲讽了……

（服务器：我的CPU不止是最快的，还是最大的）

看到这种神奇的配置，第一件事当然是……把gcc升级到8.3

BUT。。。我并不是管理员，并没有yum install的权限，无奈之下只好强行编译一只gcc

还好学过怎么编译

不就是make个j8么

设CC=gcc，CFLAGS="-pipe -O3 -march=native"

然后gmp告诉我

no such instruction: `shlx %rdx,%rax,%rax'
no such instruction: `shlx %rdx,%rax,%rax'
no such instruction: `shlx %rdx,%rax,%rax'
no such instruction: `shlx %rdx,%rax,%rax'
no such instruction: `shlx %rdx,%rax,%rax'
no such instruction: `shlx %rdx,%rax,%rax'
no such instruction: `shlx %rdx,%rax,%rax'
no such instruction: `shlx %rdx,%rax,%rax'

刷出这一屏，基本上make -j8就可以宣告

——你make的是个j8啊

……

这是我第一次遇到类似的问题，而问题的答案很简单，CentOS太老了，这个CentOS用的汇编器as并不支持哪怕是shlx这样的指令

果然，还是Gentoo保平安的好些。

---

至于这个神奇的程序……就在我觉得“这次肯定找到原因了”的时候，再一次BOOM了……而在爆炸了几天之后，我忽然发现，最开始的想法是正确的——当CPU跟硬盘都不可能出问题的时候，唯一能出问题的，只有内存了

（事实上这个BUG的确是一个内存有关的问题）