I know a lot of people met this problem and I have gone through a lot of them. But I still can't solve my problem. Here is it:
My code is c++ code. But I need to call a shared library made by c (ODE solver) for my problem.
int Reaction::cvodedensewrapper()
This is where the functions of ODE solver are called. My class is Reaction and cvodedenseWrapper() is a class function. I define a struct for data transfer from c++ to c as:
typedef struct {
unsigned int exp_ID, nExp,nRec, nComp;
realtype *Ts;
realtype *dTs;
realtype *ts;
realtype *vkp;
realtype Tr;
} User;
typedef User *UserData;
Here realtype is a name for double type. I define a UserData type variable in my code, allocate memory and initialize it as:
User data;
UserData udata=& data;
udata->Ts=new realtype[numExp];
udata->dTs=new realtype[numExp];
udata->ts=new realtype[numExp];
udata->vkp=new realtype[2*numRec];
udata->nExp=numExp;
udata->nRec=numRec;
udata->nComp=numComp;
udata->Tr=Tref;
for(i=0;i!=numExp;i++)
{
udata->Ts[i]=Tstart[i];
udata->dTs[i]=dT[i];
udata->ts[i]=tstart[i];
}
for(i=0;i!=2*numRec;i++)
{
udata->vkp[i]=vk[i];
}
The udata will be passed to ODE solver six times in a loop. At each time, I monitor the address of the dynamic arrays and its content with the code:
printf("udata->Ts: %p\n",udata->Ts);
for(i=0;i!=numExp;i++)
printf(" Ts[%d]=%.5f ",i,udata->Ts[i]);
printf("\n");
printf("udata->dTs: %p\n",udata->dTs);
for(i=0;i!=numExp;i++)
printf(" dTs[%d]=%.5f ",i,udata->dTs[i]);
printf("\n");
printf("udata->ts: %p\n",udata->ts);
for(i=0;i!=numExp;i++)
printf(" ts[%d]=%.5f ",i,udata->ts[i]);
printf("\n");
printf("udata->vkp: %p\n",udata->vkp);
for(i=0;i!=2*numRec;i++)
printf(" vkp[%d]=%.5f ",i,udata->vkp[i]);
printf("\n");
The output at the first time:
m = 0
udata->Ts: 0x62c430
Ts[0]=673.00000 Ts[1]=698.00000 Ts[2]=723.00000 Ts[3]=748.00000 Ts[4]=773.00000 Ts[5]=798.00000
udata->dTs: 0x638f30
dTs[0]=0.00000 dTs[1]=0.00000 dTs[2]=0.00000 dTs[3]=0.00000 dTs[4]=0.00000 dTs[5]=0.00000
udata->ts: 0x638f70
ts[0]=316.09846 ts[1]=283.41932 ts[2]=253.00015 ts[3]=224.61435 ts[4]=198.06463 ts[5]=173.17842
udata->vkp: 0x62e150
vkp[0]=-4.22868 vkp[1]=-2.73312 vkp[2]=2.80980 vkp[3]=-1.65961 vkp[4]=27.23419 vkp[5]=25.96871vkp[6]=16.03558 vkp[7]=10.69316
As comparison, here is the sixth output,
udata->Ts: 0x62c430
Ts[0]=673.00000 Ts[1]=698.00000 Ts[2]=723.00000 Ts[3]=748.00000 Ts[4]=773.00000 Ts[5]=798.00000
udata->dTs: 0x638f30
dTs[0]=0.00000 dTs[1]=0.00000 dTs[2]=0.00000 dTs[3]=0.00000 dTs[4]=0.00000 dTs[5]=0.00000
udata->ts: 0x638f70
ts[0]=316.09846 ts[1]=283.41932 ts[2]=253.00015 ts[3]=224.61435 ts[4]=198.06463 ts[5]=173.17842
udata->vkp: 0x62e150
vkp[0]=-4.22868 vkp[1]=-2.73312 vkp[2]=2.80980 vkp[3]=-1.65961 vkp[4]=27.23419 vkp[5]=25.96871vkp[6]=16.03558 vkp[7]=10.69316
After the loop, I just use the following syntax to release memory,
delete [] udata->Ts;
delete [] udata->dTs;
delete [] udata->ts;
delete [] udata->vkp;
udata=NULL;
It can be seen the pointer and the content pointer points to keeps the same. Nothing changed. So the problems other met is excluded. And this is the only place I use delete. ERROR information is as follows:
*** glibc detected *** /home/fenglei/mopso_oil_shale_cvode/mopso: munmap_chunk(): invalid pointer: 0x000000000062e150 ***
======= Backtrace: =========
/lib/x86_64-linux-gnu/libc.so.6(+0x7eb96)[0x7ffff7286b96]
/home/fenglei/mopso_oil_shale_cvode/mopso[0x405890]
/home/fenglei/mopso_oil_shale_cvode/mopso[0x404b28]
/home/fenglei/mopso_oil_shale_cvode/mopso[0x407482]
/home/fenglei/mopso_oil_shale_cvode/mopso[0x409da5]
/home/fenglei/mopso_oil_shale_cvode/mopso[0x40ba63]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xed)[0x7ffff722976d]
/home/fenglei/mopso_oil_shale_cvode/mopso[0x401a09]
======= Memory map: ========
00400000-00429000 r-xp 00000000 08:06 4209122 {directory}/example
00628000-00629000 r--p 00028000 08:06 4209122 {directory}/example
00629000-0062a000 rw-p 00029000 08:06 4209122 {directory}/example
0062a000-0064b000 rw-p 00000000 00:00 0 [heap]
7ffff7208000-7ffff73bd000 r-xp 00000000 08:06 131963 /lib/x86_64- linux-gnu/libc-2.15.so
7ffff73bd000-7ffff75bd000 ---p 001b5000 08:06 131963 /lib/x86_64-linux-gnu/libc-2.15.so
7ffff75bd000-7ffff75c1000 r--p 001b5000 08:06 131963 /lib/x86_64-linux-gnu/libc-2.15.so
7ffff75c1000-7ffff75c3000 rw-p 001b9000 08:06 131963 /lib/x86_64-linux-gnu/libc-2.15.so
7ffff75c3000-7ffff75c8000 rw-p 00000000 00:00 0
7ffff75c8000-7ffff75dd000 r-xp 00000000 08:06 134855 /lib/x86_64-linux-gnu/libgcc_s.so.1
7ffff75dd000-7ffff77dc000 ---p 00015000 08:06 134855 /lib/x86_64-linux-gnu/libgcc_s.so.1
7ffff77dc000-7ffff77dd000 r--p 00014000 08:06 134855 /lib/x86_64-linux-gnu/libgcc_s.so.1
7ffff77dd000-7ffff77de000 rw-p 00015000 08:06 134855 /lib/x86_64-linux-gnu/libgcc_s.so.1
7ffff77de000-7ffff78d9000 r-xp 00000000 08:06 131976 /lib/x86_64-linux-gnu/libm-2.15.so
7ffff78d9000-7ffff7ad8000 ---p 000fb000 08:06 131976 /lib/x86_64-linux-gnu/libm-2.15.so
7ffff7ad8000-7ffff7ad9000 r--p 000fa000 08:06 131976 /lib/x86_64-linux-gnu/libm-2.15.so
7ffff7ad9000-7ffff7ada000 rw-p 000fb000 08:06 131976 /lib/x86_64-linux-gnu/libm-2.15.so
7ffff7ada000-7ffff7bbc000 r-xp 00000000 08:06 5250632 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7ffff7bbc000-7ffff7dbb000 ---p 000e2000 08:06 5250632 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7ffff7dbb000-7ffff7dc3000 r--p 000e1000 08:06 5250632 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7ffff7dc3000-7ffff7dc5000 rw-p 000e9000 08:06 5250632 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7ffff7dc5000-7ffff7dda000 rw-p 00000000 00:00 0
7ffff7dda000-7ffff7dfc000 r-xp 00000000 08:06 131977 /lib/x86_64-linux-gnu/ld-2.15.so
7ffff7fdc000-7ffff7fe1000 rw-p 00000000 00:00 0
7ffff7ff6000-7ffff7ffa000 rw-p 00000000 00:00 0
7ffff7ffa000-7ffff7ffc000 r-xp 00000000 00:00 0 [vdso]
7ffff7ffc000-7ffff7ffd000 r--p 00022000 08:06 131977 /lib/x86_64-linux-gnu/ld-2.15.so
7ffff7ffd000-7ffff7fff000 rw-p 00023000 08:06 131977 /lib/x86_64-linux-gnu/ld-2.15.so
7ffffffdd000-7ffffffff000 rw-p 00000000 00:00 0 [stack]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0 [vsyscall]
Program received signal SIGABRT, Aborted.
0x00007ffff723e425 in raise () from /lib/x86_64-linux-gnu/libc.so.6
This is the first time I post my questions here, although i got a lot of stuff from this website. I hope I made my problem clear.
Code:
/*
---------------------------------
*ODE problem solver definition
---------------------------------
*/
#include <stdio.h>
/* Header files with a description of contents used */
#include <cvode/cvode.h> /* prototypes for CVODE fcts., consts. */
#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
#include <cvode/cvode_dense.h> /* prototype for CVDense */
#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
#include <sundials/sundials_types.h> /* definition of type realtype */
/* User-defined vector and matrix accessor macros: Ith, IJth */
/* These macros are defined in order to write code which exactly matches
the mathematical problem description given above.
Ith(v,i) references the ith component of the vector v, where i is in
the range [1..NEQ] and NEQ is defined below. The Ith macro is defined
using the N_VIth macro in nvector.h. N_VIth numbers the components of
a vector starting from 0.
IJth(A,i,j) references the (i,j)th element of the dense matrix A, where
i and j are in the range [1..NEQ]. The IJth macro is defined using the
DENSE_ELEM macro in dense.h. DENSE_ELEM numbers rows and columns of a
dense matrix starting from 0. */
#define Ith(v,i) NV_Ith_S(v,i-1) /* Ith numbers components 1..NEQ */
#define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */
/* Problem Constants */
#define RTOL RCONST(1.0e-4) /* scalar relative tolerance */
#define ATOL1 RCONST(1.0e-8) /* vector absolute tolerance components */
#define ATOL2 RCONST(1.0e-8)
#define ATOL3 RCONST(1.0e-8)
#define ATOL4 RCONST(1.0e-8)
#define Y(i) yt[i-1]
#define K(i) kt[i-1]
typedef struct {
unsigned int exp_ID, nExp,nRec, nComp;
realtype *Ts;
realtype *dTs;
realtype *ts;
realtype *vkp;
realtype Tr;
} User;
typedef User *UserData;
/* Functions Called by the Solver */
static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
static int Jac(long int N, realtype t,
N_Vector y, N_Vector fy, DlsMat J, void *user_data,
N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
/* Private function to check function return values */
static int check_flag(void *flagvalue, const char *funcname, int opt);
int Reaction::cvodedensewrapper()
{
int i=0,j=0,k=0,m=0,numt;
realtype reltol,t,tout,ts;
N_Vector y, abstol;
void *cvode_mem;
int iout, flag;
User data;
UserData udata=& data;
/*Allocation and initialization*/
udata->Ts=new realtype[numExp];
udata->dTs=new realtype[numExp];
udata->ts=new realtype[numExp];
udata->vkp=new realtype[2*numRec];
udata->nExp=numExp;
udata->nRec=numRec;
udata->nComp=numComp;
udata->Tr=Tref;
for(i=0;i!=numExp;i++)
{
udata->Ts[i]=Tstart[i];
udata->dTs[i]=dT[i];
udata->ts[i]=tstart[i];
}
for(i=0;i!=2*numRec;i++)
{
udata->vkp[i]=vk[i];
}
/*for debugging*/
printf("udata->Ts: %p\n",udata->Ts);
for(i=0;i!=numExp;i++)
printf(" Ts[%d]=%.5f ",i,udata->Ts[i]);
printf("\n");
printf("udata->dTs: %p\n",udata->dTs);
for(i=0;i!=numExp;i++)
printf(" dTs[%d]=%.5f ",i,udata->dTs[i]);
printf("\n");
printf("udata->ts: %p\n",udata->ts);
for(i=0;i!=numExp;i++)
printf(" ts[%d]=%.5f ",i,udata->ts[i]);
printf("\n");
printf("udata->vkp: %p\n",udata->vkp);
for(i=0;i!=2*numRec;i++)
printf(" vkp[%d]=%.5f ",i,udata->vkp[i]);
printf("\n");
/*debugging end*/
y = abstol = NULL;
cvode_mem = NULL;
/* Create serial vector of length numRec for I.C. and abstol */
y = N_VNew_Serial(numComp);
if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
abstol = N_VNew_Serial(numComp);
if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(1);
/* Initialize y */
for(i=0;i!=numComp;i++)
{
//Ith(y,i)=Yinitial[0][i-1];
NV_Ith_S(y,i)=Yinitial[0][i];
}
/*Initialize tstart*/
ts=tstart[0];
/* Set the scalar relative tolerance */
reltol = RTOL;
/* Set the vector absolute tolerance */
Ith(abstol,1)=ATOL1;
Ith(abstol,2)=ATOL2;
Ith(abstol,3)=ATOL3;
Ith(abstol,4)=ATOL3;
/* Call CVodeCreate to create the solver memory and specify the
* Backward Differentiation Formula and the use of a Newton iteration */
cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
/* Set the pointer to user-defined data*/
flag=CVodeSetUserData(cvode_mem,udata);
if(check_flag((void *) cvode_mem,"CVodeCreate",0)) return(1);
/* Call CVodeInit to initialize the integrator memory and specify the
* user's right hand side function in y'=f(t,y), the inital time T0, and
* the initial dependent variable vector y. */
flag = CVodeInit(cvode_mem, f, ts, y);
if (check_flag(&flag, "CVodeInit", 1)) return(1);
/* Call CVodeSVtolerances to specify the scalar relative tolerance
* and vector absolute tolerances */
flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
/* Call CVDense to specify the CVDENSE dense linear solver */
flag = CVDense(cvode_mem, numComp);
if (check_flag(&flag, "CVDense", 1)) return(1);
/* Set the Jacobian routine to Jac (user-supplied) */
flag = CVDlsSetDenseJacFn(cvode_mem, Jac);
if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);
/* In loop, call CVode, print results, and test for error.
Break out of loop when NOUT preset output times have been reached. */
for (m=0;m!=numExp;m++) {
numt=(int)(tend[m]-tstart[m])/dtsave[m]+1;
udata->exp_ID=m;
Yt[m].tvector[0]=tstart[m];
Yt[m].Tvector[0]=Tstart[m];
for(i=0;i!=numComp;i++)
{
Yt[m].resp[i][0]=Yinitial[m][i];
}
for(i=0;i!=numComp;i++)
{
//Ith(y,i)=Yinitial[m][i-1];
NV_Ith_S(y,i)=Yinitial[m][i];
}
ts=tstart[m];
flag = CVodeReInit(cvode_mem,ts, y);
iout = 0;
tout = tstart[m]+dtsave[m];
while(1) {
flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
if (check_flag(&flag, "CVodeInit", 1)) return(1);
if (flag == CV_SUCCESS) {
iout++;
tout=t+dtsave[m];
}
Yt[m].tvector[iout]=t;
Yt[m].Tvector[iout]=Tstart[m]+(Yt[m].tvector[iout]-Yt[m].tvector[0])*dT[m];
for(i=0;i!=numComp;i++)
{
//Yt[m].resp[i-1][iout]=Ith(y,i);
Yt[m].resp[i][iout]=NV_Ith_S(y,i);
}
if (iout == numt) break;
}
/*for debugging*/
printf("udata->Ts: %p\n",udata->Ts);
for(i=0;i!=numExp;i++)
printf(" Ts[%d]=%.5f ",i,udata->Ts[i]);
printf("\n");
printf("udata->dTs: %p\n",udata->dTs);
for(i=0;i!=numExp;i++)
printf(" dTs[%d]=%.5f ",i,udata->dTs[i]);
printf("\n");
printf("udata->ts: %p\n",udata->ts);
for(i=0;i!=numExp;i++)
printf(" ts[%d]=%.5f ",i,udata->ts[i]);
printf("\n");
printf("udata->vkp: %p\n",udata->vkp);
for(i=0;i!=2*numRec;i++)
printf(" vkp[%d]=%.5f ",i,udata->vkp[i]);
printf("\n");
/*debugging end*/
}
/*Free userdata*/
delete [] udata->Ts;
delete [] udata->dTs;
delete [] udata->ts;
delete [] udata->vkp;
udata=NULL;
/*Free ODE solver array*/
N_VDestroy_Serial(y);
N_VDestroy_Serial(abstol);
/* Free integrator memory */
CVodeFree(&cvode_mem);
return(0);
}
/*
*-------------------------------
* Functions called by the solver
*-------------------------------
*/
/*
* f routine. Compute function f(t,y).
*/
/*
* Jacobian routine. Compute J(t,y) = df/dy. *
*/
Functions with CV initils are from shared libraries.
