void N_VDiv_Serial(N_Vector x, N_Vector y, N_Vector z)
{
  long int i, N;
  realtype *xd, *yd, *zd;

  N  = NV_LENGTH_S(x);
  xd = NV_DATA_S(x);
  yd = NV_DATA_S(y);
  zd = NV_DATA_S(z);

  for (i=0; i < N; i++)
    *zd++ = (*xd++) / (*yd++);
}

static void f(realtype t, N_Vector c, N_Vector cdot, void *f_data)
{
  int i, ic, ici, idxl, idxu, idyl, idyu, iyoff, jx, jy, ns, mxns;
  realtype dcxli, dcxui, dcyli, dcyui, x, y, *cox, *coy, *fsave, dx, dy;
  realtype *cdata, *cdotdata;
  WebData wdata;
  
  wdata = (WebData) f_data;
  cdata = NV_DATA_S(c);
  cdotdata = NV_DATA_S(cdot);
  
  mxns = wdata->mxns;
  ns = wdata->ns;
  fsave = wdata->fsave;
  cox = wdata->cox;
  coy = wdata->coy;
  mxns = wdata->mxns;
  dx = wdata->dx;
  dy = wdata->dy;
  
  for (jy = 0; jy < MY; jy++) {
    y = jy*dy;
    iyoff = mxns*jy;
    idyu = (jy == MY-1) ? -mxns : mxns;
    idyl = (jy == 0) ? -mxns : mxns;
    for (jx = 0; jx < MX; jx++) {
      x = jx*dx;
      ic = iyoff + ns*jx;
      /* Get interaction rates at one point (x,y). */
      WebRates(x, y, t, cdata+ic, fsave+ic, wdata);
      idxu = (jx == MX-1) ? -ns : ns;
      idxl = (jx == 0) ? -ns : ns;
      for (i = 1; i <= ns; i++) {
        ici = ic + i-1;
        /* Do differencing in y. */
        dcyli = cdata[ici] - cdata[ici-idyl];
        dcyui = cdata[ici+idyu] - cdata[ici];
        /* Do differencing in x. */
        dcxli = cdata[ici] - cdata[ici-idxl];
        dcxui = cdata[ici+idxu] - cdata[ici];
        /* Collect terms and load cdot elements. */
        cdotdata[ici] = coy[i-1]*(dcyui - dcyli) + 
                        cox[i-1]*(dcxui - dcxli) +
                        fsave[ici];
      }
    }
  }

  /* Quadrature equation (species 1) */
  cdotdata[NEQ] = doubleIntgr(c,ISPEC,wdata);
}

  int IdasInterface::resB(double t, N_Vector xz, N_Vector xzdot, N_Vector rxz,
                                 N_Vector rxzdot, N_Vector rr, void *user_data) {
    try {
      auto m = to_mem(user_data);
      auto& s = m->self;
      m->arg[0] = NV_DATA_S(rxz);
      m->arg[1] = NV_DATA_S(rxz)+s.nrx_;
      m->arg[2] = m->rp;
      m->arg[3] = NV_DATA_S(xz);
      m->arg[4] = NV_DATA_S(xz)+s.nx_;
      m->arg[5] = m->p;
      m->arg[6] = &t;
      m->res[0] = NV_DATA_S(rr);
      m->res[1] = NV_DATA_S(rr)+s.nrx_;
      s.calc_function(m, "daeB");

      // Subtract state derivative to get residual
      casadi_axpy(s.nrx_, 1., NV_DATA_S(rxzdot), NV_DATA_S(rr));

      return 0;
    } catch(int flag) { // recoverable error
      return flag;
    } catch(exception& e) { // non-recoverable error
      userOut<true, PL_WARN>() << "resB failed: " << e.what() << endl;
      return -1;
    }
  }

realtype N_VDotProd_Serial(N_Vector x, N_Vector y)
{
  long int i, N;
  realtype sum = ZERO, *xd, *yd;

  N  = NV_LENGTH_S(x);
  xd = NV_DATA_S(x);
  yd = NV_DATA_S(y);

  for (i=0; i < N; i++)
    sum += (*xd++) * (*yd++);
  
  return(sum);
}

