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

static void PrintFinalStats(void *cvode_mem){	long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;	int flag;		flag = CVodeGetNumSteps(cvode_mem, &nst);	check_flag(&flag, "CVodeGetNumSteps", 1);	flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);	check_flag(&flag, "CVodeGetNumRhsEvals", 1);	flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);	check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);	flag = CVodeGetNumErrTestFails(cvode_mem, &netf);	check_flag(&flag, "CVodeGetNumErrTestFails", 1);	flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);	check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);	flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);	check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);		flag = CVDlsGetNumJacEvals(cvode_mem, &nje);	check_flag(&flag, "CVDlsGetNumJacEvals", 1);	flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);	check_flag(&flag, "CVDlsGetNumRhsEvals", 1);		flag = CVodeGetNumGEvals(cvode_mem, &nge);	check_flag(&flag, "CVodeGetNumGEvals", 1);		printf("/nFinal Statistics:/n");	printf("nst = %-6ld nfe  = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld/n",		   nst, nfe, nsetups, nfeLS, nje);	printf("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld/n /n",		   nni, ncfn, netf, nge);}

static void PrintFinalStats(void *cvode_mem){  long int lenrw, leniw ;  long int lenrwSPGMR, leniwSPGMR;  long int nst, nfe, nsetups, nni, ncfn, netf;  long int nli, npe, nps, ncfl, nfeSPGMR;  int flag;  flag = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);  flag = CVodeGetNumSteps(cvode_mem, &nst);  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  flag = CVSpilsGetWorkSpace(cvode_mem, &lenrwSPGMR, &leniwSPGMR);  flag = CVSpilsGetNumLinIters(cvode_mem, &nli);  flag = CVSpilsGetNumPrecEvals(cvode_mem, &npe);  flag = CVSpilsGetNumPrecSolves(cvode_mem, &nps);  flag = CVSpilsGetNumConvFails(cvode_mem, &ncfl);  flag = CVSpilsGetNumRhsEvals(cvode_mem, &nfeSPGMR);  printf("/nFinal Statistics.. /n/n");  printf("lenrw   = %6ld     leniw = %6ld/n", lenrw, leniw);  printf("llrw    = %6ld     lliw  = %6ld/n", lenrwSPGMR, leniwSPGMR);  printf("nst     = %6ld/n"                  , nst);  printf("nfe     = %6ld     nfel  = %6ld/n"  , nfe, nfeSPGMR);  printf("nni     = %6ld     nli   = %6ld/n"  , nni, nli);  printf("nsetups = %6ld     netf  = %6ld/n"  , nsetups, netf);  printf("npe     = %6ld     nps   = %6ld/n"  , npe, nps);  printf("ncfn    = %6ld     ncfl  = %6ld/n/n", ncfn, ncfl); }

static void PrintOutput(void *cvode_mem, realtype t){    long int nst, nfe, nni;    int qu, flag;    realtype hu;    flag = CVodeGetNumSteps(cvode_mem, &nst);    check_flag(&flag, "CVodeGetNumSteps", 1);    flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);    check_flag(&flag, "CVodeGetNumRhsEvals", 1);    flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);    check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);    flag = CVodeGetLastOrder(cvode_mem, &qu);    check_flag(&flag, "CVodeGetLastOrder", 1);    flag = CVodeGetLastStep(cvode_mem, &hu);    check_flag(&flag, "CVodeGetLastStep", 1);#if defined(SUNDIALS_EXTENDED_PRECISION)    printf("t = %10.2Le  nst = %ld  nfe = %ld  nni = %ld", t, nst, nfe, nni);    printf("  qu = %d  hu = %11.2Le/n/n", qu, hu);#elif defined(SUNDIALS_DOUBLE_PRECISION)    printf("t = %10.2le  nst = %ld  nfe = %ld  nni = %ld", t, nst, nfe, nni);    printf("  qu = %d  hu = %11.2le/n/n", qu, hu);#else    printf("t = %10.2e  nst = %ld  nfe = %ld  nni = %ld", t, nst, nfe, nni);    printf("  qu = %d  hu = %11.2e/n/n", qu, hu);#endif}

static void PrintOutput(void *cvode_mem, realtype t, N_Vector u){  long int nst;  int qu, flag;  realtype hu, *udata;    udata = N_VGetArrayPointer_Serial(u);  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetLastOrder(cvode_mem, &qu);  check_flag(&flag, "CVodeGetLastOrder", 1);  flag = CVodeGetLastStep(cvode_mem, &hu);  check_flag(&flag, "CVodeGetLastStep", 1);#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("%8.3Le %2d  %8.3Le %5ld/n", t, qu, hu, nst);#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("%8.3e %2d  %8.3e %5ld/n", t, qu, hu, nst);#else  printf("%8.3e %2d  %8.3e %5ld/n", t, qu, hu, nst);#endif  printf("                  Solution       ");#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("%12.4Le %12.4Le %12.4Le /n", udata[0], udata[1], udata[2]);#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("%12.4e %12.4e %12.4e /n", udata[0], udata[1], udata[2]);#else  printf("%12.4e %12.4e %12.4e /n", udata[0], udata[1], udata[2]);#endif}

static void PrintOutput(void *cvode_mem, realtype t, N_Vector u){  long int nst;  int qu, flag;  realtype hu;  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetLastOrder(cvode_mem, &qu);  check_flag(&flag, "CVodeGetLastOrder", 1);  flag = CVodeGetLastStep(cvode_mem, &hu);  check_flag(&flag, "CVodeGetLastStep", 1);#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("%8.3Le %2d  %8.3Le %5ld/n", t, qu, hu ,nst);#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("%8.3le %2d  %8.3le %5ld/n", t, qu, hu ,nst);#else  printf("%8.3e %2d  %8.3e %5ld/n", t, qu, hu ,nst);#endif  printf("                                Solution       ");#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("%12.4Le /n", N_VMaxNorm(u));#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("%12.4le /n", N_VMaxNorm(u));#else  printf("%12.4e /n", N_VMaxNorm(u));#endif}

/*======================================================================= Print the final statistics returned from the most recent call to CVODE-----------------------------------------------------------------------*/static void PrintFinalStatistics(void *cvode_mem){  long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf;  int flag;  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_flag(&flag, "CVodeGetNumRhsEvals", 1);  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_flag(&flag, "CVodeGetNumErrTestFails", 1);  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);  flag = CVDlsGetNumJacEvals(cvode_mem, &nje);  check_flag(&flag, "CVDlsGetNumJacEvals", 1);  flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);  check_flag(&flag, "CVDlsGetNumRhsEvals", 1);#pragma omp critical (statistics)#ifdef OPENMP    printf("/nFinal Statistics (thread %d):/n", omp_get_thread_num());#endif    printf("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld/n",	   nst, nfe, nsetups, nfeLS);    printf("nje = %-6ld nni = %-6ld ncfn = %-6ld    netf = %-6ld/n",	   nje, nni, ncfn, netf);}

static void PrintFinalStats(void *cvode_mem, booleantype sensi,                            booleantype err_con, int sensi_meth){  long int nst;  long int nfe, nsetups, nni, ncfn, netf;  long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;  int retval;  retval = CVodeGetNumSteps(cvode_mem, &nst);  check_retval(&retval, "CVodeGetNumSteps", 1, 0);  retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_retval(&retval, "CVodeGetNumRhsEvals", 1, 0);  retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_retval(&retval, "CVodeGetNumLinSolvSetups", 1, 0);  retval = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_retval(&retval, "CVodeGetNumErrTestFails", 1, 0);  retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1, 0);  retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1, 0);  if (sensi) {    retval = CVodeGetSensNumRhsEvals(cvode_mem, &nfSe);    check_retval(&retval, "CVodeGetSensNumRhsEvals", 1, 0);    retval = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);    check_retval(&retval, "CVodeGetNumRhsEvalsSens", 1, 0);    retval = CVodeGetSensNumLinSolvSetups(cvode_mem, &nsetupsS);    check_retval(&retval, "CVodeGetSensNumLinSolvSetups", 1, 0);    retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);    if (err_con) {      retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);      check_retval(&retval, "CVodeGetSensNumErrTestFails", 1, 0);    } else {      netfS = 0;    }    if ((sensi_meth == CV_STAGGERED) || (sensi_meth == CV_STAGGERED1)) {      retval = CVodeGetSensNumNonlinSolvIters(cvode_mem, &nniS);      check_retval(&retval, "CVodeGetSensNumNonlinSolvIters", 1, 0);      retval = CVodeGetSensNumNonlinSolvConvFails(cvode_mem, &ncfnS);      check_retval(&retval, "CVodeGetSensNumNonlinSolvConvFails", 1, 0);    } else {      nniS = 0;      ncfnS = 0;    }  }  printf("/nFinal Statistics/n/n");  printf("nst     = %5ld/n/n", nst);  printf("nfe     = %5ld/n",   nfe);  printf("netf    = %5ld    nsetups  = %5ld/n", netf, nsetups);  printf("nni     = %5ld    ncfn     = %5ld/n", nni, ncfn);  if(sensi) {    printf("/n");    printf("nfSe    = %5ld    nfeS     = %5ld/n", nfSe, nfeS);    printf("netfs   = %5ld    nsetupsS = %5ld/n", netfS, nsetupsS);    printf("nniS    = %5ld    ncfnS    = %5ld/n", nniS, ncfnS);  }}

