void Cvode::init_prepare() {
	if (init_global()) {
		if (y_) {
			N_VDestroy(y_);
			y_ = nil;
		}
		if (mem_) {
			CVodeFree(mem_);
			mem_ = nil;
		}
		if (atolnvec_) {
			N_VDestroy(atolnvec_);
			atolnvec_ = nil;
		}
		if (daspk_) {
			delete daspk_;
			daspk_ = nil;
		}
		init_eqn();
		if (neq_ > 0) {
			y_ = nvnew(neq_);
			if (use_daspk_) {
				alloc_daspk();
			}else{
				alloc_cvode();
			}
			if (maxstate_) {
				activate_maxstate(false);
				activate_maxstate(true);
			}
		}
	}
}

static void CVAckpntDelete(CkpntMem *ck_memPtr)
{
  CkpntMem tmp;
  int j;

  if (*ck_memPtr != NULL) {

    /* store head of list */
    tmp = *ck_memPtr;

    /* move head of list */
    *ck_memPtr = (*ck_memPtr)->ck_next;

    /* free N_Vectors in tmp */
    for (j=0;j<=tmp->ck_q;j++) N_VDestroy(tmp->ck_zn[j]);
    if (tmp->ck_zqm != 0) N_VDestroy(tmp->ck_zn[tmp->ck_zqm]);

    /* free N_Vectors for quadratures in tmp 
     Note that at the check point at t_initial, only znQ_[0] 
     was allocated*/
    if(tmp->ck_quadr) {
      if(tmp->ck_next != NULL) {
        for (j=0;j<=tmp->ck_q;j++) N_VDestroy(tmp->ck_znQ[j]);
        if (tmp->ck_zqm != 0) N_VDestroy(tmp->ck_znQ[tmp->ck_zqm]);
      } else {
        N_VDestroy(tmp->ck_znQ[0]);
      }
    }

    free(tmp);

  }
}

