static int cvNlsFPFunction(N_Vector ycor, N_Vector res, void* cvode_mem)
{
 CVodeMem cv_mem;
  int retval;

  if (cvode_mem == NULL) {
    cvProcessError(NULL, CV_MEM_NULL, "CVODE", "cvNlsFPFunction", MSGCV_NO_MEM);
    return(CV_MEM_NULL);
  }
  cv_mem = (CVodeMem) cvode_mem;

  /* update the state based on the current correction */
  N_VLinearSum(ONE, cv_mem->cv_zn[0], ONE, ycor, cv_mem->cv_y);

  /* evaluate the rhs function */
  retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_y, res,
                        cv_mem->cv_user_data);
  cv_mem->cv_nfe++;
  if (retval < 0) return(CV_RHSFUNC_FAIL);
  if (retval > 0) return(RHSFUNC_RECVR);

  N_VLinearSum(cv_mem->cv_h, res, -ONE, cv_mem->cv_zn[1], res);
  N_VScale(cv_mem->cv_rl1, res, res);

  return(CV_SUCCESS);
}

static int CVSpgmrDQJtimes(N_Vector v, N_Vector Jv, realtype t, 
                           N_Vector y, N_Vector fy,
                           void *jac_data, N_Vector work)
{
  CVodeMem cv_mem;
  CVSpgmrMem cvspgmr_mem;
  realtype vnrm;

  /* jac_data is cvode_mem */
  cv_mem = (CVodeMem) jac_data;
  cvspgmr_mem = (CVSpgmrMem) lmem;

  /* Evaluate norm of v */
  vnrm = N_VWrmsNorm(v, ewt);

  /* Set work = y + (1/vnrm) v */
  N_VLinearSum(ONE/vnrm, v, ONE, y, work);

  /* Set Jv = f(tn, work) */
  f(t, work, Jv, f_data); 
  nfeSG++;

  /* Replace Jv by vnrm*(Jv - fy) */
  N_VLinearSum(ONE, Jv, -ONE, fy, Jv);
  N_VScale(vnrm, Jv, Jv);

  return(0);
}

static int cvNlsFPFunctionSensStg1(N_Vector ycor, N_Vector res, void* cvode_mem)
{
  CVodeMem cv_mem;
  int retval, is;

  if (cvode_mem == NULL) {
    cvProcessError(NULL, CV_MEM_NULL, "CVODES",
                   "cvNlsFPFunctionSensStg1", MSGCV_NO_MEM);
    return(CV_MEM_NULL);
  }
  cv_mem = (CVodeMem) cvode_mem;

  /* get index of current sensitivity solve */
  is = cv_mem->sens_solve_idx;

  /* update the sensitivities based on the current correction */
  N_VLinearSum(ONE, cv_mem->cv_znS[0][is], ONE, ycor, cv_mem->cv_yS[is]);

  /* evaluate the sensitivity rhs function */
  retval = cvSensRhs1Wrapper(cv_mem, cv_mem->cv_tn,
                             cv_mem->cv_y, cv_mem->cv_ftemp,
                             is, cv_mem->cv_yS[is], res,
                             cv_mem->cv_vtemp1, cv_mem->cv_vtemp2);

  if (retval < 0) return(CV_SRHSFUNC_FAIL);
  if (retval > 0) return(SRHSFUNC_RECVR);

  /* evaluate sensitivity fixed point function */
  N_VLinearSum(cv_mem->cv_h, res, -ONE, cv_mem->cv_znS[1][is], res);
  N_VScale(cv_mem->cv_rl1, res, res);

  return(CV_SUCCESS);
}

static int IDANewyyp(IDAMem IDA_mem, realtype lambda)
{
  int retval;
  
  retval = IDA_SUCCESS;

  /* IDA_YA_YDP_INIT case: ynew  = yy0 - lambda*delta    where id_i = 0
                           ypnew = yp0 - cj*lambda*delta where id_i = 1. */
  if(icopt == IDA_YA_YDP_INIT) {

    N_VProd(id, delta, dtemp);
    N_VLinearSum(ONE, yp0, -cj*lambda, dtemp, ypnew);
    N_VLinearSum(ONE, delta, -ONE, dtemp, dtemp);
    N_VLinearSum(ONE, yy0, -lambda, dtemp, ynew);

  }else if(icopt == IDA_Y_INIT) {

    /* IDA_Y_INIT case: ynew = yy0 - lambda*delta. (ypnew = yp0 preset.) */
    N_VLinearSum(ONE, yy0, -lambda, delta, ynew);
  }

  if(sensi && (ism==IDA_SIMULTANEOUS))
    retval = IDASensNewyyp(IDA_mem, lambda);
  
  return(retval);

}

