自学教程：C++ CVode函数代码示例

int cvode_eval(solver_props *props, unsigned int modelid){  cvode_mem *mem = props->mem;  mem = &mem[modelid];  // Stop the solver if the stop time has been reached  props->running[modelid] = (props->time[modelid] + props->timestep) < props->stoptime;  if(!props->running[modelid])    return 0;  // if a positive dt is specified, then we will have this function return after it reaches the next time point,  // otherwise, it will just run one iteration and return  if(props->timestep > 0) {    // Reinitialize the function at each step    if(CVodeReInit(mem->cvmem, props->time[modelid], mem->y0) != CV_SUCCESS) {      PRINTF( "CVODE failed to reinitialize");    }    CDATAFORMAT stop_time = MIN(props->time[modelid] + props->timestep, props->stoptime);    mem->first_iteration = TRUE;    if(CVode(mem->cvmem, stop_time, mem->y0, &(props->next_time[modelid]), CV_NORMAL) != CV_SUCCESS){      PRINTF( "CVODE failed to make a fixed step in model %d./n", modelid);      return 1;    }  }  else {    mem->first_iteration = TRUE;    if(CVode(mem->cvmem, props->stoptime, mem->y0, &(props->next_time[modelid]), CV_ONE_STEP) != CV_SUCCESS){      PRINTF( "CVODE failed to make a step in model %d./n", modelid);      return 1;    }  }  return 0;}

 void CVodeInt::integrate(double tout) {   double t;   int flag;   flag = CVode(m_cvode_mem, tout, nv(m_y), &t, NORMAL);   if (flag != SUCCESS)      throw CVodeErr(" CVode error encountered."); }

 double CVodesIntegrator::step(double tout) {   double t;   int flag;   flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_ONE_STEP);   if (flag != CV_SUCCESS)      throw CVodesErr(" CVodes error encountered.");   return t; }

void CVodesIntegrator::integrate(double tout){    int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_NORMAL);    if (flag != CV_SUCCESS) {        throw CVodesErr("CVodes error encountered. Error code: " + int2str(flag) + "/n" + m_error_message +                        "/nComponents with largest weighted error estimates:/n" + getErrorInfo(10));    }    m_sens_ok = false;}

double CVodesIntegrator::step(double tout){    int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_ONE_STEP);    if (flag != CV_SUCCESS) {        throw CanteraError("CVodesIntegrator::step",            "CVodes error encountered. Error code: {}/n{}/n"            "Components with largest weighted error estimates:/n{}",            flag, m_error_message, getErrorInfo(10));    }    m_sens_ok = false;    return m_time;}

int ode_solver_solve(ode_solver* solver, const double t, double* y, double* tout){    double lTout;    /* Advance the solution */  NV_DATA_S(solver->y) = y;  int flag = CVode(solver->cvode_mem,t, solver->y, &lTout, CV_NORMAL);    if(tout)    *tout = lTout;    return flag;}

int SOLVER(cvode, eval, TARGET, SIMENGINE_STORAGE, cvode_mem *mem, unsigned int modelid) {  // Stop the solver if the stop time has been reached  mem->props->running[modelid] = mem->props->time[modelid] < mem->props->stoptime;  if(!mem->props->running[modelid])    return 0;  mem[modelid].first_iteration = TRUE;  if(CVode(mem[modelid].cvmem, mem[modelid].props->stoptime, ((N_Vector)(mem[modelid].y0)), &(mem[modelid].props->time[modelid]), CV_ONE_STEP) != CV_SUCCESS){    fprintf(stderr, "CVODE failed to make a step in model %d./n", modelid);    return 1;  }  return 0;}

real Solver::run(real tout, int &ncalls, real &rhstime){#ifdef CHECK  int msg_point = msg_stack.push("Running solver: solver::run(%e)", tout);#endif  MPI_Barrier(MPI_COMM_WORLD);  rhs_wtime = 0.0;  rhs_ncalls = 0;  pre_Wtime = 0.0;  pre_ncalls = 0.0;  int flag = CVode(cvode_mem, tout, uvec, &simtime, CV_NORMAL);    ncalls = rhs_ncalls;  rhstime = rhs_wtime;  // Copy variables  load_vars(NV_DATA_P(uvec));  // Call rhs function to get extra variables at this time  real tstart = MPI_Wtime();  (*func)(simtime);  rhs_wtime += MPI_Wtime() - tstart;  rhs_ncalls++;    if(flag < 0) {    output.write("ERROR CVODE solve failed at t = %e, flag = %d/n", simtime, flag);    return -1.0;  }  #ifdef CHECK  msg_stack.pop(msg_point);#endif  return simtime;}

int integrate(struct Integrator* integrator, double tout, double* t){  if (integrator->em->nRates > 0)  {    /* need to integrate if we have any differential equations */    int flag;    /* Make sure we don't go past the specified end time - could run into       trouble if we're almost reaching a threshold */    flag = CVodeSetStopTime(integrator->cvode_mem,(realtype)tout);    if (check_flag(&flag,"CVode",1)) return(ERR);    flag = CVode(integrator->cvode_mem,tout,integrator->y,t,CV_NORMAL);    if (check_flag(&flag,"CVode",1)) return(ERR);    /* we also need to evaluate all the other variables that are not required       to be updated during integration */    integrator->em->evaluateVariables(*t);  }  else  {	/* no differential equations so just evaluate once */	integrator->em->computeRates(tout);	integrator->em->evaluateVariables(tout);    *t = tout;  }  /*   * Now that using CV_NORMAL_TSTOP this is no longer required?   */#ifdef OLD_CODE  /* only the y array is gonna be at the desired tout, so we need to also     update the full variables array */  ud->BOUND[0] = *t;  ud->methods->ComputeVariables(ud->BOUND,ud->RATES,ud->CONSTANTS,    ud->VARIABLES);#endif    /* Make sure the outputs are up-to-date */  integrator->em->getOutputs(*t);  return(OK);}

int CVAdataStore(CVadjMem ca_mem, CkpntMem ck_mem){  CVodeMem cv_mem;  DtpntMem *dt_mem;  realtype t;  long int i;  int flag;  cv_mem = ca_mem->cv_mem;  dt_mem = ca_mem->dt_mem;  /* Initialize cv_mem with data from ck_mem */  flag = CVAckpntGet(cv_mem, ck_mem);  if (flag != CV_SUCCESS) return(CV_REIFWD_FAIL);  /* Set first structure in dt_mem[0] */  dt_mem[0]->t = t0_;  storePnt(cv_mem, dt_mem[0]);  /* Run CVode to set following structures in dt_mem[i] */  i = 1;  do {    flag = CVode(cv_mem, t1_, ytmp, &t, CV_ONE_STEP);    if (flag < 0) return(CV_FWD_FAIL);    dt_mem[i]->t = t;    storePnt(cv_mem, dt_mem[i]);    i++;  } while (t<t1_);  /* New data is now available */  ckpntData = ck_mem;  newData = TRUE;  np  = i;  return(CV_SUCCESS);}

PetscErrorCode TSStep_Sundials(TS ts){  TS_Sundials    *cvode = (TS_Sundials*)ts->data;  PetscErrorCode ierr;  PetscInt       flag;  long int       its,nsteps;  realtype       t,tout;  PetscScalar    *y_data;  void           *mem;  PetscFunctionBegin;  mem  = cvode->mem;  tout = ts->max_time;  ierr = VecGetArray(ts->vec_sol,&y_data);CHKERRQ(ierr);  N_VSetArrayPointer((realtype*)y_data,cvode->y);  ierr = VecRestoreArray(ts->vec_sol,NULL);CHKERRQ(ierr);  ierr = TSPreStep(ts);CHKERRQ(ierr);  /* We would like to call TSPreStep() when starting each step (including rejections) and TSPreStage() before each   * stage solve, but CVode does not appear to support this. */  if (cvode->monitorstep) flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);  else flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);  if (flag) { /* display error message */    switch (flag) {      case CV_ILL_INPUT:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");        break;      case CV_TOO_CLOSE:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");        break;      case CV_TOO_MUCH_WORK: {        PetscReal      tcur;        ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr);        ierr = CVodeGetCurrentTime(mem,&tcur);CHKERRQ(ierr);        SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_WORK. At t=%G, nsteps %D exceeds mxstep %D. Increase '-ts_max_steps <>' or modify TSSetDuration()",tcur,nsteps,ts->max_steps);      } break;      case CV_TOO_MUCH_ACC:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");        break;      case CV_ERR_FAILURE:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");        break;      case CV_CONV_FAILURE:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");        break;      case CV_LINIT_FAIL:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");        break;      case CV_LSETUP_FAIL:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");        break;      case CV_LSOLVE_FAIL:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");        break;      case CV_RHSFUNC_FAIL:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");        break;      case CV_FIRST_RHSFUNC_ERR:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");        break;      case CV_REPTD_RHSFUNC_ERR:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");        break;      case CV_UNREC_RHSFUNC_ERR:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");        break;      case CV_RTFUNC_FAIL:        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");        break;      default:        SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);    }  }  /* copy the solution from cvode->y to cvode->update and sol */  ierr = VecPlaceArray(cvode->w1,y_data);CHKERRQ(ierr);  ierr = VecCopy(cvode->w1,cvode->update);CHKERRQ(ierr);  ierr = VecResetArray(cvode->w1);CHKERRQ(ierr);  ierr = VecCopy(cvode->update,ts->vec_sol);CHKERRQ(ierr);  ierr = CVodeGetNumNonlinSolvIters(mem,&its);CHKERRQ(ierr);  ierr = CVSpilsGetNumLinIters(mem, &its);  ts->snes_its = its; ts->ksp_its = its;  ts->time_step = t - ts->ptime;  ts->ptime     = t;  ts->steps++;  ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr);  if (!cvode->monitorstep) ts->steps = nsteps;  PetscFunctionReturn(0);}

