On reading mecej4’s latest contribution I see yet again one of Fortran’s horror stories. The word ‘parameters’ in it may or may not involve things declared as PARAMETER, and p in f(t,y;p) looks to me when I’m a Fortran user like an argument, not a parameter, because that’s what the Fortran standards have called them ever since f66 section 8. Of course mecej4 was correctly using the mathematicians’ definition of ‘parameter’. But 1966 (or even 1957?) was the right time to replace the Fortran technical term PARAMETER by CONSTANT, which would have allowed arguments to be called parameters.
There’s no silver bullet. As always it depends on what are your requirements and what are your preferences. In Modern Fortran in Practice by @Arjen Markus provides a table comparing some of the different callback options:
In a way I think that callbacks were easier in the era of punch-cards and serial programs. You would just place your callback at the top/bottom of your stack of punch cards and declare it as external in the driver program. Still, the callback interface of an early code from Hindmarsh (epsode.f) remains just as valid today as it was in 1974:
c the user must furnish the following double precision subroutines
c and declare them external in his calling program:
c diffun(n,t,y,ydot) computes the function ydot = f(y,t),
c the right hand side of the ordinary
c differential equation system, where y
c and ydot are vectors of length n.
If you don’t like this interface you can always use an internal routine and host association to create an adaptor routine with your interface of choice. This comes at a cost of one indirection (function call), which is probably irrelevant in the grand scheme of things.
The same indirection occurs in MATLAB, but due to the dynamic typing is less intrusive. The ODE solvers in MATLAB generally expect a function of the form:
function dydt = odefun(t,y)
dydt = zeros(2,1); # size declared explicitly, or assumed from input y
dydt(1) = ...;
dydt(2) = ...;
end
If you require a parameterized function, you create a function handle to an anonymous function (here’s the indirection):
# lorenz.m
function dydt = lorenz(t,y,sigma,beta,rho)
dydt = zeros(3,1);
dydt(1) = sigma*(y(2) - y(1))
dydt(2) = y(1)*(rho - y(3)) - y(2)
dydt(3) = y(1)*y(2) - beta*y(3)
end
# solve_lorenz.m
sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,y) lorenz(t,y,sigma,beta,rho);
y0 = [1,2,1];
[tspan,y] = ode45(f,[0 100],y0);
figure
plot3(y(:,1),y(:,2),y(:,3))
It’s the same thing in Python and R. With Julia it’s slightly different in that their JIT compiler can perform some tricks to inline the callback directly into the solver. This can lead to time savings for small systems of ODEs. With Fortran you could do something similar “manually” with a preprocessor such as fypp (I don’t recommend doing this). For large (& small) systems of ODE’s I don’t think it’s worth the trouble. Just for fun, I tried to inline the function call with a solver I’m familiar with:
! lorenz.fypp
module lorenz
implicit none
public
integer, parameter :: wp = kind(1.0e0)
type :: lorenz_params
real(wp) :: sigma, beta, rho
end type
contains
#! N.B.: The argument names y, dydt, and params are fixed!
#:def lorenz()
dydt(1) = params%sigma*(y(2) - y(1))
dydt(2) = y(1)*(params%rho - y(3)) - y(2)
dydt(3) = y(1)*y(2) - params%beta*y(3)
#:enddef
#! The name jac is also fixed!
#:def lorenz_jac()
jac(1, 1) = -params%sigma
jac(2, 1) = params%rho - y(3)
jac(3, 1) = y(2)
jac(1, 2) = params%sigma
jac(2, 2) = -1
jac(3, 2) = y(1)
jac(1, 3) = 0
jac(2, 3) = -y(1)
jac(3, 3) = -params%beta
#:enddef
! subroutine solve_lorenz( ... )
#:include "stiff3_solver.fypp"
$:solve("lorenz",3,lorenz,lorenz_jac,"lorenz_params")
end module
! lorenz_main.f90
program main
use lorenz
implicit none
integer, parameter :: n = 3
real(wp) :: y(n), w(n)
real(wp) :: h0, eps, x0, x1
integer, parameter :: nprint = 1
type(lorenz_params) :: p
p = lorenz_params( &
sigma = 10.0_wp, &
beta = 8.0_wp/3.0_wp, &
rho = 28.0_wp)
! initial value
y = [1.0_wp, 1.0_wp, 1.0_wp]
! initial step size
h0 = 0.001_wp
! tolerance parameters
w = 1
eps = 1.e-4_wp
! time interval
x0 = 0.0_wp
x1 = 100.0_wp
call output(x0,y,0,0.0_wp)
call solve_lorenz(output,nprint,x0,x1,h0,eps,w,y,p)
contains
! Callback routine for output
subroutine output(x,y,iha,qa)
real(wp), intent(in) :: x
real(wp), intent(in) :: y(:)
integer, intent(in) :: iha
real(wp), intent(in) :: qa
write(*,'(4(E19.12,2X),I4,2X,E19.12)') &
x, y(1), y(2), y(3), iha, qa
end subroutine
end program
The final output:
In case anyone is interested here is the full example: lorenz.tar.gz - Google Drive (4.2 KiB).
Regarding such matters, Voltaire had this to say: “The only way to comprehend what mathematicians mean by Infinity is to contemplate the extent of human stupidity.”
2 Likes
For y = f(x) I agree that a mathematician might call x the argument, but sometimes f(x) is, when x is a complex number. Thanks to mecej4 for Voltaire’s snide comment about humans.
In addition to the approaches listed by @ivanpribec, one can use a c_ptr to pass any data structure (and unpack it in the callback function).
https://www.fortran90.org/src/best-practices.html#v-using-type-c-ptr-pointer
Scipy’s odeint has this interface:
odeint(func, y0, t, **args=()**, ...)
- args is an optional tuple of extra arguments to pass to func.
Julia’s ODEproblem has this interface:
ODEProblem(f::ODEFunction, u0, tspan, **p=NullParameters()**;kwargs...)
- p represents the parameters.
The rhs function f that is passed to CVODE’s CvodeInit function has the fom:
f(t, y, ydot, user_data) (for full details see CVRhsFn).
So, the user_data could be specified
- in the
lsoda_classas an allocatable variable or - in the rhs function as an optional, allocatable array, e.g,
subroutine rhs( self, t, y, ydot, ierr, user_data )
real(dp), intent(in) :: t
real(dp), intent(in) :: y( self%neq ) ! Assumes neq is defined in lsoda_class
real(dp), intent(out) :: ydot(size(y))
integer, intent(out) :: ierr
real(dp), intent(in), dimension(:), optional :: user_data
end subroutine
A single array of real isn’t very versatile. Any callback interface should be able to unpack arbitrary data structures through something like void * in c.
1 Like