static int IDASensNewyyp(IDAMem IDA_mem, realtype lambda)
{
  int is;

  if(icopt == IDA_YA_YDP_INIT) {

  /* IDA_YA_YDP_INIT case: 
     - ySnew  = yS0  - lambda*deltaS    where id_i = 0
     - ypSnew = ypS0 - cj*lambda*delta  where id_i = 1. */    

    for(is=0; is<Ns; is++) {
      
      /* It is ok to use dtemp as temporary vector here. */
      N_VProd(id, deltaS[is], dtemp);
      N_VLinearSum(ONE, ypS0[is], -cj*lambda, dtemp, ypS0new[is]);
      N_VLinearSum(ONE, deltaS[is], -ONE, dtemp, dtemp);
      N_VLinearSum(ONE, yyS0[is], -lambda, dtemp, yyS0new[is]);
    } /* end loop is */
  }else { 

    /* IDA_Y_INIT case: 
       - ySnew = yS0 - lambda*deltaS. (ypnew = yp0 preset.) */

    for(is=0; is<Ns; is++)
      N_VLinearSum(ONE, yyS0[is], -lambda, deltaS[is], yyS0new[is]);
  } /* end loop is */
  return(IDA_SUCCESS);
}

/*---------------------------------------------------------------
 ARKSpilsAtimes:

 This routine generates the matrix-vector product z = Av, where
 A = M - gamma*J. The product M*v is obtained either by calling 
 the mtimes routine or by just using v (if M=I).  The product 
 J*v is obtained by calling the jtimes routine. It is then scaled 
 by -gamma and added to M*v to obtain A*v. The return value is 
 the same as the values returned by jtimes and mtimes -- 
 0 if successful, nonzero otherwise.
---------------------------------------------------------------*/
int ARKSpilsAtimes(void *arkode_mem, N_Vector v, N_Vector z)
{
  ARKodeMem   ark_mem;
  ARKSpilsMem arkspils_mem;
  int jtflag, mtflag;

  ark_mem = (ARKodeMem) arkode_mem;
  arkspils_mem = (ARKSpilsMem) ark_mem->ark_lmem;

  jtflag = arkspils_mem->s_jtimes(v, z, ark_mem->ark_tn, 
				  arkspils_mem->s_ycur, 
				  arkspils_mem->s_fcur, 
				  arkspils_mem->s_j_data, 
				  arkspils_mem->s_ytemp);
  arkspils_mem->s_njtimes++;
  if (jtflag != 0) return(jtflag);

  /* Compute mass matrix vector product and add to result */
  if (ark_mem->ark_mass_matrix) {
    mtflag = ARKSpilsMtimes(arkode_mem, v, arkspils_mem->s_ytemp);
    if (mtflag != 0) return(mtflag);
    N_VLinearSum(ONE, arkspils_mem->s_ytemp, -ark_mem->ark_gamma, z, z);
  } else {
    N_VLinearSum(ONE, v, -ark_mem->ark_gamma, z, z);
  }

  return(0);
}

static int idaNlsResidual(N_Vector ycor, N_Vector res, void* ida_mem)
{
  IDAMem IDA_mem;
  int retval;

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

  /* update yy and yp based on the current correction */
  N_VLinearSum(ONE, IDA_mem->ida_yypredict, ONE, ycor, IDA_mem->ida_yy);
  N_VLinearSum(ONE, IDA_mem->ida_yppredict, IDA_mem->ida_cj, ycor, IDA_mem->ida_yp);

  /* evaluate residual */
  retval = IDA_mem->ida_res(IDA_mem->ida_tn, IDA_mem->ida_yy, IDA_mem->ida_yp,
                            res, IDA_mem->ida_user_data);

  /* increment the number of residual evaluations */
  IDA_mem->ida_nre++;

  /* save a copy of the residual vector in savres */
  N_VScale(ONE, res, IDA_mem->ida_savres);

  if (retval < 0) return(IDA_RES_FAIL);
  if (retval > 0) return(IDA_RES_RECVR);

  return(IDA_SUCCESS);
}

