C++ check_flag函数代码示例

int main()
{
    realtype abstol=ATOL, reltol=RTOL, t, tout;
    N_Vector c;
    WebData wdata;
    void *cvode_mem;
    booleantype firstrun;
    int jpre, gstype, flag;
    int ns, mxns, iout;

    c = NULL;
    wdata = NULL;
    cvode_mem = NULL;

    /* Initializations */
    c = N_VNew_Serial(NEQ);
    if(check_flag((void *)c, "N_VNew_Serial", 0)) return(1);
    wdata = AllocUserData();
    if(check_flag((void *)wdata, "AllocUserData", 2)) return(1);
    InitUserData(wdata);
    ns = wdata->ns;
    mxns = wdata->mxns;

    /* Print problem description */
    PrintIntro();

    /* Loop over jpre and gstype (four cases) */
    for (jpre = PREC_LEFT; jpre <= PREC_RIGHT; jpre++) {
        for (gstype = MODIFIED_GS; gstype <= CLASSICAL_GS; gstype++) {

            /* Initialize c and print heading */
            CInit(c, wdata);
            PrintHeader(jpre, gstype);

            /* Call CVodeInit or CVodeReInit, then CVSpgmr to set up problem */

            firstrun = (jpre == PREC_LEFT) && (gstype == MODIFIED_GS);
            if (firstrun) {
                cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
                if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

                wdata->cvode_mem = cvode_mem;

                flag = CVodeSetUserData(cvode_mem, wdata);
                if(check_flag(&flag, "CVodeSetUserData", 1)) return(1);

                flag = CVodeInit(cvode_mem, f, T0, c);
                if(check_flag(&flag, "CVodeInit", 1)) return(1);

                flag = CVodeSStolerances(cvode_mem, reltol, abstol);
                if (check_flag(&flag, "CVodeSStolerances", 1)) return(1);

                flag = CVSpgmr(cvode_mem, jpre, MAXL);
                if(check_flag(&flag, "CVSpgmr", 1)) return(1);

                flag = CVSpilsSetGSType(cvode_mem, gstype);
                if(check_flag(&flag, "CVSpilsSetGSType", 1)) return(1);

                flag = CVSpilsSetEpsLin(cvode_mem, DELT);
                if(check_flag(&flag, "CVSpilsSetEpsLin", 1)) return(1);

                flag = CVSpilsSetPreconditioner(cvode_mem, Precond, PSolve);
                if(check_flag(&flag, "CVSpilsSetPreconditioner", 1)) return(1);

            } else {

                flag = CVodeReInit(cvode_mem, T0, c);
                if(check_flag(&flag, "CVodeReInit", 1)) return(1);

                flag = CVSpilsSetPrecType(cvode_mem, jpre);
                check_flag(&flag, "CVSpilsSetPrecType", 1);
                flag = CVSpilsSetGSType(cvode_mem, gstype);
                if(check_flag(&flag, "CVSpilsSetGSType", 1)) return(1);

            }

            /* Print initial values */
            if (firstrun) PrintAllSpecies(c, ns, mxns, T0);

            /* Loop over output points, call CVode, print sample solution values. */
            tout = T1;
            for (iout = 1; iout <= NOUT; iout++) {
                flag = CVode(cvode_mem, tout, c, &t, CV_NORMAL);
                PrintOutput(cvode_mem, t);
                if (firstrun && (iout % 3 == 0)) PrintAllSpecies(c, ns, mxns, t);
                if(check_flag(&flag, "CVode", 1)) break;
                if (tout > RCONST(0.9)) tout += DTOUT;
                else tout *= TOUT_MULT;
            }

            /* Print final statistics, and loop for next case */
            PrintFinalStats(cvode_mem);

        }
    }

    /* Free all memory */
    CVodeFree(&cvode_mem);
    N_VDestroy_Serial(c);
    FreeUserData(wdata);
//.........这里部分代码省略.........

int main()
{
  realtype abstol, reltol, t, tout;
  N_Vector u;
  UserData data;
  void *cvode_mem;
  int iout, flag;

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

  /* Allocate memory, and set problem data, initial values, tolerances */ 
  u = N_VNew_Serial(NEQ);
  if(check_flag((void *)u, "N_VNew_Serial", 0)) return(1);
  data = AllocUserData();
  if(check_flag((void *)data, "AllocUserData", 2)) return(1);
  InitUserData(data);
  SetInitialProfiles(u, data->dx, data->dy);
  abstol=ATOL; 
  reltol=RTOL;

  /* Call CvodeCreate to create the solver memory 

     CV_BDF     specifies the Backward Differentiation Formula
     CV_NEWTON  specifies a Newton iteration

     A pointer to the integrator memory is returned and stored in cvode_mem. */
  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
  if(check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

  /* Set the pointer to user-defined data */
  flag = CVodeSetFdata(cvode_mem, data);
  if(check_flag(&flag, "CVodeSetFdata", 1)) return(1);

  /* Call CVodeMalloc to initialize the integrator memory: 

     f       is the user's right hand side function in u'=f(t,u)
     T0      is the initial time
     u       is the initial dependent variable vector
     CV_SS   specifies scalar relative and absolute tolerances
     reltol  is the relative tolerance
     &abstol is a pointer to the scalar absolute tolerance      */
  flag = CVodeMalloc(cvode_mem, f, T0, u, CV_SS, reltol, &abstol);
  if(check_flag(&flag, "CVodeMalloc", 1)) return(1);

  /* Call CVSpgmr to specify the linear solver CVSPGMR 
     with left preconditioning and the maximum Krylov dimension maxl */
  flag = CVSpgmr(cvode_mem, PREC_LEFT, 0);
  if(check_flag(&flag, "CVSpgmr", 1)) return(1);

  /* Set modified Gram-Schmidt orthogonalization, preconditioner 
     setup and solve routines Precond and PSolve, and the pointer 
     to the user-defined block data */
  flag = CVSpgmrSetGSType(cvode_mem, MODIFIED_GS);
  if(check_flag(&flag, "CVSpgmrSetGSType", 1)) return(1);

  flag = CVSpgmrSetPreconditioner(cvode_mem, Precond, PSolve, data);
  if(check_flag(&flag, "CVSpgmrSetPreconditioner", 1)) return(1);

  /* In loop over output points, call CVode, print results, test for error */
  printf(" \n2-species diurnal advection-diffusion problem\n\n");
  for (iout=1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {
    flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);
    PrintOutput(cvode_mem, u, t);
    if(check_flag(&flag, "CVode", 1)) break;
  }

  PrintFinalStats(cvode_mem);

  /* Free memory */
  N_VDestroy_Serial(u);
  FreeUserData(data);
  CVodeFree(cvode_mem);

  return(0);
}

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

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

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

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

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

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

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

  /* Initialize and allocate memory for KINSOL */

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

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

  /* Set optional inputs */

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

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

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

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

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

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

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

  /* Indicate exact Newton */

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

  /* Initial guess */

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

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

  /* Call KINSol to solve problem */

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

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

  /* Print final statistics and free memory */  

  PrintFinalStats(kmem);

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

int do_integrate(float t_start, float t_stop, int n_points)
{
  realtype reltol, t; 
  N_Vector y, abstol;
  void *cvode_mem;
  int flag, flagr, iout;


  y = abstol = NULL;
  cvode_mem = NULL;

  /* Create serial vector of length NEQ for I.C. and abstol */
  y = N_VNew_Serial(NEQ);
  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
  abstol = N_VNew_Serial(NEQ); 
  if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(1);

  
  //Setup the initial state values:
  setup_initial_states(y);
  setup_tolerances(abstol);

  /* Set the scalar relative tolerance */
  reltol = RTOL;


  /* Call CVodeCreate to create the solver memory and specify the 
   * Backward Differentiation Formula and the use of a Newton iteration */
  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
  if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
  
  /* Call CVodeInit to initialize the integrator memory and specify the
   * user's right hand side function in y'=f(t,y), the inital time T0, and
   * the initial dependent variable vector y. */
  flag = CVodeInit(cvode_mem, f, t_start, y);
  if (check_flag(&flag, "CVodeInit", 1)) return(1);

  /* Call CVodeSVtolerances to specify the scalar relative tolerance
   * and vector absolute tolerances */
  flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
  if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);


  /* Call CVDense to specify the CVDENSE dense linear solver */
  flag = CVDense(cvode_mem, NEQ);
  if (check_flag(&flag, "CVDense", 1)) return(1);

  /* Set the Jacobian routine to Jac (user-supplied) */
  flag = CVDlsSetDenseJacFn(cvode_mem, Jac);
  if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);

  /* In loop, call CVode, print results, and test for error.*/
  printf(" \n3-species kinetics problem\n\n");

  iout = 0;  

  int i=0;
  for(i=1;i< n_points; i++)
  {
    float t_next = t_start + (t_stop-t_start)/n_points * i;
    printf("Advancing to: %f", t_next);

    flag = CVode(cvode_mem, t_next, y, &t, CV_NORMAL);
    loop_function(t,y);


    PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));

    //printf("MH: %d %f",iout,t_next);


    if (check_flag(&flag, "CVode", 1)) break;
    assert(flag==CV_SUCCESS);
  }


  /* Print some final statistics */
  PrintFinalStats(cvode_mem);

  /* Free y and abstol vectors */
  N_VDestroy_Serial(y);
  N_VDestroy_Serial(abstol);

  /* Free integrator memory */
  CVodeFree(&cvode_mem);

  return(0);
}

