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