static void SetInitialProfiles(N_Vector cc, N_Vector cp, N_Vector id,
                               UserData webdata)
{
  long int loc, yloc, is, jx, jy, np;
  realtype xx, yy, xyfactor, fac;
  realtype *ccv, *cpv, *idv;
  
  ccv = NV_DATA_S(cc);
  cpv = NV_DATA_S(cp);
  idv = NV_DATA_S(id);
  np = webdata->np;
  
  /* Loop over grid, load cc values and id values. */
  for (jy = 0; jy < MY; jy++) {
    yy = jy * webdata->dy;
    yloc = NSMX * jy;
    for (jx = 0; jx < MX; jx++) {
      xx = jx * webdata->dx;
      xyfactor = RCONST(16.0)*xx*(ONE-xx)*yy*(ONE-yy);
      xyfactor *= xyfactor;
      loc = yloc + NUM_SPECIES*jx;
      fac = ONE + ALPHA * xx * yy + BETA * sin(FOURPI*xx) * sin(FOURPI*yy);
      
      for (is = 0; is < NUM_SPECIES; is++) {
        if (is < np) {
	    ccv[loc+is] = RCONST(10.0) + (realtype)(is+1) * xyfactor;
          idv[loc+is] = ONE;
        }
        else {
	  ccv[loc+is] = RCONST(1.0e5);
          idv[loc+is] = ZERO;
        }
      }
    }
  }
  
  /* Set c' for the prey by calling the function Fweb. */
  Fweb(ZERO, cc, cp, webdata);
  
  /* Set c' for predators to 0. */
  for (jy = 0; jy < MY; jy++) {
    yloc = NSMX * jy;
    for (jx = 0; jx < MX; jx++) {
      loc = yloc + NUM_SPECIES * jx;
      for (is = np; is < NUM_SPECIES; is++) {
        cpv[loc+is] = ZERO;
      }
    }
  }
}

// Sample residual function
int resrob(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void *rdata)
{
  realtype *yval, *ypval, *rval;

  yval = NV_DATA_S(yy); 
  ypval = NV_DATA_S(yp); 
  rval = NV_DATA_S(rr);
  
  rval[0] = -0.04*yval[0] + 1.0e4*yval[1]*yval[2] - ypval[0];
  rval[1] = 0.04*yval[0] - 1.0e4*yval[1]*yval[2] - 3.0e7*yval[1]*yval[1] - ypval[1];
  rval[2] = yval[0] + yval[1] + yval[2] - 1;

  return(0);
}

int user_fun_wrapper(CDATAFORMAT t, N_Vector y, N_Vector ydot, void *userdata){
  cvode_mem *mem = (cvode_mem*)userdata;
  solver_props *props = mem->props;
  unsigned int modelid = mem->modelid;

  model_flows(t,
	      NV_DATA_S(y) - (modelid*props->statesize), // Model flows indexes y and dydt with modelid
	      NV_DATA_S(ydot) - (modelid*props->statesize), // and this ptr arithmetic adjusts to compensate
	      props, mem->first_iteration, modelid);

  mem->first_iteration = FALSE;

  return CV_SUCCESS;
}

static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up, 
                             N_Vector id, N_Vector res)
{
  realtype xfact, yfact, *udata, *updata, *iddata;
  long int mm, mm1, i, j, offset, loc;
  
  mm = data->mm;
  mm1 = mm - 1;
  
  udata = NV_DATA_S(uu);
  updata = NV_DATA_S(up);
  iddata = NV_DATA_S(id);

  /* Initialize id to 1's. */
  N_VConst(ONE, id);

  /* Initialize uu on all grid points. */ 
  for (j = 0; j < mm; j++) {
    yfact = data->dx * j;
    offset = mm*j;
    for (i = 0;i < mm; i++) {
      xfact = data->dx * i;
      loc = offset + i;
      udata[loc] = RCONST(16.0) * xfact * (ONE - xfact) * yfact * (ONE - yfact);
    }
  }
  
  /* Initialize up vector to 0. */
  N_VConst(ZERO, up);

  /* heatres sets res to negative of ODE RHS values at interior points. */
  heatres(ZERO, uu, up, res, data);
  
  /* Copy -res into up to get correct interior initial up values. */
  N_VScale(-ONE, res, up);

  /* Finally, set values of u, up, and id at boundary points. */
  for (j = 0; j < mm; j++) {
    offset = mm*j;
    for (i = 0;i < mm; i++) {
      loc = offset + i;
      if (j == 0 || j == mm1 || i == 0 || i == mm1 ) {
        udata[loc] = BVAL; updata[loc] = ZERO; iddata[loc] = ZERO; }
    }
  }
  
  return(0);

}

