// Tries detecting a memory leak on the particular input that we have just
// executed before calling this function.
void Fuzzer::TryDetectingAMemoryLeak(uint8_t *Data, size_t Size) {
  if (!HasMoreMallocsThanFrees) return;  // mallocs==frees, a leak is unlikely.
  if (!Options.DetectLeaks) return;
  if (!&__lsan_enable || !&__lsan_disable || !__lsan_do_recoverable_leak_check)
    return;  // No lsan.
  // Run the target once again, but with lsan disabled so that if there is
  // a real leak we do not report it twice.
  __lsan_disable();
  RunOneAndUpdateCorpus(Data, Size);
  __lsan_enable();
  if (!HasMoreMallocsThanFrees) return;  // a leak is unlikely.
  if (NumberOfLeakDetectionAttempts++ > 1000) {
    Options.DetectLeaks = false;
    Printf("INFO: libFuzzer disabled leak detection after every mutation.\n"
           "      Most likely the target function accumulates allocated\n"
           "      memory in a global state w/o actually leaking it.\n"
           "      If LeakSanitizer is enabled in this process it will still\n"
           "      run on the process shutdown.\n");
    return;
  }
  // Now perform the actual lsan pass. This is expensive and we must ensure
  // we don't call it too often.
  if (__lsan_do_recoverable_leak_check()) {  // Leak is found, report it.
    CurrentUnitData = Data;
    CurrentUnitSize = Size;
    DumpCurrentUnit("leak-");
    PrintFinalStats();
    _Exit(Options.ErrorExitCode);  // not exit() to disable lsan further on.
  }
}

void Fuzzer::AlarmCallback() {
  assert(Options.UnitTimeoutSec > 0);
  if (!CurrentUnitSize)
    return; // We have not started running units yet.
  size_t Seconds =
      duration_cast<seconds>(system_clock::now() - UnitStartTime).count();
  if (Seconds == 0)
    return;
  if (Options.Verbosity >= 2)
    Printf("AlarmCallback %zd\n", Seconds);
  if (Seconds >= (size_t)Options.UnitTimeoutSec) {
    Printf("ALARM: working on the last Unit for %zd seconds\n", Seconds);
    Printf("       and the timeout value is %d (use -timeout=N to change)\n",
           Options.UnitTimeoutSec);
    if (CurrentUnitSize <= kMaxUnitSizeToPrint) {
      PrintHexArray(CurrentUnitData, CurrentUnitSize, "\n");
      PrintASCII(CurrentUnitData, CurrentUnitSize, "\n");
    }
    WriteUnitToFileWithPrefix(
        {CurrentUnitData, CurrentUnitData + CurrentUnitSize}, "timeout-");
    Printf("==%d== ERROR: libFuzzer: timeout after %d seconds\n", GetPid(),
           Seconds);
    if (__sanitizer_print_stack_trace)
      __sanitizer_print_stack_trace();
    Printf("SUMMARY: libFuzzer: timeout\n");
    PrintFinalStats();
    if (Options.AbortOnTimeout)
      abort();
    exit(Options.TimeoutExitCode);
  }
}

NO_SANITIZE_MEMORY
void Fuzzer::DeathCallback() {
  if (!CurrentUnitSize) return;
  Printf("DEATH:\n");
  DumpCurrentUnit("crash-");
  PrintFinalStats();
}

static int SolveIt(void *kmem, N_Vector u, N_Vector s, int glstr, int mset)
{
  int flag;

  printf("\n");

  if (mset==1)
    printf("Exact Newton");
  else
    printf("Modified Newton");

  if (glstr == KIN_NONE)
    printf("\n");
  else
    printf(" with line search\n");

  flag = KINSetMaxSetupCalls(kmem, mset);
  if (check_flag(&flag, "KINSetMaxSetupCalls", 1)) return(1);

  flag = KINSol(kmem, u, glstr, s, s);
  if (check_flag(&flag, "KINSol", 1)) return(1);

  printf("Solution:\n  [x1,x2] = ");
  PrintOutput(u);

  PrintFinalStats(kmem);

  return(0);

}