void ode_solver_print_stats(const ode_solver* solver, FILE* outF){  long int nst;  long int nfe, nsetups, nni, ncfn, netf;  long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;  long int nje, nfeLS;  int flag;    void* cvode_mem = solver->cvode_mem;    flag = CVodeGetNumSteps(cvode_mem, &nst);    flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);    flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);    flag = CVodeGetNumErrTestFails(cvode_mem, &netf);    flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);    flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);      if (solver->yS != 0) {    flag = CVodeGetSensNumRhsEvals(cvode_mem, &nfSe);	    flag = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);        flag = CVodeGetSensNumLinSolvSetups(cvode_mem, &nsetupsS);        flag = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);        flag = CVodeGetSensNumNonlinSolvIters(cvode_mem, &nniS);        flag = CVodeGetSensNumNonlinSolvConvFails(cvode_mem, &ncfnS);      }    flag = CVDlsGetNumJacEvals(cvode_mem, &nje);    flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);      fprintf(outF,"/n# Solver Statistics/n/n");  fprintf(outF,"# Steps            = %5ld/n/n", nst);  fprintf(outF,"# RhsEvals         = %5ld/n",   nfe);  fprintf(outF,"# ErrTestFails     = %5ld   LinSolvSetups        = %5ld/n", netf, nsetups);  fprintf(outF,"# NonlinSolvIters  = %5ld   NonlinSolvConvFails  = %5ld/n", nni, ncfn);    if(solver->yS != 0) {    fprintf(outF,"/n# Sensitivities Statistics/n");    fprintf(outF,"# SensRhsEvals     = %5ld   RhsEvals             = %5ld/n", nfSe, nfeS);    fprintf(outF,"# ErrTestFails     = %5ld   LinSolvSetups        = %5ld/n", netfS, nsetupsS);    fprintf(outF,"# NonlinSolvIters  = %5ld   NonlinSolvConvFails  = %5ld/n", nniS, ncfnS);  }    fprintf(outF,"/n# Jacobian Statistics/n");  fprintf(outF,"# JacEvals  = %5ld    RhsEvals  = %5ld/n", nje, nfeLS);  }

	int OpenSMOKE_CVODE_Sundials<T>::GetNumberOfSteps() const	{		long int nst;		int flag = CVodeGetNumSteps(cvode_mem_, &nst);		check_flag(&flag, std::string("CVodeGetNumSteps"), 1);		return int(nst);	}

static void PrintFinalStats(void *cvode_mem, booleantype sensi){  long int nst;  long int nfe, nsetups, nni, ncfn, netf;  long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;  long int nje, nfeLS;  int retval;  retval = CVodeGetNumSteps(cvode_mem, &nst);  check_retval(&retval, "CVodeGetNumSteps", 1);  retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_retval(&retval, "CVodeGetNumRhsEvals", 1);  retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);  retval = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_retval(&retval, "CVodeGetNumErrTestFails", 1);  retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);  retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);  if (sensi) {    retval = CVodeGetSensNumRhsEvals(cvode_mem, &nfSe);    check_retval(&retval, "CVodeGetSensNumRhsEvals", 1);    retval = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);    check_retval(&retval, "CVodeGetNumRhsEvalsSens", 1);    retval = CVodeGetSensNumLinSolvSetups(cvode_mem, &nsetupsS);    check_retval(&retval, "CVodeGetSensNumLinSolvSetups", 1);    retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);    check_retval(&retval, "CVodeGetSensNumErrTestFails", 1);    retval = CVodeGetSensNumNonlinSolvIters(cvode_mem, &nniS);    check_retval(&retval, "CVodeGetSensNumNonlinSolvIters", 1);    retval = CVodeGetSensNumNonlinSolvConvFails(cvode_mem, &ncfnS);    check_retval(&retval, "CVodeGetSensNumNonlinSolvConvFails", 1);  }  retval = CVDlsGetNumJacEvals(cvode_mem, &nje);  check_retval(&retval, "CVDlsGetNumJacEvals", 1);  retval = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);  check_retval(&retval, "CVDlsGetNumRhsEvals", 1);  printf("/nFinal Statistics/n/n");  printf("nst     = %5ld/n/n", nst);  printf("nfe     = %5ld/n",   nfe);  printf("netf    = %5ld    nsetups  = %5ld/n", netf, nsetups);  printf("nni     = %5ld    ncfn     = %5ld/n", nni, ncfn);  if(sensi) {    printf("/n");    printf("nfSe    = %5ld    nfeS     = %5ld/n", nfSe, nfeS);    printf("netfs   = %5ld    nsetupsS = %5ld/n", netfS, nsetupsS);    printf("nniS    = %5ld    ncfnS    = %5ld/n", nniS, ncfnS);  }  printf("/n");  printf("nje    = %5ld    nfeLS     = %5ld/n", nje, nfeLS);}

