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 AsyncExeCmd::OnTimer(wxTimerEvent &event)
{
	wxUnusedVar(event);
	if( m_stop ){
		m_proc->Terminate();
		return;
	}
	PrintOutput();
}

void Display(void)  {
	glClear(GL_COLOR_BUFFER_BIT);
	DrawAllButtons();
	if (stat.turnedON)	{
		PrintOutput();
	}
	glFlush();
	glutSwapBuffers();
}

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

void AsyncExeCmd::ProcessEnd(wxProcessEvent& event)
{
	//read all input before stopping the timer
	if( !m_stop ){
		PrintOutput();
	}
	
	m_timer->Stop();
	m_busy = false;
	m_stop = false;

	SendEndMsg(event.GetExitCode());
}

void rbBytecode::Parse()
{
	wxLogMessage(wxT(""));

	// Initialize stack
	// parse until end-script marker
	Root = new rbToken(NULL,-1,wxT("ROOT"),wxT("ROOT"));
	PushIndent();
	PrintOutput(wxT("{\n"));
	PushStack(Root);	
	while( true )
	{
		dword token = LoadToken();
		if( token == 0x53 )
			break;
	}
	PopStack();
	PopIndent();
	//ScriptText.RemoveLast(); // remove tab
	PrintOutput( wxT("}\n") );


}

