以下是解决@High Performance Mark 问题的工作示例的结果。
host system = (redacted)
compiler version = GCC version 5.1.0
compiler options = -fPIC -mmacosx-version-min=10.9.4 -mtune=core2 -Og -Wall -Wextra -Wconversion -Wpedantic -fmax-errors=5
execution command = ./a.out
Compare mesh points
1.57017982 1.57080817 1.57143652
1.57016802 1.57079625 1.57142460
Compare function values at these mesh points
1622.04211 -84420.7344 -1562.01758
1591.57471 13245402.0 -1591.65527
在您发现的演示中有一个不好的编程实践:通过添加步骤(+ 增量)而不是计算它们(k * 增量)来移动网格。这个问题很普遍,而且后果很严重(https://www.ima.umn.edu/~arnold/disasters/patriot.html)。
为了演示您的代码,网格的大小被提升到 10K 点。样本函数也从cos x
变为tan x
,我们在 pi/2 = 1.57079633 处检查奇点附近的点。虽然新手可能会发现网格值的差异微不足道,但函数值的差异是显着的。
(可以通过使用具有精确二进制表示的增量来减少网格错误,例如 2^(-13) = 1 / 8192。)
代码显示在这里。编译命令是gfortran -Wall -Wextra -Wconversion -Og -pedantic -fmax-errors=5 demo.f95
. 运行命令是./a.out
。
program demo
use iso_fortran_env
implicit none
real, parameter :: pi = acos ( -1.0 )
integer, parameter :: n = 10001
real, dimension ( 1 : n ) :: x, y, z
real :: a = 0.0, b = 2 * pi
real :: increment
integer :: k, quarter, status
character ( len = * ), parameter :: c_options = compiler_options( )
character ( len = * ), parameter :: c_version = compiler_version( )
character ( len = 255 ) :: host = " ", cmd = " "
! queries
call hostnm ( host, status )
call get_command ( cmd )
! write identifiers
write ( *, '( /, "host system = ", g0 )' ) trim ( host )
write ( *, '( "compiler version = ", g0 )' ) c_version
write ( *, '( "compiler options = ", g0 )' ) trim ( c_options )
write ( *, '( "execution command = ", g0, / )' ) trim ( cmd )
increment = ( b - a ) / ( n - 1 )
quarter = n / 4
! mesh accumulates errors
x ( 1 ) = 0.0
do k = 2, n
x ( k ) = x ( k - 1 ) + increment
end do
y = tan ( x )
print *, 'Compare mesh points'
print *, x ( quarter : quarter + 2 )
! better mesh
x ( 1 ) = 0.0
do k = 2, n
x ( k ) = ( k - 1 ) * increment
end do
z = tan ( x )
print *, x ( quarter : quarter + 2 )
print *, 'Compare function values at these mesh points'
print *, y ( quarter : quarter + 2 )
print *, z ( quarter : quarter + 2 )
end program demo