int KINSpilsDQJtimes(N_Vector v, N_Vector Jv,
                     N_Vector u, booleantype *new_u, 
                     void *data)
{
  realtype sigma, sigma_inv, sutsv, sq1norm, sign, vtv;
  KINMem kin_mem;
  KINSpilsMem kinspils_mem;
  int retval;

  /* data is kin_mem */

  kin_mem = (KINMem) data;
  kinspils_mem = (KINSpilsMem) lmem;

  /* scale the vector v and put Du*v into vtemp1 */

  N_VProd(v, uscale, vtemp1);

  /* scale u and put into Jv (used as a temporary storage) */

  N_VProd(u, uscale, Jv);

  /* compute dot product (Du*u).(Du*v) */

  sutsv = N_VDotProd(Jv, vtemp1);

  /* compute dot product (Du*v).(Du*v) */

  vtv = N_VDotProd(vtemp1, vtemp1);

  sq1norm = N_VL1Norm(vtemp1);

  sign = (sutsv >= ZERO) ? ONE : -ONE ;
 
  /*  this expression for sigma is from p. 469, Brown and Saad paper */

  sigma = sign*sqrt_relfunc*SUNMAX(SUNRabs(sutsv),sq1norm)/vtv;

  sigma_inv = ONE/sigma;

  /* compute the u-prime at which to evaluate the function func */

  N_VLinearSum(ONE, u, sigma, v, vtemp1);
 
  /* call the system function to calculate func(u+sigma*v) */

  retval = func(vtemp1, vtemp2, user_data);    
  nfes++;
  if (retval != 0) return(retval);

  /* finish the computation of the difference quotient */

  N_VLinearSum(sigma_inv, vtemp2, -sigma_inv, fval, Jv);

  return(0);
}

int IDASpilsDQJtimes(realtype tt,
                     N_Vector yy, N_Vector yp, N_Vector rr,
                     N_Vector v, N_Vector Jv,
                     realtype c_j, void *data,
                     N_Vector work1, N_Vector work2)
{
  IDAMem IDA_mem;
  IDASpilsMem idaspils_mem;
  N_Vector y_tmp, yp_tmp;
  realtype sig, siginv;
  int iter, retval;

  /* data is ida_mem */
  IDA_mem = (IDAMem) data;
  idaspils_mem = (IDASpilsMem) lmem;

  switch(ils_type) {
  case SPILS_SPGMR:
    sig = sqrtN*dqincfac;
    break;
  case SPILS_SPBCG:
    sig = dqincfac/N_VWrmsNorm(v, ewt);
    break;
  case SPILS_SPTFQMR:
    sig = dqincfac/N_VWrmsNorm(v, ewt);
    break;
  }

  /* Rename work1 and work2 for readibility */
  y_tmp  = work1;
  yp_tmp = work2;

  for (iter=0; iter<MAX_ITERS; iter++) {

    /* Set y_tmp = yy + sig*v, yp_tmp = yp + cj*sig*v. */
    N_VLinearSum(sig, v, ONE, yy, y_tmp);
    N_VLinearSum(c_j*sig, v, ONE, yp, yp_tmp);

    /* Call res for Jv = F(t, y_tmp, yp_tmp), and return if it failed. */
    retval = res(tt, y_tmp, yp_tmp, Jv, user_data);
    nres++;
    if (retval == 0) break;
    if (retval < 0)  return(-1);

    sig *= PT25;
  }

  if (retval > 0) return(+1);

  /* Set Jv to [Jv - rr]/sig and return. */
  siginv = ONE/sig;
  N_VLinearSum(siginv, Jv, -siginv, rr, Jv);

  return(0);

}

