Simd matmul 程序给出不同的数值结果

Posted

技术标签:

【中文标题】Simd matmul 程序给出不同的数值结果【英文标题】:Simd matmul program gives different numerical results 【发布时间】:2019-04-02 14:51:51 【问题描述】:

我正在尝试使用 simd 内在函数在 C 中编写矩阵乘法。我很确定我的实现,但是当我执行时,从结果矩阵系数的第 5 位开始出现一些数值错误。

REAL_T 只是一个带有 typedef 的浮点数

/* This is my matmul Version with simd, using floating simple precision*/
void matmul(int n, REAL_T *A, REAL_T *B, REAL_T *C)
  int i,j,k;
  __m256 vA, vB, vC, vRes;
  for (i=0; i<n; i++)
    for (j=0; j<n; j++)  
      for (k=0; k<n; k= k+8)
        vA = _mm256_load_ps(&A[i*n+k]);
        vB = _mm256_loadu_ps(&B[k*n+j]);
        vC = _mm256_mul_ps(vA, vB);
        vC = _mm256_hadd_ps(vC, vC);
        vC = _mm256_hadd_ps(vC, vC);
        /*To get the resulting coefficient, after doing 2 hadds,
        I have to get the first and the last element of the resulting
        Vector vC*/
        C[i*n+j] += ((float )(vC[0])) + ((float )(vC[7]));
       /* for k */
     /* for j */
   /* for i */

*/End of program
/*And this is the sequential Version*/
void matmul(int n, REAL_T *A, REAL_T *B, REAL_T *C)
  int i,j,k;
  for (i=0; i<n; i++) 
    for (j=0; j<n; j++)
      for (k=0; k<n; k++)
        C[i*n+j] +=  A[i*n+k] *  B[k*n+j];  
       /* for k */
     /* for j */
   /* for i */  

/*End of program*/
/*The matrix are initialized as follows*/
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      *(A+i*n+j) = 1 / ((REAL_T) (i+j+1));
      *(B+i*n+j) = 1.0;
      *(C+i*n+j) = 1.0;
    
/*End of initialization*/

测试的矩阵大小为 512*512。 对于顺序版本,结果矩阵的左上角给出:

+6.916512e+01  +6.916512e+01  
+5.918460e+01  +5.918460e+01  

+7.946186e+00  +7.946186e+00  
+7.936391e+00  +7.936391e+00  

但是,对于 simd 版本,正方形是:

+6.916510e+01  +6.916510e+01  
+5.918463e+01  +5.918463e+01  

+7.946147e+00  +7.946147e+00  
+7.936355e+00  +7.936355e+00 

如图所示,两个版本之间存在数字错误。 任何帮助将不胜感激!

【问题讨论】:

【参考方案1】:

这看起来很正常;以不同的顺序添加数字会在临时变量中产生不同的舍入。

FP 数学不是关联的;优化就好像它会改变结果。1Is Floating point addition and multiplication associative? / Are floating point operations in C associative?

变化量取决于数据。对于 float,仅小数点后 5 位的差异似乎是合理的。


除非您采取特殊的数字预防措施,例如先将小数字相加,否则顺序结果并不是“更正确”,它们只是有不同的错误。

事实上,对于大型列表,使用多个累加器通常会提高精度,假设您的数字都具有相似的量级。 (理想情况下,每个由多个元素组成的多个 SIMD 向量,以隐藏 FP-add 或 FMA 延迟)。 https://en.wikipedia.org/wiki/Pairwise_summation 是一种数值技术,将其提升到一个新的水平:将列表的子集相加到树中,以避免将单个数组元素添加到更大的值。例如,请参阅How to avoid less precise sum for numpy-arrays with multiple columns

使用固定数量的累加器(例如 8x __m256 = 64 float 累加器)可能会将预期误差减少 64 倍,而不是从 N 到 log N 进行完全成对求和。


脚注 1:关联性是并行化、SIMD 和多个累加器所必需的。 Associativity gives us parallelizability. But what does commutativity give?

在具有例如 4 周期延迟、每时钟 2 次吞吐量 FMA 的机器上,SIMD 宽度为 8 个浮点数,即具有 AVX2 的 Skylake 系统,潜在加速为 4*2 = 8 来自多个累加器,* 8 从 SIMD 宽度乘以内核数,与纯顺序版本相比,即使对于它可能不太准确而不仅仅是不同的问题。

大多数人认为8*8 = 64 值得! (理论上,您还可以在四核上并行化另一个可能为 4 的因子,假设对大型矩阵进行完美缩放)。

您已经在使用 float 而不是 double 来提高性能。

另请参阅Why does mulss take only 3 cycles on Haswell, different from Agner's instruction tables?,了解更多关于使用多个累加器隐藏 FMA 延迟以减少的信息,从而暴露 8 倍加速的其他因素。

另外,不要在最里面的循环中使用hadd。垂直求和,并在循环结束时使用有效的归约。 (Fastest way to do horizontal float vector sum on x86)。你真的很想避免让编译器在每一步都将你的向量提取为标量,这会破坏 SIMD 的大部分好处!除了 hadd 不值得用于 1 个向量的水平总和这一事实之外;它在所有现有 CPU 上花费 2 次随机播放 + 一个常规的 add

【讨论】:

非常感谢您的明确回答!我会改变我的 hadd 实现,因为我不知道它有那种效果......

以上是关于Simd matmul 程序给出不同的数值结果的主要内容,如果未能解决你的问题,请参考以下文章

C# 应用程序分析给出不同的结果

为啥 java.awt.Font.getStringBounds 在不同的机器上给出不同的结果?

解释 Metal 和 SIMD 中的不同类型

渲染到屏幕和 FBO 给出不同的结果

C# 和 Java 中的十六进制到字节数组给出不同的结果

如何在不影响性能的情况下抽象 SIMD 代码以处理不同的数据类型