static void PrintFinalStats(void *cvode_mem, booleantype sensi){  long int nst;  long int nfe, nsetups, nni, ncfn, netf;  long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;  long int njeD, nfeD;  int flag;  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_flag(&flag, "CVodeGetNumRhsEvals", 1);  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_flag(&flag, "CVodeGetNumErrTestFails", 1);  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);  if (sensi) {    flag = CVodeGetNumSensRhsEvals(cvode_mem, &nfSe);    check_flag(&flag, "CVodeGetNumSensRhsEvals", 1);    flag = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);    check_flag(&flag, "CVodeGetNumRhsEvalsSens", 1);    flag = CVodeGetNumSensLinSolvSetups(cvode_mem, &nsetupsS);    check_flag(&flag, "CVodeGetNumSensLinSolvSetups", 1);    flag = CVodeGetNumSensErrTestFails(cvode_mem, &netfS);    check_flag(&flag, "CVodeGetNumSensErrTestFails", 1);    flag = CVodeGetNumSensNonlinSolvIters(cvode_mem, &nniS);    check_flag(&flag, "CVodeGetNumSensNonlinSolvIters", 1);    flag = CVodeGetNumSensNonlinSolvConvFails(cvode_mem, &ncfnS);    check_flag(&flag, "CVodeGetNumSensNonlinSolvConvFails", 1);  }  flag = CVDenseGetNumJacEvals(cvode_mem, &njeD);  check_flag(&flag, "CVDenseGetNumJacEvals", 1);  flag = CVDenseGetNumRhsEvals(cvode_mem, &nfeD);  check_flag(&flag, "CVDenseGetNumRhsEvals", 1);  printf("/nFinal Statistics/n/n");  printf("nst     = %5ld/n/n", nst);  printf("nfe     = %5ld/n",   nfe);  printf("netf    = %5ld    nsetups  = %5ld/n", netf, nsetups);  printf("nni     = %5ld    ncfn     = %5ld/n", nni, ncfn);  if(sensi) {    printf("/n");    printf("nfSe    = %5ld    nfeS     = %5ld/n", nfSe, nfeS);    printf("netfs   = %5ld    nsetupsS = %5ld/n", netfS, nsetupsS);    printf("nniS    = %5ld    ncfnS    = %5ld/n", nniS, ncfnS);  }  printf("/n");  printf("njeD    = %5ld    nfeD     = %5ld/n", njeD, nfeD);}

static void PrintFinalStats(void *cvode_mem){  long int lenrw, leniw ;  long int lenrwLS, leniwLS;  long int lenrwBBDP, leniwBBDP, ngevalsBBDP;  long int nst, nfe, nsetups, nni, ncfn, netf;  long int nli, npe, nps, ncfl, nfeLS;  int flag;  flag = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);  check_flag(&flag, "CVodeGetWorkSpace", 1, 0);  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1, 0);  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_flag(&flag, "CVodeGetNumRhsEvals", 1, 0);  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_flag(&flag, "CVodeGetNumLinSolvSetups", 1, 0);  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_flag(&flag, "CVodeGetNumErrTestFails", 1, 0);  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1, 0);  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1, 0);  flag = CVSpilsGetWorkSpace(cvode_mem, &lenrwLS, &leniwLS);  check_flag(&flag, "CVSpilsGetWorkSpace", 1, 0);  flag = CVSpilsGetNumLinIters(cvode_mem, &nli);  check_flag(&flag, "CVSpilsGetNumLinIters", 1, 0);  flag = CVSpilsGetNumPrecEvals(cvode_mem, &npe);  check_flag(&flag, "CVSpilsGetNumPrecEvals", 1, 0);  flag = CVSpilsGetNumPrecSolves(cvode_mem, &nps);  check_flag(&flag, "CVSpilsGetNumPrecSolves", 1, 0);  flag = CVSpilsGetNumConvFails(cvode_mem, &ncfl);  check_flag(&flag, "CVSpilsGetNumConvFails", 1, 0);  flag = CVSpilsGetNumRhsEvals(cvode_mem, &nfeLS);  check_flag(&flag, "CVSpilsGetNumRhsEvals", 1, 0);  printf("/nFinal Statistics: /n/n");  printf("lenrw   = %5ld     leniw   = %5ld/n", lenrw, leniw);  printf("lenrwls = %5ld     leniwls = %5ld/n", lenrwLS, leniwLS);  printf("nst     = %5ld/n"                  , nst);  printf("nfe     = %5ld     nfels   = %5ld/n"  , nfe, nfeLS);  printf("nni     = %5ld     nli     = %5ld/n"  , nni, nli);  printf("nsetups = %5ld     netf    = %5ld/n"  , nsetups, netf);  printf("npe     = %5ld     nps     = %5ld/n"  , npe, nps);  printf("ncfn    = %5ld     ncfl    = %5ld/n/n", ncfn, ncfl);  flag = CVBBDPrecGetWorkSpace(cvode_mem, &lenrwBBDP, &leniwBBDP);  check_flag(&flag, "CVBBDPrecGetWorkSpace", 1, 0);  flag = CVBBDPrecGetNumGfnEvals(cvode_mem, &ngevalsBBDP);  check_flag(&flag, "CVBBDPrecGetNumGfnEvals", 1, 0);  printf("In CVBBDPRE: real/integer local work space sizes = %ld, %ld/n",	 lenrwBBDP, leniwBBDP);    printf("             no. flocal evals. = %ld/n",ngevalsBBDP);}

static void PrintFinalStats(void *cvode_mem){  long int lenrw, leniw ;  long int lenrwLS, leniwLS;  long int lenrwBP, leniwBP;  long int nst, nfe, nsetups, nni, ncfn, netf;  long int nli, npe, nps, ncfl, nfeLS;  long int nfeBP;  int flag;  flag = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);  check_flag(&flag, "CVodeGetWorkSpace", 1);  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_flag(&flag, "CVodeGetNumRhsEvals", 1);  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_flag(&flag, "CVodeGetNumErrTestFails", 1);  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);  flag = CVSpilsGetWorkSpace(cvode_mem, &lenrwLS, &leniwLS);  check_flag(&flag, "CVSpilsGetWorkSpace", 1);  flag = CVSpilsGetNumLinIters(cvode_mem, &nli);  check_flag(&flag, "CVSpilsGetNumLinIters", 1);  flag = CVSpilsGetNumPrecEvals(cvode_mem, &npe);  check_flag(&flag, "CVSpilsGetNumPrecEvals", 1);  flag = CVSpilsGetNumPrecSolves(cvode_mem, &nps);  check_flag(&flag, "CVSpilsGetNumPrecSolves", 1);  flag = CVSpilsGetNumConvFails(cvode_mem, &ncfl);  check_flag(&flag, "CVSpilsGetNumConvFails", 1);  flag = CVSpilsGetNumRhsEvals(cvode_mem, &nfeLS);  check_flag(&flag, "CVSpilsGetNumRhsEvals", 1);  flag = CVBandPrecGetWorkSpace(cvode_mem, &lenrwBP, &leniwBP);  check_flag(&flag, "CVBandPrecGetWorkSpace", 1);  flag = CVBandPrecGetNumRhsEvals(cvode_mem, &nfeBP);  check_flag(&flag, "CVBandPrecGetNumRhsEvals", 1);  printf("/nFinal Statistics.. /n/n");  printf("lenrw   = %5ld     leniw   = %5ld/n", lenrw, leniw);  printf("lenrwls = %5ld     leniwls = %5ld/n", lenrwLS, leniwLS);  printf("lenrwbp = %5ld     leniwbp = %5ld/n", lenrwBP, leniwBP);  printf("nst     = %5ld/n"                  , nst);  printf("nfe     = %5ld     nfetot  = %5ld/n"  , nfe, nfe+nfeLS+nfeBP);  printf("nfeLS   = %5ld     nfeBP   = %5ld/n"  , nfeLS, nfeBP);  printf("nni     = %5ld     nli     = %5ld/n"  , nni, nli);  printf("nsetups = %5ld     netf    = %5ld/n"  , nsetups, netf);  printf("npe     = %5ld     nps     = %5ld/n"  , npe, nps);  printf("ncfn    = %5ld     ncfl    = %5ld/n/n", ncfn, ncfl);}