void FIDA_SETVIN(char key_name[], realtype *vval, int *ier)
{
  N_Vector Vec;

  *ier = 0;

  if (!strncmp(key_name,"ID_VEC",6)) {
    Vec = NULL;
    Vec = N_VCloneEmpty(F2C_IDA_vec);
    if (Vec == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(vval, Vec);
    IDASetId(IDA_idamem, Vec);
    N_VDestroy(Vec);
  } else if (!strncmp(key_name,"CONSTR_VEC",10)) {
    Vec = NULL;
    Vec = N_VCloneEmpty(F2C_IDA_vec);
    if (Vec == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(vval, Vec);
    IDASetConstraints(IDA_idamem, Vec);
    N_VDestroy(Vec);
  } else {
    *ier = -99;
    printf("FIDASETVIN: Unrecognized key.\n\n");
  }

}

static booleantype CVAallocVectors(CVadjMem ca_mem)
{

  CVodeMem cv_mem;

  cv_mem = ca_mem->cv_mem;

  Y0 = N_VClone(tempv);
  if (Y0 == NULL) {
    return(FALSE);
  }

  Y1 = N_VClone(tempv);
  if (Y1 == NULL) {
    N_VDestroy(Y0);
    return(FALSE);
  }

  ytmp = N_VClone(tempv);
  if (ytmp == NULL) {
    N_VDestroy(Y0);
    N_VDestroy(Y1);
    return(FALSE);
  }

  return(TRUE);

}

/*---------------------------------------------------------------
 Function : PcgMalloc
 --------------------------------------------------------------*/
PcgMem PcgMalloc(int l_max, N_Vector vec_tmpl)
{
  PcgMem mem;
  N_Vector r, p, z, Ap;

  /* Check the input parameters */
  if (l_max <= 0) return(NULL);

  /* Create temporary arrays */
  r = N_VClone(vec_tmpl);
  if (r == NULL) {
    return(NULL);
  }

  p = N_VClone(vec_tmpl);
  if (p == NULL) {
    N_VDestroy(r);
    return(NULL);
  }

  z = N_VClone(vec_tmpl);
  if (z == NULL) {
    N_VDestroy(r);
    N_VDestroy(p);
    return(NULL);
  }

  Ap = N_VClone(vec_tmpl);
  if (Ap == NULL) {
    N_VDestroy(r);
    N_VDestroy(p);
    N_VDestroy(z);
    return(NULL);
  }

  /* Get memory for an PcgMemRec containing PCG vectors */
  mem = NULL;
  mem = (PcgMem) malloc(sizeof(PcgMemRec));
  if (mem == NULL) {
    N_VDestroy(r);
    N_VDestroy(p);
    N_VDestroy(z);
    N_VDestroy(Ap);
    return(NULL);
  }

  /* Set the structure fields */
  mem->l_max = l_max;
  mem->r     = r;
  mem->p     = p;
  mem->z     = z;
  mem->Ap    = Ap;

  /* Return the pointer to PCG memory */
  return(mem);
}

/*---------------------------------------------------------------
 Function : PcgFree
 --------------------------------------------------------------*/
void PcgFree(PcgMem mem)
{
  if (mem == NULL) return;

  N_VDestroy(mem->r);
  N_VDestroy(mem->p);
  N_VDestroy(mem->z);
  N_VDestroy(mem->Ap);

  free(mem); mem = NULL;
}

static void CVSpgmrFree(CVodeMem cv_mem)
{
  CVSpgmrMem cvspgmr_mem;

  cvspgmr_mem = (CVSpgmrMem) lmem;
  
  N_VDestroy(ytemp);
  N_VDestroy(x);
  SpgmrFree(spgmr_mem);
  free(cvspgmr_mem);
}

static void CVDiagFree(CVodeMem cv_mem)
{
  CVDiagMem cvdiag_mem;
  
  cvdiag_mem = (CVDiagMem) lmem;

  N_VDestroy(M);
  N_VDestroy(bit);
  N_VDestroy(bitcomp);
  free(cvdiag_mem); cvdiag_mem = NULL;
}

static void CVAhermiteFree(DtpntMem *dt_mem, long int steps)
{
  long int i;
  HermiteDataMem content;

  for (i=0; i<=steps; i++) {
    content = (HermiteDataMem) (dt_mem[i]->content);
    N_VDestroy(content->y);
    N_VDestroy(content->yd);
    free(dt_mem[i]->content);
    free(dt_mem[i]);
  }
}

static void FreeUserData(WebData wdata)
{
  int i, ngrp;

  ngrp = wdata->ngrp;
  for(i=0; i < ngrp; i++) {
    destroyMat((wdata->P)[i]);
    destroyArray((wdata->pivot)[i]);
  }
  N_VDestroy(wdata->rewt);
  N_VDestroy(wdata->vtemp);
  free(wdata);
}

static int CVDiagFree(CVodeMem cv_mem)
{
  CVDiagMem cvdiag_mem;
  
  cvdiag_mem = (CVDiagMem) lmem;

  N_VDestroy(M);
  N_VDestroy(bit);
  N_VDestroy(bitcomp);
  free(cvdiag_mem);
  cv_mem->cv_lmem = NULL;
  
  return(0);
}

static void CVSpbcgFree(CVodeMem cv_mem)
{
  CVSpilsMem cvspils_mem;
  SpbcgMem spbcg_mem;

  cvspils_mem = (CVSpilsMem) lmem;

  spbcg_mem = (SpbcgMem) spils_mem;

  N_VDestroy(ytemp);
  N_VDestroy(x);
  SpbcgFree(spbcg_mem);
  free(cvspils_mem); cvspils_mem = NULL;
}

static void cpSpbcgFree(CPodeMem cp_mem)
{
  CPSpilsMem cpspils_mem;
  SpbcgMem spbcg_mem;

  cpspils_mem = (CPSpilsMem) lmem;

  spbcg_mem = (SpbcgMem) spils_mem;

  SpbcgFree(spbcg_mem);
  N_VDestroy(x);
  N_VDestroy(ytemp);
  if (ode_type == CP_EXPL) N_VDestroy(yptemp);
  free(cpspils_mem); cpspils_mem = NULL;
}

static int IDASpgmrFree(IDAMem IDA_mem)
{
  IDASpilsMem idaspils_mem;
  SpgmrMem spgmr_mem;

  idaspils_mem = (IDASpilsMem) lmem;
  
  spgmr_mem = (SpgmrMem) spils_mem;

  N_VDestroy(ytemp);
  N_VDestroy(xx);
  SpgmrFree(spgmr_mem);
  free(lmem); lmem = NULL;

  return(0);
}

/*-----------------------------------------------------------------
  cvDlsFree
  -----------------------------------------------------------------
  This routine frees memory associates with the CVDls solver 
  interface.
  -----------------------------------------------------------------*/
int cvDlsFree(CVodeMem cv_mem)
{
  CVDlsMem cvdls_mem;

  /* Return immediately if cv_mem or cv_mem->cv_lmem are NULL */
  if (cv_mem == NULL)  return (CVDLS_SUCCESS);
  if (cv_mem->cv_lmem == NULL)  return(CVDLS_SUCCESS);
  cvdls_mem = (CVDlsMem) cv_mem->cv_lmem;

  /* Free x vector */
  if (cvdls_mem->x) {
    N_VDestroy(cvdls_mem->x);
    cvdls_mem->x = NULL;
  }

  /* Free savedJ memory */
  if (cvdls_mem->savedJ) {
    SUNMatDestroy(cvdls_mem->savedJ);
    cvdls_mem->savedJ = NULL;
  }

  /* Nullify other SUNMatrix pointer */
  cvdls_mem->A = NULL;

  /* free CVDls interface structure */
  free(cv_mem->cv_lmem);
  
  return(CVDLS_SUCCESS);
}

/* Fortran interface to C routine ARKStepFree; see farkode.h for
   further details */
void FARK_FREE() {

  ARKodeMem ark_mem;
  ark_mem = (ARKodeMem) ARK_arkodemem;

  /* free user_data structure */
  if (ark_mem->user_data)
    free(ark_mem->user_data);
  ark_mem->user_data = NULL;

  /* free main integrator memory structure (internally
     frees time step module, rootfinding, interpolation structures) */
  ARKStepFree(&ARK_arkodemem);

  /* free interface vector / matrices / linear solvers */
  N_VSetArrayPointer(NULL, F2C_ARKODE_vec);
  N_VDestroy(F2C_ARKODE_vec);
  if (F2C_ARKODE_matrix)
    SUNMatDestroy(F2C_ARKODE_matrix);
  if (F2C_ARKODE_mass_matrix)
    SUNMatDestroy(F2C_ARKODE_mass_matrix);
  if (F2C_ARKODE_linsol)
    SUNLinSolFree(F2C_ARKODE_linsol);
  if (F2C_ARKODE_mass_sol)
    SUNLinSolFree(F2C_ARKODE_mass_sol);
  return;
}

int IDASetId(void *ida_mem, N_Vector id)
{
  IDAMem IDA_mem;

  if (ida_mem==NULL) {
    IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDASetId", MSG_NO_MEM);
    return(IDA_MEM_NULL);
  }

  IDA_mem = (IDAMem) ida_mem;

  if (id == NULL) {
    if (IDA_mem->ida_idMallocDone) {
      N_VDestroy(IDA_mem->ida_id);
      lrw -= lrw1;
      liw -= liw1;
    }
    IDA_mem->ida_idMallocDone = FALSE;    
    return(IDA_SUCCESS);
  }

  if ( !(IDA_mem->ida_idMallocDone) ) {
    IDA_mem->ida_id = NULL;
    IDA_mem->ida_id = N_VClone(id);
    lrw += lrw1;
    liw += liw1;
    IDA_mem->ida_idMallocDone = TRUE;
  }

  /* Load the id vector */

  N_VScale(ONE, id, IDA_mem->ida_id);

  return(IDA_SUCCESS);
}

/*-----------------------------------------------------------------
  cvDlsDenseDQJac 
  -----------------------------------------------------------------
  This routine generates a dense difference quotient approximation 
  to the Jacobian of f(t,y). It assumes that a dense SUNMatrix is 
  stored column-wise, and that elements within each column are 
  contiguous. The address of the jth column of J is obtained via
  the accessor function SUNDenseMatrix_Column, and this pointer 
  is associated with an N_Vector using the N_VSetArrayPointer
  function.  Finally, the actual computation of the jth column of 
  the Jacobian is done with a call to N_VLinearSum.
  -----------------------------------------------------------------*/ 
int cvDlsDenseDQJac(realtype t, N_Vector y, N_Vector fy, 
                    SUNMatrix Jac, CVodeMem cv_mem, N_Vector tmp1)
{
  realtype fnorm, minInc, inc, inc_inv, yjsaved, srur;
  realtype *y_data, *ewt_data;
  N_Vector ftemp, jthCol;
  sunindextype j, N;
  int retval = 0;
  CVDlsMem cvdls_mem;

  /* access DlsMem interface structure */
  cvdls_mem = (CVDlsMem) cv_mem->cv_lmem;

  /* access matrix dimension */
  N = SUNDenseMatrix_Rows(Jac);

  /* Rename work vector for readibility */
  ftemp = tmp1;

  /* Create an empty vector for matrix column calculations */
  jthCol = N_VCloneEmpty(tmp1);

  /* Obtain pointers to the data for ewt, y */
  ewt_data = N_VGetArrayPointer(cv_mem->cv_ewt);
  y_data   = N_VGetArrayPointer(y);

  /* Set minimum increment based on uround and norm of f */
  srur = SUNRsqrt(cv_mem->cv_uround);
  fnorm = N_VWrmsNorm(fy, cv_mem->cv_ewt);
  minInc = (fnorm != ZERO) ?
    (MIN_INC_MULT * SUNRabs(cv_mem->cv_h) * cv_mem->cv_uround * N * fnorm) : ONE;

  for (j = 0; j < N; j++) {

    /* Generate the jth col of J(tn,y) */

    N_VSetArrayPointer(SUNDenseMatrix_Column(Jac,j), jthCol);

    yjsaved = y_data[j];
    inc = SUNMAX(srur*SUNRabs(yjsaved), minInc/ewt_data[j]);
    y_data[j] += inc;

    retval = cv_mem->cv_f(t, y, ftemp, cv_mem->cv_user_data);
    cvdls_mem->nfeDQ++;
    if (retval != 0) break;
    
    y_data[j] = yjsaved;

    inc_inv = ONE/inc;
    N_VLinearSum(inc_inv, ftemp, -inc_inv, fy, jthCol);

    /* DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol);   /\*UNNECESSARY?? *\/ */
  }

  /* Destroy jthCol vector */
  N_VSetArrayPointer(NULL, jthCol);  /* SHOULDN'T BE NEEDED */
  N_VDestroy(jthCol);

  return(retval);
}

void FCVBandJac(long int N, long int mupper, long int mlower,
                BandMat J, realtype t,
                N_Vector y, N_Vector fy, void *jac_data,
                N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3)
{
  N_Vector ewt;
  realtype *ydata, *fydata, *jacdata, *ewtdata, *v1data, *v2data, *v3data;
  realtype h;
  long int eband;

  ewt = N_VClone(y);

  CVodeGetErrWeights(CV_cvodemem, ewt);
  CVodeGetLastStep(CV_cvodemem, &h);

  ydata   = N_VGetArrayPointer(y);
  fydata  = N_VGetArrayPointer(fy);
  v1data  = N_VGetArrayPointer(vtemp1);
  v2data  = N_VGetArrayPointer(vtemp2);
  v3data  = N_VGetArrayPointer(vtemp3);
  ewtdata = N_VGetArrayPointer(ewt);

  eband = (J->smu) + mlower + 1;
  jacdata = BAND_COL(J,0) - mupper;


  FCV_BJAC(&N, &mupper, &mlower, &eband, 
           &t, ydata, fydata, jacdata, 
           ewtdata, &h, v1data, v2data, v3data);

  N_VDestroy(ewt);

}

/* Fortran interface routine to set residual tolerance 
   scalar/array; functions as an all-in-one interface to the C 
   routines ARKodeResStolerance and ARKodeResVtolerance; 
   see farkode.h for further details */
void FARK_SETRESTOLERANCE(int *itol, realtype *atol, int *ier) {

  N_Vector Vatol;
  realtype abstol;

  *ier = 0;

  /* Set tolerance, based on itol argument */
  abstol=1.e-9;
  switch (*itol) {
  case 1:
    if (*atol > 0.0)  abstol = *atol;
    *ier = ARKodeResStolerance(ARK_arkodemem, abstol); 
    break;
  case 2:
    Vatol = NULL;
    Vatol = N_VCloneEmpty(F2C_ARKODE_vec);
    if (Vatol == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(atol, Vatol);
    if (N_VMin(Vatol) <= 0.0)  N_VConst(abstol, Vatol);
    *ier = ARKodeResVtolerance(ARK_arkodemem, Vatol);
    N_VDestroy(Vatol);
    break;
  }

  return;
}