我的 Fortran 动态分配代码有啥问题？

Posted 2023-03-06

【中文标题】我的 Fortran 动态分配代码有啥问题？

【英文标题】：What is wrong with my Fortran dynamic allocation code?我的 Fortran 动态分配代码有什么问题？

【发布时间】：2020-02-24 21:18:24

【问题描述】：

这是我在互联网上的第一个问题！！！另外，这是一个 Fortan 问题，我刚开始学习这种语言。所以，请原谅我的无知。具体来说，如果我的示例代码很大，我深表歉意。我已尽最大努力将其尺寸缩小到最低限度，同时又不影响其清晰度。

这是问题所在：我正在尝试创建一个动态仿真模型，该模型使用一个在模型运行期间将被多次调用的函数。函数本身是一个动态集成系统，因此必须在模拟的每个时间步保留函数的值，以便在下一个时间步使用它们。

该模型包括一个简单的库存 (P)，它随着时间的推移以恒定的速率（每次 10%）减少。该函数是指数延迟。它需要 4 个参数：（输入、延迟时间、延迟顺序、输入初始值）。请注意，函数的输出具有 3 个维度。第一个维度代表程序中延迟调用的不同实例。第二个维度表示延迟的顺序。第三个维度代表模拟时间。

这里的主要挑战是动态分配将由函数创建的数组。为了解决这个问题，我尝试遵循此处提供的指南：How to increase array size on-the-fly in Fortran? 我编写的程序听起来在某些情况下运行良好，但并非总是如此。例如，如果您编译如下所示的代码，您将得到正确的结果（正确的行为是所有变量都应从 100 开始并逐渐下降）。但是如果你把 U 的延迟顺序从 5 改到 6，那么结果就会出错（U 从 100 开始，然后波动，有时会超过 100）。

我在这里缺少什么？我实现 move_alloc 的方式有什么问题吗？还是问题出在其他地方？可能是编译器（gfortran）失败了吗？ 我真的很感谢您的帮助，提前！

program model

  implicit none
  real, parameter :: start = 0, end = 10, dt = 0.25
  integer, parameter :: n = int((end - start)/dt)
  real, dimension(0:n+1) :: time, P, R, S, W, U
  integer :: i
  real, dimension(:,:,:), allocatable :: DN, temp

  P(0) = 100
  time(0) = start

  open (unit=10, file='out.txt', status='replace', action='write', position='rewind')
  write (10, *) 'Time, ', 'P, ', 'R, ', 'S, ', 'W, ', 'U'

  do i = 0, n
     P(i+1) = P(i) - .1*P(i)*dt
     R(i) = delay(P(i), 3.0, 4, P(0))
     S(i) = delay(P(i), 2.0, 4, P(0))
     W(i) = delay(P(i), 4.0, 5, P(0))
     U(i) = delay(P(i), 1.0, 5, P(0))
     time(i+1) = time(i) + dt
     write (10, *) time(i), ',', P(i), ',', R(i), ',', S(i), ',', W(i), ',', U(i)
  end do

  close (10)

contains
  function delay(input, dtime, order, init)
    real :: delay
    real, intent(in) :: input, dtime, init
    integer, save :: j, k
    integer, intent(in) :: order
    integer :: l
    data j /0/
    data k /0/
    if (j == i) then
       k = k + 1
       allocate(temp(1:k, 1:order, 0:n+1))
       if (allocated(DN)) temp(1:k-1, 1:order, 0:n+1) = DN
       call move_alloc(temp, DN)
    else
       k = 1
       j = i
    end if
    if (i == 0) then
       do l = 1, order
          DN(k, l, 0) = init * dtime / order
       end do
    end if
    DN(k, 1, i+1) = DN(k, 1, i) + (input - DN(k, 1, i) * order / dtime) * dt
    if (order > 1) then
       do l = 2, order
          DN(k, l, i+1) = DN(k, l, i) + (DN(k, l-1, i) - DN(k, l, i)) * order * dt / dtime
       end do
    end if
    delay = DN(k, order, i) * order / dtime
  end function delay

end program model

【参考方案1】：

您的数组索引不正确，应该使用编译器运行时检查来检测程序中的这些错误。

在函数delay 中，您依赖于存储在DN 中的状态，并使用order 虚拟参数引用此数组的部分。不幸的是，order 没有正确绑定到DN 的实际形状。

考虑线条

   allocate(temp(1:k, 1:order, 0:n+1))
   if (allocated(DN)) temp(1:k-1, 1:order, 0:n+1) = DN
   call move_alloc(temp, DN)

很明显，在这组行第一次执行后，DN 被分配并具有[1, 4, 42] 的形状。在第二个之后，同样清楚的是它的形状是[2 4 42]。来到函数的第三个条目，order 现在是 5，因此分配中的形状不匹配。

我发现代码非常不清楚，因此不会尝试提供有关如何更正它的建议。特别是，保存的本地状态和来自主机范围的变量的混合令人困惑。

【讨论】：

非常感谢@francescalus 花时间检查我的业余代码。我同意订单是导致问题的原因，但我想不出任何其他方法来解决它。我想我需要更多地尝试不同的编程技术来找到解决方案。