static void PrintOutput(void *cvode_mem, int my_pe, MPI_Comm comm,                        N_Vector u, realtype t){  int qu, flag;  realtype hu, *udata, tempu[2];  int npelast;  long int i0, i1, nst;  MPI_Status status;  npelast = NPEX*NPEY - 1;  udata = NV_DATA_P(u);  /* Send c1,c2 at top right mesh point to PE 0 */  if (my_pe == npelast) {    i0 = NVARS*MXSUB*MYSUB - 2;    i1 = i0 + 1;    if (npelast != 0)      MPI_Send(&udata[i0], 2, PVEC_REAL_MPI_TYPE, 0, 0, comm);    else {      tempu[0] = udata[i0];      tempu[1] = udata[i1];    }  }  /* On PE 0, receive c1,c2 at top right, then print performance data     and sampled solution values */   if (my_pe == 0) {    if (npelast != 0)      MPI_Recv(&tempu[0], 2, PVEC_REAL_MPI_TYPE, npelast, 0, comm, &status);    flag = CVodeGetNumSteps(cvode_mem, &nst);    check_flag(&flag, "CVodeGetNumSteps", 1, my_pe);    flag = CVodeGetLastOrder(cvode_mem, &qu);    check_flag(&flag, "CVodeGetLastOrder", 1, my_pe);    flag = CVodeGetLastStep(cvode_mem, &hu);    check_flag(&flag, "CVodeGetLastStep", 1, my_pe);#if defined(SUNDIALS_EXTENDED_PRECISION)    printf("t = %.2Le   no. steps = %ld   order = %d   stepsize = %.2Le/n",           t, nst, qu, hu);    printf("At bottom left:  c1, c2 = %12.3Le %12.3Le /n", udata[0], udata[1]);    printf("At top right:    c1, c2 = %12.3Le %12.3Le /n/n", tempu[0], tempu[1]);#elif defined(SUNDIALS_DOUBLE_PRECISION)    printf("t = %.2le   no. steps = %ld   order = %d   stepsize = %.2le/n",           t, nst, qu, hu);    printf("At bottom left:  c1, c2 = %12.3le %12.3le /n", udata[0], udata[1]);    printf("At top right:    c1, c2 = %12.3le %12.3le /n/n", tempu[0], tempu[1]);#else    printf("t = %.2e   no. steps = %ld   order = %d   stepsize = %.2e/n",           t, nst, qu, hu);    printf("At bottom left:  c1, c2 = %12.3e %12.3e /n", udata[0], udata[1]);    printf("At top right:    c1, c2 = %12.3e %12.3e /n/n", tempu[0], tempu[1]);#endif  }}

static void PrintFinalStats(void *cvode_mem, int linsolver){  long int lenrw, leniw ;  long int lenrwLS, leniwLS;  long int nst, nfe, nsetups, nni, ncfn, netf;  long int nli, npe, nps, ncfl, nfeLS;  int retval;  retval = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);  check_retval(&retval, "CVodeGetWorkSpace", 1);  retval = CVodeGetNumSteps(cvode_mem, &nst);  check_retval(&retval, "CVodeGetNumSteps", 1);  retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_retval(&retval, "CVodeGetNumRhsEvals", 1);  retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);  retval = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_retval(&retval, "CVodeGetNumErrTestFails", 1);  retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);  retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);  retval = CVodeGetLinWorkSpace(cvode_mem, &lenrwLS, &leniwLS);  check_retval(&retval, "CVodeGetLinWorkSpace", 1);  retval = CVodeGetNumLinIters(cvode_mem, &nli);  check_retval(&retval, "CVodeGetNumLinIters", 1);  retval = CVodeGetNumPrecEvals(cvode_mem, &npe);  check_retval(&retval, "CVodeGetNumPrecEvals", 1);  retval = CVodeGetNumPrecSolves(cvode_mem, &nps);  check_retval(&retval, "CVodeGetNumPrecSolves", 1);  retval = CVodeGetNumLinConvFails(cvode_mem, &ncfl);  check_retval(&retval, "CVodeGetNumLinConvFails", 1);  retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS);  check_retval(&retval, "CVodeGetNumLinRhsEvals", 1);  printf("/nFinal Statistics.. /n/n");  printf("lenrw   = %5ld     leniw   = %5ld/n"  , lenrw, leniw);  printf("lenrwLS = %5ld     leniwLS = %5ld/n"  , lenrwLS, leniwLS);  printf("nst     = %5ld/n"                     , nst);  printf("nfe     = %5ld     nfeLS   = %5ld/n"  , nfe, nfeLS);  printf("nni     = %5ld     nli     = %5ld/n"  , nni, nli);  printf("nsetups = %5ld     netf    = %5ld/n"  , nsetups, netf);  printf("npe     = %5ld     nps     = %5ld/n"  , npe, nps);  printf("ncfn    = %5ld     ncfl    = %5ld/n/n", ncfn, ncfl);  if (linsolver < 2)    printf("======================================================================/n/n");}

static void PrintFinalStats(void *cvode_mem){  long int nst, nfe, nni, ncfn, netf;  int retval;  retval = CVodeGetNumSteps(cvode_mem, &nst);  check_retval(&retval, "CVodeGetNumSteps", 1, 0);  retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_retval(&retval, "CVodeGetNumRhsEvals", 1, 0);  retval = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_retval(&retval, "CVodeGetNumErrTestFails", 1, 0);  retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1, 0);  retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1, 0);  printf("/nFinal Statistics: /n/n");  printf("nst = %-6ld  nfe  = %-6ld  ", nst, nfe);  printf("nni = %-6ld  ncfn = %-6ld  netf = %ld/n /n", nni, ncfn, netf);}

