\documentclass[article,nojss]{jss}
\DeclareGraphicsExtensions{.pdf,.eps}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Add-on packages and fonts
\usepackage{amsmath}
\usepackage{xspace}
\usepackage{verbatim}
\usepackage[english]{babel}
%\usepackage{mathptmx}
%\usepackage{helvet}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\newcommand{\Rmodels}{\textbf{\textsf{R models}}\xspace}
\newcommand{\DLLmodels}{\textbf{\textsf{DLL models}}\xspace}

\title{\proglang{R} Package \pkg{deSolve}, Writing Code in Compiled
Languages}

\Plaintitle{R Package deSolve, Writing Code in Compiled Languages}

\Keywords{differential equation solvers, compiled code, performance,
\proglang{FORTRAN}, \proglang{C}}

\Plainkeywords{differential equation solvers, compiled code,
performance, FORTRAN, C}


\author{
  Karline Soetaert\\
  Royal Netherlands Institute\\
  of Sea Research (NIOZ)\\
  Yerseke\\
  The Netherlands
  \And
  Thomas Petzoldt\\
  Technische Universit\"at \\
  Dresden\\
  Germany
  \And
  R. Woodrow Setzer\\
  National Center for\\ Computational Toxicology\\
  US Environmental Protection Agency
}

\Plainauthor{Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer}

