代码之家 › 专栏 › 技术社区 › Gogrik

Fortran中双精度除以整数

fortran90 division precision fortran

Gogrik · 技术社区 · 9 年前

现在我正在处理一位在我之前参与该项目的人的Fortran代码。代码相对庞大,所以我不会在这里提供。在这段代码中的很多地方都有双精度值除以整数的方法。即。:

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

这有时会导致比常规浮点数错误更严重的错误。我将为你复制其中一个错误:3.00000032899483。对他的目的来说,这是可以的,但对我来说,这个错误是可怕的。看起来,fortran正在将实数除以整数,然后将数字放大到两倍精度,添加一些随机内容。当然,如果我写:

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

它成功了!非常感谢您的帮助,如果没有您的回答,我将处理更多的问题。

2 回复 | 直到 8 年前

Gilles MrUser 9 年前

我相信你问题的最初断言,即(如果我很理解你的话) real 由 integer 导致精度低于除以 真实的 ,是错误的。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

agentp 9 年前

如前所述,您的标题声明不正确,整数除法会自动向上转换,并且是正确的。

然而,您的第一个示例确实存在问题:

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

这里有一个 单精度 1. 这个 3 提升为单精度,因此您可以使用单精度近似值 1/3 分配给您的双精度 a 。解决方法是 1.0 1.d0 .

在我看来,绕过重载运算符是个坏主意。如果需要手动修复的代码太多,请查看编译器是否有将默认文字类型更改为double的选项。