简单的程序,只产生零,错误?
Simple program, that only produce zeros, bug?
为什么这个 Fortran 程序只产生零?当我打印出来时,到处都是 -0.00000!我做错了什么?在 matlab 中它运行完美。老实说,我看不出它为什么不起作用!
似乎是分数搞砸了。如果我将 x 设置为某个小数,它就可以工作。
program main
implicit none
integer iMax, jMax
double precision, dimension(:,:), allocatable :: T
double precision x, dx,f,L2old,L2norm,y
integer i, j,n,bc
n=10
allocate(T(1:n+2, 1:n+2))
T=0.0d0
do i=2,n+1
do j=2,n+1
x=(j+1)*1/24
y=(i+1)*1/24
T(i,j)= -18*(x**2+y**2)**2
Write(*,*)'T(',i,'',j,'', T(i,j)
end do
end do
Write(*,*)'T(1,1)',T(1,1)
end program main
x=(j+1)*1/24
1/24
是一个向下舍入为 0 的整数除法。您应该能够通过使至少一个操作数浮点来强制进行浮点除法,
例如
x=(j+1)*1.0/24.0
正如 Jim Lewis 所指出的,OP 问题的答案确实是使用的整数除法。
尽管如此,我认为重要的是要指出人们应该注意浮点数是如何写下来的。正如 OP 的程序所示,x
是 DOUBLE PRECISION
类型。那么正确的结果应该是
x=(j+1)*1.0D0/24.0D0
这里的区别在于,现在您要确保除法的精度与声明的 x
相同。
以下程序演示问题::
program test
WRITE(*,'(A43)') "0.0416666666666666666666666666666666..."
WRITE(*,'(F40.34)') 1/24
WRITE(*,'(F40.34)') 1.0/24.0
WRITE(*,'(F40.34)') 1.0D0/24.0
WRITE(*,'(F40.34)') 1.0D0/24.0D0
end program test
作为输出
0.0416666666666666666666666666666666...
0.0000000000000000000000000000000000
0.0416666679084300994873046875000000
0.0416666666666666643537020320309239
0.0416666666666666643537020320309239
你清楚地看到了差异。第一行是正确的数学结果。第二行是导致零的整数除法。第三行显示除法计算为 REAL
的情况下的输出,而第四行和第五行在 DOUBLE PRECISION
中。请注意,在我的例子中,REAL
表示 32 位浮点数,DOUBLE PRECISION
表示 64 位版本。 REAL
和 DOUBLE PRECISION
的精度和表示形式取决于编译器,未在 the Standard 中定义。只要求DOUBLE PRECISION
的精度比REAL
.
高
4.4.2.3 Real type
1 The real type has values that approximate the mathematical real numbers. The processor shall provide two or more approximation methods that define sets of values for data of type real. Each such method has a representation method and is characterized by a value for the kind type parameter KIND. The kind type parameter of an approximation method is returned by the intrinsic function KIND (13.7.89).
5 If the type keyword REAL is used without a kind type parameter, the
real type with default real kind is specified and the kind value is
KIND (0.0). The type specifier DOUBLE PRECISION specifies type real
with double precision kind; the kind value is KIND (0.0D0). The
decimal precision of the double precision real approximation method
shall be greater than that of the default real method.
这实际上意味着,如果您想确保您的计算是使用 32 位、64 位或 128 位浮点表示法完成的,建议您使用内部模块中定义的正确 KIND
值 ISO_FORTRAN_ENV
.
13.8.2.21 REAL32, REAL64, and REAL128
1 The values of these default integer scalar named constants shall be
those of the kind type parameters that specify a REAL type whose
storage size expressed in bits is 32, 64, and 128 respectively. If,
for any of these constants, the processor supports more than one kind
of that size, it is processor dependent which kind value is provided.
If the processor supports no kind of a particular size, that constant
shall be equal to −2 if the processor supports kinds of a larger size
and −1 otherwise.
所以这会导致下面的代码
PROGRAM main
USE iso_fortran_env, ONLY : DP => REAL64
IMPLICIT NONE
...
REAL(DP) :: x
...
x = (j+1)*1.0_DP/24.0_DP
...
END PROGRAM main
为什么这个 Fortran 程序只产生零?当我打印出来时,到处都是 -0.00000!我做错了什么?在 matlab 中它运行完美。老实说,我看不出它为什么不起作用!
似乎是分数搞砸了。如果我将 x 设置为某个小数,它就可以工作。
program main
implicit none
integer iMax, jMax
double precision, dimension(:,:), allocatable :: T
double precision x, dx,f,L2old,L2norm,y
integer i, j,n,bc
n=10
allocate(T(1:n+2, 1:n+2))
T=0.0d0
do i=2,n+1
do j=2,n+1
x=(j+1)*1/24
y=(i+1)*1/24
T(i,j)= -18*(x**2+y**2)**2
Write(*,*)'T(',i,'',j,'', T(i,j)
end do
end do
Write(*,*)'T(1,1)',T(1,1)
end program main
x=(j+1)*1/24
1/24
是一个向下舍入为 0 的整数除法。您应该能够通过使至少一个操作数浮点来强制进行浮点除法,
例如
x=(j+1)*1.0/24.0
正如 Jim Lewis 所指出的,OP 问题的答案确实是使用的整数除法。
尽管如此,我认为重要的是要指出人们应该注意浮点数是如何写下来的。正如 OP 的程序所示,x
是 DOUBLE PRECISION
类型。那么正确的结果应该是
x=(j+1)*1.0D0/24.0D0
这里的区别在于,现在您要确保除法的精度与声明的 x
相同。
以下程序演示问题::
program test
WRITE(*,'(A43)') "0.0416666666666666666666666666666666..."
WRITE(*,'(F40.34)') 1/24
WRITE(*,'(F40.34)') 1.0/24.0
WRITE(*,'(F40.34)') 1.0D0/24.0
WRITE(*,'(F40.34)') 1.0D0/24.0D0
end program test
作为输出
0.0416666666666666666666666666666666...
0.0000000000000000000000000000000000
0.0416666679084300994873046875000000
0.0416666666666666643537020320309239
0.0416666666666666643537020320309239
你清楚地看到了差异。第一行是正确的数学结果。第二行是导致零的整数除法。第三行显示除法计算为 REAL
的情况下的输出,而第四行和第五行在 DOUBLE PRECISION
中。请注意,在我的例子中,REAL
表示 32 位浮点数,DOUBLE PRECISION
表示 64 位版本。 REAL
和 DOUBLE PRECISION
的精度和表示形式取决于编译器,未在 the Standard 中定义。只要求DOUBLE PRECISION
的精度比REAL
.
4.4.2.3 Real type
1 The real type has values that approximate the mathematical real numbers. The processor shall provide two or more approximation methods that define sets of values for data of type real. Each such method has a representation method and is characterized by a value for the kind type parameter KIND. The kind type parameter of an approximation method is returned by the intrinsic function KIND (13.7.89).
5 If the type keyword REAL is used without a kind type parameter, the real type with default real kind is specified and the kind value is KIND (0.0). The type specifier DOUBLE PRECISION specifies type real with double precision kind; the kind value is KIND (0.0D0). The decimal precision of the double precision real approximation method shall be greater than that of the default real method.
这实际上意味着,如果您想确保您的计算是使用 32 位、64 位或 128 位浮点表示法完成的,建议您使用内部模块中定义的正确 KIND
值 ISO_FORTRAN_ENV
.
13.8.2.21 REAL32, REAL64, and REAL128
1 The values of these default integer scalar named constants shall be those of the kind type parameters that specify a REAL type whose storage size expressed in bits is 32, 64, and 128 respectively. If, for any of these constants, the processor supports more than one kind of that size, it is processor dependent which kind value is provided. If the processor supports no kind of a particular size, that constant shall be equal to −2 if the processor supports kinds of a larger size and −1 otherwise.
所以这会导致下面的代码
PROGRAM main
USE iso_fortran_env, ONLY : DP => REAL64
IMPLICIT NONE
...
REAL(DP) :: x
...
x = (j+1)*1.0_DP/24.0_DP
...
END PROGRAM main