R deSolve with DLL does not handle the ip parameter

Question

In the R Package deSolve, Writing Code in Compiled Languages there is a prototype and example

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

The documentation says:

*ip points to an integer vector whose length is at least 3; the first element (ip[0]) contains the number of output values (which should be equal or larger than nout), its second element contains the length of *yout, and the third element contains the length of *ip; next are integer values, as passed by parameter ipar when calling the integrator.

But in fact there is no way to pass *ip to your C or Fortran derivs subroutine. To get code to run you have to call ode() with a named parameter nout because the routine that checks and loads the DLL requires it. However, the value passed for nout becomes ip(1) in the Fortran derivs subroutine.

It is not clear to me what is going on. The output from the example in the vignette is a matrix dim(13,7) while the values in ip() in the derivs() subroutine is always ip(1) = nout, ip(2) = nout+2, ip(3) = 4. The output is exactly the same independent of the value of nout.

If have the deSolve source code and grep searches for ip being declared or set yield NULL.

Does ip have any role in the actual calculation?

Answer 1

The *ip argument is passed to the derivs function by the solver. The solver is initialized by the initializer, e.g. N in the odec function below that calls odeparms with the number of parameters.

The number of external outputs *ip is initialized from R by the nout argument:

lsoda(y, times, func, parms, rtol = 1e-6, atol = 1e-6,
  jacfunc = NULL, jactype = "fullint", rootfunc = NULL,
  verbose = FALSE, nroot = 0, tcrit = NULL,
  hmin = 0, hmax = NULL, hini = 0, ynames = TRUE,
  maxordn = 12, maxords = 5, bandup = NULL, banddown = NULL,
  maxsteps = 5000, dllname = NULL, initfunc = dllname,
  initpar = parms, rpar = NULL, ipar = NULL, nout = 0,
  outnames = NULL, forcings = NULL, initforc = NULL,
  fcontrol = NULL, events = NULL, lags = NULL,...)

Here a C example using yout. The full source code of the example can be found in files odec.c and odedll.R on R-Forge.

#include <R.h> /* gives F77_CALL through R_ext/RS.h */

static double parms[3];

#define k1 parms[0]
#define k2 parms[1]
#define k3 parms[2]

/* initializer: same name as the dll (without extension) */
void odec(void (* odeparms)(int *, double *)) {
    int N=3;
    odeparms(&N, parms);
}

/* Derivatives */
void derivsc(int *neq, double *t, double *y, double *ydot, double *yout, int*ip) {
    if (ip[0] <3) error("nout should be at least 3");
    ydot[0] = -k1*y[0] + k2*y[1]*y[2];
    ydot[2] = k3 * y[1]*y[1];
    ydot[1] = -ydot[0]-ydot[2];
    
    yout[0] = y[0]+y[1]+y[2];
    yout[1] = y[0]*2;
    yout[2] = k3;
}