int main(){    realtype abstol=ATOL, reltol=RTOL, t, tout;    N_Vector c;    WebData wdata;    void *cvode_mem;    booleantype firstrun;    int jpre, gstype, flag;    int ns, mxns, iout;    c = NULL;    wdata = NULL;    cvode_mem = NULL;    /* Initializations */    c = N_VNew_Serial(NEQ);    if(check_flag((void *)c, "N_VNew_Serial", 0)) return(1);    wdata = AllocUserData();    if(check_flag((void *)wdata, "AllocUserData", 2)) return(1);    InitUserData(wdata);    ns = wdata->ns;    mxns = wdata->mxns;    /* Print problem description */    PrintIntro();    /* Loop over jpre and gstype (four cases) */    for (jpre = PREC_LEFT; jpre <= PREC_RIGHT; jpre++) {        for (gstype = MODIFIED_GS; gstype <= CLASSICAL_GS; gstype++) {            /* Initialize c and print heading */            CInit(c, wdata);            PrintHeader(jpre, gstype);            /* Call CVodeInit or CVodeReInit, then CVSpgmr to set up problem */            firstrun = (jpre == PREC_LEFT) && (gstype == MODIFIED_GS);            if (firstrun) {                cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);                if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);                wdata->cvode_mem = cvode_mem;                flag = CVodeSetUserData(cvode_mem, wdata);                if(check_flag(&flag, "CVodeSetUserData", 1)) return(1);                flag = CVodeInit(cvode_mem, f, T0, c);                if(check_flag(&flag, "CVodeInit", 1)) return(1);                flag = CVodeSStolerances(cvode_mem, reltol, abstol);                if (check_flag(&flag, "CVodeSStolerances", 1)) return(1);                flag = CVSpgmr(cvode_mem, jpre, MAXL);                if(check_flag(&flag, "CVSpgmr", 1)) return(1);                flag = CVSpilsSetGSType(cvode_mem, gstype);                if(check_flag(&flag, "CVSpilsSetGSType", 1)) return(1);                flag = CVSpilsSetEpsLin(cvode_mem, DELT);                if(check_flag(&flag, "CVSpilsSetEpsLin", 1)) return(1);                flag = CVSpilsSetPreconditioner(cvode_mem, Precond, PSolve);                if(check_flag(&flag, "CVSpilsSetPreconditioner", 1)) return(1);            } else {                flag = CVodeReInit(cvode_mem, T0, c);                if(check_flag(&flag, "CVodeReInit", 1)) return(1);                flag = CVSpilsSetPrecType(cvode_mem, jpre);                check_flag(&flag, "CVSpilsSetPrecType", 1);                flag = CVSpilsSetGSType(cvode_mem, gstype);                if(check_flag(&flag, "CVSpilsSetGSType", 1)) return(1);            }            /* Print initial values */            if (firstrun) PrintAllSpecies(c, ns, mxns, T0);            /* Loop over output points, call CVode, print sample solution values. */            tout = T1;            for (iout = 1; iout <= NOUT; iout++) {                flag = CVode(cvode_mem, tout, c, &t, CV_NORMAL);                PrintOutput(cvode_mem, t);                if (firstrun && (iout % 3 == 0)) PrintAllSpecies(c, ns, mxns, t);                if(check_flag(&flag, "CVode", 1)) break;                if (tout > RCONST(0.9)) tout += DTOUT;                else tout *= TOUT_MULT;            }            /* Print final statistics, and loop for next case */            PrintFinalStats(cvode_mem);        }    }    /* Free all memory */    CVodeFree(&cvode_mem);    N_VDestroy_Serial(c);    FreeUserData(wdata);//.........这里部分代码省略.........