static int cpNewtonIterationImpl(CPodeMem cp_mem)
{
  int m, retval;
  realtype del, delp, dcon;

  /* Initialize counters */
  mnewt = m = 0;

  /* Initialize delp to avoid compiler warning message */
  del = delp = ZERO;

  /* Looping point for Newton iteration */
  loop {

    /* Evaluate the residual of the nonlinear system*/
    N_VScale(-gamma, ftemp, tempv);

    /* Call the lsolve function (solution in tempv) */
    retval = lsolve(cp_mem, tempv, ewt, y, yp, ftemp);
    nni++;
    if (retval < 0) return(CP_LSOLVE_FAIL);
    if (retval > 0) return(CONV_FAIL);

    /* Get WRMS norm of correction and add correction to acor, y, and yp */
    del = N_VWrmsNorm(tempv, ewt);
    N_VLinearSum(ONE, acor, ONE,       tempv, acor);
    N_VLinearSum(ONE, y,    ONE,       tempv, y);
    N_VLinearSum(ONE, yp,   ONE/gamma, tempv, yp);

    /* Test for convergence.  If m > 0, an estimate of the convergence
       rate constant is stored in crate, and used in the test.        */
    if (m > 0) crate = MAX(NLS_CRDOWN * crate, del/delp);
    dcon = del * MIN(ONE, crate) / tq[4];    
    if (dcon <= ONE) {
      acnrm = (m==0) ? del : N_VWrmsNorm(acor, ewt);
      return(CP_SUCCESS);
    }
    
    mnewt = ++m;
    
    /* Stop at maxcor iterations or if iter. seems to be diverging. */
    if ((m == maxcor) || ((m >= 2) && (del > NLS_RDIV*delp)))  return(CONV_FAIL);
    
    /* Save norm of correction, evaluate f, and loop again */
    delp = del;
    retval = fi(tn, y, yp, ftemp, f_data);
    nfe++;
    if (retval < 0) return(CP_ODEFUNC_FAIL);
    if (retval > 0) return(ODEFUNC_RECVR);

  } /* end loop */
}

int ModifiedGS(N_Vector *v, realtype **h, int k, int p, 
               realtype *new_vk_norm)
{
  int  i, k_minus_1, i0;
  realtype new_norm_2, new_product, vk_norm, temp;
  
  vk_norm = RSqrt(N_VDotProd(v[k],v[k]));
  k_minus_1 = k - 1;
  i0 = MAX(k-p, 0);
  
  /* Perform modified Gram-Schmidt */
  
  for (i=i0; i < k; i++) {
    h[i][k_minus_1] = N_VDotProd(v[i], v[k]);
    N_VLinearSum(ONE, v[k], -h[i][k_minus_1], v[i], v[k]);
  }

  /* Compute the norm of the new vector at v[k] */

  *new_vk_norm = RSqrt(N_VDotProd(v[k], v[k]));

  /* If the norm of the new vector at v[k] is less than
     FACTOR (== 1000) times unit roundoff times the norm of the
     input vector v[k], then the vector will be reorthogonalized
     in order to ensure that nonorthogonality is not being masked
     by a very small vector length. */

  temp = FACTOR * vk_norm;
  if ((temp + (*new_vk_norm)) != temp) return(0);
  
  new_norm_2 = ZERO;

  for (i=i0; i < k; i++) {
    new_product = N_VDotProd(v[i], v[k]);
    temp = FACTOR * h[i][k_minus_1];
    if ((temp + new_product) == temp) continue;
    h[i][k_minus_1] += new_product;
    N_VLinearSum(ONE, v[k],-new_product, v[i], v[k]);
    new_norm_2 += SQR(new_product);
  }

  if (new_norm_2 != ZERO) {
    new_product = SQR(*new_vk_norm) - new_norm_2;
    *new_vk_norm = (new_product > ZERO) ? RSqrt(new_product) : ZERO;
  }

  return(0);
}

int ClassicalGS(N_Vector *v, realtype **h, int k, int p, 
                realtype *new_vk_norm, N_Vector temp, realtype *s)
{
  int  i, k_minus_1, i0;
  realtype vk_norm;

  k_minus_1 = k - 1;
  
  /* Perform Classical Gram-Schmidt */

  vk_norm = RSqrt(N_VDotProd(v[k], v[k]));

  i0 = MAX(k-p, 0);
  for (i=i0; i < k; i++) {
    h[i][k_minus_1] = N_VDotProd(v[i], v[k]);
  }

  for (i=i0; i < k; i++) {
    N_VLinearSum(ONE, v[k], -h[i][k_minus_1], v[i], v[k]);
  }

  /* Compute the norm of the new vector at v[k] */

  *new_vk_norm = RSqrt(N_VDotProd(v[k], v[k]));

  /* Reorthogonalize if necessary */

  if ((FACTOR * (*new_vk_norm)) < vk_norm) {

    for (i=i0; i < k; i++) {
      s[i] = N_VDotProd(v[i], v[k]);
    }

    if (i0 < k) {
      N_VScale(s[i0], v[i0], temp);
      h[i0][k_minus_1] += s[i0];
    }
    for (i=i0+1; i < k; i++) {
      N_VLinearSum(s[i], v[i], ONE, temp, temp);
      h[i][k_minus_1] += s[i];
    }
    N_VLinearSum(ONE, v[k], -ONE, temp, v[k]);

    *new_vk_norm = RSqrt(N_VDotProd(v[k],v[k]));
  }

  return(0);
}

