Fortran 程序没有正确更新变量的值答案

【问题标题】：Fortran program doesn't update a variable's value properlyFortran 程序没有正确更新变量的值
【发布时间】：2021-10-26 05:46:18
【问题描述】：

我正在尝试在 Fortran 中实现 Bisection 方法子例程来解决计算科学程序，而 Fortran 正在做一些奇怪的事情。因此，程序的目标是找到某个参数e0 的超越方程的解，该参数在 for 循环中的每一步都更新并传递给子程序。

问题是：

e0 未按应有的方式更新。 e0 的值在一次迭代中从 0.1799... 变为 0.8999...，而理论上它应该需要程序 9 次迭代才能到达那里。为什么会这样？
打印语句未按应有的方式打印（请参见下面的输出）。我们期望“out 1”打印，一些“in”打印（当子程序 f1 被子程序 Bisection 调用时打印），然后是“out 2”打印（使用新的 e0 值），一些“in”打印等。但我们看到第一个“out 1”打印，一些“in”打印，然后只有“out”打印。这是否意味着子程序f1 仅在“out 1”和“out 2”之间调用？（应该在每次“输出”打印之间调用）

我使用 Fortran77 进行数值求解已经有几年了，从来没有遇到过这样的事情，但是我已经快 6 个月没有编程了，所以我可能错过了一件重要的事情。

代码：

       program roots
       implicit none
       double precision A,B,eps,e0,e1
       integer i,niter
       external f1
       A = 1d-1
       B = 9d-1
       eps = 10d-6
       open(9,file='dades.dat',status='old')
       
       do i=1,9
       e0 = i*(B-A)/10 + A
c      Here is the first print. 'Out' meaning outside the subroutine
       print *, 'out', i, e0
       call Bisection(A,B,eps,f1,niter,e1,e0)
       write (9,'(2(f10.5))') e0,e1
       end do
       
       close(9)
       end program roots
       
       subroutine Bisection(A,B,eps,f,niter,xroot,e0)
       implicit none
              double precision A,B,eps,xroot,fuc,fua,e0
              integer niter,i
              niter = nint(log((B-A)/eps)/log(2.))+1
              
              do i=1,niter
                     xroot = (A+B)/2
c      Here the subroutine which uses e0 is called twice
                     call f(xroot,fuc,e0)
                     call f(A,fua,e0)
                     if (fuc .eq. 0) return
                     if (fuc*fua .lt. 0.) then
                            B = xroot
                     else
                            A = xroot
                     end if
              end do
              return
       end subroutine Bisection
       
       subroutine f1(x,f,e0)
       implicit none
       double precision x,f,e0
c      Here is the second print. 'In' meaning inside the subroutine
       print *, 'in', e0
       f = e0+1
       end subroutine f1

输出：

out           1  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 in  0.17999999999999999
 out           2  0.89999511718750003
 out           3  0.89999572753906254
 out           4  0.89999633789062505
 out           5  0.89999694824218746
 out           6  0.89999755859374997
 out           7  0.89999816894531248
 out           8  0.89999877929687500
 out           9  0.89999938964843751

【问题讨论】：

我删除了fortran77 标签。代码至少为 Fortran 90。
“这在这个例子中看不到” - 如果您要查看这个，请创建一个确实显示问题的例子。
新的 F90+ 代码中任何常见块的问题的解决方案是：不要使用常见块。没有令人信服的理由使用空白通用来解决这个编程问题，这可能会减少您在这里帮助您的人可能会看到的工作量。请考虑使用“现代”编程结构。（并且调用子程序f 和fun 等会让很多人感到困惑
在进入子程序之前，你只需要在out do-loop 中设置e0。 e0 永远不会在子例程中更新。
@JordiFuentes 在过去的 25 年中不应该在新代码中使用容易出错的公共块。如果您确实必须拥有全局数据（这不是一个好主意），请使用可以执行与 common 类似的事情但不易出错的模块。最好的方法是通过参数列表，就像您在更新中所做的那样，但是您仍然应该真正了解模块并将所有子程序合二为一——这是现代的、不易出错的方式；你现在正在做的是在没有安全带或刹车的情况下以 100 公里/小时的速度行驶。

