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

Hello Daniel,

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

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;
	 }

}

Hello Daniel,

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

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;
	 }

}