在Fortran中传递牛顿方法中的其他参数

Question

我在实现一种在Fortran程序中调用Newton方法的方法时遇到了麻烦。所以我想用牛顿的方法来解决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

我的理解是，如果我将争论对象放在求解器模块中，那么所有求解器调用都将使用相同的参数（我仍然可以使用上述方法将所有这些参数保存在模块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这些参数，而无需触及这些函数的签名。产生的效果就像改变函数的arity一样。

使用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

Answer 2

在与同事讨论后，我解决了我的问题2。

如果我们只有一个求解器模块的参数对象，那么所有的调用都将访问相同的参数，因为它们共享相同的内存空间。

为避免这种情况，我们希望将参数对象作为参数传递给求解器。因此，我们将重新编写牛顿方法，而不是使用默认的求解器子程序，因此它可以接受额外的参数。

（我之前使用过最简单的Newton子程序，因为我想保持它不受影响。）

通过这种方式，我们将定义一个参数对象数组，并在运行时传递它们。

谢谢你的意见。