fortran - 在 Fortran 中用整数除以 double

Question

现在我正在处理在我之前从事该项目的人的 Fortran 代码。代码比较庞大，这里就不提供了。在这段代码的很多地方， double 值除以整数。即:

double precision :: a,c
integer :: b
a = 1.0d0
b = 3
c = a/b

这有时会导致比常规 float 错误更严重的错误。我将为您复制粘贴其中一个错误:3.00000032899483。就他的目的而言，这是可以的，但对我来说，这个错误是可怕的。似乎 Fortran 将实数除以整数，然后将数字扩大到 double ，添加一些随机的东西。如果我写:

c = a/(real(b,8))

所以，我可以做到，但他的所有代码都需要很长时间才能搜索到。这就是为什么，我尝试以下一种方式重载 (/) 运算符:

MODULE divisionoverload
!-----------------------------------------------------
   INTERFACE OPERATOR (/)
      MODULE PROCEDURE real_divoverinteger
   END INTERFACE OPERATOR ( / )
CONTAINS
!-----------------------------------------------------

DOUBLE PRECISION FUNCTION real_divoverinteger(a,b)
   IMPLICIT NONE
   DOUBLE PRECISION, INTENT (IN) :: a 
   INTEGER, INTENT(IN) :: b 

   real_divoverinteger = a/real(b,8)
END FUNCTION real_divoverinteger

END MODULE divisionoverload

但是 gfortran4.8 给我的错误清楚地表明这与内部除法运算符冲突:

 MODULE PROCEDURE real_divoverinteger
                                          1
Error: Operator interface at (1) conflicts with intrinsic interface

那么，如果不手动处理所有整数除法，我该怎么做才能解决这个问题呢？

生成打印格式数字的代码:

 hx = (gridx(Nx1+1)-gridx(Nx1))/modx)
 do i = 2,lenx
    sub%subgridx(i)=sub%subgridx(i-1) + hx
 end do
 !*WHERE*
 !hx,gridx[1d array], sub%subgridx[1d array] - double precision
 !modx, lenx - integer

打印代码

subroutine PrintGrid(grid,N)
    integer, intent (in) :: N       
    double precision, dimension(N), intent (in) :: grid
    integer :: i
    print"(15(f20.14))", (grid(i), i=1,N)
        print*, ""
end subroutine PrintGrid

---

解决方案:我为这个问题道歉。我应该更仔细地阅读代码。我发现了一个错误。下一个问题的重现:

program main
   double precision :: a
   a = 1.0/3
   print*, a
end program main

基本上，他是在进行实数/整数除法并将这个值赋给实数*8。我用这样的代码改变了一切:

 program main
   double precision :: a
   a = 1.0d0/3
   print*, a
end program main

成功了!非常感谢您的帮助，如果没有您的回答，我会更多地处理这个问题。

Answer 1

我相信你的问题的最初断言，即(如果我理解你的话)将 real 除以 integer 导致精度低于除以 a 真实的，是错误的。在计算操作之前，Fortran 将操作的两个操作数中最弱的类型(如果需要)提升为最强类型。因此，自己做宣传是没有必要的。

例如看下面的代码:

program promotion
    implicit none
    integer, parameter :: dp = kind( 1.d0 )
    integer :: i
    real( kind = dp ) :: a

    i = 3
    a = 1
    print *, "division by a real:", a/real(i,dp), "division by an integer:", a/i
    print *, "are they equal?", a/real(i,dp) == a/i
end program promotion

结果的实际打印可能因编译器而异，但最终测试的计算结果始终为 .true。

~/tmp$ gfortran promotion.f90
~/tmp$ ./a.out 
 division by a real:  0.33333333333333331      division by an integer:  0.33333333333333331     
 are they equal? T