VkResult setup(VkDevice device, VkCommandPool command_pool, SwapChain& swap_chain, uint32_t *width, uint32_t *height)
{
    VkCommandBuffer setup_command_buffer = VK_NULL_HANDLE;  // command buffer used for setup
    VkCommandBufferAllocateInfo command_buffer_allocate_info;
    command_buffer_allocate_info.sType = VK_STRUCTURE_TYPE_COMMAND_BUFFER_ALLOCATE_INFO;
    command_buffer_allocate_info.pNext = nullptr;
    command_buffer_allocate_info.commandPool = command_pool;
    command_buffer_allocate_info.level = VK_COMMAND_BUFFER_LEVEL_PRIMARY;
    command_buffer_allocate_info.commandBufferCount = 1;
    vkAllocateCommandBuffers(device, &command_buffer_allocate_info, &setup_command_buffer);

    VkCommandBufferBeginInfo command_buffer_begin_info = {};
    command_buffer_begin_info.sType = VK_STRUCTURE_TYPE_COMMAND_BUFFER_BEGIN_INFO;
    command_buffer_begin_info.pNext = nullptr;
    command_buffer_begin_info.flags = 0;
    command_buffer_begin_info.pInheritanceInfo = nullptr;

    vkBeginCommandBuffer(setup_command_buffer, &command_buffer_begin_info);

    // SWAP CHAIN
    uint32_t present_mode_count = 0;
    vkGetPhysicalDeviceSurfacePresentModesKHR(swap_chain.gpu, swap_chain.surface, &present_mode_count, nullptr);
    std::vector<VkPresentModeKHR> present_modes(present_mode_count);
    vkGetPhysicalDeviceSurfacePresentModesKHR(swap_chain.gpu, swap_chain.surface, &present_mode_count, present_modes.data());

    // Try to use mailbox mode
    // Low latency and non-tearing
    VkPresentModeKHR swapchainPresentMode = VK_PRESENT_MODE_FIFO_KHR;
    for (uint32_t i = 0; i < present_mode_count; i++)
    {
        if (present_modes[i] == VK_PRESENT_MODE_MAILBOX_KHR)
        {
            swapchainPresentMode = VK_PRESENT_MODE_MAILBOX_KHR;
            break;
        }
        if ((swapchainPresentMode != VK_PRESENT_MODE_MAILBOX_KHR) && (present_modes[i] == VK_PRESENT_MODE_IMMEDIATE_KHR))
        {
            swapchainPresentMode = VK_PRESENT_MODE_IMMEDIATE_KHR;
        }
    }

    VkSurfaceCapabilitiesKHR surface_cap;
    vkGetPhysicalDeviceSurfaceCapabilitiesKHR(swap_chain.gpu, swap_chain.surface, &surface_cap);

    VkExtent2D swapchainExtent = {};
    // width and height are either both -1, or both not -1.
    if (surface_cap.currentExtent.width == -1)
    {
        // If the surface size is undefined, the size is set to
        // the size of the images requested.
        swapchainExtent.width = *width;
        swapchainExtent.height = *height;
    }
    else
    {
        // If the surface size is defined, the swap chain size must match
        swapchainExtent = surface_cap.currentExtent;
        *width = surface_cap.currentExtent.width;
        *height = surface_cap.currentExtent.height;
    }

    // Determine the number of images
    uint32_t desired_number_of_swapchain_images = surface_cap.minImageCount + 1;
    if ((surface_cap.maxImageCount > 0) && (desired_number_of_swapchain_images > surface_cap.maxImageCount))
    {
        desired_number_of_swapchain_images = surface_cap.maxImageCount;
    }

    VkSurfaceTransformFlagsKHR pre_transform;
    if (check_flag(surface_cap.supportedTransforms, VK_SURFACE_TRANSFORM_IDENTITY_BIT_KHR))
    {
        pre_transform = VK_SURFACE_TRANSFORM_IDENTITY_BIT_KHR;
    }
    else
    {
        pre_transform = surface_cap.currentTransform;
    }
	/*
    VkSwapchainCreateInfoKHR swapchain_creation_info = {};
    swapchain_creation_info.sType = VK_STRUCTURE_TYPE_SWAPCHAIN_CREATE_INFO_KHR;
    swapchain_creation_info.pNext = NULL;
    swapchain_creation_info.surface = surface;
    swapchain_creation_info.minImageCount = desired_number_of_swapchain_images;
    swapchain_creation_info.imageFormat = colorFormat;
    swapchain_creation_info.imageColorSpace = colorSpace;
    swapchain_creation_info.imageExtent = { swapchainExtent.width, swapchainExtent.height };
    swapchain_creation_info.imageUsage = VK_IMAGE_USAGE_COLOR_ATTACHMENT_BIT;
    swapchain_creation_info.preTransform = (VkSurfaceTransformFlagBitsKHR)preTransform;
    swapchain_creation_info.imageArrayLayers = 1;
    swapchain_creation_info.queueFamilyIndexCount = VK_SHARING_MODE_EXCLUSIVE;
    swapchain_creation_info.queueFamilyIndexCount = 0;
    swapchain_creation_info.pQueueFamilyIndices = NULL;
    swapchain_creation_info.presentMode = swapchainPresentMode;
    swapchain_creation_info.oldSwapchain = VK_NULL_HANDLE;
    swapchain_creation_info.clipped = true;
    swapchain_creation_info.compositeAlpha = VK_COMPOSITE_ALPHA_OPAQUE_BIT_KHR;

    vkCreateSwapchainKHR(device, &swapchain_creation_info, nullptr, &swap_chain.swap_chain);


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

int main(int argc, char *argv[])
{
  MPI_Comm comm;
  void *mem;
  UserData data;
  long int SystemSize, local_N, mudq, mldq, mukeep, mlkeep;
  realtype rtol, atol, t0, tout, tret;
  N_Vector uv, uvp, resid, id;
  int thispe, npes, maxl, iout, retval;

  uv = uvp = resid = id = NULL;
  data = NULL;
  mem = NULL;

  /* Set communicator, and get processor number and total number of PE's. */
  MPI_Init(&argc, &argv);
  comm = MPI_COMM_WORLD;
  MPI_Comm_rank(comm, &thispe);
  MPI_Comm_size(comm, &npes);

  if (npes != NPEX*NPEY) {
    if (thispe == 0)
      fprintf(stderr, 
              "\nMPI_ERROR(0): npes = %d not equal to NPEX*NPEY = %d\n", 
              npes, NPEX*NPEY);
    MPI_Finalize();
    return(1); 
  }
  
  /* Set local length (local_N) and global length (SystemSize). */
  local_N = MXSUB*MYSUB*NUM_SPECIES;
  SystemSize = NEQ;

  /* Set up user data block data. */
  data = (UserData) malloc(sizeof *data);

  InitUserData(data, thispe, npes, comm);
  
  /* Create needed vectors, and load initial values.
     The vector resid is used temporarily only.        */
  
  uv  = N_VNew_Parallel(comm, local_N, SystemSize);
  if(check_flag((void *)uv, "N_VNew_Parallel", 0, thispe)) MPI_Abort(comm, 1);

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

  resid = N_VNew_Parallel(comm, local_N, SystemSize);
  if(check_flag((void *)resid, "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(uv, uvp, id, resid, data);

  /* Set remaining inputs to IDAS. */
  t0 = ZERO;
  rtol = RTOL; 
  atol = ATOL;
  
  /* Call IDACreate and IDAInit to initialize solution */

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

  retval = IDASetUserData(mem, data);
  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, res, t0, uv, uvp);
  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 IDAS 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);
  
  retval = IDAGetConsistentIC(mem, uv, uvp);
  if(check_flag(&retval, "IDAGetConsistentIC", 1, thispe)) MPI_Abort(comm, 1);

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

void Ida::initialize()
{
  _properties = dynamic_cast<ISystemProperties*>(_system);
  _continuous_system = dynamic_cast<IContinuous*>(_system);
  _event_system = dynamic_cast<IEvent*>(_system);
  _mixed_system = dynamic_cast<IMixedSystem*>(_system);
  _time_system = dynamic_cast<ITime*>(_system);
  IGlobalSettings* global_settings = dynamic_cast<ISolverSettings*>(_idasettings)->getGlobalSettings();
  // Kennzeichnung, dass initialize()() (vor der Integration) aufgerufen wurde
  _idid = 5000;
  _tLastEvent = 0.0;
  _event_n = 0;
  SolverDefaultImplementation::initialize();

  _dimStates = _continuous_system->getDimContinuousStates();
  _dimZeroFunc = _event_system->getDimZeroFunc()+_event_system->getDimClock();
  _dimAE = _continuous_system->getDimAE();
   if(_dimAE>0)
		_dimSys=_dimAE+ _dimStates;
	else
		_dimSys=_dimStates;
  if (_dimStates <= 0)

  {
    _idid = -1;
    throw std::invalid_argument("Ida::initialize()");
  }
  else
  {
    // Allocate state vectors, stages and temporary arrays

   /*if (_z)
      delete[] _z;
    if (_zInit)
      delete[] _zInit;
    if (_zWrite)
      delete[] _zWrite;*/
    if (_y)
      delete[] _y;
    if (_yInit)
      delete[] _yInit;
    if (_yWrite)
      delete[] _yWrite;
    if (_ypWrite)
      delete[] _ypWrite;
    if (_yp)
      delete[] _yp;
    if (_dae_res)
      delete[] _dae_res;
    if (_zeroSign)
      delete[] _zeroSign;
    if (_absTol)
      delete[] _absTol;
    if(_delta)
      delete [] _delta;
    if(_deltaInv)
      delete [] _deltaInv;
    if(_ysave)
      delete [] _ysave;


	_y = new double[_dimSys];
	_yp = new double[_dimSys];
    _yInit = new double[_dimSys];
    _yWrite = new double[_dimSys];
	_ypWrite = new double[_dimSys];
	_dae_res = new double[_dimSys];
	/*
	_z = new double[_dimSys];
    _zInit = new double[_dimSys];
    _zWrite = new double[_dimSys];
	*/

    _zeroSign = new int[_dimZeroFunc];
    _absTol = new double[_dimSys];
    _delta =new double[_dimSys];
    _deltaInv =new double[_dimSys];
    _ysave =new double[_dimSys];

    memset(_y, 0, _dimSys * sizeof(double));
	memset(_yp, 0, _dimSys * sizeof(double));
    memset(_yInit, 0, _dimSys * sizeof(double));
    memset(_ysave, 0, _dimSys * sizeof(double));
	 std::fill_n(_absTol, _dimSys, 1.0);
    // Counter initialisieren
    _outStps = 0;

    if (_idasettings->getDenseOutput())
    {
      // Ausgabeschrittweite
      _hOut = global_settings->gethOutput();

    }

    // Allocate memory for the solver
    _idaMem = IDACreate();
    if (check_flag((void*) _idaMem, "IDACreate", 0))
    {
      _idid = -5;
      throw std::invalid_argument(/*_idid,_tCurrent,*/"Ida::initialize()");
//.........这里部分代码省略.........

void	attrib_u_char(char **str, t_var **var, int count[3])
{
    int nbr;
    int s_nbr;
    int flag;
    char c;
    int k;
    int check;
    int neg;

    k = 0;
    c = 't';
    neg = 0;
    nbr = ft_nbrlen(var[count[2]]->u_carac);
    flag = check_flag(str, count, &nbr, &c);
    check = ft_atoi_ultra(str[count[0]]);
    if (str[count[0]][0] == '.')
        nbr--;
    str[count[0]] = (char *)malloc(sizeof(char) * (nbr * 5));
    str[count[0]] = ft_ntoa_base_un(var[count[2]]->u_carac, "0123456789");
    s_nbr = ft_strlen(str[count[0]]);
    if (flag == 1000)
    {
        neg = 1;
        str[count[0]] = ft_strjoin("+", str[count[0]]);
    }
    else if (flag == 2000)
        str[count[0]] = ft_strjoin(" ", str[count[0]]);
    else if (flag == 4000)
    {
        str[count[0]] = ft_strjoin("+", str[count[0]]);
        while (k < check - s_nbr)
        {
            str[count[0]] = ft_strjoin(str[count[0]], " ");
            k++;
        }
    }
    else if (flag == 3000)
    {
        if (var[count[2]]->stars < 0)
        {
            while (k > var[count[2]]->stars + nbr - 1)
            {
                str[count[0]] = ft_strjoin(str[count[0]], " ");
                k--;
            }
        }
        else
        {
            while (k < var[count[2]]->stars - nbr + 1)
            {
                str[count[0]] = ft_strjoin(" ", str[count[0]]);
                k++;
            }
        }
    }
    else if (flag != 1 && flag != 1000 && flag != 2000 &&
             flag != 3000 && flag != 5000)
    {
        if (flag < -1 && c != '0')
        {
            while (k > flag + s_nbr + neg)
            {
                str[count[0]] = ft_strjoin(str[count[0]], " ");
                k--;
            }
        }
        else if (flag > 0 && (c == '0' || c == '.'))
        {
            while (k < flag - s_nbr - neg)
            {
                str[count[0]] = ft_strjoin("0", str[count[0]]);
                k++;
            }
        }
        else if (flag < -1 && c == '0')
        {
            while (k > flag)
            {
                str[count[0]] = ft_strjoin("0", str[count[0]]);
                k--;
            }
        }
        else
        {
            while (k < flag - s_nbr - neg)
            {
                str[count[0]] = ft_strjoin(" ", str[count[0]]);
                k++;
            }
        }
    }
    count[2]++;
}

void	attrib_x_maj_short(char **str, t_var **var, int count[3])
{
    int nbr;
    int flag;
    int k;
    int l;
    char c;

    k = 0;
    c = 't';
    nbr = ft_nbrlen(var[count[2]]->u_short);
    l = ft_atoi_spec_o(str[count[0]]);
    flag = check_flag(str, count, &nbr, &c);
    str[count[0]] = ft_ntoa_base_un(var[count[2]]->u_short, "0123456789ABCDEF");
    nbr = ft_strlen(str[count[0]]);
    if (flag == 1000)
        str[count[0]] = ft_strjoin("+", str[count[0]]);
    else if (flag == 2000)
        str[count[0]] = ft_strjoin(" ", str[count[0]]);
    else if (flag == 5000)
    {
        str[count[0]] = ft_strjoin("0x", str[count[0]]);
        flag = l;
        nbr++;
    }
    else if (flag == 3000)
    {
        if (var[count[2]]->stars < 0)
        {
            while (k > var[count[2]]->stars + nbr - 1)
            {
                str[count[0]] = ft_strjoin(str[count[0]], " ");
                k--;
            }
        }
        else
        {
            while (k < var[count[2]]->stars - nbr + 1)
            {
                str[count[0]] = ft_strjoin(" ", str[count[0]]);
                k++;
            }
        }
    }
    if (flag != -1 && flag != 1000 && flag != 2000 &&
            flag != 3000 && flag != 5000)
    {
        if (flag < -1 && c != '0')
        {
            while (k > flag + nbr)
            {
                str[count[0]] = ft_strjoin(str[count[0]], " ");
                k--;
            }
        }
        else if (flag > 0 && (c == '0' || c == '.'))
        {
            while (k < flag - nbr)
            {
                str[count[0]] = ft_strjoin("0", str[count[0]]);
                k++;
            }
        }
        else if (flag < -1 && c == '0')
        {
            while (k > flag)
            {
                str[count[0]] = ft_strjoin("0", str[count[0]]);
                k--;
            }
        }
        else
        {
            while (k < flag - nbr)
            {
                str[count[0]] = ft_strjoin(" ", str[count[0]]);
                k++;
            }
        }
    }
    if (c != '.' && c != '0')
    {
        while (str[count[0]][0] == '0' && str[count[0]][1] != '\0')
            str[count[0]] = ft_strsub(str[count[0]], 1, ft_strlen(str[count[0]]));
    }
    count[2]++;
}

int main(void)
{
  void *mem;
  UserData data;
  N_Vector uu, up, constraints, id, res;
  int ier, iout;
  long int mu, ml, netf, ncfn;
  realtype rtol, atol, t0, t1, tout, tret;
  
  mem = NULL;
  data = NULL;
  uu = up = constraints = id = res = NULL;

  /* Create vectors uu, up, res, constraints, id. */
  uu = N_VNew_Serial(NEQ);
  if(check_flag((void *)uu, "N_VNew_Serial", 0)) return(1);
  up = N_VNew_Serial(NEQ);
  if(check_flag((void *)up, "N_VNew_Serial", 0)) return(1);
  res = N_VNew_Serial(NEQ);
  if(check_flag((void *)res, "N_VNew_Serial", 0)) return(1);
  constraints = N_VNew_Serial(NEQ);
  if(check_flag((void *)constraints, "N_VNew_Serial", 0)) return(1);
  id = N_VNew_Serial(NEQ);
  if(check_flag((void *)id, "N_VNew_Serial", 0)) return(1);

  /* Create and load problem data block. */
  data = (UserData) malloc(sizeof *data);
  if(check_flag((void *)data, "malloc", 2)) return(1);
  data->mm = MGRID;
  data->dx = ONE/(MGRID - ONE);
  data->coeff = ONE/( (data->dx) * (data->dx) );

  /* Initialize uu, up, id. */
  SetInitialProfile(data, uu, up, id, res);

  /* Set constraints to all 1's for nonnegative solution values. */
  N_VConst(ONE, constraints);

  /* Set remaining input parameters. */
  t0   = ZERO;
  t1   = RCONST(0.01);
  rtol = ZERO;
  atol = RCONST(1.0e-3);

  /* Call IDACreate and IDAMalloc to initialize solution */
  mem = IDACreate();
  if(check_flag((void *)mem, "IDACreate", 0)) return(1);

  ier = IDASetUserData(mem, data);
  if(check_flag(&ier, "IDASetUserData", 1)) return(1);

  ier = IDASetId(mem, id);
  if(check_flag(&ier, "IDASetId", 1)) return(1);

  ier = IDASetConstraints(mem, constraints);
  if(check_flag(&ier, "IDASetConstraints", 1)) return(1);
  N_VDestroy_Serial(constraints);

  ier = IDAInit(mem, heatres, t0, uu, up);
  if(check_flag(&ier, "IDAInit", 1)) return(1);

  ier = IDASStolerances(mem, rtol, atol);
  if(check_flag(&ier, "IDASStolerances", 1)) return(1);

  /* Call IDABand to specify the linear solver. */
  mu = MGRID; ml = MGRID;
  ier = IDABand(mem, NEQ, mu, ml);
  if(check_flag(&ier, "IDABand", 1)) return(1);
 
  /* Call IDACalcIC to correct the initial values. */
  
  ier = IDACalcIC(mem, IDA_YA_YDP_INIT, t1);
  if(check_flag(&ier, "IDACalcIC", 1)) return(1);

  /* Print output heading. */
  PrintHeader(rtol, atol);
  
  PrintOutput(mem, t0, uu);


  /* Loop over output times, call IDASolve, and print results. */
  
  for (tout = t1, iout = 1; iout <= NOUT; iout++, tout *= TWO) {
    
    ier = IDASolve(mem, tout, &tret, uu, up, IDA_NORMAL);
    if(check_flag(&ier, "IDASolve", 1)) return(1);

    PrintOutput(mem, tret, uu);
  
  }
  
  /* Print remaining counters and free memory. */
  ier = IDAGetNumErrTestFails(mem, &netf);
  check_flag(&ier, "IDAGetNumErrTestFails", 1);
  ier = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
  check_flag(&ier, "IDAGetNumNonlinSolvConvFails", 1);
  printf("\n netf = %ld,   ncfn = %ld \n", netf, ncfn);

  IDAFree(&mem);
  N_VDestroy_Serial(uu);
//.........这里部分代码省略.........

void	attrib_d_char(char **str, t_var **var, int count[3])
{
    int nbr;
    int flag;
    char c;
    int k;
    int check;
    int check_double;
    int neg;

    k = 0;
    c = 't';
    neg = 0;
    nbr = ft_nbrlen(var[count[2]]->carac);
    flag = check_flag(str, count, &nbr, &c);
    check = ft_atoi_ultra(str[count[0]]);
    if (flag == 3500)
        check_double = ft_atoi_double(str[count[0]]);
    if (str[count[0]][0] == '.' && var[count[2]]->carac < 0)
        nbr--;
    str[count[0]] = ft_ntoa_base(var[count[2]]->carac, "0123456789");
    if (var[count[2]]->carac < 0 && str[count[0]][0] != '0')
        neg = 1;
    if (flag == 1000 && var[count[2]]->carac >= 0)
    {
        neg = 1;
        str[count[0]] = ft_strjoin("+", str[count[0]]);
    }
    else if (flag == 3500)
    {
        if (var[count[2]]->carac < 0)
        {
            str[count[0]] = ft_strsub(str[count[0]], 1, ft_nbrlen(var[count[2]]->carac) + 1);
            while (k < check_double - nbr - 1)
            {
                str[count[0]] = ft_strjoin("0", str[count[0]]);
                k++;
            }
            k = 0;
            str[count[0]] = ft_strjoin("-", str[count[0]]);
            while (k < check - check_double)
            {
                str[count[0]] = ft_strjoin(" ", str[count[0]]);
                k++;
            }
        }
        else
        {
            while (k < check_double - nbr)
            {
                str[count[0]] = ft_strjoin("0", str[count[0]]);
                k++;
            }
            k = 0;
            while (k < check - check_double)
            {
                str[count[0]] = ft_strjoin(" ", str[count[0]]);
                k++;
            }
        }
    }
    else if (flag == 2000 && var[count[2]]->carac >= 0)
        str[count[0]] = ft_strjoin(" ", str[count[0]]);
    else if (flag == 4000)
    {
        if (var[count[2]]->carac >= 0)
            str[count[0]] = ft_strjoin("+", str[count[0]]);
        while (k < check - nbr)
        {
            str[count[0]] = ft_strjoin(str[count[0]], " ");
            k++;
        }
    }
    else if (flag == 3000)
    {
        if (var[count[2]]->stars < 0)
        {
            while (k > var[count[2]]->stars + nbr - 1)
            {
                str[count[0]] = ft_strjoin(str[count[0]], " ");
                k--;
            }
        }
        else
        {
            while (k < var[count[2]]->stars - nbr + 1)
            {
                str[count[0]] = ft_strjoin(" ", str[count[0]]);
                k++;
            }
        }
    }
    else if (flag != 1 && flag != 1000 && flag != 2000 &&
             flag != 3000 && flag != 4000 && flag != 3500 && flag != 5000)
    {
        if (flag < -1 && c != '0')
        {
            while (k > flag + nbr + neg)
            {
                str[count[0]] = ft_strjoin(str[count[0]], " ");
//.........这里部分代码省略.........

static int simulate_nmr_pulse(struct bloch_sim *bs)
{
    N_Vector M = NULL;
    M = N_VNew_Serial(3 * bs->num_cells);
    if (check_flag((void *)M, "N_VNew_Serial", 0)) return(1);

    int i;
    /* Set initial (t=0) magnetization conditions */
    for(i=0; i < bs->num_cells; ++i)
    {
        X(M,i) = 0.0;
        Y(M,i) = 0.0;
        Z(M,i) = bs->cell_frequencies[i] / bs->w_avg;
    }

    realtype reltol = RCONST(1.0e-14);
    realtype abstol = RCONST(1.0e-14);

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

    int flag;
    //TODO: check if flag should be pointer;
    flag = CVodeInit(cvode_mem, bloch_equations, 0.0, M);
    if (check_flag(&flag, "CVodeInit", 1)) return(1);

    flag = CVodeSetUserData(cvode_mem, bs);
    if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);

    flag = CVodeSStolerances(cvode_mem, reltol, abstol);
    if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);

    flag = CVDense(cvode_mem, 3 * bs->num_cells);
    if (check_flag(&flag, "CVDense", 1)) return(1);

    flag = CVDlsSetDenseJacFn(cvode_mem, bloch_jacobian);
    if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);

    ///////////////////////////
    // PI/2 PULSE SIMULATION //
    ///////////////////////////
    bs->rf_on = 1;

    flag = CVodeSetStopTime(cvode_mem, bs->pi2_duration);
    if (check_flag(&flag, "CVodeSetStopTime", 1)) return 1;

    realtype time_reached;
    flag = CVode(cvode_mem, bs->pi2_duration, M, &time_reached, CV_NORMAL);

    if (flag != CV_SUCCESS)
    {
        printf("ERROR: Failed to simulate Pi/2 pulse\n");
    }

    // {{{ PI2 PULSE DEBUG STATEMENTS
    if (DEBUG)
    {
        printf("\n");
        printf("#####################################\n");
        printf("### PI/2 PULSE SIMULATION DETAILS ###\n");
        printf("#####################################\n");
        printf("\n");

        printf("TIME REACHED: %.15e\n", time_reached);
        printf("TIME REACHED - PI2_DURATION: %.15e\n", time_reached - bs->pi2_duration);

        printf("\n");
        printf("MAGNETIZATION AT END OF PI/2 PULSE:\n");
        for (i = 0; i<bs->num_cells; ++i)
        {
            printf("CELL %d: %.4e, %.4e, %.4e\n", i, X(M,i), Y(M,i), Z(M,i));
        }
    }
    // }}}

    ////////////////////
    // FID SIMULATION //
    ////////////////////
    bs->rf_on = 0;

    if (DEBUG)
    {
        printf("##############################\n");
        printf("### FID SIMULATION DETAILS ###\n");
        printf("##############################\n");
    }

    flag = CVodeReInit(cvode_mem, 0.0, M);
    if (check_flag(&flag, "CVodeReInit", 1)) return(1);
    time_reached = 0.0;

    flag = CVodeSetStopTime(cvode_mem, bs->fid_duration);
    if (check_flag(&flag, "CVodeSetStopTime", 1)) return 1;

    flag = CVodeRootInit(cvode_mem, 1, bloch_root);
    if (check_flag(&flag, "CVodeRootInit", 1)) return 1;

    realtype time_desired, M_FID_X;
    while (time_reached < bs->fid_duration)
    {
//.........这里部分代码省略.........

/*
 This routine generates the block-diagonal part of the Jacobian
 corresponding to the interaction rates, multiplies by -gamma, adds
 the identity matrix, and calls denseGETRF to do the LU decomposition of
 each diagonal block. The computation of the diagonal blocks uses
 the preset block and grouping information. One block per group is
 computed. The Jacobian elements are generated by difference
 quotients using calls to the routine fblock.

 This routine can be regarded as a prototype for the general case
 of a block-diagonal preconditioner. The blocks are of size mp, and
 there are ngrp=ngx*ngy blocks computed in the block-grouping scheme.
*/
static int Precond(realtype t, N_Vector c, N_Vector fc,
                   booleantype jok, booleantype *jcurPtr, realtype gamma,
                   void *user_data, N_Vector vtemp1, N_Vector vtemp2,
                   N_Vector vtemp3)
{
    realtype ***P;
    long int **pivot, ier;
    int i, if0, if00, ig, igx, igy, j, jj, jx, jy;
    int *jxr, *jyr, ngrp, ngx, ngy, mxmp, flag;
    long int mp;
    realtype uround, fac, r, r0, save, srur;
    realtype *f1, *fsave, *cdata, *rewtdata;
    WebData wdata;
    void *cvode_mem;
    N_Vector rewt;

    wdata = (WebData) user_data;
    cvode_mem = wdata->cvode_mem;
    cdata = N_VGetArrayPointer_Serial(c);
    rewt = wdata->rewt;
    flag = CVodeGetErrWeights(cvode_mem, rewt);
    if(check_flag(&flag, "CVodeGetErrWeights", 1)) return(1);
    rewtdata = N_VGetArrayPointer_Serial(rewt);

    uround = UNIT_ROUNDOFF;

    P = wdata->P;
    pivot = wdata->pivot;
    jxr = wdata->jxr;
    jyr = wdata->jyr;
    mp = wdata->mp;
    srur = wdata->srur;
    ngrp = wdata->ngrp;
    ngx = wdata->ngx;
    ngy = wdata->ngy;
    mxmp = wdata->mxmp;
    fsave = wdata->fsave;

    /* Make mp calls to fblock to approximate each diagonal block of Jacobian.
       Here, fsave contains the base value of the rate vector and
       r0 is a minimum increment factor for the difference quotient. */

    f1 = N_VGetArrayPointer_Serial(vtemp1);

    fac = N_VWrmsNorm (fc, rewt);
    r0 = RCONST(1000.0)*SUNRabs(gamma)*uround*NEQ*fac;
    if (r0 == ZERO) r0 = ONE;

    for (igy = 0; igy < ngy; igy++) {
        jy = jyr[igy];
        if00 = jy*mxmp;
        for (igx = 0; igx < ngx; igx++) {
            jx = jxr[igx];
            if0 = if00 + jx*mp;
            ig = igx + igy*ngx;
            /* Generate ig-th diagonal block */
            for (j = 0; j < mp; j++) {
                /* Generate the jth column as a difference quotient */
                jj = if0 + j;
                save = cdata[jj];
                r = SUNMAX(srur*SUNRabs(save),r0/rewtdata[jj]);
                cdata[jj] += r;
                fac = -gamma/r;
                fblock (t, cdata, jx, jy, f1, wdata);
                for (i = 0; i < mp; i++) {
                    P[ig][j][i] = (f1[i] - fsave[if0+i])*fac;
                }
                cdata[jj] = save;
            }
        }
    }

    /* Add identity matrix and do LU decompositions on blocks. */

    for (ig = 0; ig < ngrp; ig++) {
        denseAddIdentity(P[ig], mp);
        ier = denseGETRF(P[ig], mp, mp, pivot[ig]);
        if (ier != 0) return(1);
    }

    *jcurPtr = TRUE;
    return(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 abstol, reltol, t, tout;
  N_Vector u;
  UserData data;
  PreconData predata;
  void *cvode_mem;
  int iout, flag, my_pe, npes;
  long int neq, local_N;
  MPI_Comm comm;

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

  /* Set problem size neq */
  neq = NVARS*MX*MY;

  /* Get processor number and total number of pe's */
  MPI_Init(&argc, &argv);
  comm = MPI_COMM_WORLD;
  MPI_Comm_size(comm, &npes);
  MPI_Comm_rank(comm, &my_pe);

  if (npes != NPEX*NPEY) {
    if (my_pe == 0)
      fprintf(stderr, "\nMPI_ERROR(0): npes = %d is not equal to NPEX*NPEY = %d\n\n",
	      npes,NPEX*NPEY);
    MPI_Finalize();
    return(1);
  }

  /* Set local length */
  local_N = NVARS*MXSUB*MYSUB;

  /* Allocate and load user data block; allocate preconditioner block */
  data = (UserData) malloc(sizeof *data);
  if (check_flag((void *)data, "malloc", 2, my_pe)) MPI_Abort(comm, 1);
  InitUserData(my_pe, comm, data);
  predata = AllocPreconData (data);

  /* Allocate u, and set initial values and tolerances */ 
  u = N_VNew_Parallel(comm, local_N, neq);
  if (check_flag((void *)u, "N_VNew", 0, my_pe)) MPI_Abort(comm, 1);
  SetInitialProfiles(u, data);
  abstol = ATOL; reltol = RTOL;

  /* 
     Call CVodeCreate to create the solver memory:
     
     CV_BDF     specifies the Backward Differentiation Formula
     CV_NEWTON  specifies a Newton iteration

     A pointer to the integrator memory is returned and stored in cvode_mem.
  */
  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
  if (check_flag((void *)cvode_mem, "CVodeCreate", 0, my_pe)) MPI_Abort(comm, 1);

  /* Set the pointer to user-defined data */
  flag = CVodeSetFdata(cvode_mem, data);
  if (check_flag(&flag, "CVodeSetFdata", 1, my_pe)) MPI_Abort(comm, 1);

  /* 
     Call CVodeMalloc to initialize the integrator memory: 

     cvode_mem is the pointer to the integrator memory returned by CVodeCreate
     f       is the user's right hand side function in y'=f(t,y)
     T0      is the initial time
     u       is the initial dependent variable vector
     CV_SS   specifies scalar relative and absolute tolerances
     reltol  is the relative tolerance
     &abstol is a pointer to the scalar absolute tolerance
  */
  flag = CVodeMalloc(cvode_mem, f, T0, u, CV_SS, reltol, &abstol);
  if (check_flag(&flag, "CVodeMalloc", 1, my_pe)) MPI_Abort(comm, 1);

  /* Call CVSpgmr to specify the linear solver CVSPGMR 
     with left preconditioning and the maximum Krylov dimension maxl */
  flag = CVSpgmr(cvode_mem, PREC_LEFT, 0);
  if (check_flag(&flag, "CVSpgmr", 1, my_pe)) MPI_Abort(comm, 1);

  /* Set preconditioner setup and solve routines Precond and PSolve, 
     and the pointer to the user-defined block data */
  flag = CVSpilsSetPreconditioner(cvode_mem, Precond, PSolve, predata);
  if (check_flag(&flag, "CVSpilsSetPreconditioner", 1, my_pe)) MPI_Abort(comm, 1);

  if (my_pe == 0)
    printf("\n2-species diurnal advection-diffusion problem\n\n");

  /* In loop over output points, call CVode, print results, test for error */
  for (iout=1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {
    flag = CVode(cvode_mem, tout, u, &t, CV_NORMAL);
    if (check_flag(&flag, "CVode", 1, my_pe)) break;
    PrintOutput(cvode_mem, my_pe, comm, u, t);
  }

  /* Print final statistics */  
  if (my_pe == 0) PrintFinalStats(cvode_mem);

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