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

int CVodeCreateB(void *cvadj_mem, int lmmB, int iterB){  CVadjMem ca_mem;  void *cvode_mem;  if (cvadj_mem == NULL) return(CV_ADJMEM_NULL);  ca_mem = (CVadjMem) cvadj_mem;  cvode_mem = CVodeCreate(lmmB, iterB);  if (cvode_mem == NULL) return(CV_MEM_FAIL);  ca_mem->cvb_mem = (CVodeMem) cvode_mem;  return(CV_SUCCESS);}

cvode_mem *SOLVER(cvode, init, TARGET, SIMENGINE_STORAGE, solver_props *props){  cvode_mem *mem = (cvode_mem*) malloc(props->num_models*sizeof(cvode_mem));  unsigned int modelid;  // Only need to create this buffer in the first memory space as we are using this only for scratch  // and outside of the CVODE solver to compute the last outputs  mem[0].k1 = (CDATAFORMAT*)malloc(props->statesize*props->num_models*sizeof(CDATAFORMAT));  for(modelid=0; modelid<props->num_models; modelid++){    // Set the modelid on a per memory structure basis    mem[modelid].modelid = modelid;    // Store solver properties    mem[modelid].props = props;    // Create intial value vector    // This is done to avoid having the change the internal indexing within the flows and for the output_buffer    mem[modelid].y0 = N_VMake_Serial(props->statesize, &(props->model_states[modelid*props->statesize]));    // Create data structure for solver    mem[modelid].cvmem = CVodeCreate(CV_BDF, CV_NEWTON);    // Initialize CVODE    if(CVodeInit(mem[modelid].cvmem, user_fun_wrapper, props->starttime, ((N_Vector)(mem[modelid].y0))) != CV_SUCCESS){      fprintf(stderr, "Couldn't initialize CVODE");    }    // Set solver tolerances    if(CVodeSStolerances(mem[modelid].cvmem, props->reltol, props->abstol) != CV_SUCCESS){      fprintf(stderr, "Could not set CVODE tolerances");    }    // Set linear solver    if(CVDense(mem[modelid].cvmem, mem[modelid].props->statesize) != CV_SUCCESS){      fprintf(stderr, "Could not set CVODE linear solver");    }    // Set user data to contain pointer to memory structure for use in model_flows    if(CVodeSetUserData(mem[modelid].cvmem, &mem[modelid]) != CV_SUCCESS){      fprintf(stderr, "CVODE failed to initialize user data");    }  }  return mem;}

ode_solver*	ode_solver_alloc(ode_model* model){    ode_solver* solver = (ode_solver*) malloc( sizeof(ode_solver) );				/* alloc */  if (solver == 0){    /* TODO: write a proper error handler */    fprintf(stderr,"malloc failed to allocate memory for ode_solver/n");    return 0;  }    solver->cvode_mem = CVodeCreate(CV_BDF,CV_NEWTON);								/* alloc */  if( solver->cvode_mem == 0){    /* TODO: write a proper error handler */    fprintf(stderr,"CVodeCreate failed to allocate memory in ode_solver for cvode_mem./n");        free(solver);    return 0;  }    int P = ode_model_getP(model);  solver->params = (double*) malloc( sizeof(double) * P );						/* alloc */  if( solver->params == 0 ){    /* TODO: write a proper error handler */    fprintf(stderr,"malloc failed to allocate memory in ode_solver for params./n");    CVodeFree(solver->cvode_mem);    free(solver);    return 0;  }  ode_model_get_default_params(model, solver->params, P);	    int N = ode_model_getN(model);  solver->odeModel = model;  solver->y = N_VNewEmpty_Serial(N);												/* alloc */  NV_DATA_S(solver->y) = solver->odeModel->v;  solver->yS = 0;    return solver;}

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[]){  UserData data;  void *cvode_mem;  SUNMatrix A, AB;  SUNLinearSolver LS, LSB;  realtype dx, dy, reltol, abstol, t;  N_Vector u;  int indexB;  realtype reltolB, abstolB;  N_Vector uB;    int retval, ncheck;  data = NULL;  cvode_mem = NULL;  u = uB = NULL;  LS = LSB = NULL;  A = AB = NULL;  /* Allocate and initialize user data memory */  data = (UserData) malloc(sizeof *data);  if(check_retval((void *)data, "malloc", 2)) return(1);  dx = data->dx = XMAX/(MX+1);  dy = data->dy = YMAX/(MY+1);  data->hdcoef = ONE/(dx*dx);  data->hacoef = RCONST(1.5)/(TWO*dx);  data->vdcoef = ONE/(dy*dy);  /* Set the tolerances for the forward integration */  reltol = ZERO;  abstol = ATOL;  /* Allocate u vector */  u = N_VNew_Serial(NEQ);  if(check_retval((void *)u, "N_VNew", 0)) return(1);  /* Initialize u vector */  SetIC(u, data);  /* Create and allocate CVODES memory for forward run */  printf("/nCreate and allocate CVODES memory for forward runs/n");  cvode_mem = CVodeCreate(CV_BDF);  if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);  retval = CVodeSetUserData(cvode_mem, data);  if(check_retval(&retval, "CVodeSetUserData", 1)) return(1);  retval = CVodeInit(cvode_mem, f, T0, u);  if(check_retval(&retval, "CVodeInit", 1)) return(1);  retval = CVodeSStolerances(cvode_mem, reltol, abstol);  if(check_retval(&retval, "CVodeSStolerances", 1)) return(1);  /* Create banded SUNMatrix for the forward problem */  A = SUNBandMatrix(NEQ, MY, MY);  if(check_retval((void *)A, "SUNBandMatrix", 0)) return(1);  /* Create banded SUNLinearSolver for the forward problem */  LS = SUNLinSol_Band(u, A);  if(check_retval((void *)LS, "SUNLinSol_Band", 0)) return(1);  /* Attach the matrix and linear solver */  retval = CVodeSetLinearSolver(cvode_mem, LS, A);  if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);  /* Set the user-supplied Jacobian routine for the forward problem */  retval = CVodeSetJacFn(cvode_mem, Jac);  if(check_retval(&retval, "CVodeSetJacFn", 1)) return(1);  /* Allocate global memory */  printf("/nAllocate global memory/n");  retval = CVodeAdjInit(cvode_mem, NSTEP, CV_HERMITE);  if(check_retval(&retval, "CVodeAdjInit", 1)) return(1);  /* Perform forward run */  printf("/nForward integration/n");  retval = CVodeF(cvode_mem, TOUT, u, &t, CV_NORMAL, &ncheck);  if(check_retval(&retval, "CVodeF", 1)) return(1);  printf("/nncheck = %d/n", ncheck);  /* Set the tolerances for the backward integration */  reltolB = RTOLB;  abstolB = ATOL;  /* Allocate uB */  uB = N_VNew_Serial(NEQ);  if(check_retval((void *)uB, "N_VNew", 0)) return(1);  /* Initialize uB = 0 *///.........这里部分代码省略.........

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);//.........这里部分代码省略.........