static void PrintOutput(void *cvode_mem, realtype t, N_Vector y){    long int nst;  int qu, flag;  realtype hu;  realtype *ydata;  ydata = NV_DATA_S(y);  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetLastOrder(cvode_mem, &qu);  check_flag(&flag, "CVodeGetLastOrder", 1);  flag = CVodeGetLastStep(cvode_mem, &hu);  check_flag(&flag, "CVodeGetLastStep", 1);#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("%8.3Le %2d  %8.3Le %5ld/n", t,qu,hu,nst);#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("%8.3le %2d  %8.3le %5ld/n", t,qu,hu,nst);#else  printf("%8.3e %2d  %8.3e %5ld/n", t,qu,hu,nst);#endif  printf("                                Solution       ");#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("%12.4Le %12.4Le /n", IJKth(ydata,1,0,0), IJKth(ydata,1,MX-1,MZ-1)); #elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("%12.4le %12.4le /n", IJKth(ydata,1,0,0), IJKth(ydata,1,MX-1,MZ-1)); #else  printf("%12.4e %12.4e /n", IJKth(ydata,1,0,0), IJKth(ydata,1,MX-1,MZ-1)); #endif  printf("                                               ");#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("%12.4Le %12.4Le /n", IJKth(ydata,2,0,0), IJKth(ydata,2,MX-1,MZ-1));#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("%12.4le %12.4le /n", IJKth(ydata,2,0,0), IJKth(ydata,2,MX-1,MZ-1));#else  printf("%12.4e %12.4e /n", IJKth(ydata,2,0,0), IJKth(ydata,2,MX-1,MZ-1));#endif}

static void PrintOutput(void *cvode_mem, N_Vector u,realtype t){  long int nst;  int qu, flag;  realtype hu, *udata;  int mxh = MX/2 - 1, myh = MY/2 - 1, mx1 = MX - 1, my1 = MY - 1;  udata = N_VGetArrayPointer_Serial(u);  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetLastOrder(cvode_mem, &qu);  check_flag(&flag, "CVodeGetLastOrder", 1);  flag = CVodeGetLastStep(cvode_mem, &hu);  check_flag(&flag, "CVodeGetLastStep", 1);#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("t = %.2Le   no. steps = %ld   order = %d   stepsize = %.2Le/n",         t, nst, qu, hu);  printf("c1 (bot.left/middle/top rt.) = %12.3Le  %12.3Le  %12.3Le/n",         IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));  printf("c2 (bot.left/middle/top rt.) = %12.3Le  %12.3Le  %12.3Le/n/n",         IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("t = %.2e   no. steps = %ld   order = %d   stepsize = %.2e/n",         t, nst, qu, hu);  printf("c1 (bot.left/middle/top rt.) = %12.3e  %12.3e  %12.3e/n",         IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));  printf("c2 (bot.left/middle/top rt.) = %12.3e  %12.3e  %12.3e/n/n",         IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));#else  printf("t = %.2e   no. steps = %ld   order = %d   stepsize = %.2e/n",         t, nst, qu, hu);  printf("c1 (bot.left/middle/top rt.) = %12.3e  %12.3e  %12.3e/n",         IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));  printf("c2 (bot.left/middle/top rt.) = %12.3e  %12.3e  %12.3e/n/n",         IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));#endif}

static void PrintOutput(void *cvode_mem, int my_pe, realtype t, N_Vector u){  long int nst;  int qu, retval;  realtype hu, umax;  retval = CVodeGetNumSteps(cvode_mem, &nst);  check_retval(&retval, "CVodeGetNumSteps", 1, my_pe);  retval = CVodeGetLastOrder(cvode_mem, &qu);  check_retval(&retval, "CVodeGetLastOrder", 1, my_pe);  retval = CVodeGetLastStep(cvode_mem, &hu);  check_retval(&retval, "CVodeGetLastStep", 1, my_pe);  umax = N_VMaxNorm(u);  if (my_pe == 0) {#if defined(SUNDIALS_EXTENDED_PRECISION)    printf("%8.3Le %2d  %8.3Le %5ld/n", t,qu,hu,nst);#elif defined(SUNDIALS_DOUBLE_PRECISION)    printf("%8.3e %2d  %8.3e %5ld/n", t,qu,hu,nst);#else    printf("%8.3e %2d  %8.3e %5ld/n", t,qu,hu,nst);#endif    printf("                                Solution       ");#if defined(SUNDIALS_EXTENDED_PRECISION)    printf("%12.4Le /n", umax);#elif defined(SUNDIALS_DOUBLE_PRECISION)    printf("%12.4e /n", umax);#else    printf("%12.4e /n", umax);#endif  }  }

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

int main(int argc, char *argv[]){  realtype dx, reltol, abstol, t, tout, umax;  N_Vector u;  UserData data;  void *cvode_mem;  int iout, flag, my_pe, npes;  long int local_N, nperpe, nrem, my_base, nst;  MPI_Comm comm;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Get processor number, total number of pe's, and my_pe. */  MPI_Init(&argc, &argv);  comm = MPI_COMM_WORLD;  MPI_Comm_size(comm, &npes);  MPI_Comm_rank(comm, &my_pe);  /* Set local vector length. */  nperpe = NEQ/npes;  nrem = NEQ - npes*nperpe;  local_N = (my_pe < nrem) ? nperpe+1 : nperpe;  my_base = (my_pe < nrem) ? my_pe*local_N : my_pe*nperpe + nrem;  data = (UserData) malloc(sizeof *data);  /* Allocate data memory */  if(check_flag((void *)data, "malloc", 2, my_pe)) MPI_Abort(comm, 1);  data->comm = comm;  data->npes = npes;  data->my_pe = my_pe;  u = N_VNew_Parallel(comm, local_N, NEQ);  /* Allocate u vector */  if(check_flag((void *)u, "N_VNew", 0, my_pe)) MPI_Abort(comm, 1);  reltol = ZERO;  /* Set the tolerances */  abstol = ATOL;  dx = data->dx = XMAX/((realtype)(MX+1));  /* Set grid coefficients in data */  data->hdcoef = RCONST(1.0)/(dx*dx);  data->hacoef = RCONST(0.5)/(RCONST(2.0)*dx);  SetIC(u, dx, local_N, my_base);  /* Initialize u vector */  /*      Call CVodeCreate to create the solver memory:          CV_ADAMS   specifies the Adams Method     CV_FUNCTIONAL  specifies functional iteration     A pointer to the integrator memory is returned and stored in cvode_mem.  */  cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);  if(check_flag((void *)cvode_mem, "CVodeCreate", 0, my_pe)) MPI_Abort(comm, 1);  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);  if (my_pe == 0) PrintIntro(npes);  umax = N_VMaxNorm(u);  if (my_pe == 0) {    t = T0;    PrintData(t, umax, 0);  }  /* In loop over output points, call CVode, print results, test for error */  for (iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {    flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);    if(check_flag(&flag, "CVode", 1, my_pe)) break;    umax = N_VMaxNorm(u);    flag = CVodeGetNumSteps(cvode_mem, &nst);    check_flag(&flag, "CVodeGetNumSteps", 1, my_pe);    if (my_pe == 0) PrintData(t, umax, nst);  }  if (my_pe == 0)     PrintFinalStats(cvode_mem);  /* Print some final statistics */  N_VDestroy_Parallel(u);        /* Free the u vector *///.........这里部分代码省略.........

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

static void PrintFinalStats(void *cvode_mem){    long int lenrw, leniw ;    long int lenrwLS, leniwLS;    long int nst, nfe, nsetups, nni, ncfn, netf;    long int nli, npe, nps, ncfl, nfeLS;    int flag;    realtype avdim;    flag = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);    check_flag(&flag, "CVodeGetWorkSpace", 1);    flag = CVodeGetNumSteps(cvode_mem, &nst);    check_flag(&flag, "CVodeGetNumSteps", 1);    flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);    check_flag(&flag, "CVodeGetNumRhsEvals", 1);    flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);    check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);    flag = CVodeGetNumErrTestFails(cvode_mem, &netf);    check_flag(&flag, "CVodeGetNumErrTestFails", 1);    flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);    check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);    flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);    check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);    flag = CVSpilsGetWorkSpace(cvode_mem, &lenrwLS, &leniwLS);    check_flag(&flag, "CVSpilsGetWorkSpace", 1);    flag = CVSpilsGetNumLinIters(cvode_mem, &nli);    check_flag(&flag, "CVSpilsGetNumLinIters", 1);    flag = CVSpilsGetNumPrecEvals(cvode_mem, &npe);    check_flag(&flag, "CVSpilsGetNumPrecEvals", 1);    flag = CVSpilsGetNumPrecSolves(cvode_mem, &nps);    check_flag(&flag, "CVSpilsGetNumPrecSolves", 1);    flag = CVSpilsGetNumConvFails(cvode_mem, &ncfl);    check_flag(&flag, "CVSpilsGetNumConvFails", 1);    flag = CVSpilsGetNumRhsEvals(cvode_mem, &nfeLS);    check_flag(&flag, "CVSpilsGetNumRhsEvals", 1);    printf("/n/n Final statistics for this run:/n/n");    printf(" CVode real workspace length           = %4ld /n", lenrw);    printf(" CVode integer workspace length        = %4ld /n", leniw);    printf(" CVSPGMR real workspace length         = %4ld /n", lenrwLS);    printf(" CVSPGMR integer workspace length      = %4ld /n", leniwLS);    printf(" Number of steps                       = %4ld /n", nst);    printf(" Number of f-s                         = %4ld /n", nfe);    printf(" Number of f-s (SPGMR)                 = %4ld /n", nfeLS);    printf(" Number of f-s (TOTAL)                 = %4ld /n", nfe + nfeLS);    printf(" Number of setups                      = %4ld /n", nsetups);    printf(" Number of nonlinear iterations        = %4ld /n", nni);    printf(" Number of linear iterations           = %4ld /n", nli);    printf(" Number of preconditioner evaluations  = %4ld /n", npe);    printf(" Number of preconditioner solves       = %4ld /n", nps);    printf(" Number of error test failures         = %4ld /n", netf);    printf(" Number of nonlinear conv. failures    = %4ld /n", ncfn);    printf(" Number of linear convergence failures = %4ld /n", ncfl);    avdim = (nni > 0) ? ((realtype)nli)/((realtype)nni) : ZERO;#if defined(SUNDIALS_EXTENDED_PRECISION)    printf(" Average Krylov subspace dimension     = %.3Lf /n", avdim);#else    printf(" Average Krylov subspace dimension     = %.3f /n", avdim);#endif    printf("/n/n--------------------------------------------------------------");    printf("--------------/n");    printf(    "--------------------------------------------------------------");    printf("--------------/n");}