INT wmain(INT argc, WCHAR **argv)
{
	po::variables_map variables_map;
	printInfoStruct stPrintInfo;
	if (!ParseArguments(argc, argv, variables_map, stPrintInfo))
		return 3;

	if (variables_map.count("print-objects")) {
		PrintObjects();
		return 0;
	}

	if (variables_map.count("print-object-info")) {
		PrintObjectInfo(stPrintInfo);
		return 0;
	}

	if (QueryPerfData(stPrintInfo))
		return PrintOutput(variables_map, stPrintInfo);
	else
		return 3;
}

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

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

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

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

  /* Choose global strategy */
  globalstrategy = KIN_NONE;

  /* Allocate and initialize vectors */
  cc = N_VNew_Parallel(comm, Nlocal, NEQ);
  if (check_flag((void *)cc, "N_VNew_Parallel", 0, my_pe)) MPI_Abort(comm, 1);
  sc = N_VNew_Parallel(comm, Nlocal, NEQ);
  if (check_flag((void *)sc, "N_VNew_Parallel", 0, my_pe)) MPI_Abort(comm, 1);
  data->rates = N_VNew_Parallel(comm, Nlocal, NEQ);
  if (check_flag((void *)data->rates, "N_VNew_Parallel", 0, my_pe))
      MPI_Abort(comm, 1);
  constraints = N_VNew_Parallel(comm, Nlocal, NEQ);
  if (check_flag((void *)constraints, "N_VNew_Parallel", 0, my_pe))
      MPI_Abort(comm, 1);
  N_VConst(ZERO, constraints);
  
  SetInitialProfiles(cc, sc);

  fnormtol = FTOL; scsteptol = STOL;

  /* Call KINCreate/KINInit to initialize KINSOL: 
     nvSpec  points to machine environment data
     A pointer to KINSOL problem memory is returned and stored in kmem. */
  kmem = KINCreate();
  if (check_flag((void *)kmem, "KINCreate", 0, my_pe)) MPI_Abort(comm, 1);

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

  flag = KINSetUserData(kmem, data);
  if (check_flag(&flag, "KINSetUserData", 1, my_pe)) MPI_Abort(comm, 1);

  flag = KINSetConstraints(kmem, constraints);
  if (check_flag(&flag, "KINSetConstraints", 1, my_pe)) MPI_Abort(comm, 1);

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

  flag = KINSetFuncNormTol(kmem, fnormtol);
  if (check_flag(&flag, "KINSetFuncNormTol", 1, my_pe)) MPI_Abort(comm, 1);

  flag = KINSetScaledStepTol(kmem, scsteptol);
  if (check_flag(&flag, "KINSetScaledStepTol", 1, my_pe)) MPI_Abort(comm, 1);
  
  /* Call KINBBDPrecInit to initialize and allocate memory for the
     band-block-diagonal preconditioner, and specify the local and
     communication functions func_local and gcomm=NULL (all communication
     needed for the func_local is already done in func). */
  dq_rel_uu = ZERO;
  mudq = mldq = 2*NUM_SPECIES - 1;
  mukeep = mlkeep = NUM_SPECIES;

  /* Call KINBBDSpgmr to specify the linear solver KINSPGMR */
  maxl = 20; maxlrst = 2;
  flag = KINSpgmr(kmem, maxl);
  if (check_flag(&flag, "KINSpgmr", 1, my_pe)) MPI_Abort(comm, 1);

  /* Initialize BBD preconditioner */
  flag = KINBBDPrecInit(kmem, Nlocal, mudq, mldq, mukeep, mlkeep,
                        dq_rel_uu, func_local, NULL);
  if (check_flag(&flag, "KINBBDPrecInit", 1, my_pe)) MPI_Abort(comm, 1);


  flag = KINSpilsSetMaxRestarts(kmem, maxlrst);
  if (check_flag(&flag, "KINSpilsSetMaxRestarts", 1, my_pe)) 
    MPI_Abort(comm, 1);

  /* Print out the problem size, solution parameters, initial guess. */
  if (my_pe == 0)
    PrintHeader(globalstrategy, maxl, maxlrst, mudq, mldq, mukeep,
		mlkeep, fnormtol, scsteptol);

  /* call KINSol and print output concentration profile */
  flag = KINSol(kmem,           /* KINSol memory block */
                cc,             /* initial guesss 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, my_pe)) MPI_Abort(comm, 1);

  if (my_pe == 0) printf("\n\nComputed equilibrium species concentrations:\n");
  if (my_pe == 0 || my_pe==npelast) PrintOutput(my_pe, comm, cc);
  
  /* Print final statistics and free memory */
  if (my_pe == 0) 
    PrintFinalStats(kmem);

  N_VDestroy_Parallel(cc);
  N_VDestroy_Parallel(sc);

  KINFree(&kmem);
  FreeUserData(data);

  MPI_Finalize();

  return(0);
}

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

  /* Create and allocate CVODES memory for backward run */

  printf("\nCreate and allocate CVODES memory for backward run\n");

  retval = CVodeCreateB(cvode_mem, CV_BDF, &indexB);
  if(check_retval(&retval, "CVodeCreateB", 1)) return(1);

  retval = CVodeSetUserDataB(cvode_mem, indexB, data);
  if(check_retval(&retval, "CVodeSetUserDataB", 1)) return(1);

  retval = CVodeInitB(cvode_mem, indexB, fB, TOUT, uB);
  if(check_retval(&retval, "CVodeInitB", 1)) return(1);

  retval = CVodeSStolerancesB(cvode_mem, indexB, reltolB, abstolB);
  if(check_retval(&retval, "CVodeSStolerancesB", 1)) return(1);
 
  /* Create banded SUNMatrix for the backward problem */
  AB = SUNBandMatrix(NEQ, MY, MY);
  if(check_retval((void *)AB, "SUNBandMatrix", 0)) return(1);

  /* Create banded SUNLinearSolver for the backward problem */
  LSB = SUNLinSol_Band(uB, AB);
  if(check_retval((void *)LSB, "SUNLinSol_Band", 0)) return(1);

  /* Attach the matrix and linear solver */
  retval = CVodeSetLinearSolverB(cvode_mem, indexB, LSB, AB);
  if(check_retval(&retval, "CVodeSetLinearSolverB", 1)) return(1);

  /* Set the user-supplied Jacobian routine for the backward problem */
  retval = CVodeSetJacFnB(cvode_mem, indexB, JacB);
  if(check_retval(&retval, "CVodeSetJacFnB", 1)) return(1);

  /* Perform backward integration */
  printf("\nBackward integration\n");
  retval = CVodeB(cvode_mem, T0, CV_NORMAL);
  if(check_retval(&retval, "CVodeB", 1)) return(1);

  retval = CVodeGetB(cvode_mem, indexB, &t, uB);
  if(check_retval(&retval, "CVodeGetB", 1)) return(1);

  PrintOutput(uB, data);

  N_VDestroy(u);   /* Free the u vector                      */
  N_VDestroy(uB);  /* Free the uB vector                     */
  CVodeFree(&cvode_mem);  /* Free the CVODE problem memory          */
  SUNLinSolFree(LS);      /* Free the forward linear solver memory  */
  SUNMatDestroy(A);       /* Free the forward matrix memory         */
  SUNLinSolFree(LSB);     /* Free the backward linear solver memory */
  SUNMatDestroy(AB);      /* Free the backward matrix memory        */

  free(data);             /* Free the user data */

  return(0);
}