\Abstract{This document describes how to use the \pkg{deSolve} package
\citep{deSolve_jss}  to solve models that are written in \proglang{FORTRAN} or
\proglang{C}.}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Karline Soetaert\\
  Royal Netherlands Institute of Sea Research (NIOZ)\\
  4401 NT Yerseke, Netherlands \\
  E-mail: \email{karline.soetaert@nioz.nl}\\
  URL: \url{https://www.nioz.nl}\\
  \\
  Thomas Petzoldt\\
  Institut f\"ur Hydrobiologie\\
  Technische Universit\"at Dresden\\
  01062 Dresden, Germany\\
  E-mail: \email{thomas.petzoldt@tu-dresden.de}\\
  URL: \url{https://tu-dresden.de/Members/thomas.petzoldt/}\\
   \\
  R. Woodrow Setzer\\
  National Center for Computational Toxicology\\
  US Environmental Protection Agency
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R/Sweave specific LaTeX commands.
%% need no \usepackage{Sweave}
%\VignetteIndexEntry{R Package deSolve: Writing Code in Compiled Languages}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin of the document
\begin{document}
\SweaveOpts{engine=R,eps=FALSE}

<<preliminaries,echo=FALSE,results=hide>>=
library("deSolve")
options(prompt = "R> ")
options(width=70)
@

\maketitle

\section{Introduction}

\pkg{deSolve} \citep{deSolve_jss,deSolve}, the successor of \proglang{R} package
\pkg{odesolve} \citep{Setzer01} is a package to solve ordinary
differential equations (ODE), differential algebraic equations (DAE)
and partial differential equations (PDE).
One of the prominent features of \pkg{deSolve} is that it allows
specifying the differential equations either as:

\begin{itemize}
\item pure \proglang{R} code \citep{Rcore},
\item functions defined in lower-level languages such as
  \proglang{FORTRAN}, \proglang{C}, or \proglang{C++}, which are
  compiled into a dynamically linked library (DLL) and loaded into
  \proglang{R}.
\end{itemize}

In what follows, these implementations will be referred to as \Rmodels
and \DLLmodels respectively.
Whereas \Rmodels are easy to implement, they allow simple interactive
development, produce highly readible code and access to \proglang{R}s
high-level procedures, \DLLmodels have the benefit of increased
simulation speed. Depending on the problem, there may be a gain of up
to several orders of magnitude computing time when using compiled
code.

Here are some rules of thumb when it is worthwhile or not to switch to
\DLLmodels:

\begin{itemize}
\item As long as one makes use only of \proglang{R}s high-level
  commands, the time gain will be modest. This was demonstrated in
  \citet{deSolve_jss}, where a formulation of two interacting
  populations dispersing on a 1-dimensional or a 2-dimensional grid
  led to a time gain of a factor two only when using \DLLmodels.

\item Generally, the more statements in the model, the higher will be
  the gain of using compiled code. Thus, in the same paper
  \citep{deSolve_jss}, a very simple, 0-D, Lotka-Volterrra type of model
  describing only 2 state variables was solved 50 times faster when
  using compiled code.

\item As even \Rmodels are quite performant, the time gain induced by
  compiled code will often not be discernible when the model is only
  solved once (who can grasp the difference between a run taking 0.001
  or 0.05 seconds to finish).  However, if the model is to be applied
  multiple times, e.g.  because the model is to be fitted to data, or
  its sensitivity is to be tested, then it may be worthwhile to
  implement the model in a compiled language.
\end{itemize}

Starting from \pkg{deSolve} version 1.4, it is now also possible to use
\emph{forcing functions} in compiled code. These forcing functions are
automatically updated by the integrators. See last chapter.

\section{A simple ODE example}

Assume the following simple ODE (which is from the \code{LSODA} source
code):

\begin{align*}
 \frac{{dy_1}}{{dt}} &=  - k_1 \cdot y_1 + k_2 \cdot y_2 \cdot y_3 \\
 \frac{{dy_2}}{{dt}} &= k_1 \cdot y_1 - k_2 \cdot y_2 \cdot y_3 - k_3 \cdot y_2 \cdot y_2 \\
 \frac{{dy_3}}{{dt}} &= k_3 \cdot y_2 \cdot y_2 \\
\end{align*}

where $y_1$, $y_2$ and $y_3$ are state variables, and $k_1$, $k_2$ and
$k_3$ are parameters.

We first implement and run this model in pure \proglang{R}, then show
how to do this in \proglang{C} and in \proglang{FORTRAN}.

\subsection{ODE model implementation in R}

An ODE model implemented in \textbf{pure \proglang{R}} should be
defined as:

\begin{verbatim}
 yprime = func(t, y, parms, ...)
\end{verbatim}

where \code{t} is the current time point in the integration, \code{y}
is the current estimate of the variables in the ODE system, and
\code{parms} is a vector or list containing the parameter
values. The optional dots argument (\code{\dots}) can be used to pass any other arguments to the
function.  The return value of \code{func} should be a list, whose
first element is a vector containing the derivatives of \code{y} with
respect to time, and whose next elements contain output variables that
are required at each point in time.

The \proglang{R} implementation of the simple ODE is given below:

<<the_Rmodel>>=
model <- function(t, Y, parameters) {
  with (as.list(parameters),{

    dy1 = -k1*Y[1] + k2*Y[2]*Y[3]
    dy3 = k3*Y[2]*Y[2]
    dy2 = -dy1 - dy3

    list(c(dy1, dy2, dy3))
  })
}
@


The Jacobian ($\frac{{\partial y'}}{{\partial y}}$) associated to the
above example is:

<<Jacobian_in_R>>=
jac <- function (t, Y, parameters) {
  with (as.list(parameters),{

    PD[1,1] <- -k1
    PD[1,2] <- k2*Y[3]
    PD[1,3] <- k2*Y[2]
    PD[2,1] <- k1
    PD[2,3] <- -PD[1,3]
    PD[3,2] <- k3*Y[2]
    PD[2,2] <- -PD[1,2] - PD[3,2]

    return(PD)
  })
}
@

This model can then be run as follows:

<<Run_Rmodel>>=
parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7)
Y     <- c(1.0, 0.0, 0.0)
times <- c(0, 0.4*10^(0:11))
PD    <- matrix(nrow = 3, ncol = 3, data = 0)
out   <- ode(Y, times, model, parms = parms, jacfunc = jac)
@

\subsection{ODE model implementation in C}
\label{sec:Cexamp}

In order to create compiled models (.DLL = dynamic link libraries on
Windows or .so = shared objects on other systems) you must have a
recent version of the GNU compiler suite installed, which is quite
standard for Linux. Windows users find all the required tools on
\url{https://cran.r-project.org/bin/windows/Rtools/}.  Getting DLLs
produced by other compilers to communicate with R is much more
complicated and therefore not recommended. More details can be found
on \url{https://cran.r-project.org/doc/manuals/R-admin.html}.

The call to the derivative and Jacobian function is more complex for
compiled code compared to \proglang{R}-code, because it has to comply
with the interface needed by the integrator source codes.

Below is an implementation of this model in \proglang{C}:

\verbatiminput{mymod.c}

The implementation in \proglang{C} consists of three parts:

\begin{enumerate}
\item After defining the parameters in global \proglang{C}-variables,
  through the use of \code{\#define} statements, a function called
  \code{initmod} initialises the parameter values, passed from the
  \proglang{R}-code.

  This function has as its sole argument a pointer to
  \proglang{C}-function \code{odeparms} that fills a double array with
  double precision values, to copy the parameter values into the
  global variable.

\item Function \code{derivs} then calculates the values of the
  derivatives. The derivative function is defined as:

\begin{verbatim}
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip)
\end{verbatim}

  where \code{*neq} is the number of equations, \code{*t} is the value
  of the independent variable, \code{*y} points to a double precision
  array of length \code{*neq} that contains the current value of the
  state variables, and \code{*ydot} points to an array that will
  contain the calculated derivatives.

  \code{*yout} points to a double precision vector whose first
  \code{nout} values are other output variables (different from the
  state variables \code{y}), and the next values are double precision values
  as passed by parameter \code{rpar} when calling the integrator.  The
  key to the elements of \code{*yout} is set in \code{*ip}

  \code{*ip} points to an integer vector whose length is at least 3;
  the first element (\code{ip[0]}) contains the number of output
  values (which should be equal or larger than \code{nout}), its
  second element contains the length of \code{*yout}, and the third
  element contains the length of \code{*ip}; next are integer values,
  as passed by parameter \code{ipar} when calling the
  integrator.\footnote{Readers familiar with the source code of the
    \pkg{ODEPACK} solvers may be surprised to find the double
    precision vector \code{yout} and the integer vector \code{ip} at
    the end. Indeed none of the \pkg{ODEPACK} functions allow this,
    although it is standard in the \code{vode} and \code{daspk} codes.
    To make all integrators compatible, we have altered the
    \pkg{ODEPACK} \proglang{FORTRAN} codes to consistently pass these
    vectors.}

  Note that, in function \code{derivs}, we start by checking whether
  enough memory is allocated for the output variables (\code{if (ip[0] <
    1)}), else an error is passed to \proglang{R} and the integration
  is stopped.

\item In \proglang{C}, the call to the function that generates the
  Jacobian is as:

\begin{verbatim}
void jac(int *neq, double *t, double *y, int *ml,
           int *mu, double *pd, int *nrowpd, double *yout, int *ip)
\end{verbatim}

  where \code{*ml} and \code{*mu} are the number of non-zero bands
  below and above the diagonal of the Jacobian respectively. These
  integers are only relevant if the option of a banded Jacobian is
  selected. \code{*nrow} contains the number of rows of the Jacobian.
  Only for full Jacobian matrices, is this equal to \code{*neq}. In
  case the Jacobian is banded, the size of \code{*nrowpd} depends on the
  integrator. If the method is one of \code{lsode, lsoda, vode},
  then \code{*nrowpd}  will be equal to
  \code{*mu + 2 * *ml + 1}, where the last \code{*ml} rows should be
  filled with $0$s.

  For \code{radau}, \code{*nrowpd} will be equal to \code{*mu + *ml + 1}

  See example ``odeband'' in the directory \url{doc/examples/dynload},
  and chapter \ref{band}.
\end{enumerate}


\subsection{ODE model implementation in FORTRAN}
\label{sec:forexamp}

Models may also be defined in \proglang{FORTRAN}.

\verbatiminput{mymod.f}

In \proglang{FORTRAN}, parameters may be stored in a common block
(here called \code{myparms}).  During the initialisation, this common
block is defined to consist of a 3-valued vector (unnamed), but in the
subroutines \code{derivs} and \code{jac}, the parameters are given a
name (\code{k1}, ...).

\subsection{Running ODE models implemented in compiled code}

To run the models described above, the code in \code{mymod.f} and
\code{mymod.c} must first be compiled\footnote{This requires a
  correctly installed GNU compiler, see above.}.
This can simply be done in \proglang{R} itself, using the
\code{system} command:

<<compile_DLLmodel_F,eval=FALSE>>=
  system("R CMD SHLIB mymod.f")
@

for the \proglang{FORTRAN} code or

<<compile_DLLmodel_C,eval=FALSE>>=
  system("R CMD SHLIB mymod.c")
@

for the \proglang{C} code.

This will create file \code{mymod.dll} on windows, or \code{mymod.so} on
other platforms.

We load the DLL, in windows as:
\begin{verbatim}
dyn.load("mymod.dll")
\end{verbatim}
and in unix:
\begin{verbatim}
dyn.load("mymod.so")
\end{verbatim}
or, using a general statement:

\begin{verbatim}
dyn.load(paste("mymod", .Platform$dynlib.ext, sep = ""))
\end{verbatim}
The model can now be run as follows:
\begin{verbatim}
parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7)
Y     <- c(y1 = 1.0, y2 = 0.0, y3 = 0.0)
times <- c(0, 0.4*10^(0:11) )

out <- ode(Y, times, func = "derivs", parms = parms,
           jacfunc = "jac", dllname = "mymod",
           initfunc = "initmod", nout = 1, outnames = "Sum")
\end{verbatim}

The integration routine (here \code{ode}) recognizes that the model is
specified as a DLL due to the fact that arguments \code{func} and
\code{jacfunc} are not regular \proglang{R}-functions but character
strings. Thus, the integrator will check whether the function is
loaded in the DLL with name \code{mymod}.

Note that \code{mymod}, as specified by \code{dllname} gives the name of
the shared library \emph{without extension}. This DLL should contain
all the compiled function or subroutine definitions referred to in
\code{func}, \code{jacfunc} and \code{initfunc}.

Also, if \code{func} is specified in compiled code, then
\code{jacfunc} and \code{initfunc} (if present) should also be
specified in a compiled language.  It is not allowed to mix
\proglang{R}-functions and compiled functions.

Note also that, when invoking the integrator, we have to specify the
number of ordinary output variables, \code{nout}. This is because the
integration routine has to allocate memory to pass these output
variables back to \proglang{R}. There is no way to check for the
number of output variables in a DLL automatically.  If in the calling
of the integration routine the number of output variables is too low,
then \proglang{R} may freeze and need to be terminated! Therefore it
is advised that one checks in the code whether \code{nout} has been specified
correctly.  In the \proglang{FORTRAN} example above, the statement
\code{if (ip(1) < 1) call rexit("nout should be at least 1")} does
this.  Note that it is not an error (just a waste of memory) to set
\code{nout} to a too large value.

Finally, in order to label the output matrix, the names of the ordinary
output variables have to be passed explicitly (\code{outnames}). This is not
necessary for the state variables, as their names are known through their
initial condition (\code{y}).

\section{Alternative way of passing parameters and data in compiled code}
\label{sec:parms}

All of the solvers in \pkg{deSolve} take an argument \code{parms}
which may be an arbitrary \proglang{R} object.  In models defined in
\proglang{R} code, this argument is passed unprocessed to the various
functions that make up the model.  It is possible, as well, to pass
such R-objects to models defined in native code.

The problem is that data passed to, say, \code{ode} in the argument
\code{parms} is not visible by default to the routines that define the
model.  This is handled by a user-written initialization function, for
example \code{initmod} in the \proglang{C} and \proglang{FORTRAN}
examples from sections \ref{sec:Cexamp} and \ref{sec:forexamp}.
However, these set only the \emph{values} of the parameters.

R-objects have many attributes that may also be of interest. To have access
to these, we need to do more work, and this mode of passing parameters and
data is much more complex than what we saw in previous chapters.

In \proglang{C}, the initialization routine is declared:

\begin{verbatim}
void initmod(void (* odeparms)(int *, double *));
\end{verbatim}

That is, \code{initmod} has a single argument, a pointer to a function
that has as arguments a pointer to an \texttt{int} and a pointer to a
\texttt{double}. In \proglang{FORTRAN}, the initialization routine has
a single argument, a subroutine declared to be external.  The name of
the initialization function is passed as an argument to the
\pkg{deSolve} solver functions.

In \proglang{C}, two approaches are available for making the values
passed in \code{parms} visible to the model routines, while only the
simpler approach is available in \proglang{FORTRAN}.  The simpler
requires that \code{parms} be a numeric vector.  In \proglang{C}, the
function passed from \pkg{deSolve} to the initialization function
(called \code{odeparms} in the example) copies the values from the
parameter vector to a static array declared globally in the file where
the model is defined.  In \proglang{FORTRAN}, the values are copied
into a \code{COMMON} block.

It is possible to pass more complicated structures to \proglang{C}
functions.  Here is an example, an initializer called
\code{deltamethrin} from a model describing the pharmacokinetics of
that pesticide:

\begin{verbatim}
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Rdynload.h>
#include "deltamethrin.h"

/* initializer */
void deltamethrin(void(* odeparms)(int *, double *))
{
  int Nparms;
  DL_FUNC get_deSolve_gparms;
  SEXP gparms;
  get_deSolve_gparms = R_GetCCallable("deSolve","get_deSolve_gparms");
  gparms = get_deSolve_gparms();
  Nparms = LENGTH(gparms);
  if (Nparms != N_PARMS) {
    PROBLEM "Confusion over the length of parms"
      ERROR;
  } else {
    _RDy_deltamethrin_parms = REAL(gparms);
  }
}
\end{verbatim}

In \texttt{deltamethrin.h}, the variable
\code{\_RDy\_deltamethrin\_parms} and macro N\_PARMS are declared:

\begin{verbatim}
#define N_PARMS 63
static double *_RDy_deltamethrin_parms;
\end{verbatim}

The critical element of this method is the function
\code{R\_GetCCallable} which returns a function (called
\code{get\_deSolve\_gparms} in this implementation) that returns the
parms argument as a \code{SEXP} data type.  In this example,
\code{parms} was just a real vector, but in principle, this method can
handle arbitrarily complex objects.  For more detail on handling
\proglang{R} objects in native code, see \proglang{R} Development Core
Team (2008).

\section{deSolve integrators that support DLL models}

In the most recent version of \pkg{deSolve} all integration routines
can solve \DLLmodels. They are:

\begin{itemize}
\item all solvers of the \code{lsode} familiy: \code{lsoda},
  \code{lsode}, \code{lsodar}, \code {lsodes},
\item \code{vode}, \code{zvode},
\item \code{daspk},
\item \code{radau},
\item the Runge-Kutta integration routines (including the Euler method).
\end{itemize}

For some of these solvers the interface is slightly different
(e.g. \code{zvode, daspk}), while in others (\code{lsodar}, \code{lsodes})
different functions can be defined.  How this is implemented in a
compiled language is discussed next.

\subsection{Complex numbers, function zvode}

\code{zvode} solves ODEs that are composed of complex variables.
The program below uses \code{zvode} to solve the following system of
2 ODEs:

\begin{align*}
 \frac{{dz}}{{dt}} &=  i \cdot z\\
 \frac{{dw}}{{dt}} &= -i \cdot w \cdot w \cdot z\\
\end{align*}
where
\begin{align*}
w(0) = 1/2.1 +0i\\
z(0) = 1i
\end{align*}
on the interval t = [0, 2 $\pi$]

The example is implemented in \proglang{FORTRAN}%
\footnote{this can be found in file "zvodedll.f",
in the dynload subdirectory of the package},
\code{FEX} implements the function \code{func}:

\begin{verbatim}
      SUBROUTINE FEX (NEQ, T, Y, YDOT, RPAR, IPAR)
      INTEGER NEQ, IPAR(*)
      DOUBLE COMPLEX Y(NEQ), YDOT(NEQ), RPAR(*), CMP
      DOUBLE PRECISION T
      character(len=100) msg

c the imaginary unit i
      CMP = DCMPLX(0.0D0,1.0D0)

      YDOT(1) = CMP*Y(1)
      YDOT(2) = -CMP*Y(2)*Y(2)*Y(1)

      RETURN
      END
\end{verbatim}
\code{JEX} implements the function \code{jacfunc}
\begin{verbatim}
      SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD, RPAR, IPAR)
      INTEGER NEQ, ML, MU, NRPD, IPAR(*)
      DOUBLE COMPLEX Y(NEQ), PD(NRPD,NEQ), RPAR(*), CMP
      DOUBLE PRECISION T
c the imaginary unit i
      CMP = DCMPLX(0.0D0,1.0D0)

      PD(2,3) = -2.0D0*CMP*Y(1)*Y(2)
      PD(2,1) = -CMP*Y(2)*Y(2)
      PD(1,1) = CMP
      RETURN
      END
\end{verbatim}

Assuming this code has been compiled and is in a DLL called "zvodedll.dll",
this model is solved in R as follows:

\begin{verbatim}
dyn.load("zvodedll.dll")
outF <- zvode(func = "fex", jacfunc = "jex", y = yini, parms = NULL,
          times = times, atol = 1e-10, rtol = 1e-10, dllname = "zvodedll",
          initfunc = NULL)
\end{verbatim}

Note that in \proglang{R} names of \proglang{FORTRAN} DLL functions
(e.g. for \code{func} and \code{jacfunc}) have to be given in
lowercase letters, even if they are defined upper case in
\proglang{FORTRAN}.

Also, there is no initialiser function here (\code{initfunc = NULL}).

\subsection{DAE models, integrator daspk}

\code{daspk} is one of the integrators in the package that solve DAE
models.  In order to be used with DASPK, DAEs are specified in implicit
form: \[0 = F(t, y, y', p)\] i.e. the DAE function (passed via argument
\code{res}) specifies the ``residuals'' rather than the derivatives (as for ODEs).

Consequently the DAE function specification in a compiled language is
also different.  For code written in \proglang{C}, the calling
sequence for \code{res} must be:

\begin{verbatim}
void myres(double *t, double *y, double *ydot, double *cj,
           double *delta, int *ires, double *yout, int *ip)
\end{verbatim}

where \code{*t} is the value of the independent variable, \code{*y}
points to a double precision vector that contains the current value of
the state variables, \code{*ydot} points to an array that will contain
the derivatives, \code{*delta} points to a vector that will contain the
calculated residuals. \code{*cj} points to a scalar, which is normally
proportional to the inverse of the stepsize, while \code{*ires} points
to an integer (not used). \code{*yout} points to any other output
variables (different from the state variables y), followed by the
double precision values as passed via argument \code{rpar}; finally
\code{*ip} is an integer vector containing at least 3 elements, its
first value (\code{*ip[0]}) equals the number of output variables,
calculated in the function (and which should be equal to \code{nout}),
its second element equals the total length of \code{*yout}, its third
element equals the total length of \code{*ip}, and finally come the
integer values as passed via argument \code{ipar}.

For code written in \proglang{FORTRAN}, the calling sequence for
\code{res} must be as in the following example:

\begin{verbatim}
      subroutine myresf(t, y, ydot, cj, delta, ires, out, ip)
        integer :: ires, ip(*)
        integer, parameter :: neq = 3
        double precision :: t, y(neq), ydot(neq), delta(neq), out(*)
        double precision :: K, ka, r, prod, ra, rb
        common /myparms/K,ka,r,prod

        if(ip(1) < 1) call rexit("nout should be at least 1")
        ra = ka* y(3)
        rb = ka/K *y(1) * y(2)

!! residuals of rates of changes
        delta(3) = -ydot(3) - ra + rb + prod
        delta(1) = -ydot(1) + ra - rb
        delta(2) = -ydot(2) + ra - rb - r*y(2)
        out(1) = y(1) + y(2) + y(3)
        return
        end
\end{verbatim}

Similarly as for the ODE model discussed above, the parameters are
kept in a common block which is initialised by an initialiser
subroutine:

\begin{verbatim}
      subroutine initpar(daspkparms)

       external daspkparms
       integer, parameter :: N = 4
       double precision parms(N)
       common /myparms/parms
       call daspkparms(N, parms)
       return
      end
\end{verbatim}

See the ODE example for how to initialise parameter values in
\proglang{C}.

Similarly, the function that specifies the Jacobian in a DAE differs
from the Jacobian when the model is an ODE. The DAE Jacobian is set
with argument \code{jacres} rather than \code{jacfunc} when an ODE.

For code written in \proglang{FORTRAN}, the \code{jacres} must be as:

\begin{verbatim}
       subroutine resjacfor (t, y, dy, pd, cj, out, ipar)

        integer, parameter :: neq = 3
        integer :: ipar(*)
        double precision :: K, ka, r, prod
        double precision :: pd(neq,neq),y(neq),dy(neq),out(*)
        common /myparms/K,ka,r,prod

!res1 = -dD - ka*D + ka/K *A*B + prod
         PD(1,1) = ka/K *y(2)
         PD(1,2) = ka/K *y(1)
         PD(1,3) = -ka -cj
!res2 = -dA + ka*D - ka/K *A*B
         PD(2,1) = -ka/K *y(2) -cj
         PD(2,2) = -ka/K *y(2)
         PD(2,3) = ka
!res3 = -dB + ka*D - ka/K *A*B - r*B
         PD(3,1) = -ka/K *y(2)
         PD(3,2) = -ka/K *y(2) -r -cj
         PD(3,3) = ka
         return
        end
\end{verbatim}

\subsection{DAE models, integrator radau}
Function \code{radau} solves DAEs in linearly implicit form, i.e. in the form
$M y' = f(t, y, p)$.

The derivative function $f$ is specified in the same way as for an ODE, i.e.
\begin{verbatim}
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip)
\end{verbatim}
and
\begin{verbatim}
      subroutine derivs (neq, t, y, ydot, out, IP)
\end{verbatim}
for \proglang{C} and \proglang{FORTRAN} code respectively.

To show how it should be used, we implement the caraxis problem as in
\citep{testset}. The implementation of this index 3 DAE, comprising 8
differential, and 2 algebraic equations in R is the last example of the
\code{radau} help page. We first repeat the R implementation:

<<>>=
caraxisfun <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    yb <- r * sin(w * t)
    xb <- sqrt(L * L - yb * yb)
    Ll <- sqrt(xl^2 + yl^2)
    Lr <- sqrt((xr - xb)^2 + (yr - yb)^2)

    dxl <- ul; dyl <- vl; dxr <- ur; dyr <- vr

    dul  <- (L0-Ll) * xl/Ll      + 2 * lam2 * (xl-xr) + lam1*xb
    dvl  <- (L0-Ll) * yl/Ll      + 2 * lam2 * (yl-yr) + lam1*yb - k * g

    dur  <- (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr)
    dvr  <- (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k * g

    c1   <- xb * xl + yb * yl
    c2   <- (xl - xr)^2 + (yl - yr)^2 - L * L

    list(c(dxl, dyl, dxr, dyr, dul, dvl, dur, dvr, c1, c2))
  })
}
eps <- 0.01; M <- 10; k <- M * eps^2/2;
L <- 1; L0 <- 0.5; r <- 0.1; w <- 10; g <- 1

parameter <- c(eps = eps, M = M, k = k, L = L, L0 = L0,
               r = r, w = w, g = g)

yini <- c(xl = 0, yl = L0, xr = L, yr = L0, ul = -L0/L, vl = 0,
          ur = -L0/L, vr = 0, lam1 = 0, lam2 = 0)

# the mass matrix
Mass      <- diag(nrow = 10, 1)
Mass[5,5] <- Mass[6,6] <- Mass[7,7] <- Mass[8,8] <- M * eps * eps/2
Mass[9,9] <- Mass[10,10] <- 0
Mass

# index of the variables: 4 of index 1, 4 of index 2, 2 of index 3
index <- c(4, 4, 2)

times <- seq(0, 3, by = 0.01)
out <- radau(y = yini, mass = Mass, times = times, func = caraxisfun,
        parms = parameter, nind = index)
@

<<label=caraxis,include=FALSE>>=
plot(out, which = 1:4, type = "l", lwd = 2)
@
\setkeys{Gin}{width=0.8\textwidth}

\begin{figure}
\begin{center}
<<label=figcaraxis,fig=TRUE,echo=FALSE, width=13, height=6>>=
<<caraxis>>
@
\end{center}
\caption{Solution of the caraxis model - see text for R-code}
\label{fig:caraxis}
\end{figure}
The implementation in \proglang{FORTRAN} consists of an initialiser function
and a derivative function.
\begin{verbatim}
c----------------------------------------------------------------
c Initialiser for parameter common block
c----------------------------------------------------------------
      subroutine initcaraxis(daeparms)

      external daeparms
      integer, parameter :: N = 8
      double precision parms(N)
      common /myparms/parms

      call daeparms(N, parms)
      return
      end

c----------------------------------------------------------------
c rate of change
c----------------------------------------------------------------
      subroutine caraxis(neq, t, y, ydot, out, ip)
      implicit none
      integer          neq, IP(*)
      double precision t, y(neq), ydot(neq), out(*)
      double precision eps, M, k, L, L0, r, w, g
      common /myparms/ eps, M, k, L, L0, r, w, g

      double precision xl, yl, xr, yr, ul, vl, ur, vr, lam1, lam2
      double precision yb, xb, Ll, Lr, dxl, dyl, dxr, dyr
      double precision dul, dvl, dur, dvr, c1, c2

c expand state variables
       xl = y(1)
       yl = y(2)
       xr = y(3)
       yr = y(4)
       ul = y(5)
       vl = y(6)
       ur = y(7)
       vr = y(8)
       lam1 = y(9)
       lam2 = y(10)

       yb = r * sin(w * t)
       xb = sqrt(L * L - yb * yb)
       Ll = sqrt(xl**2 + yl**2)
       Lr = sqrt((xr - xb)**2 + (yr - yb)**2)

       dxl = ul
       dyl = vl
       dxr = ur
       dyr = vr

       dul = (L0-Ll) * xl/Ll      + 2 * lam2 * (xl-xr) + lam1*xb
       dvl = (L0-Ll) * yl/Ll      + 2 * lam2 * (yl-yr) + lam1*yb - k*g
       dur = (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr)
       dvr = (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k*g

       c1  = xb * xl + yb * yl
       c2  = (xl - xr)**2 + (yl - yr)**2 - L * L

c function values in ydot
       ydot(1)  = dxl
       ydot(2)  = dyl
       ydot(3)  = dxr
       ydot(4)  = dyr
       ydot(5)  = dul
       ydot(6)  = dvl
       ydot(7)  = dur
       ydot(8)  = dvr
       ydot(9)  = c1
       ydot(10) = c2
      return
      end
\end{verbatim}

Assuming that the code is in file ``radaudae.f'', this model is compiled,
loaded and solved in R as:
\begin{verbatim}
system("R CMD SHLIB radaudae.f")
dyn.load(paste("radaudae", .Platform$dynlib.ext, sep = ""))

outDLL <- radau(y = yini, mass = Mass, times = times, func = "caraxis",
                initfunc = "initcaraxis", parms = parameter,
                dllname = "radaudae", nind = index)

dyn.unload(paste("radaudae", .Platform$dynlib.ext, sep = ""))
\end{verbatim}

\subsection{The root function from integrators lsodar and lsode}

\code{lsodar} is an extended version of integrator \code{lsoda} that
includes a root finding function. This function is spedified via
argument \code{rootfunc}. In \code{deSolve} version 1.7, \code{lsode}
has also been extended with root finding capabilities.

Here is how to program such a function in a lower-level language.  For
code written in \proglang{C}, the calling sequence for \code{rootfunc}
must be:

\begin{verbatim}
void myroot(int *neq, double *t, double *y, int *ng, double *gout,
            double *out, int *ip )
\end{verbatim}

where \code{*neq} and \code{*ng} are the number of state variables and
root functions respectively, \code{*t} is the value of the independent
variable, \code{y} points to a double precision array that contains
the current value of the state variables, and \code{gout} points to an
array that will contain the values of the constraint function whose
root is sought.  \code{*out} and \code{*ip} are a double precision and
integer vector respectively, as described in the ODE example above.

For code written in \proglang{FORTRAN}, the calling sequence for
\code{rootfunc} must be as in following example:

\begin{verbatim}
      subroutine myroot(neq, t, y, ng, gout, out, ip)
      integer :: neq, ng, ip(*)
      double precision :: t, y(neq), gout(ng), out(*)

      gout(1) = y(1) - 1.e-4
      gout(2) = y(3) - 1e-2

      return
      end
\end{verbatim}

\subsection{jacvec, the Jacobian vector for integrator lsodes}

Finally, in integration function \code{lsodes}, not the Jacobian
\emph{matrix} is specified, but a \emph{vector}, one for each column of the
Jacobian.  This function is specified via argument \code{jacvec}.

In \proglang{FORTRAN}, the calling sequence for \code{jacvec} is:

\begin{verbatim}
SUBROUTINE JAC (NEQ, T, Y, J, IAN, JAN, PDJ, OUT, IP)
DOUBLE PRECISION T, Y(*), IAN(*), JAN(*), PDJ(*), OUT(*)
INTEGER NEQ, J, IP(*)
\end{verbatim}


\subsection{Banded jacobians in compiled code}\label{band}

In the call of the jacobian function, the number of bands below and above the
diagonal (\code{ml, mu}) and the number of rows of the Jacobian matrix,
\code{nrowPD} is specified, e.g. for \proglang{FORTRAN} code:
\begin{verbatim}
      SUBROUTINE JAC (neq, T, Y, ml, mu, PD, nrowPD, RPAR, IPAR)
\end{verbatim}
The jacobian matrix to be returned should have dimension \code{nrowPD, neq}.

In case the Jacobian is banded, the size of \code{nrowPD} depends on the
  integrator. If the method is one of \code{lsode, lsoda, vode}, or related,
  then \code{nrowPD}  will be equal to \code{mu + 2 * ml + 1}, where the
  last ml rows should be filled with $0$s.

  For \code{radau}, \code{nrowpd} will be equal to \code{mu + ml + 1}

  Thus, it is important to write the FORTRAN or C-code in such a way that
  it can be used with both types of integrators - else it is likely that
  R will freeze if the wrong integrator is used.

We implement in FORTRAN, the example of the \code{lsode} help file.
The R-code reads:
<<>>=
## the model, 5 state variables
f1 <- function  (t, y, parms) {
  ydot <- vector(len = 5)

  ydot[1] <-  0.1*y[1] -0.2*y[2]
  ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3]
  ydot[3] <-           -0.3*y[2] +0.1*y[3] -0.2*y[4]
  ydot[4] <-                     -0.3*y[3] +0.1*y[4] -0.2*y[5]
  ydot[5] <-                               -0.3*y[4] +0.1*y[5]

  return(list(ydot))
}

## the Jacobian, written in banded form
bandjac <- function  (t, y, parms) {
  jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE,
                data = c( 0  , -0.2, -0.2, -0.2, -0.2,
                          0.1,  0.1,  0.1,  0.1,  0.1,
                         -0.3, -0.3, -0.3, -0.3,    0))
  return(jac)
}

## initial conditions and output times
yini  <- 1:5
times <- 1:20

## stiff method, user-generated banded Jacobian
out  <- lsode(yini, times, f1, parms = 0, jactype = "bandusr",
              jacfunc = bandjac, bandup = 1, banddown = 1)
@

In FORTRAN, the code might look like this:
\begin{verbatim}
c Rate of change
      subroutine derivsband (neq, t, y, ydot,out,IP)
      integer           neq, IP(*)
      DOUBLE PRECISION  T, Y(5), YDOT(5), out(*)

        ydot(1) = 0.1*y(1) -0.2*y(2)
        ydot(2) = -0.3*y(1) +0.1*y(2) -0.2*y(3)
        ydot(3) =           -0.3*y(2) +0.1*y(3) -0.2*y(4)
        ydot(4) =                     -0.3*y(3) +0.1*y(4) -0.2*y(5)
        ydot(5) =                               -0.3*y(4) +0.1*y(5)
      RETURN
      END

c The banded jacobian
      subroutine jacband (neq, t, y, ml, mu, pd, nrowpd, RP, IP)
      INTEGER           neq, ml, mu, nrowpd, ip(*)
      DOUBLE PRECISION  T, Y(5), PD(nrowpd,5), rp(*)

        PD(:,:) = 0.D0

        PD(1,1) =  0.D0
        PD(1,2) = -.02D0
        PD(1,3) = -.02D0
        PD(1,4) = -.02D0
        PD(1,5) = -.02D0

        PD(2,:) = 0.1D0

        PD(3,1) = -0.3D0
        PD(3,2) = -0.3D0
        PD(3,3) = -0.3D0
        PD(3,4) = -0.3D0
        PD(3,5) = 0.D0
      RETURN
      END
\end{verbatim}
Assuming that this code is in file \code{"odeband.f"}, we compile from within
R and load the shared library (assuming the working directory holds the
source file) with:
\begin{verbatim}
system("R CMD SHLIB odeband.f")
dyn.load(paste("odeband", .Platform$dynlib.ext, sep = ""))
\end{verbatim}
To solve this problem, we write in R
\begin{verbatim}
out2 <- lsode(yini, times, "derivsband", parms = 0, jactype = "bandusr",
      jacfunc = "jacband", bandup = 1, banddown = 1, dllname = "odeband")

out2 <- radau(yini, times, "derivsband", parms = 0, jactype = "bandusr",
      jacfunc = "jacband", bandup = 1, banddown = 1, dllname = "odeband")
\end{verbatim}

This will work both for the \code{lsode} family as for \code{radau}.
In the first case, when entering subroutine \code{jacband}, \code{nrowpd} will
have the value $5$, in the second case, it will be equal to $4$.
\section{Testing functions written in compiled code}

Two utilities have been included to test the function implementation
in compiled code:

\begin{itemize}
\item \code{DLLfunc} to test the implementation of the derivative
  function as used in ODEs. This function returns the derivative
  $\frac{dy}{dt}$ and the output variables.
\item \code{DLLres} to test the implementation of the residual function
  as used in DAEs. This function returns the residual function
  $\frac{dy}{dt}-f(y,t)$ and the output variables.
\end{itemize}

These functions serve no other purpose than to test whether the
compiled code returns what it should.

\subsection{DLLfunc}

We test whether the ccl4 model, which is part of \code{deSolve}
package, returns the proper rates of changes.  (Note: see
\code{example(ccl4model)} for a more comprehensive implementation)

<<>>=
## Parameter values and initial conditions
Parms <- c(0.182, 4.0, 4.0, 0.08, 0.04, 0.74, 0.05, 0.15, 0.32,
        16.17, 281.48, 13.3, 16.17, 5.487, 153.8, 0.04321671,
        0.4027255, 1000, 0.02, 1.0, 3.8)

yini <- c( AI=21, AAM=0, AT=0, AF=0, AL=0, CLT=0,  AM=0  )

## the rate of change
DLLfunc(y = yini, dllname = "deSolve", func = "derivsccl4",
        initfunc = "initccl4", parms = Parms, times = 1,
        nout = 3, outnames = c("DOSE", "MASS", "CP")  )
@

\subsection{DLLres}

The deSolve package contains a FORTRAN implementation of the chemical
model described above (section 4.1), where the production rate is
included as a forcing function (see next section).

Here we use \code{DLLres} to test it:

<<>>=
pars <- c(K = 1, ka   = 1e6, r    = 1)

## Initial conc; D is in equilibrium with A,B
y     <- c(A = 2, B = 3, D = 2*3/pars["K"])

## Initial rate of change
dy  <- c(dA = 0, dB = 0, dD = 0)

## production increases with time
prod <- matrix(nc=2,data=c(seq(0,100,by=10),seq(0.1,0.5,len=11)))

DLLres(y=y,dy=dy,times=5,res="chemres",
          dllname="deSolve", initfunc="initparms",
          initforc="initforcs", parms=pars, forcings=prod,
          nout=2, outnames=c("CONC","Prod"))

@

\section{Using forcing functions}

Forcing functions in DLLs are implemented in a similar way as
parameters.  This means:

\begin{itemize}
\item They are initialised by means of an initialiser function. Its
  name should be passed to the solver via argument \code{initforc}.

  Similar as the parameter initialiser function, the function denoted
  by \code{initforc} has as its sole argument a pointer to the vector
  that contains the forcing funcion values in the compiled code. In
  case of \proglang{C} code, this will be a global vector; in case of
  \proglang{FORTRAN}, this will be a vector in a common block.

  The solver puts a pointer to this vector and updates the forcing
  functions in this memory area at each time step. Hence, within the
  compiled code, forcing functions can be assessed as if they are
  parameters (although, in contrast to the latter, their values will
  generally change).  No need to update the values for the current
  time step; this has been done before entering the \code{derivs}
  function.
\item The forcing function data series are passed to the integrator,
  via argument \code{forcings}; if there is only one forcing function
  data set, then a 2-columned matrix (time, value) will do; else the
  data should be passed as a list, containing (time, value) matrices
  with the individual forcing function data sets. Note that the data
  sets in this list should be \emph{in the same ordering} as the
  declaration of the forcings in the compiled code.
\end{itemize}

A number of options allow to finetune certain settings. They are in a
list called \code{fcontrol} which can be supplied as argument when
calling the solvers. The options are similar to the arguments from R
function \code{approx}, howevers the default settings are often
different.

The following options can be specified:

\begin{itemize}
\item \code{method} specifies the interpolation method to be used.
  Choices are "linear" or "constant", the default is "linear", which
  means linear interpolation (same as \code{approx})
\item \code{rule}, an integer describing how interpolation is to take
  place \emph{outside} the interval [min(times), max(times)].  If
  \code{rule} is \code{1} then an error will be triggered and the
  calculation will stop if extrapolation is necessary.  If it is
  \code{2}, the default, the value at the closest data extreme is
  used, a warning will be printed if \code{verbose} is TRUE.

  Note that the default differs from the \code{approx} default.
\item \code{f}, for method=\code{"constant"} is a number between
  \code{0} and \code{1} inclusive, indicating a compromise between
  left- and right-continuous step functions. If \code{y0} and
  \code{y1} are the values to the left and right of the point then the
  value is \code{y0*(1-f)+y1*f} so that \code{f=0} is right-continuous
  and \code{f=1} is left-continuous. The default is to have
  \code{f=0}.  For some data sets it may be more realistic to set
  \code{f=0.5}.
\item \code{ties}, the handling of tied \code{times} values. Either a
  function with a single vector argument returning a single number
  result or the string "ordered".

  Note that the default is "ordered", hence the existence of ties will
  NOT be investigated; in practice this means that, if ties exist, the
  first value will be used; if the dataset is not ordered, then
  nonsense will be produced.

  Alternative values for \code{ties} are \code{mean}, \code{min}
  etc...  which will average, or take the minimal value if multiple
  values exist at one time level.
\end{itemize}

The default settings of \code{fcontrol} are:

\code{fcontrol=list(method="linear", rule = 2, f = 0, ties =
  "ordered")}

Note that only ONE specification is allowed, even if there is more
than one forcing function data set. (may/should change in the future).

\subsection{A simple FORTRAN example}

We implement the example from chapter 3 of the book \citep{Soetaert08}
in FORTRAN.

This model describes the oxygen consumption of a (marine) sediment in
response to deposition of organic matter (the forcing function). One
state variable, the organic matter content in the sediment is modeled;
it changes as a function of the deposition \code{Flux} (forcing) and
organic matter decay (first-order decay rate \code{k}).

\[
 \frac{dC}{dt}=Flux_t-k \cdot C
\]

with initial condition $C(t=0)=C_0$; the latter is estimated as the
mean of the flux divided by the decay rate.

The FORTRAN code looks like this:

\begin{verbatim}
c Initialiser for parameter common block
      subroutine scocpar(odeparms)

      external odeparms
      integer N
      double precision parms(2)
      common /myparms/parms

      N = 1
      call odeparms(N, parms)
      return
      end

c Initialiser for forcing common block
      subroutine scocforc(odeforcs)

      external odeforcs
      integer N
      double precision forcs(1)
      common /myforcs/forcs

      N = 1
      call odeforcs(N, forcs)
      return
      end

c Rate of change and output variables

      subroutine scocder (neq, t, y, ydot,out,IP)

      integer          neq, IP(*)
      double precision t, y(neq), ydot(neq), out(*), k, depo
      common /myparms/k
      common /myforcs/depo

      if(IP(1) < 2) call rexit("nout should be at least 2")

      ydot(1) = -k*y(1) + depo

      out(1)= k*y(1)
      out(2)= depo

      return
      end

\end{verbatim}

Here the subroutine \code{scocpar} is business as usual; it
initialises the parameter common block (there is only one parameter).
Subroutine \code{odeforcs} does the same for the forcing function,
which is also positioned in a common block, called \code{myforcs}.
This common block is made available in the derivative subroutine (here
called \code{scocder}), where the forcing function is named
\code{depo}.

At each time step, the integrator updates the value of this forcing
function to the correct time point. In this way, the forcing functions
can be used as if they are (time-varying) parameters.  All that's left
to do is to pass the forcing function data set and the name of the
forcing function initialiser routine. This is how to do it in R.

First the data are inputted:

<<>>=
Flux <- matrix(ncol=2,byrow=TRUE,data=c(
  1, 0.654, 11, 0.167,   21, 0.060, 41, 0.070, 73,0.277, 83,0.186,
  93,0.140,103, 0.255,  113, 0.231,123, 0.309,133,1.127,143,1.923,
  153,1.091,163,1.001,  173, 1.691,183, 1.404,194,1.226,204,0.767,
  214, 0.893,224,0.737, 234,0.772,244, 0.726,254,0.624,264,0.439,
  274,0.168,284 ,0.280, 294,0.202,304, 0.193,315,0.286,325,0.599,
  335, 1.889,345, 0.996,355,0.681,365,1.135))
head(Flux)
@

and the parameter given a value (there is only one)

<<>>=
parms <- 0.01
@

The initial condition \code{Yini} is estimated as the annual mean of
the Flux and divided by the decay rate (parameter).

<<>>=
meanDepo <- mean(approx(Flux[,1],Flux[,2], xout=seq(1,365,by=1))$y)
Yini <- c(y=meanDepo/parms)
@

After defining the output times, the model is run, using integration routine
\code{ode}.

The \emph{name} of the derivate function \code{"scocder"}, of the dll
\code{"deSolve"}\footnote{this example is made part of the deSolve
  package, hence the name of the dll is "deSolve"} and of the
initialiser function \code{"scocpar"} are passed, as in previous
examples.

In addition, the forcing function data set is also passed
(\code{forcings=Flux}) as is the name of the forcing initialisation
function (\code{initforc="scocforc"}).

<<>>=
times <- 1:365
out <- ode(y=Yini, times, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc="scocforc",  forcings=Flux,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation","Depo"))
head(out)
@

Now, the way the forcing functions are interpolated are changed:
Rather than linear interpolation, constant (block, step) interpolation is used.

<<>>=
fcontrol <- list(method="constant")
out2 <- ode(y=Yini, times, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc="scocforc",  forcings=Flux, fcontrol=fcontrol,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation","Depo"))
@

Finally, the results are plotted:

<<label=scoc,include=FALSE>>=
par (mfrow=c(1,2))
plot(out, which = "Depo", col="red",
  xlab="days", ylab="mmol C/m2/ d", main="method='linear'")
lines(out[,"time"], out[,"Mineralisation"], lwd=2, col="blue")
legend("topleft",lwd=1:2,col=c("red","blue"), c("Flux","Mineralisation"))

plot(out, which = "Depo", col="red",
  xlab="days", ylab="mmol C/m2/ d", main="method='constant'")
lines(out2[,"time"], out2[,"Mineralisation"], lwd=2, col="blue")
@
\setkeys{Gin}{width=0.8\textwidth}

\begin{figure}
\begin{center}
<<label=figscoc,fig=TRUE,echo=FALSE, width=13, height=6>>=
<<scoc>>
@
\end{center}
\caption{Solution of the SCOC model, implemented in compiled code, and
including a forcing function - see text for R-code}
\label{fig:scoc}
\end{figure}

\subsection{An example in C}

Consider the following R-code which implements a
resource-producer-consumer Lotka-Volterra type of model in R (it is a
modified version of the example of function \code{ode}):

<<>>=
SPCmod <- function(t, x, parms, input)  {
  with(as.list(c(parms, x)), {
    import <- input(t)
    dS <- import - b*S*P + g*C    # substrate
    dP <- c*S*P  - d*C*P          # producer
    dC <- e*P*C  - f*C            # consumer
    res <- c(dS, dP, dC)
    list(res, signal = import)
  })
}

## The parameters
parms  <- c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)

## vector of timesteps
times  <- seq(0, 100, by=0.1)

## external signal with several rectangle impulses
signal <- as.data.frame(list(times = times,
                 import = rep(0, length(times))))

signal$import <- ifelse((trunc(signal$times) %% 2 == 0), 0, 1)
sigimp <- approxfun(signal$times, signal$import, rule = 2)

## Start values for steady state
xstart <- c(S = 1, P = 1, C = 1)

## Solve model
print (system.time(
out <- ode(y = xstart,times = times,
           func = SPCmod, parms, input = sigimp)
))
@

All output is printed at once:

<<label=lv,include=FALSE>>=
plot(out)
@

\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
<<label=figlv,fig=TRUE,echo=FALSE>>=
<<lv>>
@
\end{center}
\caption{Solution of the Lotka-Volterra resource (S)-producer (P) -
consumer (C) model with time-variable input (signal) - see text for R-code}
\label{fig:lv}
\end{figure}

The C-code, in file \url{Forcing\_lv.c}, can be found in the packages
\url{/doc/examples/dynload} subdirectory\footnote{this can be opened by typing
\code{browseURL(paste(system.file(package = "deSolve"), "/doc/examples/dynload", sep = ""))}}.
It can be compiled, from within R by

\begin{verbatim}
system("R CMD SHLIB Forcing_lv.c")
\end{verbatim}

After defining the parameter and forcing vectors, and giving them
comprehensible names, the parameter and forcing initialiser functions
are defined (\code{parmsc} and \code{forcc} respectively). Next is the
derivative function, \code{derivsc}.

\begin{verbatim}
#include <R.h>

static double parms[6];
static double forc[1];

/* A trick to keep up with the parameters and forcings */
#define b parms[0]
#define c parms[1]
#define d parms[2]
#define e parms[3]
#define f parms[4]
#define g parms[5]

#define import forc[0]

/* initializers: */
void odec(void (* odeparms)(int *, double *))
{
    int N=6;
    odeparms(&N, parms);
}

void forcc(void (* odeforcs)(int *, double *))
{
    int N=1;
    odeforcs(&N, forc);
}

/* derivative function */
void derivsc(int *neq, double *t, double *y, double *ydot,
  double *yout, int*ip)
{
    if (ip[0] <2) error("nout should be at least 2");
    ydot[0] = import - b*y[0]*y[1] + g*y[2];
    ydot[1] =          c*y[0]*y[1] - d*y[2]*y[1];
    ydot[2] =          e*y[1]*y[2] - f*y[2];

    yout[0] = y[0] + y[1] + y[2];
    yout[1] = import;
}

\end{verbatim}

After defining the forcing function time series, which is to be
interpolated by the integration routine, and loading the DLL, the
model is run:

\begin{verbatim}
Sigimp <- approx(signal$times, signal$import, xout=ftime,rule = 2)$y
forcings <- cbind(ftime,Sigimp)

dyn.load("Forcing_lv.dll")
out <- ode(y=xstart, times, func = "derivsc",
   parms = parms, dllname = "Forcing_lv",initforc = "forcc",
   forcings=forcings, initfunc = "parmsc", nout = 2,
   outnames = c("Sum","signal"), method = rkMethod("rk34f"))
dyn.unload("Forcing_lv.dll")
\end{verbatim}

This code executes about 30 times faster than the \proglang{R}-code.

With a longer simulation time, the difference becomes more pronounced,
e.g. with times till 800 days, the DLL code executes 200 times faster%
\footnote{this is due to the sequential update of the forcing
  functions by the solvers, compared to the bisectioning approach used
  by approxfun}.

\section{Implementing events in compiled code}
  An \code{event} occurs when the value of a state variable is suddenly changed,
  e.g. a certain amount is added, or part is removed. The integration
  routines cannot deal easily with such state variable changes. Typically these
  events occur only at specific times.
  In \code{deSolve}, events can be imposed by means of an input file that
  specifies at which time a certain state variable is altered, or via
  an event function.

  Both types of events combine with compiled code.

  Take the previous example, the Lotka-Volterra SPC model. Suppose that every
  10 days, half of the consumer is removed.

  We first implement these events as a \code{data.frame}

<<>>=
eventdata <- data.frame(var=rep("C",10),time=seq(10,100,10),value=rep(0.5,10),
  method=rep("multiply",10))
eventdata
@

  This model is solved, and plotted as:
\begin{verbatim}
dyn.load("Forcing_lv.dll")
  out2 <-  ode(y = y, times, func = "derivsc",
            parms = parms, dllname = "Forcing_lv", initforc="forcc",
            forcings = forcings, initfunc = "parmsc", nout = 2,
            outnames = c("Sum", "signal"), events=list(data=eventdata))
dyn.unload("Forcing_lv.dll")
plot(out2, which = c("S","P","C"), type = "l")

\end{verbatim}

The event can also be implemented in \proglang{C} as:
\begin{verbatim}
void event(int *n, double *t, double *y) {
   y[2] = y[2]*0.5;
}
\end{verbatim}
Here n is the length of the state variable vector \code{y}.
and is then solved as:
\begin{verbatim}
dyn.load("Forcing_lv.dll")
  out3 <-  ode(y = y, times, func = "derivsc",
            parms = parms, dllname = "Forcing_lv", initforc="forcc",
            forcings = forcings, initfunc = "parmsc", nout = 2,
            outnames = c("Sum", "signal"),
            events = list(func="event",time=seq(10,90,10)))
dyn.unload("Forcing_lv.dll")
\end{verbatim}

\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
\includegraphics{comp-event}
\end{center}
\caption{Solution of the Lotka-Volterra resource (S)~-- producer (P)~--
consumer (C) model with time-variable input (signal) and with half of
the consumer removed every 10 days - see text for R-code}
\label{fig:lv2}
\end{figure}

\section{Delay differential equations}

It is now also very simple to implement delay 
differential equations in compiled code and solve them with \code{dede}.
In order to do so, you need to get access to the R-functions \code{lagvalue} and \code{lagderiv}
that will give you the past value of the state variable or its derivative 
respectively.

\subsection{Delays implemented in Fortran}

If you use \proglang{Fortran}, then the easiest way is to link your code with a file called 
\code{dedeUtils.c} that you will find in the packages subdirectory \code{inst/doc/dynload-dede}. 
This file contains Fortran-callable interfaces to the delay-differential utility
functions from package \pkg{deSolve}, and that are written in \proglang{C}.
Its content is: 
\begin{verbatim}
void F77_SUB(lagvalue)(double *T, int *nr, int *N, double *ytau) {
  static void(*fun)(double, int*, int, double*) = NULL;
  if (fun == NULL)
    fun =  (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagvalue");
  fun(*T, nr, *N, ytau);
  return;
}

void F77_SUB(lagderiv)(double *T, int *nr, int *N, double *ytau) {
  static void(*fun)(double, int*, int, double*) = NULL;
  if (fun == NULL)
    fun =  (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagderiv");
  fun(*T, nr, *N, ytau);
  return;
}
\end{verbatim}

Here \code{T} is the time at which the value needs to be retrieved, 
\code{nr} is an integer that defines the number of the state variable or its
derivative whose delay we want, \code{N} is the total number of state variabes 
and \code{ytau} will have the result.

We start with an example, a Lotka-Volterra system with delay, that we will implement
 in \proglang{Fortran} (you will find this example in the package directory
\code{inst/doc/dynload-dede}, in file \code{dede_lvF.f}

The R-code would be:
<<>>=
derivs <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    if (t < tau)
      ytau <- c(1, 1)
    else
      ytau <- lagvalue(t - tau, c(1, 2))

    dN <- f * N - g * N * P
    dP <- e * g * ytau[1] * ytau[2] - m * P
    list(c(dN, dP), tau=ytau[1], tau=ytau[2])
  })
}

yinit <- c(N = 1, P = 1)
times <- seq(0, 500)
parms <- c(f = 0.1, g = 0.2, e = 0.1, m = 0.1, tau = .2)
yout <- dede(y = yinit, times = times, func = derivs, parms = parms)
head(yout)
@

In Fortran the code looks like this:
\begin{verbatim}
!  file dede_lfF.f
!  Initializer for parameter common block
      subroutine initmod(odeparms)
      external odeparms
      double precision parms(5)
      common /myparms/parms

      call odeparms(5, parms)
      return
      end

!  Derivatives and one output variable
      subroutine derivs(neq, t, y, ydot, yout, ip)
      integer neq, ip(*)
      double precision t, y(neq), ydot(neq), yout(*)

      double precision N, P, ytau(2), tlag
      integer nr(2)

      double precision f, g, e, m, tau
      common /myparms/f, g, e, m, tau

      if (ip(1) < 2) call rexit("nout should be at least 2")

      N = y(1)
      P = y(2)
      nr(1) = 0
      nr(2) = 1
      ytau(1) = 1.0
      ytau(2) = 1.0
 
      tlag = t - tau
      if (tlag .GT. 0.0) call lagvalue(tlag, nr, 2, ytau)

      ydot(1) = f * N - g * N * P
      ydot(2) = e * g * ytau(1) * ytau(2) - m * P

      yout(1) = ytau(1)
      yout(2) = ytau(2)

      return
      end
\end{verbatim}

During compilation, we need to also compile the file \code{dedeUtils.c}. Assuming that the
above \proglang{Fortran} code is in file \code{dede_lvF.f}, which is found in the working directory that 
also contains file \code{dedeUtils.c}, the problem is compiled and run as:
\begin{verbatim}
system("R CMD SHLIB dede_lvF.f dedeUtils.c")
dyn.load(paste("dede_lvF", .Platform$dynlib.ext, sep=""))
yout3 <- dede(yinit, times = times, func = "derivs", parms = parms,
    dllname = "dede_lvF", initfunc = "initmod", nout = 2)
\end{verbatim}

\subsection{Delays implemented in C}

We now give the same example in \proglang{C}-code (you will find this in directory \code{inst/doc/dynload-dede/dede_lv.c}).

\begin{verbatim}
#include <R.h>
#include <Rinternals.h>
#include <Rdefines.h>
#include <R_ext/Rdynload.h>

static double parms[5];
#define f parms[0]
#define g parms[1]
#define e parms[2]
#define m parms[3]
#define tau parms[4]

/* Interface to dede utility functions in package deSolve */
void lagvalue(double T, int *nr, int N, double *ytau) {
  static void(*fun)(double, int*, int, double*) = NULL;
  if (fun == NULL)
    fun =  (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagvalue");
  return fun(T, nr, N, ytau);
}

void lagderiv(double T, int *nr, int N, double *ytau) {
  static void(*fun)(double, int*, int, double*) = NULL;
  if (fun == NULL)
    fun =  (void(*)(double, int*, int, double*))R_GetCCallable("deSolve", "lagderiv");
  return fun(T, nr, N, ytau);
}

/* Initializer  */
void initmod(void (* odeparms)(int *, double *)) {
  int N = 5;
  odeparms(&N, parms);
}

/* Derivatives */
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip) {

  if (ip[0] < 2) error("nout should be at least 1");

  double N = y[0];
  double P = y[1];

  int Nout  = 2;                  // number of returned lags  ( <= n_eq !!)
  int nr[2] = {0, 1};             // which lags are needed? 
                                  // numbering starts from zero !
  double ytau[2] = {1.0, 1.0};    // array; initialize with default values !

  double T = *t - tau;
  
  if (*t > tau) {
    lagvalue(T, nr, Nout, ytau);
  }

  ydot[0] = f * N - g * N * P;
  ydot[1] = e * g * ytau[0] * ytau[1] - m * P;

  yout[0] = ytau[0];
  yout[1] = ytau[1];

}
\end{verbatim}

Assuming this code is in a file called \code{dede_lv.c}, which is in the working
directory, this file is then compiled and run as:
\begin{verbatim}
system("R CMD SHLIB dede_lv.c")
dyn.load(paste("dede_lv", .Platform$dynlib.ext, sep=""))
yout2 <- dede(yinit, times = times, func = "derivs", parms = parms,
    dllname = "dede_lv", initfunc = "initmod", nout = 2)
dyn.unload(paste("dede_lv", .Platform$dynlib.ext, sep=""))
\end{verbatim}


\section{Difference equations in compiled code}
There is one special-purpose solver, triggered with \code{method = "iteration"}
which can be used in cases where the new values of the state variables are
estimated by the user, and need not be found by integration.

This is for instance useful when the model consists of difference equations,
or for 1-D models when transport is implemented by an implicit or a
semi-implicit method.

An example of a discrete time model, represented by a difference
equation is given in the help file of solver \code{ode}. It consists of the
host-parasitoid model described as from \citet[p283]{Soetaert08}.

We first give the R-code, and how it is solved:
\begin{verbatim}
Parasite <- function (t, y, ks) {
  P <- y[1]
  H <- y[2]
  f    <- A * P / (ks +H)
  Pnew <- H* (1-exp(-f))
  Hnew <- H * exp(rH*(1.-H) - f)

  list (c(Pnew, Hnew))
}
rH <- 2.82  # rate of increase
A  <- 100   # attack rate
ks <- 15.   # half-saturation density

out <- ode (func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50,
            parms = ks, method = "iteration")
\end{verbatim}
Note that the function returns the updated value of the state variables rather
than the rate of change (derivative). The method ``iteration'' does not perform
any integration.

The implementation in \proglang{FORTRAN} consists of an initialisation
function to pass the parameter values (\code{initparms}) and the "update"
function that returns the new values of the state variables (\code{parasite}):
\begin{verbatim}
      subroutine initparms(odeparms)
        external odeparms
        double precision parms(3)
        common /myparms/parms
        call odeparms(3, parms)
        return
      end

      subroutine parasite (neq, t, y, ynew, out, iout)
        integer          neq, iout(*)
        double precision t, y(neq), ynew(neq), out(*), rH, A, ks
        common /myparms/ rH, A, ks
	      double precision P, H, f

    	  P = y(1)
    	  H = y(2)
        f = A * P / (ks + H)

        ynew(1) = H * (1.d0 - exp(-f))
        ynew(2) = H * exp (rH * (1.d0 - H) - f)
        return

      end
\end{verbatim}

The model is compiled, loaded and executed in R as:
\begin{verbatim}
system("R CMD SHLIB difference.f")
dyn.load(paste("difference", .Platform$dynlib.ext, sep = ""))

require(deSolve)
rH <- 2.82 # rate of increase
A  <- 100  # attack rate
ks <- 15.  # half-saturation density

parms <- c(rH = rH, A = A, ks = ks)
out <- ode (func = "parasite", y = c(P = 0.5, H = 0.5), times = 0:50,
            initfunc = "initparms", dllname = "difference", parms = parms,
            method = "iteration")
\end{verbatim}

\section{Final remark}

Detailed information about communication between \proglang{C},
\proglang{FORTRAN} and \proglang{R} can be found in \citet{Rexts2009}.

Notwithstanding the speed gain when using compiled code, one should
not carelessly decide to always resort to this type of modelling.

Because the code needs to be formally compiled and linked to
\proglang{R} much of the elegance when using pure \proglang{R} models
is lost. Moreover, mistakes are easily made and paid harder in
compiled code: often a programming error will terminate
\proglang{R}. In addition, these errors may not be simple to trace.

\clearpage
%% this adds References to the PDF-Index without adding an obsolete section
\phantomsection
\addcontentsline{toc}{section}{References}
\bibliography{integration}

\end{document}