int main(int argc, char *argv[]){  realtype dx, reltol, abstol, t, tout, umax;  N_Vector u;  UserData data;  void *cvode_mem;  int iout, flag, my_pe, npes;  long int nst;  HYPRE_Int local_N, nperpe, nrem, my_base;  HYPRE_ParVector Upar; /* Declare HYPRE parallel vector */  HYPRE_IJVector  Uij;  /* Declare "IJ" interface to HYPRE vector */  MPI_Comm comm;  u = NULL;  data = NULL;  cvode_mem = NULL;  /* Get processor number, total number of pe's, and my_pe. */  MPI_Init(&argc, &argv);  comm = MPI_COMM_WORLD;  MPI_Comm_size(comm, &npes);  MPI_Comm_rank(comm, &my_pe);  /* Set partitioning. */  nperpe = NEQ/npes;  nrem = NEQ - npes*nperpe;  local_N = (my_pe < nrem) ? nperpe+1 : nperpe;  my_base = (my_pe < nrem) ? my_pe*local_N : my_pe*nperpe + nrem;    /* Allocate hypre vector */  HYPRE_IJVectorCreate(comm, my_base, my_base + local_N - 1, &Uij);  HYPRE_IJVectorSetObjectType(Uij, HYPRE_PARCSR);  HYPRE_IJVectorInitialize(Uij);  /* Allocate user defined data */  data = (UserData) malloc(sizeof *data);  /* Allocate data memory */  if(check_flag((void *)data, "malloc", 2, my_pe)) MPI_Abort(comm, 1);  data->comm = comm;  data->npes = npes;  data->my_pe = my_pe;  reltol = ZERO;  /* Set the tolerances */  abstol = ATOL;  dx = data->dx = XMAX/((realtype)(MX+1));  /* Set grid coefficients in data */  data->hdcoef = RCONST(1.0)/(dx*dx);  data->hacoef = RCONST(0.5)/(RCONST(2.0)*dx);  /* Initialize solutin vector. */  SetIC(Uij, dx, local_N, my_base);  HYPRE_IJVectorAssemble(Uij);  HYPRE_IJVectorGetObject(Uij, (void**) &Upar);  u = N_VMake_ParHyp(Upar);  /* Create wrapper u around hypre vector */  if(check_flag((void *)u, "N_VNew", 0, my_pe)) MPI_Abort(comm, 1);    /* Call CVodeCreate to create the solver memory and specify the    * Adams-Moulton LMM and the use of a functional iteration */  cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);  if(check_flag((void *)cvode_mem, "CVodeCreate", 0, my_pe)) MPI_Abort(comm, 1);  flag = CVodeSetUserData(cvode_mem, data);  if(check_flag(&flag, "CVodeSetUserData", 1, my_pe)) MPI_Abort(comm, 1);  /* Call CVodeInit to initialize the integrator memory and specify the   * user's right hand side function in u'=f(t,u), the inital time T0, and   * the initial dependent variable vector u. */  flag = CVodeInit(cvode_mem, f, T0, u);  if(check_flag(&flag, "CVodeInit", 1, my_pe)) 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, my_pe)) return(1);  if (my_pe == 0) PrintIntro(npes);  umax = N_VMaxNorm(u);  if (my_pe == 0) {    t = T0;    PrintData(t, umax, 0);  }  /* In loop over output points, call CVode, print results, test for error */  for (iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {    flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);    if(check_flag(&flag, "CVode", 1, my_pe)) break;    umax = N_VMaxNorm(u);    flag = CVodeGetNumSteps(cvode_mem, &nst);    check_flag(&flag, "CVodeGetNumSteps", 1, my_pe);    if (my_pe == 0) PrintData(t, umax, nst);  }  if (my_pe == 0)     PrintFinalStats(cvode_mem);  /* Print some final statistics *///.........这里部分代码省略.........