/*-----------------------------------------------------------------
  cvLsATimes

  This routine generates the matrix-vector product z = Mv, where
  M = I - gamma*J. The product J*v is obtained by calling the jtimes
  routine. It is then scaled by -gamma and added to v to obtain M*v.
  The return value is the same as the value returned by jtimes --
  0 if successful, nonzero otherwise.
  -----------------------------------------------------------------*/
int cvLsATimes(void *cvode_mem, N_Vector v, N_Vector z)
{
  CVodeMem cv_mem;
  CVLsMem  cvls_mem;
  int      retval;

  /* access CVLsMem structure */
  retval = cvLs_AccessLMem(cvode_mem, "cvLsATimes",
                           &cv_mem, &cvls_mem);
  if (retval != CVLS_SUCCESS)  return(retval);

  /* call Jacobian-times-vector product routine
     (either user-supplied or internal DQ) */
  retval = cvls_mem->jtimes(v, z, cv_mem->cv_tn,
                            cvls_mem->ycur,
                            cvls_mem->fcur,
                            cvls_mem->jt_data,
                            cvls_mem->ytemp);
  cvls_mem->njtimes++;
  if (retval != 0) return(retval);

  /* add contribution from identity matrix */
  N_VLinearSum(ONE, v, -cv_mem->cv_gamma, z, z);

  return(0);
}

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

static int KINDenseDQJac(long int n, DenseMat J,
                         N_Vector u, N_Vector fu, void *jac_data,
                         N_Vector tmp1, N_Vector tmp2)
{
  realtype inc, inc_inv, ujsaved, ujscale, sign;
  realtype *tmp2_data, *u_data, *uscale_data;
  N_Vector ftemp, jthCol;
  long int j;
  int retval;

  KINMem kin_mem;
  KINDenseMem  kindense_mem;

  /* jac_data points to kin_mem */
  kin_mem = (KINMem) jac_data;
  kindense_mem = (KINDenseMem) lmem;

  /* Save pointer to the array in tmp2 */
  tmp2_data = N_VGetArrayPointer(tmp2);

  /* Rename work vectors for readibility */
  ftemp = tmp1; 
  jthCol = tmp2;

  /* Obtain pointers to the data for u and uscale */
  u_data   = N_VGetArrayPointer(u);
  uscale_data = N_VGetArrayPointer(uscale);

  /* This is the only for loop for 0..N-1 in KINSOL */

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

    /* Generate the jth col of J(u) */

    N_VSetArrayPointer(DENSE_COL(J,j), jthCol);

    ujsaved = u_data[j];
    ujscale = ONE/uscale_data[j];
    sign = (ujsaved >= ZERO) ? ONE : -ONE;
    inc = sqrt_relfunc*MAX(ABS(ujsaved), ujscale)*sign;
    u_data[j] += inc;

    retval = func(u, ftemp, f_data);
    if (retval != 0) return(-1); 

    u_data[j] = ujsaved;

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

  }

  /* Restore original array pointer in tmp2 */
  N_VSetArrayPointer(tmp2_data, tmp2);

  /* Increment counter nfeD */
  nfeD += n;

  return(0);
}

static int IDANewy(IDAMem IDA_mem)
{
  
  /* IDA_YA_YDP_INIT case: ynew = yy0 - delta    where id_i = 0. */
  if(icopt == IDA_YA_YDP_INIT) {
    N_VProd(id, delta, dtemp);
    N_VLinearSum(ONE, delta, -ONE, dtemp, dtemp);
    N_VLinearSum(ONE, yy0, -ONE, dtemp, ynew);
    return(IDA_SUCCESS);
  }

  /* IDA_Y_INIT case: ynew = yy0 - delta. */
  N_VLinearSum(ONE, yy0, -ONE, delta, ynew);
  return(IDA_SUCCESS);

}

