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

[Daniel Kaschek]

Hi Karline,

On Mo, 2014-07-07 at 08:09 +0000, Karline Soetaert wrote:
> Hello Daniel,
> 
> Sorry for answering this late; I was at the useR conference.

I am glad that you respond :-)

> 
> 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).

I had a look at the documentation of sparsetype="sparsejan". I don't
understand why "ian" should have n+1 elements. Let's assume the jacobian
matrix

0 1 0
0 0 1
1 1 1

Then, I have

jan = c(3, 1, 3, 2, 3) #third element in column 1, first and third
element in column 2, ...

ian = c(1, 2, 4) #column1 starts at element 1 of jan, column2 starts at
element 2 of jan, column3 starts at element 4 of jan.

In case, I did not understand the documentation correctly. Could you
tell me "jan" and "ian" for the above example?

> 
> 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.

I tried to run my exmample with verbose=TRUE. However, when specifying
jacvec, I get the following error:

-------------------- 
Integration method 
--------------------
 
Error in cat(message, "\n") : object 'txt' not found

> 
> For large problems, the explicit specification should be faster.

I have one more question about the correct arguments of jacvec in the C
code. In the fortran code (odepack.f from
http://people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.f) I find
the specifications

SUBROUTINE JAC (NEQ, T, Y, J, IAN, JAN, PDJ)
DOUBLE PRECISION T, Y(*), IAN(*), JAN(*), PDJ(*)

In contrast, in the deSolve package in the file "call_lsoda.c", I find

static void C_jac_vec (int *neq, double *t, double *y, int *j,
            int *ian, int *jan, double *pdj, double *yout, int *iout)

My point is, that IAN and JAN are double in fortran and integers in C.
Is this a contradiction?

Sorry for so many questions.

Best regards,
Daniel



> 
> 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
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

-- 
Daniel Kaschek
Institute of Physics, Freiburg
Office: 210
Phone: +49-761-203-8531