int main(int argc, char *argv[]){  void *cvode_mem;  UserData data;  realtype t, tout;  N_Vector y;  int iout, flag, nthreads, nnz;  realtype pbar[NS];  int is;   N_Vector *yS;  booleantype sensi, err_con;  int sensi_meth;  cvode_mem = NULL;  data      = NULL;  y         =  NULL;  yS        = NULL;  /* Process arguments */  ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);  /* User data structure */  data = (UserData) malloc(sizeof *data);  if (check_flag((void *)data, "malloc", 2)) return(1);  data->p[0] = RCONST(0.04);  data->p[1] = RCONST(1.0e4);  data->p[2] = RCONST(3.0e7);  /* Initial conditions */  y = N_VNew_Serial(NEQ);  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);  Ith(y,1) = Y1;  Ith(y,2) = Y2;  Ith(y,3) = Y3;  /* 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);  /* Call CVodeInit to initialize the integrator memory and specify the     user's right hand side function in y'=f(t,y), the initial time T0, and     the initial dependent variable vector y. */  flag = CVodeInit(cvode_mem, f, T0, y);  if (check_flag(&flag, "CVodeInit", 1)) return(1);  /* Call CVodeWFtolerances to specify a user-supplied function ewt that sets     the multiplicative error weights W_i for use in the weighted RMS norm */  flag = CVodeWFtolerances(cvode_mem, ewt);  if (check_flag(&flag, "CVodeSetEwtFn", 1)) return(1);  /* Attach user data */  flag = CVodeSetUserData(cvode_mem, data);  if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);  /* Call CVKLU to specify the CVKLU sparse direct linear solver */  nthreads = 1;                 /* no. of threads to use when factoring the system*/  nnz = NEQ * NEQ;              /* max no. of nonzeros entries in the Jac */  flag = CVSuperLUMT(cvode_mem, nthreads, NEQ, nnz);  if (check_flag(&flag, "CVSuperLUMT", 1)) return(1);  /* Set the Jacobian routine to Jac (user-supplied) */  flag = CVSlsSetSparseJacFn(cvode_mem, Jac);  if (check_flag(&flag, "CVSlsSetSparseJacFn", 1)) return(1);  printf("/n3-species chemical kinetics problem/n");  /* Sensitivity-related settings */  if (sensi) {    /* Set parameter scaling factor */    pbar[0] = data->p[0];    pbar[1] = data->p[1];    pbar[2] = data->p[2];    /* Set sensitivity initial conditions */    yS = N_VCloneVectorArray_Serial(NS, y);    if (check_flag((void *)yS, "N_VCloneVectorArray_Serial", 0)) return(1);    for (is=0;is<NS;is++) N_VConst(ZERO, yS[is]);    /* Call CVodeSensInit1 to activate forward sensitivity computations       and allocate internal memory for COVEDS related to sensitivity       calculations. Computes the right-hand sides of the sensitivity       ODE, one at a time */    flag = CVodeSensInit1(cvode_mem, NS, sensi_meth, fS, yS);    if(check_flag(&flag, "CVodeSensInit", 1)) return(1);    /* Call CVodeSensEEtolerances to estimate tolerances for sensitivity        variables based on the rolerances supplied for states variables and        the scaling factor pbar */    flag = CVodeSensEEtolerances(cvode_mem);    if(check_flag(&flag, "CVodeSensEEtolerances", 1)) return(1);    /* Set sensitivity analysis optional inputs */    /* Call CVodeSetSensErrCon to specify the error control strategy for        sensitivity variables */    flag = CVodeSetSensErrCon(cvode_mem, err_con);    if (check_flag(&flag, "CVodeSetSensErrCon", 1)) return(1);//.........这里部分代码省略.........

int run_rate_state_sim(std::vector<std::vector<realtype> > &results, RSParams &params) {	realtype		long_term_reltol, event_reltol, t, tout, tbase=0;	N_Vector		y, long_term_abstol, event_abstol;	unsigned int	i, n;	int				flag, err_code;	void			*long_term_cvode, *event_cvode, *current_cvode;		// Create serial vector of length NEQ for I.C. and abstol	y = N_VNew_Serial(params.num_eqs()*params.num_blocks());	if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);	long_term_abstol = N_VNew_Serial(params.num_eqs()*params.num_blocks()); 	if (check_flag((void *)long_term_abstol, "N_VNew_Serial", 0)) return(1);	event_abstol = N_VNew_Serial(params.num_eqs()*params.num_blocks()); 	if (check_flag((void *)event_abstol, "N_VNew_Serial", 0)) return(1);		// Initialize y	for (i=0;i<params.num_blocks();++i) {		NV_Ith_S(y,i*params.num_eqs()+EQ_X) = params.init_val(i, EQ_X);		NV_Ith_S(y,i*params.num_eqs()+EQ_V) = params.init_val(i, EQ_V);		NV_Ith_S(y,i*params.num_eqs()+EQ_H) = params.init_val(i, EQ_H);	}		/* Initialize interactions */	/*interaction = new realtype[NBLOCKS*NBLOCKS];	double int_level = 1e-2;	double dropoff = 1.1;	for (i=0;i<NBLOCKS;++i) {		for (n=0;n<NBLOCKS;++n) {			interaction[i*NBLOCKS+n] = (i==n?(1.0-int_level):int_level);		}	}*/		/* Set the scalar relative tolerance */	long_term_reltol = RCONST(1.0e-12);	event_reltol = RCONST(1.0e-12);	/* Set the vector absolute tolerance */	for (i=0;i<params.num_blocks();++i) {		Xth(long_term_abstol,i) = RCONST(1.0e-12);		Vth(long_term_abstol,i) = RCONST(1.0e-12);		Hth(long_term_abstol,i) = RCONST(1.0e-12);		Xth(event_abstol,i) = RCONST(1.0e-12);		Vth(event_abstol,i) = RCONST(1.0e-12);		Hth(event_abstol,i) = RCONST(1.0e-12);	}		/* Call CVodeCreate to create the solver memory and specify the 	 * Backward Differentiation Formula and the use of a Newton iteration */	long_term_cvode = CVodeCreate(CV_BDF, CV_NEWTON);	if (check_flag((void *)long_term_cvode, "CVodeCreate", 0)) return(1);	event_cvode = CVodeCreate(CV_BDF, CV_NEWTON);	if (check_flag((void *)event_cvode, "CVodeCreate", 0)) return(1);		// Turn off error messages	//CVodeSetErrFile(long_term_cvode, NULL);	//CVodeSetErrFile(event_cvode, NULL);		/* 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(long_term_cvode, func, T0, y);	if (check_flag(&flag, "CVodeInit", 1)) return(1);	flag = CVodeInit(event_cvode, func, T0, y);	if (check_flag(&flag, "CVodeInit", 1)) return(1);		/* Call CVodeSVtolerances to specify the scalar relative tolerance	 * and vector absolute tolerances */	flag = CVodeSVtolerances(long_term_cvode, long_term_reltol, long_term_abstol);	if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);	flag = CVodeSVtolerances(event_cvode, event_reltol, event_abstol);	if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);		/* Set the root finding function */	//flag = CVodeRootInit(long_term_cvode, params.num_blocks(), vel_switch_finder);	//flag = CVodeRootInit(event_cvode, params.num_blocks(), vel_switch_finder);	//if (check_flag(&flag, "CVodeRootInit", 1)) return(1);		/* Call CVDense to specify the CVDENSE dense linear solver */	//flag = CVSpbcg(cvode_mem, PREC_NONE, 0);	//if (check_flag(&flag, "CVSpbcg", 1)) return(1);	//flag = CVSpgmr(cvode_mem, PREC_NONE, 0);	//if (check_flag(&flag, "CVSpgmr", 1)) return(1);	flag = CVDense(long_term_cvode, params.num_eqs()*params.num_blocks());	if (check_flag(&flag, "CVDense", 1)) return(1);	flag = CVDense(event_cvode, params.num_eqs()*params.num_blocks());	if (check_flag(&flag, "CVDense", 1)) return(1);		flag = CVodeSetUserData(long_term_cvode, &params);	if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);	flag = CVodeSetUserData(event_cvode, &params);	if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);		CVodeSetMaxNumSteps(long_term_cvode, 100000);	CVodeSetMaxNumSteps(event_cvode, 100000);		/* Set the Jacobian routine to Jac (user-supplied) */	flag = CVDlsSetDenseJacFn(long_term_cvode, Jac);	if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);	flag = CVDlsSetDenseJacFn(event_cvode, Jac);	if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);	//.........这里部分代码省略.........

int main(int narg, char **args){    realtype reltol, t, tout;    N_Vector state, abstol;    void *cvode_mem;    int flag, flagr;    int rootsfound[NRF];    int rootdir[] = {1,};    FILE *pout;    if(!(pout = fopen("results/iaf_v.dat", "w"))){        fprintf(stderr, "Cannot open file results/iaf_v.dat. Are you trying to write to a non-existent directory? Exiting.../n");        exit(1);    }    state = abstol = NULL;    cvode_mem = NULL;    state = N_VNew_Serial(NEQ);    if (check_flag((void *)state, "N_VNew_Serial", 0)) return(1);    abstol = N_VNew_Serial(NEQ);     if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(1);        realtype reset = -0.07;    realtype C = 3.2e-12;    realtype thresh = -0.055;    realtype gleak = 2e-10;    realtype eleak = -0.053;    realtype p[] = {reset, C, thresh, gleak, eleak, };    realtype v = reset;    NV_Ith_S(state, 0) = reset;    reltol = RTOL;    NV_Ith_S(abstol,0) = ATOL0;     /* Allocations and initializations */    cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);    if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);        flag = CVodeInit(cvode_mem, dstate_dt, T0, state);    if (check_flag(&flag, "CVodeInit", 1)) return(1);       flag = CVodeSetUserData(cvode_mem, p);    if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);       flag = CVodeSVtolerances(cvode_mem, reltol, abstol);    if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);           flag = CVodeRootInit(cvode_mem, NRF, root_functions);    if (check_flag(&flag, "CVodeRootInit", 1)) return(1);    CVodeSetRootDirection(cvode_mem, rootdir);    if (check_flag(&flag, "CVodeSetRootDirection", 1)) return(1);               flag = CVDense(cvode_mem, NEQ);    if (check_flag(&flag, "CVDense", 1)) return(1);    printf(" /n Integrating iaf /n/n");    printf("#t v, /n");    PrintOutput(pout, t, state);       tout = DT;    while(1) {        flag = CVode(cvode_mem, tout, state, &t, CV_NORMAL);                if(flag == CV_ROOT_RETURN) {            /* Event detected */            flagr = CVodeGetRootInfo(cvode_mem, rootsfound);            if (check_flag(&flagr, "CVodeGetRootInfo", 1)) return(1);            PrintRootInfo(t, state, rootsfound);                      if(rootsfound[0]){                //condition_0                v = NV_Ith_S(state, 0);                NV_Ith_S(state, 0) = reset;        }                    /* Restart integration with event-corrected state */            flag = CVodeSetUserData(cvode_mem, p);            if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);        CVodeReInit(cvode_mem, t, state);        //PrintRootInfo(t, state, rootsfound);    }        else                {            PrintOutput(pout, t, state);            if(check_flag(&flag, "CVode", 1)) break;            if(flag == CV_SUCCESS) {                tout += DT;            }            if (t >= T1) break;        }    }//.........这里部分代码省略.........

int dynamixMain (int argc, char * argv[]) {  //// DECLARING VARIABLES  // Struct of parameters  PARAMETERS p;  // CVode variables  void * cvode_mem = NULL;			// pointer to block of CVode memory  N_Vector y, yout;			// arrays of populations  // arrays for energetic parameters  realtype ** V = NULL;				// pointer to k-c coupling constants  realtype * Vbridge = NULL;			// pointer to array of bridge coupling constants.  // first element [0] is Vkb1, last [Nb] is VcbN  realtype * Vnobridge = NULL;			// coupling constant when there is no bridge  //// Setting defaults for parameters to be read from input  //// done setting defaults  int flag;  realtype * k_pops = NULL;				// pointers to arrays of populations  realtype * l_pops = NULL;  realtype * c_pops = NULL;  realtype * b_pops = NULL;  realtype * ydata = NULL;				// pointer to ydata (contains all populations)  realtype * wavefunction = NULL;			// (initial) wavefunction  realtype * dm = NULL;					// density matrix  realtype * dmt = NULL;				// density matrix in time  realtype * wfnt = NULL;				// density matrix in time  realtype * k_energies = NULL;				// pointers to arrays of energies  realtype * c_energies = NULL;  realtype * b_energies = NULL;  realtype * l_energies = NULL;  realtype t0 = 0.0;				// initial time  realtype t = 0;  realtype tret = 0;					// time returned by the solver  time_t startRun;				// time at start of log  time_t endRun;					// time at end of log  struct tm * currentTime = NULL;			// time structure for localtime#ifdef DEBUG  FILE * realImaginary;				// file containing real and imaginary parts of the wavefunction#endif  FILE * log;					// log file with run times  realtype * tkprob = NULL; 				// total probability in k, l, c, b states at each timestep  realtype * tlprob = NULL;  realtype * tcprob = NULL;  realtype * tbprob = NULL;  double ** allprob = NULL;				// populations in all states at all times  realtype * times = NULL;  realtype * qd_est = NULL;  realtype * qd_est_diag = NULL;  std::string inputFile = "ins/parameters.in";			// name of input file  std::string cEnergiesInput = "ins/c_energies.in";  std::string cPopsInput = "ins/c_pops.in";  std::string bEnergiesInput = "ins/b_energies.in";  std::string VNoBridgeInput = "ins/Vnobridge.in";  std::string VBridgeInput = "ins/Vbridge.in";  std::map<const std::string, bool> outs;	// map of output file names to bool  // default output directory  p.outputDir = "outs/";  double summ = 0;			// sum variable  // ---- process command line flags ---- //  opterr = 0;  int c;  std::string insDir;  /* process command line options */  while ((c = getopt(argc, argv, "i:o:")) != -1) {    switch (c) {      case 'i':	// check that it ends in a slash	std::cerr << "[dynamix]: assigning input directory" << std::endl;	insDir = optarg;	if (strcmp(&(insDir.at(insDir.length() - 1)), "/")) {	  std::cerr << "ERROR: option -i requires argument ("	            << insDir << ") to have a trailing slash (/)." << std::endl;	  return 1;	}	else {	  // ---- assign input files ---- //	  inputFile = insDir + "parameters.in";	  cEnergiesInput = insDir + "c_energies.in";	  cPopsInput = insDir + "c_pops.in";	  bEnergiesInput = insDir + "b_energies.in";	  VNoBridgeInput = insDir + "Vnobridge.in";	  VBridgeInput = insDir + "Vbridge.in";	}	break;      case 'o':	std::cerr << "[dynamix]: assigning output directory" << std::endl;	p.outputDir = optarg;	break;      case '?':	if (optopt == 'i') {	  fprintf(stderr, "Option -%c requires a directory argument./n", optopt);//.........这里部分代码省略.........

int main(int argc, char *argv[]){  void *cvode_mem;  UserData data;  realtype t, tout;  N_Vector y;  int iout, flag;  realtype pbar[NS];  int is;   N_Vector *yS;  booleantype sensi, err_con;  int sensi_meth;  cvode_mem = NULL;  data      = NULL;  y         =  NULL;  yS        = NULL;  /* Process arguments */  ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);  /* User data structure */  data = (UserData) malloc(sizeof *data);  if (check_flag((void *)data, "malloc", 2)) return(1);  data->p[0] = RCONST(0.04);  data->p[1] = RCONST(1.0e4);  data->p[2] = RCONST(3.0e7);  /* Initial conditions */  y = N_VNew_Serial(NEQ);  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);  Ith(y,1) = Y1;  Ith(y,2) = Y2;  Ith(y,3) = Y3;  /* Create CVODES object */  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);  if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);  /* Allocate space for CVODES */  flag = CVodeMalloc(cvode_mem, f, T0, y, CV_WF, 0.0, NULL);  if (check_flag(&flag, "CVodeMalloc", 1)) return(1);  /* Use private function to compute error weights */  flag = CVodeSetEwtFn(cvode_mem, ewt, NULL);  if (check_flag(&flag, "CVodeSetEwtFn", 1)) return(1);  /* Attach user data */  flag = CVodeSetFdata(cvode_mem, data);  if (check_flag(&flag, "CVodeSetFdata", 1)) return(1);  /* Attach linear solver */  flag = CVDense(cvode_mem, NEQ);  if (check_flag(&flag, "CVDense", 1)) return(1);  flag = CVDenseSetJacFn(cvode_mem, Jac, data);  if (check_flag(&flag, "CVDenseSetJacFn", 1)) return(1);  printf("/n3-species chemical kinetics problem/n");  /* Sensitivity-related settings */  if (sensi) {    pbar[0] = data->p[0];    pbar[1] = data->p[1];    pbar[2] = data->p[2];    yS = N_VNewVectorArray_Serial(NS, NEQ);    if (check_flag((void *)yS, "N_VNewVectorArray_Serial", 0)) return(1);    for (is=0;is<NS;is++) N_VConst(ZERO, yS[is]);    flag = CVodeSensMalloc(cvode_mem, NS, sensi_meth, yS);    if(check_flag(&flag, "CVodeSensMalloc", 1)) return(1);    flag = CVodeSetSensRhs1Fn(cvode_mem, fS);    if (check_flag(&flag, "CVodeSetSensRhs1Fn", 1)) return(1);    flag = CVodeSetSensErrCon(cvode_mem, err_con);    if (check_flag(&flag, "CVodeSetSensFdata", 1)) return(1);    flag = CVodeSetSensFdata(cvode_mem, data);    if (check_flag(&flag, "CVodeSetSensFdata", 1)) return(1);    flag = CVodeSetSensParams(cvode_mem, NULL, pbar, NULL);    if (check_flag(&flag, "CVodeSetSensParams", 1)) return(1);    printf("Sensitivity: YES ");    if(sensi_meth == CV_SIMULTANEOUS)         printf("( SIMULTANEOUS +");    else       if(sensi_meth == CV_STAGGERED) printf("( STAGGERED +");      else                           printf("( STAGGERED1 +");       if(err_con) printf(" FULL ERROR CONTROL )");    else        printf(" PARTIAL ERROR CONTROL )");  } else {    printf("Sensitivity: NO ");  }  //.........这里部分代码省略.........

int main(int argc, char *argv[]){  realtype dx, reltol, abstol, t, tout, umax;  N_Vector u;  UserData data;  void *cvode_mem;  int iout, retval, my_pe, npes;  sunindextype local_N, nperpe, nrem, my_base;  long int nst;  MPI_Comm comm;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Get processor number, total number of pe's, and my_pe. */  MPI_Init(&argc, &argv);  comm = MPI_COMM_WORLD;  MPI_Comm_size(comm, &npes);  MPI_Comm_rank(comm, &my_pe);  /* Set local vector length. */  nperpe = NEQ/npes;  nrem = NEQ - npes*nperpe;  local_N = (my_pe < nrem) ? nperpe+1 : nperpe;  my_base = (my_pe < nrem) ? my_pe*local_N : my_pe*nperpe + nrem;  data = (UserData) malloc(sizeof *data);  /* Allocate data memory */  if(check_retval((void *)data, "malloc", 2, my_pe)) MPI_Abort(comm, 1);  data->comm = comm;  data->npes = npes;  data->my_pe = my_pe;  u = N_VNew_Parallel(comm, local_N, NEQ);  /* Allocate u vector */  if(check_retval((void *)u, "N_VNew", 0, my_pe)) MPI_Abort(comm, 1);  reltol = ZERO;  /* Set the tolerances */  abstol = ATOL;  dx = data->dx = XMAX/((realtype)(MX+1));  /* Set grid coefficients in data */  data->hdcoef = RCONST(1.0)/(dx*dx);  data->hacoef = RCONST(0.5)/(RCONST(2.0)*dx);  SetIC(u, dx, local_N, my_base);  /* Initialize u vector */  /* Call CVodeCreate to create the solver memory and specify the   * Adams-Moulton LMM */  cvode_mem = CVodeCreate(CV_ADAMS);  if(check_retval((void *)cvode_mem, "CVodeCreate", 0, my_pe)) MPI_Abort(comm, 1);  retval = CVodeSetUserData(cvode_mem, data);  if(check_retval(&retval, "CVodeSetUserData", 1, my_pe)) MPI_Abort(comm, 1);  /* Call CVodeInit to initialize the integrator memory and specify the   * user's right hand side function in u'=f(t,u), the inital time T0, and   * the initial dependent variable vector u. */  retval = CVodeInit(cvode_mem, f, T0, u);  if(check_retval(&retval, "CVodeInit", 1, my_pe)) return(1);  /* Call CVodeSStolerances to specify the scalar relative tolerance   * and scalar absolute tolerances */  retval = CVodeSStolerances(cvode_mem, reltol, abstol);  if (check_retval(&retval, "CVodeSStolerances", 1, my_pe)) return(1);  /* Call CVDiag to create and attach CVODE-specific diagonal linear solver */  retval = CVDiag(cvode_mem);  if(check_retval(&retval, "CVDiag", 1, my_pe)) return(1);  if (my_pe == 0) PrintIntro(npes);  umax = N_VMaxNorm(u);  if (my_pe == 0) {    t = T0;    PrintData(t, umax, 0);  }  /* In loop over output points, call CVode, print results, test for error */  for (iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {    retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL);    if(check_retval(&retval, "CVode", 1, my_pe)) break;    umax = N_VMaxNorm(u);    retval = CVodeGetNumSteps(cvode_mem, &nst);    check_retval(&retval, "CVodeGetNumSteps", 1, my_pe);    if (my_pe == 0) PrintData(t, umax, nst);  }  if (my_pe == 0)    PrintFinalStats(cvode_mem);  /* Print some final statistics */  N_VDestroy_Parallel(u);        /* Free the u vector */  CVodeFree(&cvode_mem);         /* Free the integrator memory */  free(data);                    /* Free user data */  MPI_Finalize();  return(0);//.........这里部分代码省略.........

static int Problem2(void){  realtype reltol=RTOL, abstol=ATOL, t, tout, er, erm, ero;  int miter, flag, temp_flag, nerr=0;  N_Vector y;  void *cvode_mem;  booleantype firstrun;  int qu, iout;  realtype hu;  y = NULL;  cvode_mem = NULL;  y = N_VNew_Serial(P2_NEQ);  if(check_flag((void *)y, "N_VNew", 0)) return(1);  PrintIntro2();  cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);  for (miter=FUNC; miter <= BAND_DQ; miter++) {    if ((miter==DENSE_USER) || (miter==DENSE_DQ)) continue;    ero = ZERO;    N_VConst(ZERO, y);    NV_Ith_S(y,0) = ONE;          firstrun = (miter==FUNC);    if (firstrun) {      flag = CVodeMalloc(cvode_mem, f2, P2_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeMalloc", 1)) return(1);    } else {      flag = CVodeSetIterType(cvode_mem, CV_NEWTON);      if(check_flag(&flag, "CVodeSetIterType", 1)) ++nerr;      flag = CVodeReInit(cvode_mem, f2, P2_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeReInit", 1)) return(1);    }          flag = PrepareNextRun(cvode_mem, CV_ADAMS, miter, P2_MU, P2_ML);    if(check_flag(&flag, "PrepareNextRun", 1)) return(1);    PrintHeader2();    for(iout=1, tout=P2_T1; iout <= P2_NOUT; iout++, tout*=P2_TOUT_MULT) {      flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);      check_flag(&flag, "CVode", 1);      erm = MaxError(y, t);      temp_flag = CVodeGetLastOrder(cvode_mem, &qu);      if(check_flag(&temp_flag, "CVodeGetLastOrder", 1)) ++nerr;      temp_flag = CVodeGetLastStep(cvode_mem, &hu);      if(check_flag(&temp_flag, "CVodeGetLastStep", 1)) ++nerr;      PrintOutput2(t, erm, qu, hu);      if (flag != CV_SUCCESS) {        nerr++;        break;      }      er = erm / abstol;        if (er > ero) ero = er;        if (er > P2_TOL_FACTOR) {          nerr++;          PrintErrOutput(P2_TOL_FACTOR);        }    }        PrintFinalStats(cvode_mem, miter, ero);  }  CVodeFree(cvode_mem);  cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);  for (miter=FUNC; miter <= BAND_DQ; miter++) {    if ((miter==DENSE_USER) || (miter==DENSE_DQ)) continue;    ero = ZERO;    N_VConst(ZERO, y);    NV_Ith_S(y,0) = ONE;          firstrun = (miter==FUNC);    if (firstrun) {      flag = CVodeMalloc(cvode_mem, f2, P2_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeMalloc", 1)) return(1);    } else {      flag = CVodeSetIterType(cvode_mem, CV_NEWTON);      if(check_flag(&flag, "CVodeSetIterType", 1)) ++nerr;      flag = CVodeReInit(cvode_mem, f2, P2_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeReInit", 1)) return(1);    }    flag = PrepareNextRun(cvode_mem, CV_BDF, miter, P2_MU, P2_ML);    if(check_flag(&flag, "PrepareNextRun", 1)) return(1);    PrintHeader2();          for(iout=1, tout=P2_T1; iout <= P2_NOUT; iout++, tout*=P2_TOUT_MULT) {      flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);      check_flag(&flag, "CVode", 1);      erm = MaxError(y, t);      temp_flag = CVodeGetLastOrder(cvode_mem, &qu);      if(check_flag(&temp_flag, "CVodeGetLastOrder", 1)) ++nerr;//.........这里部分代码省略.........