标签： fortran

【解决方案1】：

在深入查看代码后，我终于找到了问题所在。问题是变量A 和B 既是输入又是输出，因为它们是在子例程中修改的。因此，当第一次调用子程序时，它会更改A 和B 的值（与代码开头设置的不同）。其余的都是由此引发的问题。
解决方案是在每次迭代开始时重置值，如下所示：

       program roots
       implicit none
       double precision A,B,eps,e0,e1
       integer i,niter
       external f1
       eps = 10d-6
       open(9,file='dades.dat',status='old')
       
       do i=1,9
       A = 1d-1
       B = 9d-1
       e0 = i*(B-A)/10 + A
c      Here is the first print. 'Out' meaning outside the subroutine
       print *, 'out', i, e0
       call Bisection(A,B,eps,f1,niter,e1,e0)
       write (9,'(2(f10.5))') e0,e1
       end do
       
       close(9)
       end program roots
       
       subroutine Bisection(A,B,eps,f,niter,xroot,e0)
       implicit none
              double precision A,B,eps,xroot,fuc,fua,e0
              integer niter,i
              niter = nint(log((B-A)/eps)/log(2.))+1
              
              do i=1,niter
                     xroot = (A+B)/2
c      Here the subroutine which uses e0 is called twice
                     call f(xroot,fuc,e0)
                     call f(A,fua,e0)
                     if (fuc .eq. 0) return
                     if (fuc*fua .lt. 0.) then
                            B = xroot
                     else
                            A = xroot
                     end if
              end do
              return
       end subroutine Bisection
       
       subroutine f1(x,f,e0)
       implicit none
       double precision x,f,e0
c      Here is the second print. 'In' meaning inside the subroutine
       print *, 'in', e0
       f = e0+1
       end subroutine f1

【讨论】：

如果你真的想保留A 和B 的值，为什么不让Bisection 自己通过使用它们的本地副本来处理呢？此外，我会在Bisection 中添加f 的接口，从而避免需要external。代码中还有几件事使它更像是 Fortran 77，而不是其他任何东西。这有什么特殊原因吗？
@Pap - 制作本地副本的要点。这是一个我没有想到的简单概念。当您说“在Bisection 中添加f 的接口”时，您到底是什么意思？用子程序的实际计算替换call f(xroot,fuc,e0)？原因是我主要是教Fortran77，因为我只使用Fortran进行零星的数值求解，所以我没有花更多的时间学习Fortran90。
您的代码有几个问题阻止我发布完整的答案。仅举几例：(1) 子程序 f1 得到 x 并且根本没有真正使用它。 (2) 在Bisection 中，您将实数与零进行比较：(fuc .eq. 0)。这是一个很大的“禁忌”（事实上，如果你这样做，大多数编译器至少会发出警告）。请改用(abs(fuc) <= epsilon(1.)) 之类的东西。 (3) f1 和 Bisection 中的所有变量都是在没有任何 intent 属性的情况下定义的，这使得任何人，可能你也很难阅读代码并弄清楚过程得到什么以及它返回什么。
@Pap - 是的，你是对的。 (1) 为了简单起见，我简化了很多 f1 的计算（与实际代码相比），因为使用 e0 + 1 会产生相同的错误。在实际功能中，它确实使用了x。 (2) 是的，我通常这样做，但由于这是我开始 Fortran 时的旧子程序，所以我完全忘记了更改它。也谢谢你的建议！ (3) 我不使用 intent 属性，因为（我相信）它不在 Fortran77 中，这是他们教给我的，但也许我应该花一些时间来学习它。