static int func(N_Vector u, N_Vector f, void *user_data)
{
  realtype dx, dy, hdiff, vdiff;
  realtype hdc, vdc;
  realtype uij, udn, uup, ult, urt;
  realtype *udata, *fdata;
  realtype x,y;

  int i, j;

  dx = ONE/(NX+1);  
  dy = ONE/(NY+1);
  hdc = ONE/(dx*dx);
  vdc = ONE/(dy*dy);

  udata = NV_DATA_S(u);
  fdata = NV_DATA_S(f);

  for (j=1; j <= NY; j++) {

    y = j*dy;

    for (i=1; i <= NX; i++) {

      x = i*dx;

      /* Extract u at x_i, y_j and four neighboring points */

      uij = IJth(udata, i, j);
      udn = (j == 1)  ? ZERO : IJth(udata, i, j-1);
      uup = (j == NY) ? ZERO : IJth(udata, i, j+1);
      ult = (i == 1)  ? ZERO : IJth(udata, i-1, j);
      urt = (i == NX) ? ZERO : IJth(udata, i+1, j);

      /* Evaluate diffusion components */

      hdiff = hdc*(ult - TWO*uij + urt);
      vdiff = vdc*(uup - TWO*uij + udn);

      /* Set residual at x_i, y_j */

      IJth(fdata, i, j) = hdiff + vdiff + uij - uij*uij*uij + 2.0;

    }

  }

  return(0);
}