static int Problem1(void){  realtype reltol=RTOL, abstol=ATOL, t, tout, ero, er;  int miter, flag, temp_flag, iout, nerr=0;  N_Vector y;  void *cvode_mem;  booleantype firstrun;  int qu;  realtype hu;  y = NULL;  cvode_mem = NULL;  y = N_VNew_Serial(P1_NEQ);  if(check_flag((void *)y, "N_VNew_Serial", 0)) return(1);  PrintIntro1();  cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);  for (miter=FUNC; miter <= DIAG; miter++) {    ero = ZERO;    NV_Ith_S(y,0) = TWO;    NV_Ith_S(y,1) = ZERO;    firstrun = (miter==FUNC);    if (firstrun) {      flag = CVodeMalloc(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeMalloc", 1)) return(1);    } else {      flag = CVodeSetIterType(cvode_mem, CV_NEWTON);      if(check_flag(&flag, "CVodeSetIterType", 1)) ++nerr;      flag = CVodeReInit(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeReInit", 1)) return(1);    }          flag = PrepareNextRun(cvode_mem, CV_ADAMS, miter, 0, 0);    if(check_flag(&flag, "PrepareNextRun", 1)) return(1);    PrintHeader1();    for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {      flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);      check_flag(&flag, "CVode", 1);      temp_flag = CVodeGetLastOrder(cvode_mem, &qu);      if(check_flag(&temp_flag, "CVodeGetLastOrder", 1)) ++nerr;      temp_flag = CVodeGetLastStep(cvode_mem, &hu);      if(check_flag(&temp_flag, "CVodeGetLastStep", 1)) ++nerr;      PrintOutput1(t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);      if (flag != CV_SUCCESS) {        nerr++;        break;      }      if (iout%2 == 0) {        er = ABS(NV_Ith_S(y,0)) / abstol;        if (er > ero) ero = er;        if (er > P1_TOL_FACTOR) {          nerr++;	  PrintErrOutput(P1_TOL_FACTOR);        }      }    }        PrintFinalStats(cvode_mem, miter, ero);  }  CVodeFree(cvode_mem);  cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);  for (miter=FUNC; miter <= DIAG; miter++) {    ero = ZERO;    NV_Ith_S(y,0) = TWO;    NV_Ith_S(y,1) = ZERO;          firstrun = (miter==FUNC);    if (firstrun) {      flag = CVodeMalloc(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeMalloc", 1)) return(1);    } else {      flag = CVodeSetIterType(cvode_mem, CV_NEWTON);      if(check_flag(&flag, "CVodeSetIterType", 1)) ++nerr;      flag = CVodeReInit(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);      if(check_flag(&flag, "CVodeReInit", 1)) return(1);    }          flag = PrepareNextRun(cvode_mem, CV_BDF, miter, 0, 0);         if(check_flag(&flag, "PrepareNextRun", 1)) return(1);    PrintHeader1();          for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {      flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);      check_flag(&flag, "CVode", 1);      temp_flag = CVodeGetLastOrder(cvode_mem, &qu);      if(check_flag(&temp_flag, "CVodeGetLastOrder", 1)) ++nerr;      temp_flag = CVodeGetLastStep(cvode_mem, &hu);      if(check_flag(&temp_flag, "CVodeGetLastStep", 1)) ++nerr;      PrintOutput1(t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);//.........这里部分代码省略.........

int main(){  realtype abstol, reltol, t, tout;  N_Vector u;  UserData data;  void *cvode_mem;  int flag, iout, jpre;  long int ml, mu;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Allocate and initialize u, and set problem data and tolerances */   u = N_VNew_Serial(NEQ);  if(check_flag((void *)u, "N_VNew_Serial", 0)) return(1);  data = (UserData) malloc(sizeof *data);  if(check_flag((void *)data, "malloc", 2)) return(1);  InitUserData(data);  SetInitialProfiles(u, data->dx, data->dy);  abstol = ATOL;   reltol = RTOL;  /* 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, data);  if(check_flag(&flag, "CVodeSetUserData", 1)) return(1);  /* Call CVodeInit to initialize the integrator memory and specify the   * user's right hand side function in u'=f(t,u), the inital time T0, and   * the initial dependent variable vector u. */  flag = CVodeInit(cvode_mem, f, T0, u);  if(check_flag(&flag, "CVodeInit", 1)) return(1);  /* Call CVodeSStolerances to specify the scalar relative tolerance   * and scalar absolute tolerances */  flag = CVodeSStolerances(cvode_mem, reltol, abstol);  if (check_flag(&flag, "CVodeSStolerances", 1)) return(1);  /* Call CVSpgmr to specify the linear solver CVSPGMR      with left preconditioning and the maximum Krylov dimension maxl */  flag = CVSpgmr(cvode_mem, PREC_LEFT, 0);  if(check_flag(&flag, "CVSpgmr", 1)) return(1);  /* Call CVBandPreInit to initialize band preconditioner */  ml = mu = 2;  flag = CVBandPrecInit(cvode_mem, NEQ, mu, ml);  if(check_flag(&flag, "CVBandPrecInit", 0)) return(1);  PrintIntro(mu, ml);  /* Loop over jpre (= PREC_LEFT, PREC_RIGHT), and solve the problem */  for (jpre = PREC_LEFT; jpre <= PREC_RIGHT; jpre++) {        /* On second run, re-initialize u, the solver, and CVSPGMR */        if (jpre == PREC_RIGHT) {            SetInitialProfiles(u, data->dx, data->dy);            flag = CVodeReInit(cvode_mem, T0, u);      if(check_flag(&flag, "CVodeReInit", 1)) return(1);      flag = CVSpilsSetPrecType(cvode_mem, PREC_RIGHT);      check_flag(&flag, "CVSpilsSetPrecType", 1);            printf("/n/n-------------------------------------------------------");      printf("------------/n");    }        printf("/n/nPreconditioner type is:  jpre = %s/n/n",           (jpre == PREC_LEFT) ? "PREC_LEFT" : "PREC_RIGHT");        /* In loop over output points, call CVode, print results, test for error */        for (iout = 1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {      flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);      check_flag(&flag, "CVode", 1);      PrintOutput(cvode_mem, u, t);      if (flag != CV_SUCCESS) {        break;      }    }        /* Print final statistics */        PrintFinalStats(cvode_mem);      } /* End of jpre loop */  /* Free memory */  N_VDestroy_Serial(u);  free(data);  CVodeFree(&cvode_mem);//.........这里部分代码省略.........

int main(int argc, char *argv[]){  void *cvode_mem;  UserData data;  realtype abstol, reltol, t, tout;  N_Vector y;  int iout, flag;  realtype *pbar;  int is, *plist;  N_Vector *uS;  booleantype sensi, err_con;  int sensi_meth;  pbar = NULL;  plist = NULL;  uS = NULL;  y = NULL;  data = NULL;  cvode_mem = NULL;  /* Process arguments */  ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);  /* Problem parameters */  data = AllocUserData();  if(check_flag((void *)data, "AllocUserData", 2)) return(1);  InitUserData(data);  /* Initial states */  y = N_VNew_Serial(NEQ);  if(check_flag((void *)y, "N_VNew_Serial", 0)) return(1);  SetInitialProfiles(y, data->dx, data->dz);    /* Tolerances */  abstol=ATOL;   reltol=RTOL;  /* Create CVODES object */  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);  flag = CVodeSetFdata(cvode_mem, data);  if(check_flag(&flag, "CVodeSetFdata", 1)) return(1);  flag = CVodeSetMaxNumSteps(cvode_mem, 2000);  if(check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return(1);  /* Allocate CVODES memory */  flag = CVodeMalloc(cvode_mem, f, T0, y, CV_SS, reltol, &abstol);  if(check_flag(&flag, "CVodeMalloc", 1)) return(1);  /* Attach CVSPGMR linear solver */  flag = CVSpgmr(cvode_mem, PREC_LEFT, 0);  if(check_flag(&flag, "CVSpgmr", 1)) return(1);  flag = CVSpilsSetPreconditioner(cvode_mem, Precond, PSolve, data);  if(check_flag(&flag, "CVSpilsSetPreconditioner", 1)) return(1);  printf("/n2-species diurnal advection-diffusion problem/n");  /* Forward sensitivity analysis */  if(sensi) {    plist = (int *) malloc(NS * sizeof(int));    if(check_flag((void *)plist, "malloc", 2)) return(1);    for(is=0; is<NS; is++) plist[is] = is;    pbar = (realtype *) malloc(NS * sizeof(realtype));    if(check_flag((void *)pbar, "malloc", 2)) return(1);    for(is=0; is<NS; is++) pbar[is] = data->p[plist[is]];    uS = N_VCloneVectorArray_Serial(NS, y);    if(check_flag((void *)uS, "N_VCloneVectorArray_Serial", 0)) return(1);    for(is=0;is<NS;is++)      N_VConst(ZERO,uS[is]);    flag = CVodeSensMalloc(cvode_mem, NS, sensi_meth, uS);    if(check_flag(&flag, "CVodeSensMalloc", 1)) return(1);    flag = CVodeSetSensErrCon(cvode_mem, err_con);    if(check_flag(&flag, "CVodeSetSensErrCon", 1)) return(1);    flag = CVodeSetSensRho(cvode_mem, ZERO);    if(check_flag(&flag, "CVodeSetSensRho", 1)) return(1);    flag = CVodeSetSensParams(cvode_mem, data->p, pbar, plist);    if(check_flag(&flag, "CVodeSetSensParams", 1)) return(1);    printf("Sensitivity: YES ");    if(sensi_meth == CV_SIMULTANEOUS)         printf("( SIMULTANEOUS +");    else       if(sensi_meth == CV_STAGGERED) printf("( STAGGERED +");      else                           printf("( STAGGERED1 +");       if(err_con) printf(" FULL ERROR CONTROL )");    else        printf(" PARTIAL ERROR CONTROL )");      } else {//.........这里部分代码省略.........

int main(void){  realtype dx, dy, reltol, abstol, t, tout, umax;  N_Vector u;  UserData data;  void *cvode_mem;  int iout, flag;  long int nst;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Create a serial vector */  u = N_VNew_Serial(NEQ);  /* Allocate u vector */  if(check_flag((void*)u, "N_VNew_Serial", 0)) return(1);  reltol = ZERO;  /* Set the tolerances */  abstol = ATOL;  data = (UserData) malloc(sizeof *data);  /* Allocate data memory */  if(check_flag((void *)data, "malloc", 2)) return(1);  dx = data->dx = XMAX/(MX+1);  /* Set grid coefficients in data */  dy = data->dy = YMAX/(MY+1);  data->hdcoef = ONE/(dx*dx);  data->hacoef = HALF/(TWO*dx);  data->vdcoef = ONE/(dy*dy);  SetIC(u, data);  /* Initialize u vector */  /* 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);  /* Call CVodeInit to initialize the integrator memory and specify the   * user's right hand side function in u'=f(t,u), the inital time T0, and   * the initial dependent variable vector u. */  flag = CVodeInit(cvode_mem, f, T0, u);  if(check_flag(&flag, "CVodeInit", 1)) return(1);  /* Call CVodeSStolerances to specify the scalar relative tolerance   * and scalar absolute tolerance */  flag = CVodeSStolerances(cvode_mem, reltol, abstol);  if (check_flag(&flag, "CVodeSStolerances", 1)) return(1);  /* Set the pointer to user-defined data */  flag = CVodeSetUserData(cvode_mem, data);  if(check_flag(&flag, "CVodeSetUserData", 1)) return(1);  /* Call CVLapackBand to specify the CVBAND band linear solver */  flag = CVLapackBand(cvode_mem, NEQ, MY, MY);  if(check_flag(&flag, "CVLapackBand", 1)) return(1);  /* Set the user-supplied Jacobian routine Jac */  flag = CVDlsSetBandJacFn(cvode_mem, Jac);  if(check_flag(&flag, "CVDlsSetBandJacFn", 1)) return(1);  /* In loop over output points: call CVode, print results, test for errors */  umax = N_VMaxNorm(u);  PrintHeader(reltol, abstol, umax);  for(iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {    flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);    if(check_flag(&flag, "CVode", 1)) break;    umax = N_VMaxNorm(u);    flag = CVodeGetNumSteps(cvode_mem, &nst);    check_flag(&flag, "CVodeGetNumSteps", 1);    PrintOutput(t, umax, nst);  }  PrintFinalStats(cvode_mem);  /* Print some final statistics   */  N_VDestroy_Serial(u);   /* Free the u vector */  CVodeFree(&cvode_mem);  /* Free the integrator memory */  free(data);             /* Free the user data */  return(0);}