[R-sig-dyn-mod] Possible bug specifying root functions for lsodes

[McGrath, Justin M]

There doesn't seem to be any developers' Web site or a place to submit
bugs, so I thought I'd try here.

I am using deSolve version 1.10-4, R version 2.15.2. I am unable to
figure out how to use root functions with the lsodes solver.
Regardless of what the function returns, I always receive the
following error:

DLSODES- MF (=I1) illegal.
In above message, I =
[1] 1413567809
Error in lsodes(func = LVmod0D, y = yini, parms = pars, times = times,  :
  illegal input detected before taking any integration steps - see
written message

This happens with models I specify myself and also with the examples
given the in the vignette.

LVmod0D <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
IngestC <- rI * P * C
GrowthP <- rG * P * (1 - P/K)
MortC <- rM * C
dP <- GrowthP - IngestC
dC <- IngestC * AE - MortC
return(list(c(dP, dC)))
})}

pars <- c(rI = 0.2, rG = 1.0, rM = 0.2, AE = 0.5, K = 10)
yini <- c(P = 1, C = 2)
times <- seq(0, 200, by = 1)
print(system.time(out <- ode(func = LVmod0D, y = yini, parms = pars,
times = times)))

rootfun <- function(Time, State, Pars) {
dstate <- unlist(LVmod0D(Time, State, Pars))
sum(abs(dstate)) - 1e-4
}

# These produce the same output. Presumably because calling lsoda with
rootfun just calls lsodar eventually anyway.
outar <- lsodar(func = LVmod0D, y = yini, parms = pars, times = times,
rootfun = rootfun)
outa <- lsoda(func = LVmod0D, y = yini, parms = pars, times = times,
rootfun = rootfun)

# Using lsodes without a root function works as expected.
outs <- lsodes(func = LVmod0D, y = yini, parms = pars, times = times)

# However, adding a root function breaks it.
outsr <- lsodes(func = LVmod0D, y = yini, parms = pars, times = times,
rootfun = rootfun)

# If you use a constant function, it still doesn't work.
rootfun <- function(Time, State, Pars) {
return(1)
}
outsr2 <- lsodes(func = LVmod0D, y = yini, parms = pars, times =
times, rootfun = rootfun)

I am not familiar with the code behind deSolve, but it seems strange
that "rwork" appears twice in the call to dlsodesr in call_lsoda.c
between lines 472 and 478 (copied below). However, the same thing
occurs in the call to dlsodes, and that seems to work fine. In any
case, I cannot figure out why I can't use root functions with lsodes.

     } else if (solver == 7) {
        F77_CALL(dlsodesr) (deriv_func, &n_eq, xytmp, &tin, &tout,
               &itol, Rtol, Atol, &itask, &istate, &iopt, rwork,
               &lrw, iwork, &liw, rwork, jac_vec, &jt, root_func,
&nroot, jroot, /*rwork: iwk in fortran*/
        out, ipar);
 lyh = iwork[21];
      }

If this is a bug, could it please be fixed? If it is not a bug, could
you explain what I am doing wrong with lsodes?

Best wishes,
Justin