5

首先,抱歉标题不好:/

我试图重现一篇论文关于计算三对角对称矩阵的特征值的结果。我分别通过四舍五入到正无穷和负无穷来确定一些值“上限和下限”。

我不是每次都更改舍入模式,而是使用“技巧”:fl⁻(y) = -fl⁺(-y),其中 fl⁻(y) 是使用负无穷大舍入模式时 y 的值,并且fl⁺(y) 是使用舍入模式到正无穷大时的 y 值。因此,我在 C 中有以下代码:

fesetround(FE_UPWARD);
first =  - (-d[i] + x);
second =  ( - (( e[i-1]*e[i-1] ) / a_inf ));
a_inf = first + second;

first = d[i] - x;
second = - ( ( e[i-1]*e[i-1] ) / a_sup );
a_sup = first + second;

它工作正常,除了一个例子 a_inf 给了我正确的结果,但 a_sup 给出了错误的结果,尽管第一个和第二个变量似乎具有相同的值。

但是,如果我这样做:

  fesetround(FE_UPWARD);
  first =  - (-d[i] + x);
  second =  ( - (( e[i-1]*e[i-1] ) / a_inf ));

  fesetround(FE_DOWNWARD);
  first =  - (-d[i] + x);
  second =  ( - (( e[i-1]*e[i-1] ) / a_sup ));

我得到了正确的结果。因此,如果我使用技巧 fl⁻(y) = -fl⁺(-y),我会得到正确的结果,如果我更改舍入模式并使用原始表达式,我会得到错误的结果。知道为什么吗?

在这两种情况下,变量 first 和 second 的值如下:

first  1.031250000000000e+07,  second -1.031250000000000e+07
first  1.031250000000000e+07,  second -1.031250000000000e+07

a_inf 和 a_sup 的正确值分别是 -1.862645149230957e-09 和 +1.862645149230957e-09,但在第一种情况下 a_sup = 0,这是错误的

我猜它正在发生的是某种灾难性的取消,但我不知道在这种情况下如何解决它......

提前致谢!

4

1 回答 1

1

您遇到的第一个问题是您仅使用“技巧”来first向 -inf 舍入,而不是second,因此它仍然向 +inf 舍入。

第二个问题是 C 没有给你任何关于它如何进行浮点计算的保证。特别是,编译器可以自由地重新排序和重新关联操作,因为它认为适合性能,即使这种重新排列可能会改变程序的舍入行为。所以当你说:

first =  - (-d[i] + x);

编译器可能会重新排列它并执行一次减法,而不是求反、加法、求反,这会反转舍入的方向。现在,您有时可以通过禁用所有优化来使事情按您期望的方式工作,但即使这样,也不能保证。

于 2012-06-17T01:43:07.240 回答