//.........这里部分代码省略.........  retval = CVodeQuadInit(cvode_mem, fQ, q);  if (check_retval(&retval, "CVodeQuadInit", 1)) return(1);  /* Call CVodeSetQuadErrCon to specify whether or not the quadrature variables     are to be used in the step size control mechanism within CVODES. Call     CVodeQuadSStolerances or CVodeQuadSVtolerances to specify the integration     tolerances for the quadrature variables. */  retval = CVodeSetQuadErrCon(cvode_mem, SUNTRUE);  if (check_retval(&retval, "CVodeSetQuadErrCon", 1)) return(1);  /* Call CVodeQuadSStolerances to specify scalar relative and absolute     tolerances. */  retval = CVodeQuadSStolerances(cvode_mem, reltolQ, abstolQ);  if (check_retval(&retval, "CVodeQuadSStolerances", 1)) return(1);  /* Allocate global memory */  /* Call CVodeAdjInit to update CVODES memory block by allocting the internal      memory needed for backward integration.*/  steps = STEPS; /* no. of integration steps between two consecutive ckeckpoints*/  retval = CVodeAdjInit(cvode_mem, steps, CV_HERMITE);  /*  retval = CVodeAdjInit(cvode_mem, steps, CV_POLYNOMIAL);  */  if (check_retval(&retval, "CVodeAdjInit", 1)) return(1);  /* Perform forward run */  printf("Forward integration ... ");  /* Call CVodeF to integrate the forward problem over an interval in time and     saves checkpointing data */  retval = CVodeF(cvode_mem, TOUT, y, &time, CV_NORMAL, &ncheck);  if (check_retval(&retval, "CVodeF", 1)) return(1);  retval = CVodeGetNumSteps(cvode_mem, &nst);  if (check_retval(&retval, "CVodeGetNumSteps", 1)) return(1);  printf("done ( nst = %ld )/n",nst);  printf("/nncheck = %d/n/n", ncheck);  retval = CVodeGetQuad(cvode_mem, &time, q);  if (check_retval(&retval, "CVodeGetQuad", 1)) return(1);  printf("--------------------------------------------------------/n");#if defined(SUNDIALS_EXTENDED_PRECISION)  printf("G:          %12.4Le /n",Ith(q,1));#elif defined(SUNDIALS_DOUBLE_PRECISION)  printf("G:          %12.4e /n",Ith(q,1));#else  printf("G:          %12.4e /n",Ith(q,1));#endif  printf("--------------------------------------------------------/n/n");  /* Test check point linked list      (uncomment next block to print check point information) */    /*  {    int i;        printf("/nList of Check Points (ncheck = %d)/n/n", ncheck);    ckpnt = (CVadjCheckPointRec *) malloc ( (ncheck+1)*sizeof(CVadjCheckPointRec));    CVodeGetAdjCheckPointsInfo(cvode_mem, ckpnt);    for (i=0;i<=ncheck;i++) {      printf("Address:       %p/n",ckpnt[i].my_addr);      printf("Next:          %p/n",ckpnt[i].next_addr);      printf("Time interval: %le  %le/n",ckpnt[i].t0, ckpnt[i].t1);

/*  * Get and print some final statistics */void PrintFinalStats(struct Integrator* integrator){  void* cvode_mem = integrator->cvode_mem;  long int lenrw = -1, leniw = -1, nst = -1, nfe = -1, nsetups = -1, nni = -1, ncfn = -1, netf = -1;  long int lenrwLS = -1, leniwLS = -1, nje = -1, nfeLS = -1,npe = -1,nps = -1,ncfl = -1,nli = -1;  int flag;  flag = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);  check_flag(&flag, "CVodeGetWorkSpace", 1);  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_flag(&flag, "CVodeGetNumRhsEvals", 1);  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_flag(&flag, "CVodeGetNumErrTestFails", 1);  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);  printf("/n Final integrator statistics for this run:/n");  printf(" (MM: %s; IM: %s; LS: %s; max-step: %0.4le)/n",    multistepMethodToString(      simulationGetMultistepMethod(integrator->simulation)),    iterationMethodToString(      simulationGetIterationMethod(integrator->simulation)),    linearSolverToString(simulationGetLinearSolver(integrator->simulation)),    simulationGetBvarMaxStep(integrator->simulation));  printf(" CVode real workspace length              = %4ld /n", lenrw);  printf(" CVode integer workspace length           = %4ld /n", leniw);  printf(" Number of steps                          = %4ld /n",  nst);  printf(" Number of f-s                            = %4ld /n",  nfe);  printf(" Number of setups                         = %4ld /n",  nsetups);  printf(" Number of nonlinear iterations           = %4ld /n",  nni);  printf(" Number of nonlinear convergence failures = %4ld /n",  ncfn);  printf(" Number of error test failures            = %4ld /n/n",netf);  if (simulationGetIterationMethod(integrator->simulation) == NEWTON)  {    enum LinearSolver solver =      simulationGetLinearSolver(integrator->simulation);    switch(solver)    {      case DENSE:      {        //flag = CVDenseGetNumJacEvals(cvode_mem, &nje);        //check_flag(&flag, "CVDenseGetNumJacEvals", 1);        //flag = CVDenseGetNumRhsEvals(cvode_mem, &nfeLS);        //check_flag(&flag, "CVDenseGetNumRhsEvals", 1);        //flag = CVDenseGetWorkSpace(cvode_mem, &lenrwLS, &leniwLS);        //check_flag(&flag, "CVDenseGetWorkSpace", 1);      } break;      case BAND:      {        //flag = CVBandGetNumJacEvals(cvode_mem, &nje);        //check_flag(&flag, "CVBandGetNumJacEvals", 1);        //flag = CVBandGetNumRhsEvals(cvode_mem, &nfeLS);        //check_flag(&flag, "CVBandGetNumRhsEvals", 1);        //flag = CVBandGetWorkSpace(cvode_mem, &lenrwLS, &leniwLS);        //check_flag(&flag, "CVBandGetWorkSpace", 1);      } break;      case DIAG:      {        nje = nsetups;        flag = CVDiagGetNumRhsEvals(cvode_mem, &nfeLS);        check_flag(&flag, "CVDiagGetNumRhsEvals", 1);        flag = CVDiagGetWorkSpace(cvode_mem, &lenrwLS, &leniwLS);        check_flag(&flag, "CVDiagGetWorkSpace", 1);      } break;      case SPGMR:      case SPBCG:      case SPTFQMR:      {        nje = nsetups;        flag = CVSpilsGetWorkSpace(cvode_mem,&lenrwLS,&leniwLS);        check_flag(&flag, "CVSpilsGetWorkSpace", 1);        flag = CVSpilsGetNumRhsEvals(cvode_mem, &nfeLS);        check_flag(&flag, "CVSpilsGetNumRhsEvals", 1);                flag = CVSpilsGetNumLinIters(cvode_mem, &nli);        check_flag(&flag, "CVSpilsGetNumLinIters", 1);        flag = CVSpilsGetNumPrecEvals(cvode_mem, &npe);        check_flag(&flag, "CVSpilsGetNumPrecEvals", 1);        flag = CVSpilsGetNumPrecSolves(cvode_mem, &nps);        check_flag(&flag, "CVSpilsGetNumPrecSolves", 1);        flag = CVSpilsGetNumConvFails(cvode_mem, &ncfl);        check_flag(&flag, "CVSpilsGetNumConvFails", 1);      } break;      default:      {        nje = -1;        nfeLS = -1;        lenrwLS = -1;        leniwLS = -1;      }//.........这里部分代码省略.........

static void PrintFinalStats(void *cvode_mem, int miter, realtype ero){  long int lenrw, leniw, nst, nfe, nsetups, nni, ncfn, netf;  long int lenrwL, leniwL, nje, nfeL;  int flag;  flag = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);  check_flag(&flag, "CVodeGetWorkSpace", 1);  flag = CVodeGetNumSteps(cvode_mem, &nst);  check_flag(&flag, "CVodeGetNumSteps", 1);  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);  check_flag(&flag, "CVodeGetNumRhsEvals", 1);  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);  check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);  check_flag(&flag, "CVodeGetNumErrTestFails", 1);  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);  check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);  check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);  printf("/n Final statistics for this run:/n/n");  printf(" CVode real workspace length              = %4ld /n", lenrw);  printf(" CVode integer workspace length           = %4ld /n", leniw);  printf(" Number of steps                          = %4ld /n",  nst);  printf(" Number of f-s                            = %4ld /n",  nfe);  printf(" Number of setups                         = %4ld /n",  nsetups);  printf(" Number of nonlinear iterations           = %4ld /n",  nni);  printf(" Number of nonlinear convergence failures = %4ld /n",  ncfn);  printf(" Number of error test failures            = %4ld /n/n",netf);    if (miter != FUNC) {    switch(miter) {    case DENSE_USER :    case DENSE_DQ   :      flag = CVDenseGetNumJacEvals(cvode_mem, &nje);      check_flag(&flag, "CVDenseGetNumJacEvals", 1);      flag = CVDenseGetNumRhsEvals(cvode_mem, &nfeL);      check_flag(&flag, "CVDenseGetNumRhsEvals", 1);      flag = CVDenseGetWorkSpace(cvode_mem, &lenrwL, &leniwL);      check_flag(&flag, "CVDenseGetWorkSpace", 1);      break;    case BAND_USER  :    case BAND_DQ    :      flag = CVBandGetNumJacEvals(cvode_mem, &nje);      check_flag(&flag, "CVBandGetNumJacEvals", 1);      flag = CVBandGetNumRhsEvals(cvode_mem, &nfeL);      check_flag(&flag, "CVBandGetNumRhsEvals", 1);      flag = CVBandGetWorkSpace(cvode_mem, &lenrwL, &leniwL);      check_flag(&flag, "CVBandGetWorkSpace", 1);      break;      case DIAG       :      nje = nsetups;      flag = CVDiagGetNumRhsEvals(cvode_mem, &nfeL);      check_flag(&flag, "CVDiagGetNumRhsEvals", 1);      flag = CVDiagGetWorkSpace(cvode_mem, &lenrwL, &leniwL);      check_flag(&flag, "CVDiagGetWorkSpace", 1);      break;    }    printf(" Linear solver real workspace length      = %4ld /n", lenrwL);    printf(" Linear solver integer workspace length   = %4ld /n", leniwL);    printf(" Number of Jacobian evaluations           = %4ld  /n", nje);    printf(" Number of f-s evaluations                = %4ld /n/n", nfeL);  }  #if defined(SUNDIALS_EXTENDED_PRECISION)  printf(" Error overrun = %.3Lf /n", ero);#else  printf(" Error overrun = %.3f /n", ero);#endif}

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

