[R-sig-dyn-mod] lsodes with jacobian
Karline Soetaert
Karline.Soetaert at nioz.nl
Mon Jul 7 19:16:46 CEST 2014
Daniel,
1. ian should be c(1, 2, 4, 6) where the last element equals the number of nonzeros plus 1. Else the solver does not know where jan stops (:-)
2. I think the description in the LSODES fortran code is wrong, and it should be an integer (as in the example in the fortran code).
However, it does not matter, as ian and jan can be ignored in the jacobian function (they should be known, as passed to the solver).
3 I think your pd function, for (*j == 0) is wrong:
You write:
pdj[0] = rGrow*(1-y[0]/K)-rGrow*y[0]*(1/K)-rIng*y[1];
I think it should be:
pdj[0] = rGrow*(1-2*y[0]/K)-rGrow*y[0]*(1/K)-rIng*y[1];
( as d(y^2)/dy = 2*y)
I wonder if all this effort is worth it for such a small model? But maybe you are doing this for larger problems?
4. You may also want to have a look at the file deSolve_utils.c, where I have implemented the determination of the sparsity (ian and jan) for 1D, 2D and 3dimensional problems.
Hope all 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: maandag 7 juli 2014 16:38
To: r-sig-dynamic-models at r-project.org
Subject: Re: [R-sig-dyn-mod] lsodes with jacobian
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
_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
More information about the R-sig-dynamic-models
mailing list