00001 #include <math.h>
00002 #include <iostream>
00003 #include "../nr.h"
00004 #include "../def.h"
00005 #include "../model.h"
00006 #include <string.h>
00007 #include "../libCVODES/sundialstypes.h"
00008 #include "../libCVODES/cvodes.h"
00009 #include "../libCVODES/cvdense.h"
00010 #include "../libCVODES/nvector_serial.h"
00011 #include "../libCVODES/dense.h"
00012 #include "../libCVODES/sundialsmath.h"
00013
00014 using namespace std;
00015
00016
00017
00018 #define Ith(v,i) NV_Ith_S(v,i-1) // i-th vector component
00019 #define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) // (i,j)-th matrix component
00020
00021 #define RTOL RCONST(1e-8)
00022 #define ATOL RCONST(1e-8)
00023
00024
00025
00026 typedef struct {
00027 double *p;
00028 long nvar;
00029 int sensi;
00030 } *UserData;
00031
00032
00033
00034 void f(realtype t, N_Vector y, N_Vector ydot, void *f_data);
00035
00036 static int ewt(N_Vector y, N_Vector w, void *e_data);
00037
00038 void sensDerivs (double t, double *Yt, double *Ft, double *p);
00039
00040 void integrateCVODES(Glob *globs,double ystart[], double p[], long nvar, double x1,
00041 double x2,double eps,int sensi)
00042 {
00043
00044 void *cvode_mem;
00045 UserData data;
00046 realtype t,t0,t1;
00047 N_Vector y;
00048 long i;
00049
00050 cvode_mem = NULL;
00051 data = NULL;
00052 y = NULL;
00053
00054
00055 data = (UserData) malloc(sizeof *data);
00056
00057
00058 data->p=p;
00059 data->nvar=nvar;
00060 data->sensi=sensi;
00061
00062
00063 y = N_VNew_Serial(nvar);
00064 for(i=1;i<=nvar;i++)
00065 Ith(y,i)=ystart[i];
00066
00067 t0= RCONST(x1);
00068 t1= RCONST(x2);
00069 t = RCONST(t1);
00070
00071
00072 cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
00073
00074 realtype atol=ATOL;
00075 CVodeMalloc(cvode_mem, f, t0, y, CV_SS, RTOL, &atol);
00076
00077 CVodeSetFdata(cvode_mem, data);
00078
00079 CVDense(cvode_mem, nvar);
00080
00081
00082 CVode(cvode_mem, t1, y, &t, CV_NORMAL);
00083
00084
00085 for(i=1;i<=nvar;i++)
00086 ystart[i]=Ith(y,i);
00087
00088
00089 N_VDestroy_Serial(y);
00090 free(data);
00091 CVodeFree(cvode_mem);
00092 }
00093
00094
00095
00096
00097
00098 void f(realtype t, N_Vector y, N_Vector ydot, void *f_data)
00099 {
00100 UserData data;
00101 data = (UserData) f_data;
00102 long i,nvar=data->nvar;
00103 double *ydot_nr=(double*)malloc((nvar+1)*sizeof(double));
00104 double *y_nr=(double*)malloc((nvar+1)*sizeof(double));
00105
00106 for(i=1;i<=nvar;i++)
00107 y_nr[i]=Ith(y,i);
00108
00109 if(data->sensi==TRUE)
00110 sensDerivs((double)t,y_nr,ydot_nr,data->p);
00111 else
00112 ode((double)t,y_nr,ydot_nr,data->p);
00113
00114 for(i=1;i<=nvar;i++)
00115 Ith(ydot,i)=ydot_nr[i];
00116
00117 free(y_nr);
00118 free(ydot_nr);
00119
00120 }
00121
00122