与其他语言相比,Fortran 在数值准确性方面是否存在固有的局限性?
在进行简单的编程练习时,我生成了一个 while 循环(Fortran 中的 DO 循环),该循环的目的是在实际变量达到精确值时退出。
我注意到,由于使用的精度,从未满足相等性并且循环变得无限。当然,这并非闻所未闻,建议不要比较两个数字是否相等,最好看看两个数字之间的绝对差是否小于设定的阈值。
我发现令人失望的是我必须将这个阈值设置得有多低,即使变量是双精度的,我的循环才能正确退出。此外,当我用 Perl 重写这个循环的“精炼”版本时,我在数值精度方面没有任何问题,并且循环退出得很好。
由于在 Perl 和 Fortran 中产生问题的代码都很小,我想在这里重现它,以防我掩盖了一个重要的细节:
Fortran 代码
PROGRAM precision_test
IMPLICIT NONE
! Data Dictionary
INTEGER :: count = 0 ! Number of times the loop has iterated
REAL(KIND=8) :: velocity
REAL(KIND=8), PARAMETER :: MACH_2_METERS_PER_SEC = 340.0
velocity = 0.5 * MACH_2_METERS_PER_SEC ! Initial Velocity
DO
WRITE (*, 300) velocity
300 FORMAT (F20.8)
IF (count == 50) EXIT
IF (velocity == 5.0 * MACH_2_METERS_PER_SEC) EXIT
! IF (abs(velocity - (5.0 * MACH_2_METERS_PER_SEC)) < 1E-4) EXIT
velocity = velocity + 0.1 * MACH_2_METERS_PER_SEC
count = count + 1
END DO
END PROGRAM precision_test
Perl 代码< /strong>
#! /usr/bin/perl -w
use strict;
my $mach_2_meters_per_sec = 340.0;
my $velocity = 0.5 * $mach_2_meters_per_sec;
while (1) {
printf "%20.8f\n", $velocity;
exit if ($velocity == 5.0 * $mach_2_meters_per_sec);
$velocity = $velocity + 0.1 * $mach_2_meters_per_sec;
}
Fortran 中的注释行是我需要使用的循环正常退出的行。请注意,阈值设置为 1E-4,我觉得这很可悲。
变量的名称来自我正在执行的基于自学的编程练习,没有任何相关性。
目的是当速度变量达到 1700 时循环停止。
以下是截断的输出:
Perl 输出
170.00000000
204.00000000
238.00000000
272.00000000
306.00000000
340.00000000
...
1564.00000000
1598.00000000
1632.00000000
1666.00000000
1700.00000000
Fortran 输出
170.00000000
204.00000051
238.00000101
272.00000152
306.00000203
340.00000253
...
1564.00002077
1598.00002128
1632.00002179
1666.00002229
1700.00002280
Fortran 的速度和速度有什么好处?如果并行化的准确性很糟糕,那么并行化是否容易?让我想起做事的三种方法:
正确的方法
错误的方法
错误的最大功率方法
” “这不是方法不对吗
“是啊!不过要快一些!”
开个玩笑,我一定做错了什么。
与其他语言相比,Fortran 在数值准确性方面是否存在固有的局限性,或者我(很可能)是有问题的人?
我的编译器是 gfortran(gcc 版本 4.1.2)、Perl v5.12.1,在双核 AMD Opteron @ 1 GHZ 上。
While working on a simple programming exercise, I produced a while loop (DO loop in Fortran) that was meant to exit when a real variable had reached a precise value.
I noticed that due to the precision being used, the equality was never met and the loop became infinite. This is, of course, not unheard of and one is advised that, rather than comparing two numbers for equality, it is best see if the absolute difference between two numbers is less than a set threshold.
What I found disappointing was how low I had to set this threshold, even with variables at double precision, for my loop to exit properly. Furthermore, when I rewrote a "distilled" version of this loop in Perl, I had no problems with numerical accuracy and the loop exited fine.
Since the code to produce the problem is so small, in both Perl and Fortran, I'd like to reproduce it here in case I am glossing over an important detail:
Fortran Code
PROGRAM precision_test
IMPLICIT NONE
! Data Dictionary
INTEGER :: count = 0 ! Number of times the loop has iterated
REAL(KIND=8) :: velocity
REAL(KIND=8), PARAMETER :: MACH_2_METERS_PER_SEC = 340.0
velocity = 0.5 * MACH_2_METERS_PER_SEC ! Initial Velocity
DO
WRITE (*, 300) velocity
300 FORMAT (F20.8)
IF (count == 50) EXIT
IF (velocity == 5.0 * MACH_2_METERS_PER_SEC) EXIT
! IF (abs(velocity - (5.0 * MACH_2_METERS_PER_SEC)) < 1E-4) EXIT
velocity = velocity + 0.1 * MACH_2_METERS_PER_SEC
count = count + 1
END DO
END PROGRAM precision_test
Perl Code
#! /usr/bin/perl -w
use strict;
my $mach_2_meters_per_sec = 340.0;
my $velocity = 0.5 * $mach_2_meters_per_sec;
while (1) {
printf "%20.8f\n", $velocity;
exit if ($velocity == 5.0 * $mach_2_meters_per_sec);
$velocity = $velocity + 0.1 * $mach_2_meters_per_sec;
}
The commented-out line in Fortran is what I would need to use for the loop to exit normally. Notice that the threshold is set to 1E-4, which I feel is quite pathetic.
The names of the variables come from the self-study-based programming exercise I was performing and don't have any relevance.
The intent is that the loop stops when the velocity variable reaches 1700.
Here are the truncated outputs:
Perl Output
170.00000000
204.00000000
238.00000000
272.00000000
306.00000000
340.00000000
...
1564.00000000
1598.00000000
1632.00000000
1666.00000000
1700.00000000
Fortran Output
170.00000000
204.00000051
238.00000101
272.00000152
306.00000203
340.00000253
...
1564.00002077
1598.00002128
1632.00002179
1666.00002229
1700.00002280
What good is Fortran's speed and ease of parallelization if its accuracy stinks? Reminds me of the three ways to do things:
The Right Way
The Wrong Way
The Max Power Way
"Isn't that just the wrong way?"
"Yeah! But faster!"
All kidding aside, I must be doing something wrong.
Does Fortran have inherent limitations on numerical accuracy compared to other languages, or am I (quite likely) the one at fault?
My compiler is gfortran (gcc version 4.1.2), Perl v5.12.1, on a Dual Core AMD Opteron @ 1 GHZ.
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(2)
您的分配意外地将值转换为单精度,然后再转换回双精度。
尝试将您的
0.1 *
更改为0.1D0 *
,您应该会看到问题得到解决。Your assignment is accidentally converting the value to single precision and then back to double.
Try making your
0.1 *
be0.1D0 *
and you should see your problem fixed.正如已经回答的那样,Fortran 中的“普通”浮点常量将默认为默认实数类型,这可能是单精度的。这几乎是一个经典的错误。
另外,使用“kind=8”是不可移植的——它会给你gfortran的双精度,但不会给其他一些编译器。在 Fortran >= 90 中指定变量和常量精度的安全、可移植的方法是使用内部函数,并请求所需的精度。然后在精度很重要的常量上指定“种类”。一种方便的方法是定义您自己的符号。例如:
这对于整数也很重要,例如,如果需要大值。有关整数的相关问题,请参阅 Fortran:整数*4 与整数(4) 与整数( kind=4) 和 Fortran 中的长整型
As already answered, "plain" floating point constants in Fortran will default to the default real type, which will likely be single-precision. This is an almost classic mistake.
Also, using "kind=8" is not portable -- it will give you double precision with gfortran, but not with some other compilers. The safe, portable way to specify precisions for both variables and constants in Fortran >= 90 is to use the intrinsic functions, and request the precision that you need. Then specify "kinds" on the constants where precision is important. A convenient method is to define your own symbols. For example:
This can also be important for integers, e.g., if large values are needed. For related questions for integers, see Fortran: integer*4 vs integer(4) vs integer(kind=4) and Long ints in Fortran