# 'Varargin' in Fortran

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

• What's not working about your code above? It seems like the best approach is to use an assumed shape array as you have done. – Doug Lipinski May 23 '14 at 15:51
• I agree with @doug-lipinski, I think you have a workable solution here although I am uncertain if you need to really make the copy here. If you want to go beyond what you have done here (a more complex argument for instance) you can pass in arbitrary data through functions without them needing to know what the data looks like by leveraging the ISO C bindings. If you are interested in this solution I can write a more in depth answer. – Kyle Mandli May 23 '14 at 20:35
• @DougLipinski The code is working. But I like to be able to pass a generic list as vararg_in. That means I would like to pass a REAL*8, a COMPLEX*16(:,:) and other stuff all bundled in vararg_in. – DaPhil May 24 '14 at 6:19
• @KyleMandli Yes I am interested in such a solution and would be grateful if you can supply more details. – DaPhil May 24 '14 at 6:21
• Oh. Got it. As long as the number of types you need to support is small, it's simplest to just create multiple assumed shape arrays, one for each type. If desired, you could call the routine using keyword arguments too. Otherwise, Fortran is not the right solution if you have to have this capability. You'll have to use hacks or cross-language solutions like ISO_C_BINDING or just move to another language like C/C++ that supports variadic arguments if you have to do it this way. – Doug Lipinski May 24 '14 at 13:06

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