void        handleError( void *cvode_mem, N_Vector y0, int flag, mxArray *plhs[], int nrhs, int sensitivity, N_Vector *yS0, realtype *sensitivities, struct mData *data ) {        #ifdef DEBUG        long int temp;        realtype tempreal;    #endif    	#ifdef DEBUG        printf( "<<< DEBUG OUTPUT >>>/n" );        printf( "PARAMETERS: /n" );        for ( temp = 0; temp < N_PARAMS; temp++ ) {            printf( "P(%d) = %f/n", temp, data->p[ temp ] );        }        CVodeGetNumSteps(cvode_mem, &temp);        printf( "Number of steps taken by the solver: %d/n", temp );        CVodeGetNumErrTestFails(cvode_mem, &temp);        printf( "Number of local error test failures: %d/n", temp );        CVodeGetLastStep(cvode_mem, &tempreal);        printf( "Last stepsize: %f/n", tempreal );        CVodeGetCurrentStep(cvode_mem, &tempreal);        printf( "Last step: %f/n", tempreal );        	#endif                                    if ( sensitivity == 1 ) {        if ( nrhs < 7 ) free( sensitivities );		free( data->p );        N_VDestroyVectorArray_Serial( yS0, N_STATES + N_PARAMS );		if ( plhs[2] != NULL ) mxDestroyArray(plhs[2]);    }        if ( plhs[1] != NULL ) mxDestroyArray(plhs[1]);    if ( plhs[0] != NULL ) mxDestroyArray(plhs[0]);    N_VDestroy_Serial( y0 );    /*printf( "Freeing..." );*/    /*CVodeFree( &cvode_mem );*/    /*printf( "Success!/n" );*/    	switch( flag ) {        case CV_MEM_NULL:            printf( "ERROR: No memory was allocated for cvode_mem/n" );            break;        case CV_NO_MALLOC:          printf( "ERROR: Forgot or failed CVodeInit/n" );          break;        case CV_ILL_INPUT:          printf( "ERROR: Input for CVode was illegal/n" );          break;        case CV_TOO_CLOSE:          printf( "ERROR: Initial time too close to final time/n" );          break;        case CV_TOO_MUCH_WORK:          printf( "ERROR: Solver took maximum number of internal steps, but hasn't reached t_out/n" );          break;        case CV_TOO_MUCH_ACC:          printf( "ERROR: Could not attain desired accuracy/n" );          break;        case CV_ERR_FAILURE:          printf( "ERROR: Error tests failed too many times/n" );          break;                              case CV_CONV_FAILURE:          printf( "ERROR: Convergence failure in solving the linear system/n" );          break;        case CV_LINIT_FAIL:          printf( "ERROR: Linear solver failed to initialize/n" );          break;        case CV_LSETUP_FAIL:          printf( "ERROR: Linear solver setup failed/n" );          break;        case CV_RHSFUNC_FAIL:          printf( "ERROR: Right hand side failed in an unrecoverable manner/n" );          break;        case CV_REPTD_RHSFUNC_ERR:          printf( "ERROR: Convergence test failures occured too many times in RHS/n" );          break;        case CV_UNREC_RHSFUNC_ERR:          printf( "ERROR: Unrecoverable error in the RHS/n" );          break;        case CV_RTFUNC_FAIL:          printf( "ERROR: Rootfinding function failed!/n" );          break;        default:          printf( "ERROR: I have no idea what's going on :(/n" );          break;    }        mexErrMsgTxt( "Aborting" );}