十行代码块中高达 20% 的数字错误或错误

Question

重写应用程序的一小段代码已经产生了相当大的性能改进。代码是 100% 顺序的，因此存储在内存中的值应该没有隐藏的扰动。

仔细检查计算后的结果是否相同，结果显示相对误差高达 20%，所以问题是这是否可以用算法中的数值错误来解释，或者存在错误在算法本身和两个块是不等价的。

从概念上讲，主要的修改是替换这个

for m=0,...
  result += gradient * temp[m]

有了这个

for m=0,...
  sum_temp += temp[m]
result = gradient * sum_temp

同时在处理定义新数组和初始化数组的过程中

编辑：（为清楚起见）下面代码中所谓的 'result' 是 SoGrSca 和 SoGrSca 的元素

编辑：（最终）结果的典型值为 +-1e-8

实际代码如下（C++，怪异命名的嵌套数组为了清晰起见一直保持不变）

// Three relevant macros for tidyness
//
#define FOR_I3 for (int i = 0; i < 3; ++i)
#define LOOP_n_N(NUM)  for (int n = 0; n < NUM; n++)
#define LOOP_m_N(NUM)  for (int m = 0; m < NUM; m++)

// Toy definitions for eye-friendlyness
// on Whosebug
//
const int NDIM = 6; 
const int ORDER = 20;

// More definitions...
//
// 1) Arrays of doubles
//
// SoGrSca, SoGrVec, PrimSca, Jaco_Sol, PrimVec, mat_der_sol
//
// 2) Arrays of integers
//
// mat_co1, mat_co2



// Version 1 of the code
//
FOR_I3 SoGrSca[0][nk + l][i]           = 0.0;
FOR_I3 FOR_J3 SoGrVec[0][nk + l][j][i] = 0.0;
//
LOOP_n_N(NDIM) LOOP_m_N(ORDER) {

  int idof1 = ORDER*mat_co1[NDIM*l + n] + m;
  int idof2 = mat_Gr[ndof*n +   ORDER*mat_co2[NDIM*l + n] + m           ];

  double grp_temp = mat_der_sol[idof1]*GsGm1*PrimSca[1][nk + idof2]/PrimSca[0][nk + idof2];
  FOR_I3 SoGrSca[0][nk + l][i]      += Jaco_Sol[nk + l][n][i]*grp_temp;
  FOR_J3 {
    double grp_temp = mat_der_sol[idof1]*PrimVec[0][nk + idof2][j];
    FOR_I3 SoGrVec[0][nk + l][j][i] += Jaco_Sol[nk + l][n][i]*grp_temp;
  }
}

// Store SoGrSca and SoGrVec values ...

// Version 2 of the code
//
FOR_I3 SoGrSca[0][nk + l][i]           = 0.0;
FOR_I3 FOR_J3 SoGrVec[0][nk + l][j][i] = 0.0;
//
double grp_temp_A_v[NDIM];
double grp_temp_B_v[NDIM][3];
LOOP_n_N(NDIM){                    
        grp_temp_A_v[n] = 0.0;                  
        FOR_I3  grp_temp_B_v[n][i] = 0.0;       
}
LOOP_n_N(NDIM) LOOP_m_N(ORDER) {
  int idof1 = ORDER * mat_co1[NDIM*l + n] + m;
  int idof2 = mat_Gr[ndof*n +   ORDER*mat_co2[NDIM*l + n] + m           ];
  grp_temp_A_v[n] += mat_der_sol[idof1] * GsGm1*PrimSca[1][nk + idof2]/PrimSca[0][nk + idof2];
  FOR_J3 grp_temp_B_v[n][j] += mat_der_sol[idof1] * PrimVec[0][nk + idof2][j];
}

LOOP_n_N(NDIM) {
  FOR_I3 SoGrSca[0][nk + l][i]      += Jaco_Sol[nk + l][n][i] * grp_temp_A_v[n];
  FOR_J3 {
    FOR_I3 SoGrVec[0][nk + l][j][i] += Jaco_Sol[nk + l][n][i] * grp_temp_B_v[n][j];
  }
}


// Compare SoGrSca and SoGrVec... 
//
// ERROR: they are different

比较算法是，给定y_old和y_new，

( abs( y_old - y_new ) / ( abs(y_old) + z ) ) < epsilon // 'true' means equality

其中 z 防止零除法错误，

z=1e-10

并且通过所需的 epsilon 很大，即 epsilon=0.2 虽然它应该更接近例如1e-7

有趣的编辑：这只是加入带有已经完善的宏的项目的另一个例子哈哈

Answer 1

浮点数运算不是distributive。这意味着在某些情况下以下陈述是正确的

x * (a + b) != x * a + x * b;

如果您将乘法移出循环并期望得到完全相同的结果，则您假设分配律始终为真，但事实并非如此。例如，我想出了 this godbolt snippet.

十行代码块中高达 20% 的数字错误或错误

up to 20% Numerical error or Bug in ten line code block

c

c++

rounding-error

numerical-methods