int main(void)
{
  void *mem;
  N_Vector yy, yp, avtol;
  realtype rtol, *yval, *ypval, *atval;
  realtype t0, tout1, tout, tret;
  int iout, retval, retvalr;
  int rootsfound[2];
  SUNMatrix A;
  SUNLinearSolver LS;
  sunindextype nnz;

  mem = NULL;
  yy = yp = avtol = NULL;
  yval = ypval = atval = NULL;
  A = NULL;
  LS = NULL;

  /* Allocate N-vectors. */
  yy = N_VNew_Serial(NEQ);
  if(check_retval((void *)yy, "N_VNew_Serial", 0)) return(1);
  yp = N_VNew_Serial(NEQ);
  if(check_retval((void *)yp, "N_VNew_Serial", 0)) return(1);
  avtol = N_VNew_Serial(NEQ);
  if(check_retval((void *)avtol, "N_VNew_Serial", 0)) return(1);

  /* Create and initialize  y, y', and absolute tolerance vectors. */
  yval  = N_VGetArrayPointer(yy);
  yval[0] = ONE;
  yval[1] = ZERO;
  yval[2] = ZERO;

  ypval = N_VGetArrayPointer(yp);
  ypval[0]  = RCONST(-0.04);
  ypval[1]  = RCONST(0.04);
  ypval[2]  = ZERO;

  rtol = RCONST(1.0e-4);

  atval = N_VGetArrayPointer(avtol);
  atval[0] = RCONST(1.0e-8);
  atval[1] = RCONST(1.0e-6);
  atval[2] = RCONST(1.0e-6);

  /* Integration limits */
  t0 = ZERO;
  tout1 = RCONST(0.4);

  PrintHeader(rtol, avtol, yy);

  /* Call IDACreate and IDAInit to initialize IDA memory */
  mem = IDACreate();
  if(check_retval((void *)mem, "IDACreate", 0)) return(1);
  retval = IDAInit(mem, resrob, t0, yy, yp);
  if(check_retval(&retval, "IDAInit", 1)) return(1);
  /* Call IDASVtolerances to set tolerances */
  retval = IDASVtolerances(mem, rtol, avtol);
  if(check_retval(&retval, "IDASVtolerances", 1)) return(1);

  /* Free avtol */
  N_VDestroy(avtol);

  /* Call IDARootInit to specify the root function grob with 2 components */
  retval = IDARootInit(mem, 2, grob);
  if (check_retval(&retval, "IDARootInit", 1)) return(1);

  /* Create sparse SUNMatrix for use in linear solves */
  nnz = NEQ * NEQ;
  A = SUNSparseMatrix(NEQ, NEQ, nnz, CSC_MAT);
  if(check_retval((void *)A, "SUNSparseMatrix", 0)) return(1);

  /* Create SuperLUMT SUNLinearSolver object (one thread) */
  LS = SUNLinSol_SuperLUMT(yy, A, 1);
  if(check_retval((void *)LS, "SUNLinSol_SuperLUMT", 0)) return(1);

  /* Attach the matrix and linear solver */
  retval = IDASetLinearSolver(mem, LS, A);
  if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);

  /* Set the user-supplied Jacobian routine */
  retval = IDASetJacFn(mem, jacrob);
  if(check_retval(&retval, "IDASetJacFn", 1)) return(1);

  /* In loop, call IDASolve, print results, and test for error.
     Break out of loop when NOUT preset output times have been reached. */

  iout = 0; tout = tout1;
  while(1) {

    retval = IDASolve(mem, tout, &tret, yy, yp, IDA_NORMAL);

    PrintOutput(mem,tret,yy);

    if(check_retval(&retval, "IDASolve", 1)) return(1);

    if (retval == IDA_ROOT_RETURN) {
      retvalr = IDAGetRootInfo(mem, rootsfound);
      check_retval(&retvalr, "IDAGetRootInfo", 1);
      PrintRootInfo(rootsfound[0],rootsfound[1]);
    }
//.........这里部分代码省略.........

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

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

    printf("Sensitivity: NO ");

  }

  /* In loop over output points, call CVode, print results, test for error */

  printf("\n\n");
  printf("========================================================================\n");
  printf("     T     Q       H      NST                    Bottom left  Top right \n");
  printf("========================================================================\n");

  for (iout=1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {
    flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
    if(check_flag(&flag, "CVode", 1)) break;
    PrintOutput(cvode_mem, t, y);
    if (sensi) {
      flag = CVodeGetSens(cvode_mem, t, uS);
      if(check_flag(&flag, "CVodeGetSens", 1)) break;
      PrintOutputS(uS);
    }
    
    printf("------------------------------------------------------------------------\n");

  }

  /* Print final statistics */
  PrintFinalStats(cvode_mem, sensi);

  /* Free memory */
  N_VDestroy_Serial(y);
  if (sensi) {
    N_VDestroyVectorArray_Serial(uS, NS);
    free(pbar);
    free(plist);
  }
  FreeUserData(data);
  CVodeFree(&cvode_mem);

  return(0);
}

