'Varargin' in Fortran

Question

Consider the following: I would like to solve a differential equation in Fortran (with a Runge-Kutta routine), e.g. x'=-a*x, where a is a constant. I then want to solve this for different a. I now would like to write the code such that it can be applied to other problems. But in general I do not know what additional input my differential equation may need (and what size and type they have). So my question is if it is possible in Fortran to pass all additional parameters to my Runge-Kutta method (at best with a single property) so that it passes it on to my differential equation so that it can be used there? With a known set of parameters this is easy but what about an unknown set?

My approach so far is this:

PROGRAM main
  USE RKModule

  IMPLICIT NONE

  INTEGER :: l

  REAL*8, PARAMETER :: num_t = 1e3
  REAL*8, PARAMETER :: min_t = 0.0d0
  REAL*8, PARAMETER :: max_t = 10.0d0

  REAL*8, DIMENSION(1), PARAMETER :: a = 1.8d0

  REAL*8, DIMENSION(num_t) :: dGrid_t, grid_t
  COMPLEX*16, DIMENSION(num_t) :: sol
  COMPLEX*16, DIMENSION(1) :: x_in, x_out

  DO l = 1, num_t
    dGrid_t(l) = (max_t - min_t) / DBLE(num_t)
    grid_t(l) = min_t + DBLE(l-0.5d0)*dGrid_t(l)
  END DO

  x_in(1) = dcmplx(1.0d0, 0.0d0)
  x_out(1) = dcmplx(0.0d0, 0.0d0)

  DO l = 1, num_t
    sol(l) = x_in(1)
    CALL rk45(x_out, DGL, x_in, grid_t(l), dGrid_t(l), a)
    x_in = x_out
  END DO

  OPEN(UNIT = 1000, FILE = 'sol.dat', ACTION = 'WRITE')
  WRITE(1000, FMT = '(e14.7e3, 1x, e14.7e3)', ADVANCE = "NO") sol
  CLOSE(1000)

  CONTAINS

  SUBROUTINE DGL(t,f,df,vararg_in)
    REAL*8, INTENT(IN) :: t
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in

    REAL*8, DIMENSION(1) :: b

    b = 1.0d0
    IF (PRESENT(vararg_in)) THEN
      b = vararg_in
    END IF

    df = -b*f

  END SUBROUTINE DGL

END PROGRAM main

And the Runge Kutta Routine:

MODULE RKModule
  IMPLICIT NONE

  CONTAINS

  SUBROUTINE rk45(x_out, RHS, x_in, t, dt, vararg_in)
    INTERFACE
      SUBROUTINE RHS(t,f,df,a)
    REAL*8, INTENT(IN) :: t
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: a
      END SUBROUTINE
    END INTERFACE
    REAL*8, INTENT(IN) :: t, dt 
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: x_in
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: x_out
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in

    REAL*8 :: tHlp
    COMPLEX*16, DIMENSION(SIZE(x_in)) :: k1 ,k2, k3, k4, dx, x_hlp
    REAL*8, DIMENSION(:), ALLOCATABLE :: vararg_in_copy

    IF (PRESENT(vararg_in)) THEN
      ALLOCATE(vararg_in_copy(SIZE(vararg_in)))
      vararg_in_copy = vararg_in

      tHlp = t
      x_hlp = x_in
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k1 = dt * dx

      tHlp = t + 0.5d0 * dt
      x_hlp = x_in + 0.5d0 * k1
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k2 = dt * dx

      x_hlp = x_in + 0.5d0 * k2
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k3 = dt * dx

      tHlp = t + dt
      x_hlp = x_in + k3
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k4 = dt * dx

      x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0

      DEALLOCATE(vararg_in_copy)
    ELSE

      tHlp = t
      x_hlp = x_in
      CALL RHS(tHlp, x_hlp, dx)
      k1 = dt * dx

      tHlp = t + 0.5d0 * dt
      x_hlp = x_in + 0.5d0 * k1
      CALL RHS(tHlp, x_hlp, dx)
      k2 = dt * dx

      x_hlp = x_in + 0.5d0 * k2
      CALL RHS(tHlp, x_hlp, dx)
      k3 = dt * dx

      tHlp = t + dt
      x_hlp = x_in + k3
      CALL RHS(tHlp, x_hlp, dx)
      k4 = dt * dx

      x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0

    END IF

  END SUBROUTINE rk45

END MODULE RKModule

Answer 1

This question better suits to StackOverflow.

You can use class(*) or type(c_ptr). Both come from Fortran 2003. The latter is widely supported, the former only since gfortran 4.8 and other recent versions of other compilers. Note type(c_ptr) is very useful even inside Fortran, without calling any C.

The first approach:

subroutine RK4(normal arguments, vararg_in)
  class(*), intent(in) :: vararg_in(:)  !for scalar delete (:)

  !if you need to preserve a copy (why?)
  allocate(vararg_copy, source=varrarg_in)

  call RHS(...,vararg_copy) !or directly vararg_in
end subroutine


subroutine DGL(...,vararg_in)
  class(*), intent(in) :: vararg_in(:)  !for scalar delete (:)

  select type (vararg_in)
    select type (expected_type)
     !do something
    class default
     !some exception
  end select
end subroutine

You can pass anything to class(*), possibly a derived type, if you need variables of differing types.

The other approach:

subroutine RK4(normal arguments, vararg_in)
  use iso_c_binding
  type(c_ptr), intent(in) :: vararg_in_ptr

  !Doing a copy of the pointer doesn't make much sense

  call RHS(...,vararg_in_ptr)
end subroutine


subroutine DGL(...,vararg_in_ptr)
  use iso_c_binding
  type(c_ptr), intent(in) :: vararg_in_ptr
  type(expected_type), pointer :: vararg_in(:) !for scalar delete (:)

  !for an array
  call c_f_pointer(vararg_in_ptr, vararg_in, expected_shape)
  !for a scalar just
  call c_f_pointer(vararg_in_ptr, vararg_in)

  !do something
end subroutine