static void CVDenseDQJac(long int N, DenseMat J, realtype t,
                         N_Vector y, N_Vector fy, void *jac_data,
                         N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
    realtype fnorm, minInc, inc, inc_inv, yjsaved, srur;
    realtype *tmp2_data, *y_data, *ewt_data;
    N_Vector ftemp, jthCol;
    long int j;

    CVodeMem cv_mem;
    CVDenseMem  cvdense_mem;

    /* jac_data points to cvode_mem */
    cv_mem = (CVodeMem) jac_data;
    cvdense_mem = (CVDenseMem) lmem;

    /* Save pointer to the array in tmp2 */
    tmp2_data = N_VGetArrayPointer(tmp2);

    /* Rename work vectors for readibility */
    ftemp = tmp1;
    jthCol = tmp2;

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

    /* Set minimum increment based on uround and norm of f */
    srur = RSqrt(uround);
    fnorm = N_VWrmsNorm(fy, ewt);
    minInc = (fnorm != ZERO) ?
             (MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;

    /* This is the only for loop for 0..N-1 in CVODE */

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

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

        N_VSetArrayPointer(DENSE_COL(J,j), jthCol);

        yjsaved = y_data[j];
        inc = MAX(srur*ABS(yjsaved), minInc/ewt_data[j]);
        y_data[j] += inc;
        f(tn, y, ftemp, f_data);
        y_data[j] = yjsaved;

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

        DENSE_COL(J,j) = N_VGetArrayPointer(jthCol);
    }

    /* Restore original array pointer in tmp2 */
    N_VSetArrayPointer(tmp2_data, tmp2);

    /* Increment counter nfeD */
    nfeD += N;
}

void GetSol(void *cpode_mem, N_Vector yy0, realtype tol, 
            realtype tout, booleantype proj, N_Vector yref)
{
  N_Vector yy, yp;
  realtype t, x, y, xd, yd, g;
  int flag;
  long int nst, nfe, nsetups, nje, nfeLS, ncfn, netf;

  if (proj) {
    printf(" YES   ");
    CPodeSetProjFrequency(cpode_mem, 1);
  } else {
    CPodeSetProjFrequency(cpode_mem, 0);
    printf(" NO    ");
  }

  yy = N_VNew_Serial(4);
  yp = N_VNew_Serial(4);

  flag = CPodeReInit(cpode_mem, (void *)f, NULL, 0.0, yy0, NULL, CP_SS, tol, &tol);

  flag = CPode(cpode_mem, tout, &t, yy, yp, CP_NORMAL_TSTOP);

  x  = Ith(yy,1);
  y  = Ith(yy,2);
  g = ABS(x*x + y*y - 1.0);

  N_VLinearSum(1.0, yy, -1.0, yref, yy);

  N_VAbs(yy, yy);

  x  = Ith(yy,1);
  y  = Ith(yy,2);  
  xd = Ith(yy,3);
  yd = Ith(yy,4);


  printf("%9.2e  %9.2e  %9.2e  %9.2e  |  %9.2e  |",  
         Ith(yy,1),Ith(yy,2),Ith(yy,3),Ith(yy,4),g);

  CPodeGetNumSteps(cpode_mem, &nst);
  CPodeGetNumFctEvals(cpode_mem, &nfe);
  CPodeGetNumLinSolvSetups(cpode_mem, &nsetups);
  CPodeGetNumErrTestFails(cpode_mem, &netf);
  CPodeGetNumNonlinSolvConvFails(cpode_mem, &ncfn);

  CPDlsGetNumJacEvals(cpode_mem, &nje);
  CPDlsGetNumFctEvals(cpode_mem, &nfeLS);


  printf(" %6ld   %6ld+%-4ld  %4ld (%3ld)  |  %3ld  %3ld\n",
         nst, nfe, nfeLS, nsetups, nje, ncfn, netf);

  N_VDestroy_Serial(yy);
  N_VDestroy_Serial(yp);

  return;
}

void N_VLinearSum_SensWrapper(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z)
{
  int i;

  for (i=0; i < NV_NVECS_SW(x); i++)
    N_VLinearSum(a, NV_VEC_SW(x,i), b, NV_VEC_SW(y,i), NV_VEC_SW(z,i));

  return;
}

static int IDANewyyp(IDAMem IDA_mem, realtype lambda)
{
  
  /* IDA_YA_YDP_INIT case: ynew  = yy0 - lambda*delta    where id_i = 0
                           ypnew = yp0 - cj*lambda*delta where id_i = 1. */
  if(icopt == IDA_YA_YDP_INIT) {
    N_VProd(id, delta, dtemp);
    N_VLinearSum(ONE, yp0, -cj*lambda, dtemp, ypnew);
    N_VLinearSum(ONE, delta, -ONE, dtemp, dtemp);
    N_VLinearSum(ONE, yy0, -lambda, dtemp, ynew);
    return(IDA_SUCCESS);
  }

  /* IDA_Y_INIT case: ynew = yy0 - lambda*delta. (ypnew = yp0 preset.) */
  N_VLinearSum(ONE, yy0, -lambda, delta, ynew);
  return(IDA_SUCCESS);

}