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

Question

我正在尝试使用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 

如图所示，这两个版本之间存在数值误差。任何帮助将非常感激 ！

Answer 1

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

FP数学不是关联的；好像是进行优化都会改变结果。¹ 浮点加法和乘法关联吗？ / C中的浮点运算是否关联？

更改量取决于数据。对于，仅在小数点后第五位的差异似乎是合理的float。

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

实际上，假设所有数字的大小相似，使用多个累加器通常可以提高大型列表的精度。（理想情况下，每个SIMD向量都由多个元素组成，以隐藏FP-add或FMA延迟）。 https://en.wikipedia.org/wiki/Pairwise_summation是一种数字技术，可将其带入一个新的水平：将树中列表的子集求和，以避免将单个数组元素添加到更大的值。例如，请参见如何避免多列numpy数组的总和不那么精确

使用固定数量的累加器（例如8x __m256= 64个float累加器）可能会将预期误差降低64 倍，而不是从N到对数N进行完整的成对求和。

脚注1：对于并行化，SIMD和多个累加器，关联性是必需的。 关联性使我们具有可并行性。但是可交换性给什么呢？

在具有4个周期的延迟，每时钟2个吞吐量的FMA的计算机上，SIMD宽度为8个浮点数，即使用AVX2的Skylake系统，多个累加器的潜在加速为4 * 2 = 8，SIMD为* 8宽度，内核数乘以纯顺序版本，即使对于精度可能较低而并非完全不同的问题。

大多数人认为8*8 = 64值得这样做！（并且理论上，您也可以在四核上并行化另一个可能为4的因数，假设对大型矩阵进行完美缩放）。

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

另请参阅为什么mulss在Haswell上仅需要3个周期，而与Agner的指令表不同？详细了解如何使用多个累加器来隐藏FMA延迟，从而减少了其他8倍加速。

另外，千万不能使用hadd的最内循环中。垂直求和，并在循环结束时使用有效的归约。（在x86上执行水平浮点向量总和的最快方法）。您真的很想避免编译器在每一步都将您的向量提取为标量，这使SIMD的大部分优势都无法兑现！除了hadd不值得用于1个向量的水平和的事实之外，add在所有现有的CPU上，它需要2次改组+定期执行。