fortran - 在 Fortran 牛顿法中传递附加参数

Question

我在实现在 Fortran 程序中调用牛顿法的方法时遇到问题。所以我想用牛顿法来求解the link下面的方程

但是，我的程序与上面的示例略有不同。就我而言，该方程需要一些在运行时生成的附加信息。

subroutine solve(f, fp, x0, x, iters, debug)

这意味着 f 不仅基于 x 计算，还基于其他一些变量(但 x 是未知的)。

我有一个解决方案，仅适用于简单的情况:我使用了一个模块来包含牛顿求解器。我还定义了一个派生数据类型来保存模块内的所有参数。现在效果很好。

我的问题是:我需要多次调用牛顿方法，并且每次的参数都不同。我应该如何设计模块的结构？或者我应该使用不同的解决方案？

我在下面提供了一个简单的示例:

module solver
  type argu
    integer :: m
  end type argu
  type(argu):: aArgu_test  !should I put here?
  contains
    subroutine solve(f, fp, x0, x, iters, debug)
       ...
       !m is used inside here
    end subroutine solve
    subroutine set_parameter(m_in)
       aArgu%m = m_in
    end subroutine set_parameter()

end module solver

调用模块是:

!only one set of argument, but used many times
module A
   use module solver

   do i = 1, 4, 1
     set_parameter(i) 
     !call Newtow method
     ...
   enddo
end module A

!can I use an array for argu type if possible?
module B
   use module solver
   type(argu), dimension(:), allocable :: aArgu ! or should I put here or inside solver module?

end module B

我的理解是，如果我将 argu 对象放入求解器模块中，那么所有求解器调用都将使用相同的参数(我仍然可以使用上述方法将它们全部保存在模块 A 中)。在这种情况下，我必须在每个 for 循环期间更新参数？

因为程序使用 MPI/OpenMP 运行，所以我想确保线程之间没有覆盖。谢谢。

Answer 1

现代 Fortran 中有一个常见的模式可以解决您面临的问题 ( partial function application )。与其他语言不同，Fortran 没有函数闭包，因此为函数创建词法作用域有点“复杂”并且受到限制。

您确实应该考虑重新访问 @VladmirF 在评论中分享的所有链接，其中大多数直接适用于您的案例。我会给你一个解决方案的例子。

这是一个不使用包装类型的解决方案。我将使用 Fortran 2008 标准中包含的一项功能:将内部过程作为参数传递。它与最新的 gfortran、Intel 等兼容。如果您无法访问具有此功能的编译器或者您更喜欢具有派生类型的解决方案，您可以引用this回答。

module without_custom_type
  use, intrinsic :: iso_fortran_env, only: r8 => real64
  use :: solver

contains
  subroutine solve_quad(a, b, c, x0, x, iters, debug)
    integer, intent(in) :: a, b, c
    real(r8), intent(in) :: x0
    real(r8), intent(out) :: x
    integer, intent(out) :: iters
    logical, intent(in) :: debug

    call solve(f, fp, x0, x, iters, debug)

  contains
    real(r8) function f(x)
      real(r8),intent(in) :: x
      f = a * x * x + b * x + c
    end

    real(r8) function fp(x)
      real(r8),intent(in) :: x
      fp = 2 * a * x + b
    end
  end
end

此代码的基本原理是:由于 f 和 fp 位于 solve_quad 过程内部，因此它们可以访问参数 a、b 和 c 通过主机关联，而不触及这些函数的签名。产生的效果就像改变函数的数量。

使用 gfortran 8.0 和您共享的链接中的 solver 实现对其进行测试，我得到了以下结果:

program test
  use, intrinsic :: iso_fortran_env, only: r8 => real64
  use :: without_custom_type
  implicit none

  real(r8) :: x, x0
  integer :: iters
  integer :: a = 1, b = -5, c = 4

  x0 = 0
  call solve_quad(a, b, c, x0, x, iters, .false.)
  print *, x, iters
  ! output: 1.0000000000000000, 5

  x0 = 7
  call solve_quad(a, b, c, x0, x, iters, .false.)
  print *, x, iters
  ! output: 4.0000000000000000, 6    
end