static int simulate_nmr_pulse(struct bloch_sim *bs){    N_Vector M = NULL;    M = N_VNew_Serial(3 * bs->num_cells);    if (check_flag((void *)M, "N_VNew_Serial", 0)) return(1);    int i;    /* Set initial (t=0) magnetization conditions */    for(i=0; i < bs->num_cells; ++i)    {        X(M,i) = 0.0;        Y(M,i) = 0.0;        Z(M,i) = bs->cell_frequencies[i] / bs->w_avg;    }    realtype reltol = RCONST(1.0e-14);    realtype abstol = RCONST(1.0e-14);    void *cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);    if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);    int flag;    //TODO: check if flag should be pointer;    flag = CVodeInit(cvode_mem, bloch_equations, 0.0, M);    if (check_flag(&flag, "CVodeInit", 1)) return(1);    flag = CVodeSetUserData(cvode_mem, bs);    if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);    flag = CVodeSStolerances(cvode_mem, reltol, abstol);    if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);    flag = CVDense(cvode_mem, 3 * bs->num_cells);    if (check_flag(&flag, "CVDense", 1)) return(1);    flag = CVDlsSetDenseJacFn(cvode_mem, bloch_jacobian);    if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);    ///////////////////////////    // PI/2 PULSE SIMULATION //    ///////////////////////////    bs->rf_on = 1;    flag = CVodeSetStopTime(cvode_mem, bs->pi2_duration);    if (check_flag(&flag, "CVodeSetStopTime", 1)) return 1;    realtype time_reached;    flag = CVode(cvode_mem, bs->pi2_duration, M, &time_reached, CV_NORMAL);    if (flag != CV_SUCCESS)    {        printf("ERROR: Failed to simulate Pi/2 pulse/n");    }    // {{{ PI2 PULSE DEBUG STATEMENTS    if (DEBUG)    {        printf("/n");        printf("#####################################/n");        printf("### PI/2 PULSE SIMULATION DETAILS ###/n");        printf("#####################################/n");        printf("/n");        printf("TIME REACHED: %.15e/n", time_reached);        printf("TIME REACHED - PI2_DURATION: %.15e/n", time_reached - bs->pi2_duration);        printf("/n");        printf("MAGNETIZATION AT END OF PI/2 PULSE:/n");        for (i = 0; i<bs->num_cells; ++i)        {            printf("CELL %d: %.4e, %.4e, %.4e/n", i, X(M,i), Y(M,i), Z(M,i));        }    }    // }}}    ////////////////////    // FID SIMULATION //    ////////////////////    bs->rf_on = 0;    if (DEBUG)    {        printf("##############################/n");        printf("### FID SIMULATION DETAILS ###/n");        printf("##############################/n");    }    flag = CVodeReInit(cvode_mem, 0.0, M);    if (check_flag(&flag, "CVodeReInit", 1)) return(1);    time_reached = 0.0;    flag = CVodeSetStopTime(cvode_mem, bs->fid_duration);    if (check_flag(&flag, "CVodeSetStopTime", 1)) return 1;    flag = CVodeRootInit(cvode_mem, 1, bloch_root);    if (check_flag(&flag, "CVodeRootInit", 1)) return 1;    realtype time_desired, M_FID_X;    while (time_reached < bs->fid_duration)    {//.........这里部分代码省略.........

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[]){  ProblemData d;  MPI_Comm comm;  int npes, npes_needed;  int myId;   long int neq, l_neq;  void *cvode_mem;  N_Vector y, q;  realtype abstol, reltol, abstolQ, reltolQ;  long int mudq, mldq, mukeep, mlkeep;  int indexB;  N_Vector yB, qB;  realtype abstolB, reltolB, abstolQB, reltolQB;  long int mudqB, mldqB, mukeepB, mlkeepB;  realtype tret, *qdata, G;  int ncheckpnt, flag;  booleantype output;  /* Initialize MPI and set Ids */  MPI_Init(&argc, &argv);  comm = MPI_COMM_WORLD;  MPI_Comm_rank(comm, &myId);  /* Check number of processes */  npes_needed = NPX * NPY;#ifdef USE3D  npes_needed *= NPZ;#endif  MPI_Comm_size(comm, &npes);  if (npes_needed != npes) {    if (myId == 0)      fprintf(stderr,"I need %d processes but I only got %d/n",              npes_needed, npes);    MPI_Abort(comm, EXIT_FAILURE);  }  /* Test if matlab output is requested */  if (argc > 1) output = TRUE;  else          output = FALSE;  /* Allocate and set problem data structure */  d = (ProblemData) malloc(sizeof *d);  SetData(d, comm, npes, myId, &neq, &l_neq);    if (myId == 0) PrintHeader();  /*--------------------------     Forward integration phase    --------------------------*/  /* Allocate space for y and set it with the I.C. */  y = N_VNew_Parallel(comm, l_neq, neq);  N_VConst(ZERO, y);    /* Allocate and initialize qB (local contribution to cost) */  q = N_VNew_Parallel(comm, 1, npes);   N_VConst(ZERO, q);  /* Create CVODES object, attach user data, and allocate space */  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);  flag = CVodeSetUserData(cvode_mem, d);  flag = CVodeInit(cvode_mem, f, ti, y);  abstol = ATOL;    reltol = RTOL;     flag = CVodeSStolerances(cvode_mem, reltol, abstol);  /* attach linear solver */  flag = CVSpgmr(cvode_mem, PREC_LEFT, 0);    /* Attach preconditioner and linear solver modules */  mudq = mldq = d->l_m[0]+1;  mukeep = mlkeep = 2;    flag = CVBBDPrecInit(cvode_mem, l_neq, mudq, mldq,                        mukeep, mlkeep, ZERO,                       f_local, NULL);    /* Initialize quadrature calculations */  abstolQ = ATOL_Q;  reltolQ = RTOL_Q;  flag = CVodeQuadInit(cvode_mem, fQ, q);  flag = CVodeQuadSStolerances(cvode_mem, reltolQ, abstolQ);  flag = CVodeSetQuadErrCon(cvode_mem, TRUE);  /* Allocate space for the adjoint calculation */  flag = CVodeAdjInit(cvode_mem, STEPS, CV_HERMITE);  /* Integrate forward in time while storing check points */  if (myId == 0) printf("Begin forward integration... ");  flag = CVodeF(cvode_mem, tf, y, &tret, CV_NORMAL, &ncheckpnt);  if (myId == 0) printf("done. ");   /* Extract quadratures *///.........这里部分代码省略.........

//.........这里部分代码省略.........    }    case BDF:    {      mlmm = CV_BDF;      break;    }    default:    {      ERROR("CreateIntegrator","Invalid multistep method choice/n");      DestroyIntegrator(&integrator);      return(NULL);    }  }  switch (simulationGetIterationMethod(integrator->simulation))  {    case FUNCTIONAL:    {      miter = CV_FUNCTIONAL;      break;    }    case NEWTON:    {      miter = CV_NEWTON;      break;    }    default:    {      ERROR("CreateIntegrator","Invalid iteration method choice/n");      DestroyIntegrator(&integrator);      return(NULL);    }  }  /*      Call CVodeCreate to create the solver memory:          A pointer to the integrator problem memory is returned and     stored in cvode_mem.  */  integrator->cvode_mem = CVodeCreate(mlmm,miter);  if (check_flag((void *)(integrator->cvode_mem),"CVodeCreate",0))  {    DestroyIntegrator(&integrator);    return(NULL);  }    /*      Call CVodeMalloc to initialize the integrator memory:           cvode_mem is the pointer to the integrator memory returned by CVodeCreate     f         is the user's right hand side function in y'=f(t,y)     tStart    is the initial time     y         is the initial dependent variable vector     CV_SS     specifies scalar relative and absolute tolerances     reltol    the scalar relative tolerance     &abstol   is the absolute tolerance vector  */  flag = CVodeInit(integrator->cvode_mem,f,    simulationGetBvarStart(integrator->simulation),integrator->y);  if (check_flag(&flag,"CVodeMalloc",1))  {	  DestroyIntegrator(&integrator);      return(NULL);  }  double* atol = simulationGetATol(integrator->simulation);  if (simulationGetATolLength(integrator->simulation) == 1)  {

int main(){  realtype abstol, reltol, t, tout;  N_Vector u;  UserData data;  void *cvode_mem;  int iout, flag;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Allocate memory, and set problem data, initial values, tolerances */   u = N_VNew_Serial(NEQ);  if(check_flag((void *)u, "N_VNew_Serial", 0)) return(1);  data = AllocUserData();  if(check_flag((void *)data, "AllocUserData", 2)) return(1);  InitUserData(data);  SetInitialProfiles(u, data->dx, data->dy);  abstol=ATOL;   reltol=RTOL;  /* Call CvodeCreate to create the solver memory      CV_BDF     specifies the Backward Differentiation Formula     CV_NEWTON  specifies a Newton iteration     A pointer to the integrator memory is returned and stored in cvode_mem. */  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 = CVodeSetFdata(cvode_mem, data);  if(check_flag(&flag, "CVodeSetFdata", 1)) return(1);  /* Call CVodeMalloc to initialize the integrator memory:      f       is the user's right hand side function in u'=f(t,u)     T0      is the initial time     u       is the initial dependent variable vector     CV_SS   specifies scalar relative and absolute tolerances     reltol  is the relative tolerance     &abstol is a pointer to the scalar absolute tolerance      */  flag = CVodeMalloc(cvode_mem, f, T0, u, CV_SS, reltol, &abstol);  if(check_flag(&flag, "CVodeMalloc", 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);  /* Set modified Gram-Schmidt orthogonalization, preconditioner      setup and solve routines Precond and PSolve, and the pointer      to the user-defined block data */  flag = CVSpgmrSetGSType(cvode_mem, MODIFIED_GS);  if(check_flag(&flag, "CVSpgmrSetGSType", 1)) return(1);  flag = CVSpgmrSetPreconditioner(cvode_mem, Precond, PSolve, data);  if(check_flag(&flag, "CVSpgmrSetPreconditioner", 1)) return(1);  /* In loop over output points, call CVode, print results, test for error */  printf(" /n2-species diurnal advection-diffusion problem/n/n");  for (iout=1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {    flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);    PrintOutput(cvode_mem, u, t);    if(check_flag(&flag, "CVode", 1)) break;  }  PrintFinalStats(cvode_mem);  /* Free memory */  N_VDestroy_Serial(u);  FreeUserData(data);  CVodeFree(cvode_mem);  return(0);}

//.........这里部分代码省略.........  if (p.wavefunction) {    y = N_VMake_Serial(2*p.NEQ, wavefunction);  }  else {    y = N_VMake_Serial(2*p.NEQ2, dm);  }  // put in t = 0 information  if (! p.wavefunction) {    updateDM(y, dmt, 0, &p);  }  else {    updateWfn(y, wfnt, 0, &p);  }  // the vector yout has the same dimensions as y  yout = N_VClone(y);#ifdef DEBUG  realImaginary = fopen("real_imaginary.out", "w");#endif  // Make plot files  makePlots(outs, &p);  // only do propagation if not just making plots  if (! p.justPlots) {    // Make outputs independent of time propagation    computeGeneralOutputs(outs, &p);    // create CVode object    // this is a stiff problem, I guess?#ifdef DEBUG    std::cout << "/nCreating cvode_mem object./n";#endif    cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);    flag = CVodeSetUserData(cvode_mem, (void *) &p);#ifdef DEBUG    std::cout << "/nInitializing CVode solver./n";#endif    // initialize CVode solver //    if (p.wavefunction) {      //flag = CVodeInit(cvode_mem, &RHS_WFN, t0, y);      flag = CVodeInit(cvode_mem, &RHS_WFN_SPARSE, t0, y);    }    else {      if (p.kinetic) {	flag = CVodeInit(cvode_mem, &RHS_DM_RELAX, t0, y);      }      else if (p.rta) {	flag = CVodeInit(cvode_mem, &RHS_DM_RTA, t0, y);	//flag = CVodeInit(cvode_mem, &RHS_DM_RTA_BLAS, t0, y);      }      else if (p.dephasing) {	flag = CVodeInit(cvode_mem, &RHS_DM_dephasing, t0, y);      }      else {	//flag = CVodeInit(cvode_mem, &RHS_DM, t0, y);	flag = CVodeInit(cvode_mem, &RHS_DM_BLAS, t0, y);      }    }#ifdef DEBUG    std::cout << "/nSpecifying integration tolerances./n";#endif    // specify integration tolerances //

int do_integrate(float t_start, float t_stop, int n_points){  realtype reltol, t;   N_Vector y, abstol;  void *cvode_mem;  int flag, flagr, iout;  y = abstol = NULL;  cvode_mem = NULL;  /* Create serial vector of length NEQ for I.C. and abstol */  y = N_VNew_Serial(NEQ);  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);  abstol = N_VNew_Serial(NEQ);   if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(1);    //Setup the initial state values:  setup_initial_states(y);  setup_tolerances(abstol);  /* Set the scalar relative tolerance */  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);    /* 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, t_start, 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, NEQ);  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.*/  printf(" /n3-species kinetics problem/n/n");  iout = 0;    int i=0;  for(i=1;i< n_points; i++)  {    float t_next = t_start + (t_stop-t_start)/n_points * i;    printf("Advancing to: %f", t_next);    flag = CVode(cvode_mem, t_next, y, &t, CV_NORMAL);    loop_function(t,y);    PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));    //printf("MH: %d %f",iout,t_next);    if (check_flag(&flag, "CVode", 1)) break;    assert(flag==CV_SUCCESS);  }  /* Print some final statistics */  PrintFinalStats(cvode_mem);  /* Free y and abstol vectors */  N_VDestroy_Serial(y);  N_VDestroy_Serial(abstol);  /* Free integrator memory */  CVodeFree(&cvode_mem);  return(0);}

int main(int argc, char *argv[]){  realtype abstol, reltol, t, tout;  N_Vector u;  UserData data;  PreconData predata;  void *cvode_mem;  int iout, flag, my_pe, npes;  long int neq, local_N;  MPI_Comm comm;  u = NULL;  data = NULL;  predata = NULL;  cvode_mem = NULL;  /* Set problem size neq */  neq = NVARS*MX*MY;  /* Get processor number and total number of pe's */  MPI_Init(&argc, &argv);  comm = MPI_COMM_WORLD;  MPI_Comm_size(comm, &npes);  MPI_Comm_rank(comm, &my_pe);  if (npes != NPEX*NPEY) {    if (my_pe == 0)      fprintf(stderr, "/nMPI_ERROR(0): npes = %d is not equal to NPEX*NPEY = %d/n/n",	      npes,NPEX*NPEY);    MPI_Finalize();    return(1);  }  /* Set local length */  local_N = NVARS*MXSUB*MYSUB;  /* Allocate and load user data block; allocate preconditioner block */  data = (UserData) malloc(sizeof *data);  if (check_flag((void *)data, "malloc", 2, my_pe)) MPI_Abort(comm, 1);  InitUserData(my_pe, comm, data);  predata = AllocPreconData (data);  /* Allocate u, and set initial values and tolerances */   u = N_VNew_Parallel(comm, local_N, neq);  if (check_flag((void *)u, "N_VNew", 0, my_pe)) MPI_Abort(comm, 1);  SetInitialProfiles(u, data);  abstol = ATOL; reltol = RTOL;  /*      Call CVodeCreate to create the solver memory:          CV_BDF     specifies the Backward Differentiation Formula     CV_NEWTON  specifies a Newton iteration     A pointer to the integrator memory is returned and stored in cvode_mem.  */  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);  if (check_flag((void *)cvode_mem, "CVodeCreate", 0, my_pe)) MPI_Abort(comm, 1);  /* Set the pointer to user-defined data */  flag = CVodeSetFdata(cvode_mem, data);  if (check_flag(&flag, "CVodeSetFdata", 1, my_pe)) MPI_Abort(comm, 1);  /*      Call CVodeMalloc to initialize the integrator memory:      cvode_mem is the pointer to the integrator memory returned by CVodeCreate     f       is the user's right hand side function in y'=f(t,y)     T0      is the initial time     u       is the initial dependent variable vector     CV_SS   specifies scalar relative and absolute tolerances     reltol  is the relative tolerance     &abstol is a pointer to the scalar absolute tolerance  */  flag = CVodeMalloc(cvode_mem, f, T0, u, CV_SS, reltol, &abstol);  if (check_flag(&flag, "CVodeMalloc", 1, my_pe)) MPI_Abort(comm, 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, my_pe)) MPI_Abort(comm, 1);  /* Set preconditioner setup and solve routines Precond and PSolve,      and the pointer to the user-defined block data */  flag = CVSpilsSetPreconditioner(cvode_mem, Precond, PSolve, predata);  if (check_flag(&flag, "CVSpilsSetPreconditioner", 1, my_pe)) MPI_Abort(comm, 1);  if (my_pe == 0)    printf("/n2-species diurnal advection-diffusion problem/n/n");  /* 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);    if (check_flag(&flag, "CVode", 1, my_pe)) break;    PrintOutput(cvode_mem, my_pe, comm, u, t);  }  /* Print final statistics */    if (my_pe == 0) PrintFinalStats(cvode_mem);//.........这里部分代码省略.........

int main(){  realtype abstol, reltol, t, tout;  N_Vector u;  UserData data;  void *cvode_mem;  int iout, flag;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Allocate memory, and set problem data, initial values, tolerances */   u = N_VNew_Serial(NEQ);  if(check_flag((void *)u, "N_VNew_Serial", 0)) return(1);  data = AllocUserData();  if(check_flag((void *)data, "AllocUserData", 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);  /* set the JAcobian-times-vector function */  flag = CVSpilsSetJacTimesVecFn(cvode_mem, jtv);  if(check_flag(&flag, "CVSpilsSetJacTimesVecFn", 1)) return(1);  /* Set the preconditioner solve and setup functions */  flag = CVSpilsSetPreconditioner(cvode_mem, Precond, PSolve);  if(check_flag(&flag, "CVSpilsSetPreconditioner", 1)) return(1);  /* In loop over output points, call CVode, print results, test for error */  printf(" /n2-species diurnal advection-diffusion problem/n/n");  for (iout=1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {    flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);    PrintOutput(cvode_mem, u, t);    if(check_flag(&flag, "CVode", 1)) break;  }  PrintFinalStats(cvode_mem);  /* Free memory */  N_VDestroy_Serial(u);  FreeUserData(data);  CVodeFree(&cvode_mem);  return(0);}

//.........这里部分代码省略.........            if (absoluteTolerance < 0) {                emit error(QObject::tr("the 'absolute tolerance' property must have a value greater than or equal to 0"));                return;            }        } else {            emit error(QObject::tr("the 'absolute tolerance' property value could not be retrieved"));            return;        }        if (mProperties.contains(InterpolateSolutionId)) {            mInterpolateSolution = mProperties.value(InterpolateSolutionId).toBool();        } else {            emit error(QObject::tr("the 'interpolate solution' property value could not be retrieved"));            return;        }        // Initialise the ODE solver itself        OpenCOR::Solver::OdeSolver::initialize(pVoiStart, pRatesStatesCount,                                               pConstants, pRates, pStates,                                               pAlgebraic, pComputeRates);        // Create the states vector        mStatesVector = N_VMake_Serial(pRatesStatesCount, pStates);        // Create the CVODE solver        bool newtonIteration = !iterationType.compare(NewtonIteration);        mSolver = CVodeCreate(!integrationMethod.compare(BdfMethod)?CV_BDF:CV_ADAMS,                              newtonIteration?CV_NEWTON:CV_FUNCTIONAL);        // Use our own error handler        CVodeSetErrHandlerFn(mSolver, errorHandler, this);        // Initialise the CVODE solver        CVodeInit(mSolver, rhsFunction, pVoiStart, mStatesVector);        // Set some user data        mUserData = new CvodeSolverUserData(pConstants, pAlgebraic,                                            pComputeRates);        CVodeSetUserData(mSolver, mUserData);        // Set the maximum step        CVodeSetMaxStep(mSolver, maximumStep);        // Set the maximum number of steps        CVodeSetMaxNumSteps(mSolver, maximumNumberOfSteps);        // Set the linear solver, if needed        if (newtonIteration) {            if (!linearSolver.compare(DenseLinearSolver)) {                CVDense(mSolver, pRatesStatesCount);            } else if (!linearSolver.compare(BandedLinearSolver)) {                CVBand(mSolver, pRatesStatesCount, upperHalfBandwidth, lowerHalfBandwidth);

int main(void){  realtype abstol, reltol, t, tout;  N_Vector u;  UserData data;  void *cvode_mem;  int linsolver, iout, flag;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Allocate memory, and set problem data, initial values, tolerances */   u = N_VNew_Serial(NEQ);  if(check_flag((void *)u, "N_VNew_Serial", 0)) return(1);  data = AllocUserData();  if(check_flag((void *)data, "AllocUserData", 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);  /* START: Loop through SPGMR, SPBCG and SPTFQMR linear solver modules */  for (linsolver = 0; linsolver < 3; ++linsolver) {    if (linsolver != 0) {      /* Re-initialize user data */      InitUserData(data);      SetInitialProfiles(u, data->dx, data->dy);    /* Re-initialize CVode for the solution of the same problem, but       using a different linear solver module */      flag = CVodeReInit(cvode_mem, T0, u);      if (check_flag(&flag, "CVodeReInit", 1)) return(1);    }    /* Attach a linear solver module */    switch(linsolver) {    /* (a) SPGMR */    case(USE_SPGMR):      /* Print header */      printf(" -------");      printf(" /n| SPGMR |/n");      printf(" -------/n");      /* 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);      /* Set modified Gram-Schmidt orthogonalization, preconditioner 	 setup and solve routines Precond and PSolve, and the pointer 	 to the user-defined block data */      flag = CVSpilsSetGSType(cvode_mem, MODIFIED_GS);      if(check_flag(&flag, "CVSpilsSetGSType", 1)) return(1);      break;    /* (b) SPBCG */    case(USE_SPBCG):      /* Print header */      printf(" -------");      printf(" /n| SPBCG |/n");      printf(" -------/n");      /* Call CVSpbcg to specify the linear solver CVSPBCG 	 with left preconditioning and the maximum Krylov dimension maxl */      flag = CVSpbcg(cvode_mem, PREC_LEFT, 0);      if(check_flag(&flag, "CVSpbcg", 1)) return(1);      break;    /* (c) SPTFQMR */    case(USE_SPTFQMR)://.........这里部分代码省略.........

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);}

int main(){  realtype t, tout;  N_Vector y;  void *cvode_mem;  int flag, flagr, iout;  int rootsfound[2];  y = NULL;  cvode_mem = NULL;  /* Create serial vector of length NEQ for I.C. */  y = N_VNew_Serial(NEQ);  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);  /* Initialize y */  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 inital time T0, and   * the initial dependent variable vector y. */  flag = CVodeInit(cvode_mem, f, T0, y);  if (check_flag(&flag, "CVodeInit", 1)) return(1);  /* Use private function to compute error weights */  flag = CVodeWFtolerances(cvode_mem, ewt);  if (check_flag(&flag, "CVodeSetEwtFn", 1)) return(1);  /* Call CVodeRootInit to specify the root function g with 2 components */  flag = CVodeRootInit(cvode_mem, 2, g);  if (check_flag(&flag, "CVodeRootInit", 1)) return(1);  /* Call CVDense to specify the CVDENSE dense linear solver */  flag = CVDense(cvode_mem, NEQ);  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.  */  printf(" /n3-species kinetics problem/n/n");  iout = 0;  tout = T1;  while(1) {    flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);    PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));    if (flag == CV_ROOT_RETURN) {      flagr = CVodeGetRootInfo(cvode_mem, rootsfound);      check_flag(&flagr, "CVodeGetRootInfo", 1);      PrintRootInfo(rootsfound[0],rootsfound[1]);    }    if (check_flag(&flag, "CVode", 1)) break;    if (flag == CV_SUCCESS) {      iout++;      tout *= TMULT;    }    if (iout == NOUT) break;  }  /* Print some final statistics */  PrintFinalStats(cvode_mem);  /* Free y vector */  N_VDestroy_Serial(y);  /* Free integrator memory */  CVodeFree(&cvode_mem);  return(0);}

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 main(int argc, char *argv[]){  UserData data;  SUNMatrix A, AB;  SUNLinearSolver LS, LSB;  void *cvode_mem;  realtype reltolQ, abstolQ;  N_Vector y, q, constraints;  int steps;  int indexB;  realtype reltolB, abstolB, abstolQB;  N_Vector yB, qB, constraintsB;  realtype time;  int retval, ncheck;  long int nst, nstB;  CVadjCheckPointRec *ckpnt;  data = NULL;  A = AB = NULL;  LS = LSB = NULL;  cvode_mem = NULL;  ckpnt = NULL;  y = yB = qB = NULL;  constraints = NULL;  constraintsB = NULL;  /* Print problem description */  printf("/nAdjoint Sensitivity Example for Chemical Kinetics/n");  printf("-------------------------------------------------/n/n");  printf("ODE: dy1/dt = -p1*y1 + p2*y2*y3/n");  printf("     dy2/dt =  p1*y1 - p2*y2*y3 - p3*(y2)^2/n");  printf("     dy3/dt =  p3*(y2)^2/n/n");  printf("Find dG/dp for/n");  printf("     G = int_t0^tB0 g(t,p,y) dt/n");  printf("     g(t,p,y) = y3/n/n/n");  /* User data structure */  data = (UserData) malloc(sizeof *data);  if (check_retval((void *)data, "malloc", 2)) return(1);  data->p[0] = RCONST(0.04);  data->p[1] = RCONST(1.0e4);  data->p[2] = RCONST(3.0e7);  /* Initialize y */  y = N_VNew_Serial(NEQ);  if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);  Ith(y,1) = RCONST(1.0);  Ith(y,2) = ZERO;  Ith(y,3) = ZERO;  /* Set constraints to all 1's for nonnegative solution values. */  constraints = N_VNew_Serial(NEQ);  if(check_retval((void *)constraints, "N_VNew_Serial", 0)) return(1);  N_VConst(ONE, constraints);  /* Initialize q */  q = N_VNew_Serial(1);  if (check_retval((void *)q, "N_VNew_Serial", 0)) return(1);  Ith(q,1) = ZERO;  /* Set the scalar realtive and absolute tolerances reltolQ and abstolQ */  reltolQ = RTOL;  abstolQ = ATOLq;  /* Create and allocate CVODES memory for forward run */  printf("Create and allocate CVODES memory for forward runs/n");  /* Call CVodeCreate to create the solver memory and specify the      Backward Differentiation Formula */  cvode_mem = CVodeCreate(CV_BDF);  if (check_retval((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. */  retval = CVodeInit(cvode_mem, f, T0, y);  if (check_retval(&retval, "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 */  retval = CVodeWFtolerances(cvode_mem, ewt);  if (check_retval(&retval, "CVodeWFtolerances", 1)) return(1);  /* Attach user data */  retval = CVodeSetUserData(cvode_mem, data);  if (check_retval(&retval, "CVodeSetUserData", 1)) return(1);  /* Call CVodeSetConstraints to initialize constraints */  retval = CVodeSetConstraints(cvode_mem, constraints);  if (check_retval(&retval, "CVODESetConstraints", 1)) return(1);  N_VDestroy(constraints);//.........这里部分代码省略.........

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(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);//.........这里部分代码省略.........

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);//.........这里部分代码省略.........

PetscErrorCode TSSetUp_Sundials(TS ts){  TS_Sundials    *cvode = (TS_Sundials*)ts->data;  PetscErrorCode ierr;  PetscInt       glosize,locsize,i,flag;  PetscScalar    *y_data,*parray;  void           *mem;  PC             pc;  PCType         pctype;  PetscBool      pcnone;  PetscFunctionBegin;  /* get the vector size */  ierr = VecGetSize(ts->vec_sol,&glosize);CHKERRQ(ierr);  ierr = VecGetLocalSize(ts->vec_sol,&locsize);CHKERRQ(ierr);  /* allocate the memory for N_Vec y */  cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);  if (!cvode->y) SETERRQ(PETSC_COMM_SELF,1,"cvode->y is not allocated");  /* initialize N_Vec y: copy ts->vec_sol to cvode->y */  ierr   = VecGetArray(ts->vec_sol,&parray);CHKERRQ(ierr);  y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);  for (i = 0; i < locsize; i++) y_data[i] = parray[i];  ierr = VecRestoreArray(ts->vec_sol,NULL);CHKERRQ(ierr);  ierr = VecDuplicate(ts->vec_sol,&cvode->update);CHKERRQ(ierr);  ierr = VecDuplicate(ts->vec_sol,&cvode->ydot);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);CHKERRQ(ierr);  /*    Create work vectors for the TSPSolve_Sundials() routine. Note these are    allocated with zero space arrays because the actual array space is provided    by Sundials and set using VecPlaceArray().  */  ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w1);CHKERRQ(ierr);  ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w2);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);CHKERRQ(ierr);  /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */  mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);  if (!mem) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"CVodeCreate() fails");  cvode->mem = mem;  /* Set the pointer to user-defined data */  flag = CVodeSetUserData(mem, ts);  if (flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetUserData() fails");  /* Sundials may choose to use a smaller initial step, but will never use a larger step. */  flag = CVodeSetInitStep(mem,(realtype)ts->time_step);  if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetInitStep() failed");  if (cvode->mindt > 0) {    flag = CVodeSetMinStep(mem,(realtype)cvode->mindt);    if (flag) {      if (flag == CV_MEM_NULL) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL");      else if (flag == CV_ILL_INPUT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, hmin is nonpositive or it exceeds the maximum allowable step size");      else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed");    }  }  if (cvode->maxdt > 0) {    flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);    if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMaxStep() failed");  }  /* 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 cvode->y */  flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y);  if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag);  /* specifies scalar relative and absolute tolerances */  flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol);  if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag);  /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */  flag = CVodeSetMaxNumSteps(mem,ts->max_steps);  /* call CVSpgmr to use GMRES as the linear solver.        */  /* setup the ode integrator with the given preconditioner */  ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr);  ierr = PCGetType(pc,&pctype);CHKERRQ(ierr);  ierr = PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);CHKERRQ(ierr);  if (pcnone) {    flag = CVSpgmr(mem,PREC_NONE,0);    if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);  } else {    flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl);    if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);    /* Set preconditioner and solve routines Precond and PSolve,     and the pointer to the user-defined block data */    flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);    if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag);  }  flag = CVSpilsSetGSType(mem, MODIFIED_GS);  if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);  PetscFunctionReturn(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;//.........这里部分代码省略.........

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);//.........这里部分代码省略.........

  void CVodesIntegrator::initialize(double t0, FuncEval& func)   {    m_neq = func.neq();    m_t0  = t0;    if (m_y) {      N_VDestroy_Serial(nv(m_y));    // free solution vector if already allocated    }    m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq));   // allocate solution vector    for (int i=0; i<m_neq; i++) {      NV_Ith_S(nv(m_y), i) = 0.0;    }    // check abs tolerance array size    if (m_itol == CV_SV && m_nabs < m_neq)       throw CVodesErr("not enough absolute tolerance values specified.");    func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));    if (m_cvode_mem) CVodeFree(&m_cvode_mem);    /*     *  Specify the method and the iteration type:     *      Cantera Defaults:     *         CV_BDF  - Use BDF methods      *         CV_NEWTON - use newton's method      */    m_cvode_mem = CVodeCreate(m_method, m_iter);    if (!m_cvode_mem) throw CVodesErr("CVodeCreate failed.");    int flag = 0;#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)    if (m_itol == CV_SV) {      // vector atol      flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,			 m_reltol, nv(m_abstol));    }    else {      // scalar atol      flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,			 m_reltol, &m_abstols);    }    if (flag != CV_SUCCESS) {      if (flag == CV_MEM_FAIL) {	throw CVodesErr("Memory allocation failed.");      }      else if (flag == CV_ILL_INPUT) {	throw CVodesErr("Illegal value for CVodeMalloc input argument.");      }      else 	throw CVodesErr("CVodeMalloc failed.");    }#elif defined(SUNDIALS_VERSION_24)    flag = CVodeInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y));    if (flag != CV_SUCCESS) {      if (flag == CV_MEM_FAIL) {	throw CVodesErr("Memory allocation failed.");      } else if (flag == CV_ILL_INPUT) {	throw CVodesErr("Illegal value for CVodeInit input argument.");      } else {	throw CVodesErr("CVodeInit failed.");      }    }    if (m_itol == CV_SV) {      flag = CVodeSVtolerances(m_cvode_mem, m_reltol, nv(m_abstol));    } else {      flag = CVodeSStolerances(m_cvode_mem, m_reltol, m_abstols);    }    if (flag != CV_SUCCESS) {      if (flag == CV_MEM_FAIL) {	throw CVodesErr("Memory allocation failed.");      } else if (flag == CV_ILL_INPUT) {	throw CVodesErr("Illegal value for CVodeInit input argument.");      } else {	throw CVodesErr("CVodeInit failed.");      }    }#else    printf("unknown sundials verson/n");    exit(-1);#endif    if (m_type == DENSE + NOJAC) {      long int N = m_neq;      CVDense(m_cvode_mem, N);    }    else if (m_type == DIAG) {      CVDiag(m_cvode_mem);    }    else if (m_type == GMRES) {      CVSpgmr(m_cvode_mem, PREC_NONE, 0);    }    else if (m_type == BAND + NOJAC) {      long int N = m_neq;      long int nu = m_mupper;      long int nl = m_mlower;      CVBand(m_cvode_mem, N, nu, nl);//.........这里部分代码省略.........