void Fuzzer::DeathCallback() {
  if (!CurrentUnitSize) return;
  Printf("DEATH:\n");
  if (CurrentUnitSize <= kMaxUnitSizeToPrint) {
    PrintHexArray(CurrentUnitData, CurrentUnitSize, "\n");
    PrintASCII(CurrentUnitData, CurrentUnitSize, "\n");
  }
  WriteUnitToFileWithPrefix(
      {CurrentUnitData, CurrentUnitData + CurrentUnitSize}, "crash-");
  PrintFinalStats();
}

int main(void)
{
  realtype dx, dy, reltol, abstol, t, tout, umax;
  N_Vector u, up;
  UserData data;
  void *cvode_mem;
  int iout, flag;
  long int nst;

  u = NULL;
  data = NULL;
  cvode_mem = NULL;

  u  = N_VNew_Serial(NEQ);
  up = N_VNew_Serial(NEQ);

  reltol = ZERO;
  abstol = ATOL;

  data = (UserData) malloc(sizeof *data);
  dx = data->dx = XMAX/(MX+1);
  dy = data->dy = YMAX/(MY+1);
  data->hdcoef = ONE/(dx*dx);
  data->hacoef = HALF/(TWO*dx);
  data->vdcoef = ONE/(dy*dy);

  SetIC(u, data);

  cvode_mem = CPodeCreate(CP_EXPL, CP_BDF, CP_NEWTON);
  flag = CPodeInit(cvode_mem, (void *)f, data, T0, u, NULL, CP_SS, reltol, &abstol);

  flag = CPLapackBand(cvode_mem, NEQ, MY, MY);
  flag = CPDlsSetJacFn(cvode_mem, (void *)Jac, data);

  /* In loop over output points: call CPode, print results, test for errors */
  umax = N_VMaxNorm(u);
  PrintHeader(reltol, abstol, umax);
  for(iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {
    flag = CPode(cvode_mem, tout, &t, u, up, CP_NORMAL);
    umax = N_VMaxNorm(u);
    flag = CPodeGetNumSteps(cvode_mem, &nst);
    PrintOutput(t, umax, nst);
  }

  PrintFinalStats(cvode_mem);

  N_VDestroy_Serial(u);
  CPodeFree(&cvode_mem);
  free(data);

  return(0);
}

NO_SANITIZE_MEMORY
void Fuzzer::AlarmCallback() {
  assert(Options.UnitTimeoutSec > 0);
  if (InOOMState) {
    Printf("==%d== ERROR: libFuzzer: out-of-memory (used: %zdMb; limit: %zdMb)\n",
           GetPid(), GetPeakRSSMb(), Options.RssLimitMb);
    Printf("   To change the out-of-memory limit use -rss_limit_mb=<N>\n");
    if (CurrentUnitSize && CurrentUnitData) {
      DumpCurrentUnit("oom-");
      if (__sanitizer_print_stack_trace)
        __sanitizer_print_stack_trace();
    }
    Printf("SUMMARY: libFuzzer: out-of-memory\n");
    PrintFinalStats();
    _Exit(Options.ErrorExitCode); // Stop right now.
  }

  if (!CurrentUnitSize)
    return; // We have not started running units yet.
  size_t Seconds =
      duration_cast<seconds>(system_clock::now() - UnitStartTime).count();
  if (Seconds == 0)
    return;
  if (Options.Verbosity >= 2)
    Printf("AlarmCallback %zd\n", Seconds);
  if (Seconds >= (size_t)Options.UnitTimeoutSec) {
    Printf("ALARM: working on the last Unit for %zd seconds\n", Seconds);
    Printf("       and the timeout value is %d (use -timeout=N to change)\n",
           Options.UnitTimeoutSec);
    DumpCurrentUnit("timeout-");
    Printf("==%d== ERROR: libFuzzer: timeout after %d seconds\n", GetPid(),
           Seconds);
    if (__sanitizer_print_stack_trace)
      __sanitizer_print_stack_trace();
    Printf("SUMMARY: libFuzzer: timeout\n");
    PrintFinalStats();
    _Exit(Options.TimeoutExitCode); // Stop right now.
  }
}

static void Problem2(void)
{
  void *fct;
  void *cpode_mem;
  N_Vector y, yp;
  realtype reltol=RTOL, abstol=ATOL, t, tout, erm, hu;
  int flag, iout, qu;

  y = NULL;
  yp = NULL;
  cpode_mem = NULL;

  y = N_VNew_Serial(P2_NEQ);
  N_VConst(ZERO, y);
  NV_Ith_S(y,0) = ONE;

  yp = N_VNew_Serial(P2_NEQ);

  if (ODE == CP_EXPL) {
    fct = (void *)f2;
  } else {
    fct = (void *)res2;
    f2(P2_T0, y, yp, NULL);
  }

  cpode_mem = CPodeCreate(ODE, CP_ADAMS, CP_FUNCTIONAL);
  /*  flag = CPodeSetInitStep(cpode_mem, 2.0e-9);*/
  flag = CPodeInit(cpode_mem, fct, NULL, P2_T0, y, yp, CP_SS, reltol, &abstol);

  printf("\n      t        max.err      qu     hu \n");
  for(iout=1, tout=P2_T1; iout <= P2_NOUT; iout++, tout*=P2_TOUT_MULT) {
    flag = CPode(cpode_mem, tout, &t, y, yp, CP_NORMAL);
    if (flag != CP_SUCCESS) break;
    erm = MaxError(y, t);
    flag = CPodeGetLastOrder(cpode_mem, &qu);
    flag = CPodeGetLastStep(cpode_mem, &hu);
    printf("%10.3f  %12.4le   %2d   %12.4le\n", t, erm, qu, hu);
  }

  PrintFinalStats(cpode_mem);

  CPodeFree(&cpode_mem);
  N_VDestroy_Serial(y);
  N_VDestroy_Serial(yp);

  return;
}

int main()
{
  void *fct;
  void *cpode_mem;
  N_Vector y, yp;
  realtype reltol=RTOL, abstol=ATOL, t, tout, hu;
  int flag, iout, qu;

  y = NULL;
  yp = NULL;
  cpode_mem = NULL;

  y = N_VNew_Serial(P1_NEQ);
  NV_Ith_S(y,0) = TWO;
  NV_Ith_S(y,1) = ZERO;

  yp = N_VNew_Serial(P1_NEQ);
  
  if (ODE == CP_EXPL) {
    fct = (void *)f;
  } else {
    fct = (void *)res;
    f(P1_T0, y, yp, NULL);
  }

  cpode_mem = CPodeCreate(ODE, CP_ADAMS, CP_FUNCTIONAL);
  /*  flag = CPodeSetInitStep(cpode_mem, 4.0e-9);*/
  flag = CPodeInit(cpode_mem, fct, NULL, P1_T0, y, yp, CP_SS, reltol, &abstol);

  printf("\n     t           x              xdot         qu     hu \n");
  for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {
    flag = CPode(cpode_mem, tout, &t, y, yp, CP_NORMAL);
    if (flag != CP_SUCCESS)  break;
    flag = CPodeGetLastOrder(cpode_mem, &qu);
    flag = CPodeGetLastStep(cpode_mem, &hu);    
    printf("%10.5f    %12.5le   %12.5le   %2d    %6.4le\n", t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);
  }
  
  PrintFinalStats(cpode_mem);

  CPodeFree(&cpode_mem);
  N_VDestroy_Serial(y);
  N_VDestroy_Serial(yp);

  return 0;
}

int main()
{
  realtype fnormtol, scsteptol;
  N_Vector y, scale, constraints;
  int mset, flag, i;
  void *kmem;
  SUNMatrix J;
  SUNLinearSolver LS;

  y = scale = constraints = NULL;
  kmem = NULL;
  J = NULL;
  LS = NULL;

  printf("\nRobot Kinematics Example\n");
  printf("8 variables; -1 <= x_i <= 1\n");
  printf("KINSOL problem size: 8 + 2*8 = 24 \n\n");

  /* Create vectors for solution, scales, and constraints */

  y = N_VNew_Serial(NEQ);
  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);

  scale = N_VNew_Serial(NEQ);
  if (check_flag((void *)scale, "N_VNew_Serial", 0)) return(1);

  constraints = N_VNew_Serial(NEQ);
  if (check_flag((void *)constraints, "N_VNew_Serial", 0)) return(1);

  /* Initialize and allocate memory for KINSOL */

  kmem = KINCreate();
  if (check_flag((void *)kmem, "KINCreate", 0)) return(1);

  flag = KINInit(kmem, func, y); /* y passed as a template */
  if (check_flag(&flag, "KINInit", 1)) return(1);

  /* Set optional inputs */

  N_VConst_Serial(ZERO,constraints);
  for (i = NVAR+1; i <= NEQ; i++) Ith(constraints, i) = ONE;
  
  flag = KINSetConstraints(kmem, constraints);
  if (check_flag(&flag, "KINSetConstraints", 1)) return(1);

  fnormtol  = FTOL; 
  flag = KINSetFuncNormTol(kmem, fnormtol);
  if (check_flag(&flag, "KINSetFuncNormTol", 1)) return(1);

  scsteptol = STOL;
  flag = KINSetScaledStepTol(kmem, scsteptol);
  if (check_flag(&flag, "KINSetScaledStepTol", 1)) return(1);

  /* Create dense SUNMatrix */
  J = SUNDenseMatrix(NEQ, NEQ);
  if(check_flag((void *)J, "SUNDenseMatrix", 0)) return(1);

  /* Create dense SUNLinearSolver object */
  LS = SUNLinSol_Dense(y, J);
  if(check_flag((void *)LS, "SUNLinSol_Dense", 0)) return(1);

  /* Attach the matrix and linear solver to KINSOL */
  flag = KINSetLinearSolver(kmem, LS, J);
  if(check_flag(&flag, "KINSetLinearSolver", 1)) return(1);

  /* Set the Jacobian function */
  flag = KINSetJacFn(kmem, jac);
  if (check_flag(&flag, "KINSetJacFn", 1)) return(1);

  /* Indicate exact Newton */

  mset = 1;
  flag = KINSetMaxSetupCalls(kmem, mset);
  if (check_flag(&flag, "KINSetMaxSetupCalls", 1)) return(1);

  /* Initial guess */

  N_VConst_Serial(ONE, y);
  for(i = 1; i <= NVAR; i++) Ith(y,i) = SUNRsqrt(TWO)/TWO;

  printf("Initial guess:\n");
  PrintOutput(y);

  /* Call KINSol to solve problem */

  N_VConst_Serial(ONE,scale);
  flag = KINSol(kmem,           /* KINSol memory block */
                y,              /* initial guess on input; solution vector */
                KIN_LINESEARCH, /* global strategy choice */
                scale,          /* scaling vector, for the variable cc */
                scale);         /* scaling vector for function values fval */
  if (check_flag(&flag, "KINSol", 1)) return(1);

  printf("\nComputed solution:\n");
  PrintOutput(y);

  /* Print final statistics and free memory */  

  PrintFinalStats(kmem);

//.........这里部分代码省略.........

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

static int Problem1(void)
{
  realtype reltol=RTOL, abstol=ATOL, t, tout, ero, er;
  int miter, flag, temp_flag, iout, nerr=0;
  N_Vector y;
  void *cvode_mem;
  booleantype firstrun;
  int qu;
  realtype hu;

  y = NULL;
  cvode_mem = NULL;

  y = N_VNew_Serial(P1_NEQ);
  if(check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
  PrintIntro1();

  cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

  for (miter=FUNC; miter <= DIAG; miter++) {
    ero = ZERO;
    NV_Ith_S(y,0) = TWO;
    NV_Ith_S(y,1) = ZERO;

    firstrun = (miter==FUNC);
    if (firstrun) {
      flag = CVodeMalloc(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);
      if(check_flag(&flag, "CVodeMalloc", 1)) return(1);
    } else {
      flag = CVodeSetIterType(cvode_mem, CV_NEWTON);
      if(check_flag(&flag, "CVodeSetIterType", 1)) ++nerr;
      flag = CVodeReInit(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);
      if(check_flag(&flag, "CVodeReInit", 1)) return(1);
    }
      
    flag = PrepareNextRun(cvode_mem, CV_ADAMS, miter, 0, 0);
    if(check_flag(&flag, "PrepareNextRun", 1)) return(1);

    PrintHeader1();

    for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {
      flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
      check_flag(&flag, "CVode", 1);
      temp_flag = CVodeGetLastOrder(cvode_mem, &qu);
      if(check_flag(&temp_flag, "CVodeGetLastOrder", 1)) ++nerr;
      temp_flag = CVodeGetLastStep(cvode_mem, &hu);
      if(check_flag(&temp_flag, "CVodeGetLastStep", 1)) ++nerr;
      PrintOutput1(t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);
      if (flag != CV_SUCCESS) {
        nerr++;
        break;
      }
      if (iout%2 == 0) {
        er = ABS(NV_Ith_S(y,0)) / abstol;
        if (er > ero) ero = er;
        if (er > P1_TOL_FACTOR) {
          nerr++;
	  PrintErrOutput(P1_TOL_FACTOR);
        }
      }
    }
    
    PrintFinalStats(cvode_mem, miter, ero);
  }

  CVodeFree(cvode_mem);

  cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);
  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

  for (miter=FUNC; miter <= DIAG; miter++) {
    ero = ZERO;
    NV_Ith_S(y,0) = TWO;
    NV_Ith_S(y,1) = ZERO;
      
    firstrun = (miter==FUNC);
    if (firstrun) {
      flag = CVodeMalloc(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);
      if(check_flag(&flag, "CVodeMalloc", 1)) return(1);
    } else {
      flag = CVodeSetIterType(cvode_mem, CV_NEWTON);
      if(check_flag(&flag, "CVodeSetIterType", 1)) ++nerr;
      flag = CVodeReInit(cvode_mem, f1, P1_T0, y, CV_SS, reltol, &abstol);
      if(check_flag(&flag, "CVodeReInit", 1)) return(1);
    }
      
    flag = PrepareNextRun(cvode_mem, CV_BDF, miter, 0, 0);     
    if(check_flag(&flag, "PrepareNextRun", 1)) return(1);

    PrintHeader1();
      
    for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {
      flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
      check_flag(&flag, "CVode", 1);
      temp_flag = CVodeGetLastOrder(cvode_mem, &qu);
      if(check_flag(&temp_flag, "CVodeGetLastOrder", 1)) ++nerr;
      temp_flag = CVodeGetLastStep(cvode_mem, &hu);
      if(check_flag(&temp_flag, "CVodeGetLastStep", 1)) ++nerr;
      PrintOutput1(t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);
//.........这里部分代码省略.........

int main(void)
{
  int globalstrategy;
  realtype fnormtol, scsteptol;
  N_Vector cc, sc, constraints;
  UserData data;
  int flag, maxl, maxlrst;
  void *kmem;

  cc = sc = constraints = NULL;
  kmem = NULL;
  data = NULL;

  /* Allocate memory, and set problem data, initial values, tolerances */ 
  globalstrategy = KIN_NONE;

  data = AllocUserData();
  if (check_flag((void *)data, "AllocUserData", 2)) return(1);
  InitUserData(data);

  /* Create serial vectors of length NEQ */
  cc = N_VNew_Serial(NEQ);
  if (check_flag((void *)cc, "N_VNew_Serial", 0)) return(1);
  sc = N_VNew_Serial(NEQ);
  if (check_flag((void *)sc, "N_VNew_Serial", 0)) return(1);
  data->rates = N_VNew_Serial(NEQ);
  if (check_flag((void *)data->rates, "N_VNew_Serial", 0)) return(1);

  constraints = N_VNew_Serial(NEQ);
  if (check_flag((void *)constraints, "N_VNew_Serial", 0)) return(1);
  N_VConst(TWO, constraints);

  SetInitialProfiles(cc, sc);

  fnormtol=FTOL; scsteptol=STOL;

  /* Call KINCreate/KINInit to initialize KINSOL.
     A pointer to KINSOL problem memory is returned and stored in kmem. */
  kmem = KINCreate();
  if (check_flag((void *)kmem, "KINCreate", 0)) return(1);

  /* Vector cc passed as template vector. */
  flag = KINInit(kmem, func, cc);
  if (check_flag(&flag, "KINInit", 1)) return(1);

  flag = KINSetUserData(kmem, data);
  if (check_flag(&flag, "KINSetUserData", 1)) return(1);
  flag = KINSetConstraints(kmem, constraints);
  if (check_flag(&flag, "KINSetConstraints", 1)) return(1);
  flag = KINSetFuncNormTol(kmem, fnormtol);
  if (check_flag(&flag, "KINSetFuncNormTol", 1)) return(1);
  flag = KINSetScaledStepTol(kmem, scsteptol);
  if (check_flag(&flag, "KINSetScaledStepTol", 1)) return(1);

  /* We no longer need the constraints vector since KINSetConstraints
     creates a private copy for KINSOL to use. */
  N_VDestroy_Serial(constraints);

  /* Call KINSpgmr to specify the linear solver KINSPGMR with preconditioner
     routines PrecSetupBD and PrecSolveBD. */
  maxl = 15; 
  maxlrst = 2;
  flag = KINSpgmr(kmem, maxl);
  if (check_flag(&flag, "KINSpgmr", 1)) return(1);

  flag = KINSpilsSetMaxRestarts(kmem, maxlrst);
  if (check_flag(&flag, "KINSpilsSetMaxRestarts", 1)) return(1);
  flag = KINSpilsSetPreconditioner(kmem,
				   PrecSetupBD,
				   PrecSolveBD);
  if (check_flag(&flag, "KINSpilsSetPreconditioner", 1)) return(1);

  /* Print out the problem size, solution parameters, initial guess. */
  PrintHeader(globalstrategy, maxl, maxlrst, fnormtol, scsteptol);

  /* Call KINSol and print output concentration profile */
  flag = KINSol(kmem,           /* KINSol memory block */
                cc,             /* initial guess on input; solution vector */
                globalstrategy, /* global stragegy choice */
                sc,             /* scaling vector, for the variable cc */
                sc);            /* scaling vector for function values fval */
  if (check_flag(&flag, "KINSol", 1)) return(1);

  printf("\n\nComputed equilibrium species concentrations:\n");
  PrintOutput(cc);

  /* Print final statistics and free memory */  
  PrintFinalStats(kmem);

  N_VDestroy_Serial(cc);
  N_VDestroy_Serial(sc);
  KINFree(&kmem);
  FreeUserData(data);

  return(0);
}

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

//.........这里部分代码省略.........
  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 */
  flag = CVodeGetQuad(cvode_mem, &tret, q);
  qdata = NV_DATA_P(q);
  MPI_Allreduce(&qdata[0], &G, 1, PVEC_REAL_MPI_TYPE, MPI_SUM, comm);
#if defined(SUNDIALS_EXTENDED_PRECISION)
  if (myId == 0) printf("  G = %Le\n",G);
#elif defined(SUNDIALS_DOUBLE_PRECISION)
  if (myId == 0) printf("  G = %e\n",G);
#else
  if (myId == 0) printf("  G = %e\n",G);
#endif

  /* Print statistics for forward run */
  if (myId == 0) PrintFinalStats(cvode_mem);

  /*-------------------------- 
    Backward integration phase
    --------------------------*/
 
  /* Allocate and initialize yB */
  yB = N_VNew_Parallel(comm, l_neq, neq); 
  N_VConst(ZERO, yB);

  /* Allocate and initialize qB (gradient) */
  qB = N_VNew_Parallel(comm, l_neq, neq); 
  N_VConst(ZERO, qB);

  /* Create and allocate backward CVODE memory */
  flag = CVodeCreateB(cvode_mem, CV_BDF, CV_NEWTON, &indexB);
  flag = CVodeSetUserDataB(cvode_mem, indexB, d);
  flag = CVodeInitB(cvode_mem, indexB, fB, tf, yB);
  abstolB = ATOL_B;  
  reltolB = RTOL_B; 
  flag = CVodeSStolerancesB(cvode_mem, indexB, reltolB, abstolB);

  /* Attach preconditioner and linear solver modules */
  flag = CVSpgmrB(cvode_mem, indexB, PREC_LEFT, 0); 
  mudqB = mldqB = d->l_m[0]+1;
  mukeepB = mlkeepB = 2;  
  flag = CVBBDPrecInitB(cvode_mem, indexB, l_neq, mudqB, mldqB, 
                        mukeepB, mlkeepB, ZERO, fB_local, NULL);

  /* Initialize quadrature calculations */
  abstolQB = ATOL_QB;
  reltolQB = RTOL_QB;
  flag = CVodeQuadInitB(cvode_mem, indexB, fQB, qB);

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

    }

    PrintFinalStats(cvode_mem);

    N_VDestroy_Serial(state);
    N_VDestroy_Serial(abstol);

    CVodeFree(&cvode_mem);

    fclose(pout);
    return(0);
}

//.........这里部分代码省略.........

  cp  = N_VNew_Parallel(comm, local_N, SystemSize);
  if(check_flag((void *)cp, "N_VNew_Parallel", 0, thispe)) MPI_Abort(comm, 1);

  res = N_VNew_Parallel(comm, local_N, SystemSize);
  if(check_flag((void *)res, "N_VNew_Parallel", 0, thispe)) MPI_Abort(comm, 1);

  id  = N_VNew_Parallel(comm, local_N, SystemSize);
  if(check_flag((void *)id, "N_VNew_Parallel", 0, thispe)) MPI_Abort(comm, 1);
  
  SetInitialProfiles(cc, cp, id, res, webdata);
  
  N_VDestroy_Parallel(res);
  
  /* Set remaining inputs to IDAMalloc. */
  
  t0 = ZERO;
  rtol = RTOL; 
  atol = ATOL;
  
  /* Call IDACreate and IDAMalloc to initialize solution */

  mem = IDACreate();
  if(check_flag((void *)mem, "IDACreate", 0, thispe)) MPI_Abort(comm, 1);

  retval = IDASetUserData(mem, webdata);
  if(check_flag(&retval, "IDASetUserData", 1, thispe)) MPI_Abort(comm, 1);

  retval = IDASetId(mem, id);
  if(check_flag(&retval, "IDASetId", 1, thispe)) MPI_Abort(comm, 1);

  retval = IDAInit(mem, resweb, t0, cc, cp);
  if(check_flag(&retval, "IDAInit", 1, thispe)) MPI_Abort(comm, 1);
  
  retval = IDASStolerances(mem, rtol, atol);
  if(check_flag(&retval, "IDASStolerances", 1, thispe)) MPI_Abort(comm, 1);

  /* Call IDASpgmr to specify the IDA linear solver IDASPGMR */

  maxl = 16;
  retval = IDASpgmr(mem, maxl);
  if(check_flag(&retval, "IDASpgmr", 1, thispe)) MPI_Abort(comm, 1);

  /* Call IDABBDPrecInit to initialize the band-block-diagonal preconditioner.
     The half-bandwidths for the difference quotient evaluation are exact
     for the system Jacobian, but only a 5-diagonal band matrix is retained. */
  
  mudq = mldq = NSMXSUB;
  mukeep = mlkeep = 2;
  retval = IDABBDPrecInit(mem, local_N, mudq, mldq, mukeep, mlkeep, 
                          ZERO, reslocal, NULL);
  if(check_flag(&retval, "IDABBDPrecInit", 1, thispe)) MPI_Abort(comm, 1);
  
  /* Call IDACalcIC (with default options) to correct the initial values. */
  
  tout = RCONST(0.001);
  retval = IDACalcIC(mem, IDA_YA_YDP_INIT, tout);
  if(check_flag(&retval, "IDACalcIC", 1, thispe)) MPI_Abort(comm, 1);
  
  /* On PE 0, print heading, basic parameters, initial values. */
 
  if (thispe == 0) PrintHeader(SystemSize, maxl, 
                               mudq, mldq, mukeep, mlkeep,
                               rtol, atol);
  PrintOutput(mem, cc, t0, webdata, comm);

  /* Call IDA in tout loop, normal mode, and print selected output. */
  
  for (iout = 1; iout <= NOUT; iout++) {
    
    retval = IDASolve(mem, tout, &tret, cc, cp, IDA_NORMAL);
    if(check_flag(&retval, "IDASolve", 1, thispe)) MPI_Abort(comm, 1);
    
    PrintOutput(mem, cc, tret, webdata, comm);
    
    if (iout < 3) tout *= TMULT; 
    else          tout += TADD;

  }
  
  /* On PE 0, print final set of statistics. */
  
  if (thispe == 0)  PrintFinalStats(mem);

  /* Free memory. */

  N_VDestroy_Parallel(cc);
  N_VDestroy_Parallel(cp);
  N_VDestroy_Parallel(id);

  IDAFree(&mem);

  destroyMat(webdata->acoef);
  N_VDestroy_Parallel(webdata->rates);
  free(webdata);

  MPI_Finalize();

  return(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_