int resrob(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
{
  realtype *yval, *ypval, *rval;

  yval = NV_DATA_S(yy); 
  ypval = NV_DATA_S(yp); 
  rval = NV_DATA_S(rr);

  rval[0]  = RCONST(-0.04)*yval[0] + RCONST(1.0e4)*yval[1]*yval[2];
  rval[1]  = -rval[0] - RCONST(3.0e7)*yval[1]*yval[1] - ypval[1];
  rval[0] -=  ypval[0];
  rval[2]  =  yval[0] + yval[1] + yval[2] - ONE;

  return(0);
}

static int f(realtype t,N_Vector y,N_Vector ydot,void *f_data)
{
  ExecutableModel* em = (ExecutableModel*)f_data;
  realtype* yD = NV_DATA_S(y);
  realtype* ydotD = NV_DATA_S(ydot);
  long int len = NV_LENGTH_S(y);
  long int i;

  for (i=0;i<len;i++) em->states[i]=(double)yD[i];
  /* while integrating we only need to compute the rates */
  em->computeRates(t);
  for (i=0;i<len;i++) ydotD[i] = (realtype)(em->rates[i]);

  return(0);
}

static void PrintAllSpecies(N_Vector c, int ns, int mxns, realtype t)
{
    int i, jx ,jy;
    realtype *cdata;

    cdata = NV_DATA_S(c);
#if defined(SUNDIALS_EXTENDED_PRECISION)
    printf("c values at t = %Lg:\n\n", t);
#elif defined(SUNDIALS_DOUBLE_PRECISION)
    printf("c values at t = %lg:\n\n", t);
#else
    printf("c values at t = %g:\n\n", t);
#endif
    for (i=1; i <= ns; i++) {
        printf("Species %d\n", i);
        for (jy=MY-1; jy >= 0; jy--) {
            for (jx=0; jx < MX; jx++) {
#if defined(SUNDIALS_EXTENDED_PRECISION)
                printf("%-10.6Lg", cdata[(i-1) + jx*ns + jy*mxns]);
#elif defined(SUNDIALS_DOUBLE_PRECISION)
                printf("%-10.6lg", cdata[(i-1) + jx*ns + jy*mxns]);
#else
                printf("%-10.6g", cdata[(i-1) + jx*ns + jy*mxns]);
#endif
            }
            printf("\n");
        }
        printf("\n");
    }
}

/* This routine computes and loads the vector of initial values. */
static void CInit(N_Vector c, WebData wdata)
{
    int jx, jy, ns, mxns, ioff, iyoff, i, ici;
    realtype argx, argy, x, y, dx, dy, x_factor, y_factor, *cdata;

    cdata = NV_DATA_S(c);
    ns = wdata->ns;
    mxns = wdata->mxns;
    dx = wdata->dx;
    dy = wdata->dy;

    x_factor = RCONST(4.0)/SQR(AX);
    y_factor = RCONST(4.0)/SQR(AY);
    for (jy = 0; jy < MY; jy++) {
        y = jy*dy;
        argy = SQR(y_factor*y*(AY-y));
        iyoff = mxns*jy;
        for (jx = 0; jx < MX; jx++) {
            x = jx*dx;
            argx = SQR(x_factor*x*(AX-x));
            ioff = iyoff + ns*jx;
            for (i = 1; i <= ns; i++) {
                ici = ioff + i-1;
                cdata[ici] = RCONST(10.0) + i*argx*argy;
            }
        }
    }
}

static void SetInitialProfiles(N_Vector u, realtype dx, realtype dy)
{
  int jx, jy;
  realtype x, y, cx, cy;
  realtype *udata;

  /* Set pointer to data array in vector u. */

  udata = NV_DATA_S(u);

  /* Load initial profiles of c1 and c2 into u vector */

  for (jy = 0; jy < MY; jy++) {
    y = YMIN + jy*dy;
    cy = SQR(RCONST(0.1)*(y - YMID));
    cy = ONE - cy + RCONST(0.5)*SQR(cy);
    for (jx = 0; jx < MX; jx++) {
      x = XMIN + jx*dx;
      cx = SQR(RCONST(0.1)*(x - XMID));
      cx = ONE - cx + RCONST(0.5)*SQR(cx);
      IJKth(udata,1,jx,jy) = C1_SCALE*cx*cy; 
      IJKth(udata,2,jx,jy) = C2_SCALE*cx*cy;
    }
  }
}

  void IdasInterface::reset(IntegratorMemory* mem, double t, const double* _x,
                            const double* _z, const double* _p) const {
    log("IdasInterface::reset", "begin");
    auto m = to_mem(mem);

    // Reset the base classes
    SundialsInterface::reset(mem, t, _x, _z, _p);

    // Re-initialize
    copy(init_xdot_.begin(), init_xdot_.end(), NV_DATA_S(m->xzdot));
    THROWING(IDAReInit, m->mem, grid_.front(), m->xz, m->xzdot);

    // Re-initialize quadratures
    if (nq_>0) THROWING(IDAQuadReInit, m->mem, m->q);

    // Correct initial conditions, if necessary
    if (calc_ic_) {
      THROWING(IDACalcIC, m->mem, IDA_YA_YDP_INIT , first_time_);
      THROWING(IDAGetConsistentIC, m->mem, m->xz, m->xzdot);
    }

    // Re-initialize backward integration
    if (nrx_>0) THROWING(IDAAdjReInit, m->mem);

    // Set the stop time of the integration -- don't integrate past this point
    if (stop_at_end_) setStopTime(m, grid_.back());

    log("IdasInterface::reset", "end");
  }

int PsetupHeat(realtype tt, 
               N_Vector uu, N_Vector up, N_Vector rr, 
               realtype c_j, void *prec_data, 
               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
  
  long int i, j, offset, loc, mm;
  realtype *ppv, pelinv;
  UserData data;
  
  data = (UserData) prec_data;
  ppv = NV_DATA_S(data->pp);
  mm = data->mm;

  /* Initialize the entire vector to 1., then set the interior points to the
     correct value for preconditioning. */
  N_VConst(ONE,data->pp);
  
  /* Compute the inverse of the preconditioner diagonal elements. */
  pelinv = ONE/(c_j + FOUR*data->coeff); 
  
  for (j = 1; j < mm-1; j++) {
    offset = mm * j;
    for (i = 1; i < mm-1; i++) {
      loc = offset + i;
      ppv[loc] = pelinv;
    }
  }
  
  return(0);  
}

static void SetInitialProfiles(N_Vector y, realtype dx, realtype dz)
{
  int jx, jz;
  realtype x, z, cx, cz;
  realtype *ydata;

  /* Set pointer to data array in vector y. */

  ydata = NV_DATA_S(y);

  /* Load initial profiles of c1 and c2 into y vector */

  for (jz=0; jz < MZ; jz++) {
    z = ZMIN + jz*dz;
    cz = SQR(RCONST(0.1)*(z - ZMID));
    cz = ONE - cz + RCONST(0.5)*SQR(cz);
    for (jx=0; jx < MX; jx++) {
      x = XMIN + jx*dx;
      cx = SQR(RCONST(0.1)*(x - XMID));
      cx = ONE - cx + RCONST(0.5)*SQR(cx);
      IJKth(ydata,1,jx,jz) = C1_SCALE*cx*cz; 
      IJKth(ydata,2,jx,jz) = C2_SCALE*cx*cz;
    }
  }
}

void ode_solver_init(ode_solver* solver, const double t0,  double* y0, int lenY, double* p, int lenP ){
  int i,flag;
  /* Get parameters */
  
  if (p != 0){
    if (lenP != solver->odeModel->P) {
      fprintf(stderr,"ode_solver_init: lenP must be equal %d, the number of parameters in the ode model.\n",solver->odeModel->P);
      return ;
    }
    
    for(i = 0; i < solver->odeModel->P; i++)
      solver->params[i] = p[i];
  }
  
  /* Get initial conditions */
  if(y0 != 0){
    if( lenY != solver->odeModel->N ){
      fprintf(stderr,"ode_solver_init: lenY must be equal %d, the number of variables in the ode model.\n",solver->odeModel->N);
      return ;
    }
    NV_DATA_S(solver->y) = y0;
  }
  
  /* initialise */
  flag = CVodeInit(solver->cvode_mem, solver->odeModel->vf_eval, t0, solver->y);
  flag = CVodeSetUserData(solver->cvode_mem, solver->params);
  flag = CVDense(solver->cvode_mem, solver->odeModel->N);
  flag = CVDlsSetDenseJacFn(solver->cvode_mem, solver->odeModel->vf_jac);
  flag = CVodeSStolerances(solver->cvode_mem, ODE_SOLVER_REL_ERR, ODE_SOLVER_ABS_ERR);
  flag = CVodeSetMaxNumSteps(solver->cvode_mem, ODE_SOLVER_MX_STEPS);
  
}

static void PrintOutput(void *mem, realtype t, N_Vector y)
{
  realtype *yval;
  int retval, kused;
  long int nst;
  realtype hused;

  yval  = NV_DATA_S(y);

  retval = IDAGetLastOrder(mem, &kused);
  check_flag(&retval, "IDAGetLastOrder", 1);
  retval = IDAGetNumSteps(mem, &nst);
  check_flag(&retval, "IDAGetNumSteps", 1);
  retval = IDAGetLastStep(mem, &hused);
  check_flag(&retval, "IDAGetLastStep", 1);
#if defined(SUNDIALS_EXTENDED_PRECISION)
  printf("%10.4Le %12.4Le %12.4Le %12.4Le | %3ld  %1d %12.4Le\n", 
         t, yval[0], yval[1], yval[2], nst, kused, hused);
#elif defined(SUNDIALS_DOUBLE_PRECISION)
  printf("%10.4e %12.4e %12.4e %12.4e | %3ld  %1d %12.4e\n", 
         t, yval[0], yval[1], yval[2], nst, kused, hused);
#else
  printf("%10.4e %12.4e %12.4e %12.4e | %3ld  %1d %12.4e\n", 
         t, yval[0], yval[1], yval[2], nst, kused, hused);
#endif
}

static int PSolve(realtype tn, N_Vector y, N_Vector fy,
                  N_Vector r, N_Vector z,
                  realtype gamma, realtype delta,
                  int lr, void *P_data, N_Vector vtemp)
{
  realtype **(*P)[MZ];
  long int *(*pivot)[MZ];
  int jx, jz;
  realtype *zdata, *v;
  UserData data;

  /* Extract the P and pivot arrays from P_data. */

  data = (UserData) P_data;
  P = data->P;
  pivot = data->pivot;
  zdata = NV_DATA_S(z);

  N_VScale(ONE, r, z);

  /* Solve the block-diagonal system Px = r using LU factors stored
     in P and pivot data in pivot, and return the solution in z. */

  for (jx=0; jx < MX; jx++) {
    for (jz=0; jz < MZ; jz++) {
      v = &(IJKth(zdata, 1, jx, jz));
      gesl(P[jx][jz], NUM_SPECIES, pivot[jx][jz], v);
    }
  }

  return(0);
}