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

Question

我正在尝试使用 simd intrinsics 在 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 数学不是关联的；像这样优化会改变结果。¹ Is Floating point addition and multiplication associative? / Are floating point operations in C associative?

变化量取决于数据。对于 float.

，仅小数点后 5 位的差异似乎是合理的

除非您采取特殊的数值预防措施，例如先将小数相加，否则顺序结果不是 "more correct"，它们只是有不同的错误。

事实上，使用多个累加器通常提高大型列表的精度，假设您的数字都具有相似的大小。（理想情况下，多个 SIMD 向量每个都由多个元素组成，以隐藏 FP-add 或 FMA 延迟）。 https://en.wikipedia.org/wiki/Pairwise_summation is a numerical technique that takes this to the next level: summing subsets of the list in a tree, to avoid adding single array elements to a much larger value. See for example,

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

脚注 1：关联性对于并行化、SIMD 和多个累加器是必需的。

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

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

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

另请参阅了解更多关于使用多个累加器隐藏 FMA 延迟的减少，暴露 8 倍加速的其他因素。

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

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

Simd matmul program gives different numerical results

c

floating-point

simd

vectorization

avx