[R-sig-dyn-mod] lsodes with jacobian

[Karline Soetaert]

Hello Daniel,

Sorry for answering this late; I was at the useR conference.

This is a very technical setting that is determined by the original source code. 

I think specifiying inz is enough. You have to set sparsetype = "sparsejan", and then specify the vector inz such that it has the necessary information about where to find the nonzero elements in the jacobian. It is to contain two vectors in the original FORTRAN code of lsodes: ian and jan. How to do this is explained in the "details" section of lsodes ( use ?lsodes).

As your model with the explicit jacobian runs longer than in the other case, this probably means that the jacobian itself is not correctly specified. 
One way to see this is to ask for the diagnostics of both runs. This will tell you the number of steps it took; it should be comparable.

For large problems, the explicit specification should be faster.

Hope this helps,


Karline

-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Daniel Kaschek
Sent: zaterdag 28 juni 2014 16:31
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] lsodes with jacobian

Hi again,

by looking at call_losda.c in the source code of the deSolve package, I found the necessary call. You can find my fist example at the end of the mail.

Next, I tried this out for an ODE with 600 equations. Unfortunately, it takes 60 times longer with jacobian than without :-( I suppose, this has to do with the sparsity specifications. Or is it the "if"-lines in the jacobian? 

* What would be the correct specification of sparsetype, nnz and inz in case of an explicitly given jacvec function?
* Knowing the zeros of my jacobian, does it make sens to write down the jacvec function explicitly or is it sufficient to specify nnz and inz?

Cheers,
Daniel

#include <R.h>
 #include <math.h> 

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

#define rGrow parms[0]
 #define K parms[1]
 #define rIng parms[2]
 #define rMort parms[3]
 #define y0_0 parms[4]
 #define y1_0 parms[5]
 #define assEff forc[0]
 #define assEffDot forc[1] 

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

void initforc(void (* odeforcs)(int *, double *)) {
	 int N=2;
	 odeforcs(&N, forc);
}

/** Derivatives (ODE system) **/
void derivs (int *n, double *t, double *y, double *ydot, double *RPAR, int *IPAR) {


	 ydot[0] = rGrow*y[0]*(1-y[0]/K)-rIng*y[0]*y[1];
 	 ydot[1] = rIng*y[0]*y[1]*assEff-rMort*y[1];

}

/** Jacobian vector of the ODE system **/ void jacvec (int *neq, double *t, double *y, int *j, int *ian, int *jan, double *pdj, double *yout, int *iout) {

	 int i;
for(i=0; i<*neq; i++) pdj[i] = 0.;
	 if(*j == 0 ) {
	 pdj[0] = rGrow*(1-y[0]/K)-rGrow*y[0]*(1/K)-rIng*y[1];
 	 pdj[1] = rIng*y[1]*assEff;
	 }
	 else if(*j == 1 ) {
	 pdj[0] = -(rIng*y[0]);
 	 pdj[1] = rIng*y[0]*assEff-rMort;
	 }

}

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models