Passing additional arguments in Newton’s method in Fortran

Question

I am having trouble in implementing an approach to call Newton's method in a Fortran program. So I want to use Newton's method to solve an equation following the link

However, my program is slightly different with the example above. In my case, the equation requires some additional information which are produced during runtime.

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

which means the f is calculated not only based on x, but also a few other variables (but x is the unknown).

I have a solution, which only works for a simple case: I used a module to include Newton's solver. I also defined a derived data type to hold all the arguments inside the module. It works good now.

My question is: I need to call the Newton's method many times, and each time the arguments are different. How should I design the structure of the modules? Or should I use a different solution?

I provided a simple example below:

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

And the calling module is:

!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

My understanding is that if I put the argu object inside the solver module, then all solver calling will use the same arguments (I can still save all of them inside module A using the above method). In that case, I have to update the arguments during each for loop?

Because the program runs using MPI/OpenMP, I want to make sure there is no overwritten among threads. Thank you.

Answer 1

There is a common pattern in modern Fortran for the problem you are facing (partial function application). Unlike other languages, Fortran doesn't have function closures, so making a lexical scope for a function is a little "convoluted" and kind of limited.

You should really consider revisiting all the links @VladmirF shared on the comment, most of them apply straightforwardly to your case. I will give you an example of a solution.

This is a solution without using a wrapper type. I will use a feature included in Fortran 2008 standard: passing an internal procedure as an argument. It is compatible with the latest gfortran, Intel and many others. If you can't access a compiler with this feature or if you prefer a solution with a derived type, you can refer to this answer.

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

The rationale of this code is: as f and fp lay inside of the solve_quad procedure, they have access to the arguments a, b and c by host association, without touching those function's signatures. The resulting effect is like changing the arity of the function.

Testing it with gfortran 8.0 and the solver implementation from the link you shared, I got this:

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