简单的程序，只产生零，错误？ [复制]答案

【问题标题】：Simple program, that only produce zeros, bug? [duplicate]简单的程序，只产生零，错误？ [复制]
【发布时间】：2016-10-29 00:07:51
【问题描述】：

为什么这个 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

【问题讨论】：

我不记得 any Fortran，但 *1/24 看起来很像整数除法。

标签： fortran fortran90

【解决方案1】：

x=(j+1)*1/24

1/24 是一个向下舍入为 0 的整数除法。您应该能够通过使至少一个操作数为浮点数来强制浮点除法，例如

x=(j+1)*1.0/24.0

【讨论】：

【解决方案2】：

正如 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 实型

1 实数类型具有近似于数学实数的值。处理器应提供两个或多个近似方法，为实数类型的数据定义一组值。每个这样的方法都有一个表示方法，并由种类类型参数 KIND 的值来表征。近似方法的 kind 类型参数由内部函数 KIND (13.7.89) 返回。

5 如果使用类型关键字 REAL 而没有种类类型参数，则指定了具有默认真实种类的真实类型，并且种类值为种类 (0.0)。类型说明符 DOUBLE PRECISION 指定类型 real 双精度类； kind 值为 KIND (0.0D0)。这双精度实数逼近法的小数精度应大于默认的实方法。

这实际上意味着，如果您想确保使用 32 位、64 位或 128 位浮点表示来完成计算，建议您使用内部模块 ISO_FORTRAN_ENV 中定义的正确 KIND 值。

13.8.2.21 REAL32、REAL64 和 REAL128

1 这些默认整数标量命名常量的值应为那些指定 REAL 类型的类型参数以位表示的存储大小分别为 32、64 和 128。如果，对于这些常量中的任何一个，处理器支持不止一种对于那个大小，提供哪种类型的值取决于处理器。如果处理器不支持某种特定大小，则该常量如果处理器支持更大尺寸的种类，则应等于 -2 否则为 -1。

所以这会导致下面的代码

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

【讨论】：