int main()
{
  UserData data;
  realtype fnormtol, scsteptol;
  N_Vector u1, u2, u, s, c;
  int glstr, mset, flag;
  void *kmem;

  u1 = u2 = u = NULL;
  s = c = NULL;
  kmem = NULL;
  data = NULL;

  /* User data */

  data = (UserData)malloc(sizeof *data);
  data->lb[0] = PT25;       data->ub[0] = ONE;
  data->lb[1] = ONEPT5;     data->ub[1] = TWO*PI;

  /* Create serial vectors of length NEQ */
  u1 = N_VNew_Serial(NEQ);
  if (check_flag((void *)u1, "N_VNew_Serial", 0)) return(1);

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

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

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

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

  SetInitialGuess1(u1,data);
  SetInitialGuess2(u2,data);

  N_VConst_Serial(ONE,s); /* no scaling */

  Ith(c,1) =  ZERO;   /* no constraint on x1 */
  Ith(c,2) =  ZERO;   /* no constraint on x2 */
  Ith(c,3) =  ONE;    /* l1 = x1 - x1_min >= 0 */
  Ith(c,4) = -ONE;    /* L1 = x1 - x1_max <= 0 */
  Ith(c,5) =  ONE;    /* l2 = x2 - x2_min >= 0 */
  Ith(c,6) = -ONE;    /* L2 = x2 - x22_min <= 0 */
  
  fnormtol=FTOL; scsteptol=STOL;


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

  flag = KINSetUserData(kmem, data);
  if (check_flag(&flag, "KINSetUserData", 1)) return(1);
  flag = KINSetConstraints(kmem, c);
  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);

  flag = KINInit(kmem, func, u);
  if (check_flag(&flag, "KINInit", 1)) return(1);

  /* Call KINDense to specify the linear solver */

  flag = KINDense(kmem, NEQ);
  if (check_flag(&flag, "KINDense", 1)) return(1);

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

  /* --------------------------- */

  printf("\n------------------------------------------\n");
  printf("\nInitial guess on lower bounds\n");
  printf("  [x1,x2] = ");
  PrintOutput(u1);

  N_VScale_Serial(ONE,u1,u);
  glstr = KIN_NONE;
  mset = 1;
  SolveIt(kmem, u, s, glstr, mset);

  /* --------------------------- */

  N_VScale_Serial(ONE,u1,u);
  glstr = KIN_LINESEARCH;
  mset = 1;
  SolveIt(kmem, u, s, glstr, mset);

  /* --------------------------- */

  N_VScale_Serial(ONE,u1,u);
  glstr = KIN_NONE;
  mset = 0;
  SolveIt(kmem, u, s, glstr, mset);

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

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 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()
{
  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 sta