diff --git a/include/sundials/sundials_config.in b/include/sundials/sundials_config.in
index 732fe0df3..813e7953a 100644
--- a/include/sundials/sundials_config.in
+++ b/include/sundials/sundials_config.in
@@ -17,6 +17,9 @@
 #ifndef _SUNDIALS_CONFIG_H
 #define _SUNDIALS_CONFIG_H
 
+#define NVECTOR_MACRO(fname) N_V##fname##_Serial
+//#define NVECTOR_MACRO(fname) N_V##fname
+
 #include "sundials/sundials_export.h"
 
 #ifndef SUNDIALS_DEPRECATED_MSG
diff --git a/src/cvode/cvode.c b/src/cvode/cvode.c
index 91859026f..df5f9c580 100644
--- a/src/cvode/cvode.c
+++ b/src/cvode/cvode.c
@@ -432,7 +432,7 @@ int CVodeInit(void *cvode_mem, CVRhsFn f, realtype t0, N_Vector y0)
   /* Set space requirements for one N_Vector */
 
   if (y0->ops->nvspace != NULL) {
-    N_VSpace(y0, &lrw1, &liw1);
+    NVECTOR_MACRO(Space)(y0, &lrw1, &liw1);
   } else {
     lrw1 = 0;
     liw1 = 0;
@@ -456,7 +456,7 @@ int CVodeInit(void *cvode_mem, CVRhsFn f, realtype t0, N_Vector y0)
   cv_mem->cv_tn = t0;
 
   /* Initialize zn[0] in the history array */
-  N_VScale(ONE, y0, cv_mem->cv_zn[0]);
+  NVECTOR_MACRO(Scale)(ONE, y0, cv_mem->cv_zn[0]);
 
   /* create a Newton nonlinear solver object by default */
   NLS = SUNNonlinSol_Newton(y0, cv_mem->cv_sunctx);
@@ -611,7 +611,7 @@ int CVodeReInit(void *cvode_mem, realtype t0, N_Vector y0)
 
   /* Initialize zn[0] in the history array */
 
-  N_VScale(ONE, y0, cv_mem->cv_zn[0]);
+  NVECTOR_MACRO(Scale)(ONE, y0, cv_mem->cv_zn[0]);
 
   /* Initialize all the counters */
 
@@ -744,7 +744,7 @@ int CVodeSVtolerances(void *cvode_mem, realtype reltol, N_Vector abstol)
                    "Missing N_VMin routine from N_Vector");
     return(CV_ILL_INPUT);
   }
-  atolmin = N_VMin(abstol);
+  atolmin = NVECTOR_MACRO(Min)(abstol);
   if (atolmin < ZERO) {
     cvProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeSVtolerances",
                    MSGCV_BAD_ABSTOL);
@@ -754,14 +754,14 @@ int CVodeSVtolerances(void *cvode_mem, realtype reltol, N_Vector abstol)
   /* Copy tolerances into memory */
 
   if ( !(cv_mem->cv_VabstolMallocDone) ) {
-    cv_mem->cv_Vabstol = N_VClone(cv_mem->cv_ewt);
+    cv_mem->cv_Vabstol = NVECTOR_MACRO(Clone)(cv_mem->cv_ewt);
     cv_mem->cv_lrw += cv_mem->cv_lrw1;
     cv_mem->cv_liw += cv_mem->cv_liw1;
     cv_mem->cv_VabstolMallocDone = SUNTRUE;
   }
 
   cv_mem->cv_reltol = reltol;
-  N_VScale(ONE, abstol, cv_mem->cv_Vabstol);
+  NVECTOR_MACRO(Scale)(ONE, abstol, cv_mem->cv_Vabstol);
   cv_mem->cv_atolmin0 = (atolmin == ZERO);
 
   cv_mem->cv_itol = CV_SV;
@@ -1137,7 +1137,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     cv_mem->cv_h0u    = cv_mem->cv_h;
     cv_mem->cv_hprime = cv_mem->cv_h;
 
-    N_VScale(cv_mem->cv_h, cv_mem->cv_zn[1], cv_mem->cv_zn[1]);
+    NVECTOR_MACRO(Scale)(cv_mem->cv_h, cv_mem->cv_zn[1], cv_mem->cv_zn[1]);
 
     /* Check for zeros of root function g at and near t0. */
 
@@ -1211,7 +1211,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
           cv_mem->cv_irfnd = 0;
           if ((irfndp == 1) && (itask == CV_ONE_STEP)) {
             cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;
-            N_VScale(ONE, cv_mem->cv_zn[0], yout);
+            NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], yout);
             SUNDIALS_MARK_FUNCTION_END(CV_PROFILER);
             return(CV_SUCCESS);
           }
@@ -1249,7 +1249,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     if ( itask == CV_ONE_STEP &&
          SUNRabs(cv_mem->cv_tn - cv_mem->cv_tretlast) > troundoff ) {
       cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;
-      N_VScale(ONE, cv_mem->cv_zn[0], yout);
+      NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], yout);
       SUNDIALS_MARK_FUNCTION_END(CV_PROFILER);
       return(CV_SUCCESS);
     }
@@ -1319,7 +1319,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
 
         istate = CV_ILL_INPUT;
         cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;
-        N_VScale(ONE, cv_mem->cv_zn[0], yout);
+        NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], yout);
         break;
       }
     }
@@ -1330,19 +1330,19 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
                      MSGCV_MAX_STEPS, cv_mem->cv_tn);
       istate = CV_TOO_MUCH_WORK;
       cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;
-      N_VScale(ONE, cv_mem->cv_zn[0], yout);
+      NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], yout);
       break;
     }
 
     /* Check for too much accuracy requested */
-    nrm = N_VWrmsNorm(cv_mem->cv_zn[0], cv_mem->cv_ewt);
+    nrm = NVECTOR_MACRO(WrmsNorm)(cv_mem->cv_zn[0], cv_mem->cv_ewt);
     cv_mem->cv_tolsf = cv_mem->cv_uround * nrm;
     if (cv_mem->cv_tolsf > ONE) {
       cvProcessError(cv_mem, CV_TOO_MUCH_ACC, "CVODE", "CVode",
                      MSGCV_TOO_MUCH_ACC, cv_mem->cv_tn);
       istate = CV_TOO_MUCH_ACC;
       cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;
-      N_VScale(ONE, cv_mem->cv_zn[0], yout);
+      NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], yout);
       cv_mem->cv_tolsf *= TWO;
       break;
     } else {
@@ -1366,7 +1366,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     if (kflag != CV_SUCCESS) {
       istate = cvHandleFailure(cv_mem, kflag);
       cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;
-      N_VScale(ONE, cv_mem->cv_zn[0], yout);
+      NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], yout);
       break;
     }
 
@@ -1452,7 +1452,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     if (itask == CV_ONE_STEP) {
       istate = CV_SUCCESS;
       cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;
-      N_VScale(ONE, cv_mem->cv_zn[0], yout);
+      NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], yout);
       cv_mem->cv_next_q = cv_mem->cv_qprime;
       cv_mem->cv_next_h = cv_mem->cv_hprime;
       break;
@@ -1542,7 +1542,7 @@ int CVodeGetDky(void *cvode_mem, realtype t, int k, N_Vector dky)
     cv_mem->cv_Xvecs[nvec] = cv_mem->cv_zn[j];
     nvec += 1;
   }
-  ier = N_VLinearCombination(nvec, cv_mem->cv_cvals, cv_mem->cv_Xvecs, dky);
+  ier = NVECTOR_MACRO(LinearCombination)(nvec, cv_mem->cv_cvals, cv_mem->cv_Xvecs, dky);
   if (ier != CV_SUCCESS) {
     SUNDIALS_MARK_FUNCTION_END(CV_PROFILER);
     return(CV_VECTOROP_ERR);
@@ -1553,7 +1553,7 @@ int CVodeGetDky(void *cvode_mem, realtype t, int k, N_Vector dky)
     return(CV_SUCCESS);
   }
   r = SUNRpowerI(cv_mem->cv_h, -k);
-  N_VScale(r, dky, dky);
+  NVECTOR_MACRO(Scale)(r, dky, dky);
 
   SUNDIALS_MARK_FUNCTION_END(CV_PROFILER);
   return(CV_SUCCESS);
@@ -1579,7 +1579,7 @@ int CVodeComputeState(void *cvode_mem, N_Vector ycor, N_Vector y)
 
   SUNDIALS_MARK_FUNCTION_BEGIN(CV_PROFILER);
 
-  N_VLinearSum(ONE, cv_mem->cv_zn[0], ONE, ycor, y);
+  NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[0], ONE, ycor, y);
 
   SUNDIALS_MARK_FUNCTION_END(CV_PROFILER);
   return(CV_SUCCESS);
@@ -1681,73 +1681,73 @@ static booleantype cvAllocVectors(CVodeMem cv_mem, N_Vector tmpl)
 
   /* Allocate ewt, acor, tempv, ftemp */
 
-  cv_mem->cv_ewt = N_VClone(tmpl);
+  cv_mem->cv_ewt = NVECTOR_MACRO(Clone)(tmpl);
   if (cv_mem->cv_ewt == NULL) return(SUNFALSE);
 
-  cv_mem->cv_acor = N_VClone(tmpl);
+  cv_mem->cv_acor = NVECTOR_MACRO(Clone)(tmpl);
   if (cv_mem->cv_acor == NULL) {
-    N_VDestroy(cv_mem->cv_ewt);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
     return(SUNFALSE);
   }
 
-  cv_mem->cv_tempv = N_VClone(tmpl);
+  cv_mem->cv_tempv = NVECTOR_MACRO(Clone)(tmpl);
   if (cv_mem->cv_tempv == NULL) {
-    N_VDestroy(cv_mem->cv_ewt);
-    N_VDestroy(cv_mem->cv_acor);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_acor);
     return(SUNFALSE);
   }
 
-  cv_mem->cv_ftemp = N_VClone(tmpl);
+  cv_mem->cv_ftemp = NVECTOR_MACRO(Clone)(tmpl);
   if (cv_mem->cv_ftemp == NULL) {
-    N_VDestroy(cv_mem->cv_tempv);
-    N_VDestroy(cv_mem->cv_ewt);
-    N_VDestroy(cv_mem->cv_acor);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_tempv);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_acor);
     return(SUNFALSE);
   }
 
-  cv_mem->cv_vtemp1 = N_VClone(tmpl);
+  cv_mem->cv_vtemp1 = NVECTOR_MACRO(Clone)(tmpl);
   if (cv_mem->cv_vtemp1 == NULL) {
-    N_VDestroy(cv_mem->cv_ftemp);
-    N_VDestroy(cv_mem->cv_tempv);
-    N_VDestroy(cv_mem->cv_ewt);
-    N_VDestroy(cv_mem->cv_acor);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ftemp);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_tempv);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_acor);
     return(SUNFALSE);
   }
 
-  cv_mem->cv_vtemp2 = N_VClone(tmpl);
+  cv_mem->cv_vtemp2 = NVECTOR_MACRO(Clone)(tmpl);
   if (cv_mem->cv_vtemp2 == NULL) {
-    N_VDestroy(cv_mem->cv_vtemp1);
-    N_VDestroy(cv_mem->cv_ftemp);
-    N_VDestroy(cv_mem->cv_tempv);
-    N_VDestroy(cv_mem->cv_ewt);
-    N_VDestroy(cv_mem->cv_acor);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp1);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ftemp);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_tempv);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_acor);
     return(SUNFALSE);
   }
 
-  cv_mem->cv_vtemp3 = N_VClone(tmpl);
+  cv_mem->cv_vtemp3 = NVECTOR_MACRO(Clone)(tmpl);
   if (cv_mem->cv_vtemp3 == NULL) {
-    N_VDestroy(cv_mem->cv_vtemp2);
-    N_VDestroy(cv_mem->cv_vtemp1);
-    N_VDestroy(cv_mem->cv_ftemp);
-    N_VDestroy(cv_mem->cv_tempv);
-    N_VDestroy(cv_mem->cv_ewt);
-    N_VDestroy(cv_mem->cv_acor);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp2);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp1);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ftemp);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_tempv);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_acor);
     return(SUNFALSE);
   }
 
   /* Allocate zn[0] ... zn[qmax] */
 
   for (j=0; j <= cv_mem->cv_qmax; j++) {
-    cv_mem->cv_zn[j] = N_VClone(tmpl);
+    cv_mem->cv_zn[j] = NVECTOR_MACRO(Clone)(tmpl);
     if (cv_mem->cv_zn[j] == NULL) {
-      N_VDestroy(cv_mem->cv_ewt);
-      N_VDestroy(cv_mem->cv_acor);
-      N_VDestroy(cv_mem->cv_tempv);
-      N_VDestroy(cv_mem->cv_ftemp);
-      N_VDestroy(cv_mem->cv_vtemp1);
-      N_VDestroy(cv_mem->cv_vtemp2);
-      N_VDestroy(cv_mem->cv_vtemp3);
-      for (i=0; i < j; i++) N_VDestroy(cv_mem->cv_zn[i]);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_acor);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_tempv);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_ftemp);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp1);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp2);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp3);
+      for (i=0; i < j; i++) NVECTOR_MACRO(Destroy)(cv_mem->cv_zn[i]);
       return(SUNFALSE);
     }
   }
@@ -1774,26 +1774,26 @@ static void cvFreeVectors(CVodeMem cv_mem)
 
   maxord = cv_mem->cv_qmax_alloc;
 
-  N_VDestroy(cv_mem->cv_ewt);
-  N_VDestroy(cv_mem->cv_acor);
-  N_VDestroy(cv_mem->cv_tempv);
-  N_VDestroy(cv_mem->cv_ftemp);
-  N_VDestroy(cv_mem->cv_vtemp1);
-  N_VDestroy(cv_mem->cv_vtemp2);
-  N_VDestroy(cv_mem->cv_vtemp3);
-  for (j=0; j <= maxord; j++) N_VDestroy(cv_mem->cv_zn[j]);
+  NVECTOR_MACRO(Destroy)(cv_mem->cv_ewt);
+  NVECTOR_MACRO(Destroy)(cv_mem->cv_acor);
+  NVECTOR_MACRO(Destroy)(cv_mem->cv_tempv);
+  NVECTOR_MACRO(Destroy)(cv_mem->cv_ftemp);
+  NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp1);
+  NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp2);
+  NVECTOR_MACRO(Destroy)(cv_mem->cv_vtemp3);
+  for (j=0; j <= maxord; j++) NVECTOR_MACRO(Destroy)(cv_mem->cv_zn[j]);
 
   cv_mem->cv_lrw -= (maxord + 8)*cv_mem->cv_lrw1;
   cv_mem->cv_liw -= (maxord + 8)*cv_mem->cv_liw1;
 
   if (cv_mem->cv_VabstolMallocDone) {
-    N_VDestroy(cv_mem->cv_Vabstol);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_Vabstol);
     cv_mem->cv_lrw -= cv_mem->cv_lrw1;
     cv_mem->cv_liw -= cv_mem->cv_liw1;
   }
 
   if (cv_mem->cv_constraintsMallocDone) {
-    N_VDestroy(cv_mem->cv_constraints);
+    NVECTOR_MACRO(Destroy)(cv_mem->cv_constraints);
     cv_mem->cv_lrw -= cv_mem->cv_lrw1;
     cv_mem->cv_liw -= cv_mem->cv_liw1;
   }
@@ -1841,7 +1841,7 @@ static int cvInitialSetup(CVodeMem cv_mem)
 
   /* Check to see if y0 satisfies constraints */
   if (cv_mem->cv_constraintsSet) {
-    conOK = N_VConstrMask(cv_mem->cv_constraints, cv_mem->cv_zn[0], cv_mem->cv_tempv);
+    conOK = NVECTOR_MACRO(ConstrMask)(cv_mem->cv_constraints, cv_mem->cv_zn[0], cv_mem->cv_tempv);
     if (!conOK) {
       cvProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "cvInitialSetup",
                      MSGCV_Y0_FAIL_CONSTR);
@@ -2065,15 +2065,15 @@ static realtype cvUpperBoundH0(CVodeMem cv_mem, realtype tdist)
   temp1 = cv_mem->cv_tempv;
   temp2 = cv_mem->cv_acor;
 
-  N_VAbs(cv_mem->cv_zn[0], temp2);
+  NVECTOR_MACRO(Abs)(cv_mem->cv_zn[0], temp2);
   cv_mem->cv_efun(cv_mem->cv_zn[0], temp1, cv_mem->cv_e_data);
-  N_VInv(temp1, temp1);
-  N_VLinearSum(HUB_FACTOR, temp2, ONE, temp1, temp1);
+  NVECTOR_MACRO(Inv)(temp1, temp1);
+  NVECTOR_MACRO(LinearSum)(HUB_FACTOR, temp2, ONE, temp1, temp1);
 
-  N_VAbs(cv_mem->cv_zn[1], temp2);
+  NVECTOR_MACRO(Abs)(cv_mem->cv_zn[1], temp2);
 
-  N_VDiv(temp2, temp1, temp1);
-  hub_inv = N_VMaxNorm(temp1);
+  NVECTOR_MACRO(Div)(temp2, temp1, temp1);
+  hub_inv = NVECTOR_MACRO(MaxNorm)(temp1);
 
   /*
    * bound based on tdist -- allow at most a step of magnitude
@@ -2100,16 +2100,16 @@ static int cvYddNorm(CVodeMem cv_mem, realtype hg, realtype *yddnrm)
 {
   int retval;
 
-  N_VLinearSum(hg, cv_mem->cv_zn[1], ONE, cv_mem->cv_zn[0], cv_mem->cv_y);
+  NVECTOR_MACRO(LinearSum)(hg, cv_mem->cv_zn[1], ONE, cv_mem->cv_zn[0], cv_mem->cv_y);
   retval = cv_mem->cv_f(cv_mem->cv_tn+hg, cv_mem->cv_y,
                         cv_mem->cv_tempv, cv_mem->cv_user_data);
   cv_mem->cv_nfe++;
   if (retval < 0) return(CV_RHSFUNC_FAIL);
   if (retval > 0) return(RHSFUNC_RECVR);
 
-  N_VLinearSum(ONE/hg, cv_mem->cv_tempv, -ONE/hg, cv_mem->cv_zn[1], cv_mem->cv_tempv);
+  NVECTOR_MACRO(LinearSum)(ONE/hg, cv_mem->cv_tempv, -ONE/hg, cv_mem->cv_zn[1], cv_mem->cv_tempv);
 
-  *yddnrm = N_VWrmsNorm(cv_mem->cv_tempv, cv_mem->cv_ewt);
+  *yddnrm = NVECTOR_MACRO(WrmsNorm)(cv_mem->cv_tempv, cv_mem->cv_ewt);
 
   return(CV_SUCCESS);
 }
@@ -2241,7 +2241,7 @@ static int cvStep(CVodeMem cv_mem)
   /*  Finally, we rescale the acor array to be the
       estimated local error vector. */
 
-  N_VScale(cv_mem->cv_tq[2], cv_mem->cv_acor, cv_mem->cv_acor);
+  NVECTOR_MACRO(Scale)(cv_mem->cv_tq[2], cv_mem->cv_acor, cv_mem->cv_acor);
 
   return(CV_SUCCESS);
 }
@@ -2312,7 +2312,7 @@ static void cvAdjustAdams(CVodeMem cv_mem, int deltaq)
   /* On an order increase, set new column of zn to zero and return */
 
   if (deltaq==1) {
-    N_VConst(ZERO, cv_mem->cv_zn[cv_mem->cv_L]);
+    NVECTOR_MACRO(Const)(ZERO, cv_mem->cv_zn[cv_mem->cv_L]);
     return;
   }
 
@@ -2342,7 +2342,7 @@ static void cvAdjustAdams(CVodeMem cv_mem, int deltaq)
     cv_mem->cv_cvals[j-2] = -cv_mem->cv_l[j];
 
   if (cv_mem->cv_q > 2)
-    (void) N_VScaleAddMulti(cv_mem->cv_q-2, cv_mem->cv_cvals,
+    (void) NVECTOR_MACRO(ScaleAddMulti)(cv_mem->cv_q-2, cv_mem->cv_cvals,
                             cv_mem->cv_zn[cv_mem->cv_q],
                             cv_mem->cv_zn+2, cv_mem->cv_zn+2);
 }
@@ -2402,12 +2402,12 @@ static void cvIncreaseBDF(CVodeMem cv_mem)
     }
   }
   A1 = (-alpha0 - alpha1) / prod;
-  N_VScale(A1, cv_mem->cv_zn[cv_mem->cv_indx_acor],
+  NVECTOR_MACRO(Scale)(A1, cv_mem->cv_zn[cv_mem->cv_indx_acor],
            cv_mem->cv_zn[cv_mem->cv_L]);
 
   /* for (j=2; j <= cv_mem->cv_q; j++) */
   if (cv_mem->cv_q > 1)
-    (void) N_VScaleAddMulti(cv_mem->cv_q-1, cv_mem->cv_l+2,
+    (void) NVECTOR_MACRO(ScaleAddMulti)(cv_mem->cv_q-1, cv_mem->cv_l+2,
                             cv_mem->cv_zn[cv_mem->cv_L],
                             cv_mem->cv_zn+2, cv_mem->cv_zn+2);
 }
@@ -2441,7 +2441,7 @@ static void cvDecreaseBDF(CVodeMem cv_mem)
     cv_mem->cv_cvals[j-2] = -cv_mem->cv_l[j];
 
   if (cv_mem->cv_q > 2)
-    (void) N_VScaleAddMulti(cv_mem->cv_q-2, cv_mem->cv_cvals,
+    (void) NVECTOR_MACRO(ScaleAddMulti)(cv_mem->cv_q-2, cv_mem->cv_cvals,
                             cv_mem->cv_zn[cv_mem->cv_q],
                             cv_mem->cv_zn+2, cv_mem->cv_zn+2);
 }
@@ -2463,7 +2463,7 @@ void cvRescale(CVodeMem cv_mem)
   for (j=1; j <= cv_mem->cv_q; j++)
     cv_mem->cv_cvals[j] = cv_mem->cv_eta * cv_mem->cv_cvals[j-1];
 
-  (void) N_VScaleVectorArray(cv_mem->cv_q, cv_mem->cv_cvals,
+  (void) NVECTOR_MACRO(ScaleVectorArray)(cv_mem->cv_q, cv_mem->cv_cvals,
                              cv_mem->cv_zn+1, cv_mem->cv_zn+1);
 
   cv_mem->cv_h = cv_mem->cv_hscale * cv_mem->cv_eta;
@@ -2494,13 +2494,13 @@ static void cvPredict(CVodeMem cv_mem)
 
   for (k = 1; k <= cv_mem->cv_q; k++)
     for (j = cv_mem->cv_q; j >= k; j--)
-      N_VLinearSum(ONE, cv_mem->cv_zn[j-1], ONE,
+      NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[j-1], ONE,
                    cv_mem->cv_zn[j], cv_mem->cv_zn[j-1]);
 
 #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG
   SUNLogger_QueueMsg(CV_LOGGER, SUN_LOGLEVEL_DEBUG,
     "CVODE::cvPredict", "return", "predictor =", "");
-  N_VPrintFile(cv_mem->cv_zn[0], CV_LOGGER->debug_fp);
+  NVECTOR_MACRO(PrintFile)(cv_mem->cv_zn[0], CV_LOGGER->debug_fp);
 #endif
 }
 
@@ -2801,7 +2801,7 @@ static int cvNls(CVodeMem cv_mem, int nflag)
   }
 
   /* initial guess for the correction to the predictor */
-  N_VConst(ZERO, cv_mem->cv_acor);
+  NVECTOR_MACRO(Const)(ZERO, cv_mem->cv_acor);
 
   /* call nonlinear solver setup if it exists */
   if ((cv_mem->NLS)->ops->setup) {
@@ -2827,11 +2827,11 @@ static int cvNls(CVodeMem cv_mem, int nflag)
   /* solve successful */
 
   /* update the state based on the final correction from the nonlinear solver */
-  N_VLinearSum(ONE, cv_mem->cv_zn[0], ONE, cv_mem->cv_acor, cv_mem->cv_y);
+  NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[0], ONE, cv_mem->cv_acor, cv_mem->cv_y);
 
   /* compute acnrm if is was not already done by the nonlinear solver */
   if (!cv_mem->cv_acnrmcur)
-    cv_mem->cv_acnrm = N_VWrmsNorm(cv_mem->cv_acor, cv_mem->cv_ewt);
+    cv_mem->cv_acnrm = NVECTOR_MACRO(WrmsNorm)(cv_mem->cv_acor, cv_mem->cv_ewt);
 
   /* update Jacobian status */
   cv_mem->cv_jcur = SUNFALSE;
@@ -2866,7 +2866,7 @@ static int cvCheckConstraints(CVodeMem cv_mem)
   N_Vector tmp = cv_mem->cv_tempv;
 
   /* Get mask vector mm, set where constraints failed */
-  constraintsPassed = N_VConstrMask(cv_mem->cv_constraints, cv_mem->cv_y, mm);
+  constraintsPassed = NVECTOR_MACRO(ConstrMask)(cv_mem->cv_constraints, cv_mem->cv_y, mm);
   if (constraintsPassed) return(CV_SUCCESS);
 
   /* Constraints not met */
@@ -2884,19 +2884,19 @@ static int cvCheckConstraints(CVodeMem cv_mem)
   else
 #endif
   {
-    N_VCompare(ONEPT5, cv_mem->cv_constraints, tmp); /* a[i]=1 when |c[i]|=2  */
-    N_VProd(tmp, cv_mem->cv_constraints, tmp);       /* a * c                 */
-    N_VDiv(tmp, cv_mem->cv_ewt, tmp);                /* a * c * wt            */
-    N_VLinearSum(ONE, cv_mem->cv_y, -PT1, tmp, tmp); /* y - 0.1 * a * c * wt  */
-    N_VProd(tmp, mm, tmp);                           /* v = mm*(y-0.1*a*c*wt) */
+    NVECTOR_MACRO(Compare)(ONEPT5, cv_mem->cv_constraints, tmp); /* a[i]=1 when |c[i]|=2  */
+    NVECTOR_MACRO(Prod)(tmp, cv_mem->cv_constraints, tmp);       /* a * c                 */
+    NVECTOR_MACRO(Div)(tmp, cv_mem->cv_ewt, tmp);                /* a * c * wt            */
+    NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_y, -PT1, tmp, tmp); /* y - 0.1 * a * c * wt  */
+    NVECTOR_MACRO(Prod)(tmp, mm, tmp);                           /* v = mm*(y-0.1*a*c*wt) */
   }
 
-  vnorm = N_VWrmsNorm(tmp, cv_mem->cv_ewt);        /* ||v|| */
+  vnorm = NVECTOR_MACRO(WrmsNorm)(tmp, cv_mem->cv_ewt);        /* ||v|| */
 
   /* If vector v of constraint corrections is small in norm, correct and
      accept this step */
   if (vnorm <= cv_mem->cv_tq[4]) {
-    N_VLinearSum(ONE, cv_mem->cv_acor,
+    NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_acor,
                  -ONE, tmp, cv_mem->cv_acor);      /* acor <- acor - v */
     return(CV_SUCCESS);
   }
@@ -2905,9 +2905,9 @@ static int cvCheckConstraints(CVodeMem cv_mem)
   if (SUNRabs(cv_mem->cv_h) <= cv_mem->cv_hmin*ONEPSM) return(CV_CONSTR_FAIL);
 
   /* Constraint correction is too large, reduce h by computing eta = h'/h */
-  N_VLinearSum(ONE, cv_mem->cv_zn[0], -ONE, cv_mem->cv_y, tmp);
-  N_VProd(mm, tmp, tmp);
-  cv_mem->cv_eta = PT9*N_VMinQuotient(cv_mem->cv_zn[0], tmp);
+  NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[0], -ONE, cv_mem->cv_y, tmp);
+  NVECTOR_MACRO(Prod)(mm, tmp, tmp);
+  cv_mem->cv_eta = PT9*NVECTOR_MACRO(MinQuotient)(cv_mem->cv_zn[0], tmp);
   cv_mem->cv_eta = SUNMAX(cv_mem->cv_eta, PT1);
   cv_mem->cv_eta = SUNMAX(cv_mem->cv_eta,
                           cv_mem->cv_hmin / SUNRabs(cv_mem->cv_h));
@@ -3012,7 +3012,7 @@ void cvRestore(CVodeMem cv_mem, realtype saved_t)
   cv_mem->cv_tn = saved_t;
   for (k = 1; k <= cv_mem->cv_q; k++)
     for (j = cv_mem->cv_q; j >= k; j--)
-      N_VLinearSum(ONE, cv_mem->cv_zn[j-1], -ONE,
+      NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[j-1], -ONE,
                    cv_mem->cv_zn[j], cv_mem->cv_zn[j-1]);
 }
 
@@ -3130,7 +3130,7 @@ static int cvDoErrorTest(CVodeMem cv_mem, int *nflagPtr, realtype saved_t,
   if (retval < 0) return(CV_RHSFUNC_FAIL);
   if (retval > 0) return(CV_UNREC_RHSFUNC_ERR);
 
-  N_VScale(cv_mem->cv_h, cv_mem->cv_tempv, cv_mem->cv_zn[1]);
+  NVECTOR_MACRO(Scale)(cv_mem->cv_h, cv_mem->cv_tempv, cv_mem->cv_zn[1]);
 
 #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO
   SUNLogger_QueueMsg(CV_LOGGER, SUN_LOGLEVEL_INFO,
@@ -3174,19 +3174,19 @@ static void cvCompleteStep(CVodeMem cv_mem)
   cv_mem->cv_tau[1] = cv_mem->cv_h;
 
   /* Apply correction to column j of zn: l_j * Delta_n */
-  (void) N_VScaleAddMulti(cv_mem->cv_q+1, cv_mem->cv_l, cv_mem->cv_acor,
+  (void) NVECTOR_MACRO(ScaleAddMulti)(cv_mem->cv_q+1, cv_mem->cv_l, cv_mem->cv_acor,
                           cv_mem->cv_zn, cv_mem->cv_zn);
 
   /* Apply the projection correction to column j of zn: p_j * Delta_n */
   if (cv_mem->proj_applied) {
-    (void) N_VScaleAddMulti(cv_mem->cv_q+1,
+    (void) NVECTOR_MACRO(ScaleAddMulti)(cv_mem->cv_q+1,
                             cv_mem->proj_p, cv_mem->cv_tempv, /* tempv = acorP */
                             cv_mem->cv_zn, cv_mem->cv_zn);
   }
 
   cv_mem->cv_qwait--;
   if ((cv_mem->cv_qwait == 1) && (cv_mem->cv_q != cv_mem->cv_qmax)) {
-    N_VScale(ONE, cv_mem->cv_acor, cv_mem->cv_zn[cv_mem->cv_qmax]);
+    NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_acor, cv_mem->cv_zn[cv_mem->cv_qmax]);
     cv_mem->cv_saved_tq5 = cv_mem->cv_tq[5];
     cv_mem->cv_indx_acor = cv_mem->cv_qmax;
   }
@@ -3302,7 +3302,7 @@ static realtype cvComputeEtaqm1(CVodeMem cv_mem)
 
   cv_mem->cv_etaqm1 = ZERO;
   if (cv_mem->cv_q > 1) {
-    ddn = N_VWrmsNorm(cv_mem->cv_zn[cv_mem->cv_q], cv_mem->cv_ewt) * cv_mem->cv_tq[1];
+    ddn = NVECTOR_MACRO(WrmsNorm)(cv_mem->cv_zn[cv_mem->cv_q], cv_mem->cv_ewt) * cv_mem->cv_tq[1];
     cv_mem->cv_etaqm1 = ONE/(SUNRpowerR(BIAS1*ddn, ONE/cv_mem->cv_q) + ADDON);
   }
   return(cv_mem->cv_etaqm1);
@@ -3324,9 +3324,9 @@ static realtype cvComputeEtaqp1(CVodeMem cv_mem)
     if (cv_mem->cv_saved_tq5 == ZERO) return(cv_mem->cv_etaqp1);
     cquot = (cv_mem->cv_tq[5] / cv_mem->cv_saved_tq5) *
       SUNRpowerI(cv_mem->cv_h/cv_mem->cv_tau[2], cv_mem->cv_L);
-    N_VLinearSum(-cquot, cv_mem->cv_zn[cv_mem->cv_qmax], ONE,
+    NVECTOR_MACRO(LinearSum)(-cquot, cv_mem->cv_zn[cv_mem->cv_qmax], ONE,
                  cv_mem->cv_acor, cv_mem->cv_tempv);
-    dup = N_VWrmsNorm(cv_mem->cv_tempv, cv_mem->cv_ewt) * cv_mem->cv_tq[3];
+    dup = NVECTOR_MACRO(WrmsNorm)(cv_mem->cv_tempv, cv_mem->cv_ewt) * cv_mem->cv_tq[3];
     cv_mem->cv_etaqp1 = ONE / (SUNRpowerR(BIAS3*dup, ONE/(cv_mem->cv_L+1)) + ADDON);
   }
   return(cv_mem->cv_etaqp1);
@@ -3381,7 +3381,7 @@ static void cvChooseEta(CVodeMem cv_mem)
          * to order q+1, so it represents Delta_n in the ELTE at q+1
          */
 
-        N_VScale(ONE, cv_mem->cv_acor, cv_mem->cv_zn[cv_mem->cv_qmax]);
+        NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_acor, cv_mem->cv_zn[cv_mem->cv_qmax]);
       }
     }
   }
@@ -3406,8 +3406,8 @@ static int cvHandleFailure(CVodeMem cv_mem, int flag)
 
   /* Set vector of  absolute weighted local errors */
   /*
-  N_VProd(acor, ewt, tempv);
-  N_VAbs(tempv, tempv);
+  NVECTOR_MACRO(Prod)(acor, ewt, tempv);
+  NVECTOR_MACRO(Abs)(tempv, tempv);
   */
 
   /* Depending on flag, print error message and return error flag */
@@ -3525,9 +3525,9 @@ static void cvBDFStab(CVodeMem cv_mem)
     sq = factorial * cv_mem->cv_q * (cv_mem->cv_q+1) *
       cv_mem->cv_acnrm / SUNMAX(cv_mem->cv_tq[5],TINY);
     sqm1 = factorial * cv_mem->cv_q *
-      N_VWrmsNorm(cv_mem->cv_zn[cv_mem->cv_q], cv_mem->cv_ewt);
+      NVECTOR_MACRO(WrmsNorm)(cv_mem->cv_zn[cv_mem->cv_q], cv_mem->cv_ewt);
     sqm2 = factorial *
-      N_VWrmsNorm(cv_mem->cv_zn[cv_mem->cv_q-1], cv_mem->cv_ewt);
+      NVECTOR_MACRO(WrmsNorm)(cv_mem->cv_zn[cv_mem->cv_q-1], cv_mem->cv_ewt);
     cv_mem->cv_ssdat[1][1] = sqm2*sqm2;
     cv_mem->cv_ssdat[1][2] = sqm1*sqm1;
     cv_mem->cv_ssdat[1][3] = sq*sq;
@@ -3919,7 +3919,7 @@ static int cvRcheck1(CVodeMem cv_mem)
   hratio = SUNMAX(cv_mem->cv_ttol/SUNRabs(cv_mem->cv_h), PT1);
   smallh = hratio*cv_mem->cv_h;
   tplus = cv_mem->cv_tlo + smallh;
-  N_VLinearSum(ONE, cv_mem->cv_zn[0], hratio, cv_mem->cv_zn[1], cv_mem->cv_y);
+  NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[0], hratio, cv_mem->cv_zn[1], cv_mem->cv_y);
   retval = cv_mem->cv_gfun(tplus, cv_mem->cv_y,
                            cv_mem->cv_ghi, cv_mem->cv_user_data);
   cv_mem->cv_nge++;
@@ -3989,7 +3989,7 @@ static int cvRcheck2(CVodeMem cv_mem)
   tplus = cv_mem->cv_tlo + smallh;
   if ( (tplus - cv_mem->cv_tn)*cv_mem->cv_h >= ZERO) {
     hratio = smallh/cv_mem->cv_h;
-    N_VLinearSum(ONE, cv_mem->cv_y, hratio, cv_mem->cv_zn[1], cv_mem->cv_y);
+    NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_y, hratio, cv_mem->cv_zn[1], cv_mem->cv_y);
   } else {
     (void) CVodeGetDky(cv_mem, tplus, 0, cv_mem->cv_y);
   }
@@ -4036,12 +4036,12 @@ static int cvRcheck3(CVodeMem cv_mem)
   /* Set thi = tn or tout, whichever comes first; set y = y(thi). */
   if (cv_mem->cv_taskc == CV_ONE_STEP) {
     cv_mem->cv_thi = cv_mem->cv_tn;
-    N_VScale(ONE, cv_mem->cv_zn[0], cv_mem->cv_y);
+    NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], cv_mem->cv_y);
   }
   if (cv_mem->cv_taskc == CV_NORMAL) {
     if ( (cv_mem->cv_toutc - cv_mem->cv_tn)*cv_mem->cv_h >= ZERO) {
       cv_mem->cv_thi = cv_mem->cv_tn;
-      N_VScale(ONE, cv_mem->cv_zn[0], cv_mem->cv_y);
+      NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_zn[0], cv_mem->cv_y);
     } else {
       cv_mem->cv_thi = cv_mem->cv_toutc;
       (void) CVodeGetDky(cv_mem, cv_mem->cv_thi, 0, cv_mem->cv_y);
@@ -4382,19 +4382,19 @@ static int cvEwtSetSS(CVodeMem cv_mem, N_Vector ycur, N_Vector weight)
                      cv_mem->cv_Sabstol, ycur, cv_mem->cv_tempv,
                      weight);
     if (cv_mem->cv_atolmin0) {
-      if (N_VMin(cv_mem->cv_tempv) <= ZERO) return(-1);
+      if (NVECTOR_MACRO(Min)(cv_mem->cv_tempv) <= ZERO) return(-1);
     }
   }
   else
 #endif
   {
-    N_VAbs(ycur, cv_mem->cv_tempv);
-    N_VScale(cv_mem->cv_reltol, cv_mem->cv_tempv, cv_mem->cv_tempv);
-    N_VAddConst(cv_mem->cv_tempv, cv_mem->cv_Sabstol, cv_mem->cv_tempv);
+    NVECTOR_MACRO(Abs)(ycur, cv_mem->cv_tempv);
+    NVECTOR_MACRO(Scale)(cv_mem->cv_reltol, cv_mem->cv_tempv, cv_mem->cv_tempv);
+    NVECTOR_MACRO(AddConst)(cv_mem->cv_tempv, cv_mem->cv_Sabstol, cv_mem->cv_tempv);
     if (cv_mem->cv_atolmin0) {
-      if (N_VMin(cv_mem->cv_tempv) <= ZERO) return(-1);
+      if (NVECTOR_MACRO(Min)(cv_mem->cv_tempv) <= ZERO) return(-1);
     }
-    N_VInv(cv_mem->cv_tempv, weight);
+    NVECTOR_MACRO(Inv)(cv_mem->cv_tempv, weight);
   }
 
   return(0);
@@ -4421,19 +4421,19 @@ static int cvEwtSetSV(CVodeMem cv_mem, N_Vector ycur, N_Vector weight)
                      cv_mem->cv_Vabstol, ycur, cv_mem->cv_tempv,
                      weight);
     if (cv_mem->cv_atolmin0) {
-      if (N_VMin(cv_mem->cv_tempv) <= ZERO) return(-1);
+      if (NVECTOR_MACRO(Min)(cv_mem->cv_tempv) <= ZERO) return(-1);
     }
   }
   else
 #endif
   {
-    N_VAbs(ycur, cv_mem->cv_tempv);
-    N_VLinearSum(cv_mem->cv_reltol, cv_mem->cv_tempv, ONE,
+    NVECTOR_MACRO(Abs)(ycur, cv_mem->cv_tempv);
+    NVECTOR_MACRO(LinearSum)(cv_mem->cv_reltol, cv_mem->cv_tempv, ONE,
                  cv_mem->cv_Vabstol, cv_mem->cv_tempv);
     if (cv_mem->cv_atolmin0) {
-      if (N_VMin(cv_mem->cv_tempv) <= ZERO) return(-1);
+      if (NVECTOR_MACRO(Min)(cv_mem->cv_tempv) <= ZERO) return(-1);
     }
-    N_VInv(cv_mem->cv_tempv, weight);
+    NVECTOR_MACRO(Inv)(cv_mem->cv_tempv, weight);
   }
 
   return(0);
diff --git a/src/cvode/cvode_bandpre.c b/src/cvode/cvode_bandpre.c
index d48ba5c31..f3c6ecff2 100644
--- a/src/cvode/cvode_bandpre.c
+++ b/src/cvode/cvode_bandpre.c
@@ -143,7 +143,7 @@ int CVBandPrecInit(void *cvode_mem, sunindextype N,
 
   /* allocate memory for temporary N_Vectors */
   pdata->tmp1 = NULL;
-  pdata->tmp1 = N_VClone(cv_mem->cv_tempv);
+  pdata->tmp1 = NVECTOR_MACRO(Clone)(cv_mem->cv_tempv);
   if (pdata->tmp1 == NULL) {
     SUNLinSolFree(pdata->LS);
     SUNMatDestroy(pdata->savedP);
@@ -154,12 +154,12 @@ int CVBandPrecInit(void *cvode_mem, sunindextype N,
     return(CVLS_MEM_FAIL);
   }
   pdata->tmp2 = NULL;
-  pdata->tmp2 = N_VClone(cv_mem->cv_tempv);
+  pdata->tmp2 = NVECTOR_MACRO(Clone)(cv_mem->cv_tempv);
   if (pdata->tmp2 == NULL) {
     SUNLinSolFree(pdata->LS);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
-    N_VDestroy(pdata->tmp1);
+    NVECTOR_MACRO(Destroy)(pdata->tmp1);
     free(pdata); pdata = NULL;
     cvProcessError(cv_mem, CVLS_MEM_FAIL, "CVBANDPRE",
                     "CVBandPrecInit", MSGBP_MEM_FAIL);
@@ -172,8 +172,8 @@ int CVBandPrecInit(void *cvode_mem, sunindextype N,
     SUNLinSolFree(pdata->LS);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
-    N_VDestroy(pdata->tmp1);
-    N_VDestroy(pdata->tmp2);
+    NVECTOR_MACRO(Destroy)(pdata->tmp1);
+    NVECTOR_MACRO(Destroy)(pdata->tmp2);
     free(pdata); pdata = NULL;
     cvProcessError(cv_mem, CVLS_SUNLS_FAIL, "CVBANDPRE",
                     "CVBandPrecInit", MSGBP_SUNLS_FAIL);
@@ -233,7 +233,7 @@ int CVBandPrecGetWorkSpace(void *cvode_mem, long int *lenrwBP,
   *leniwBP = 4;
   *lenrwBP = 0;
   if (cv_mem->cv_tempv->ops->nvspace) {
-    N_VSpace(cv_mem->cv_tempv, &lrw1, &liw1);
+    NVECTOR_MACRO(Space)(cv_mem->cv_tempv, &lrw1, &liw1);
     *leniwBP += 2*liw1;
     *lenrwBP += 2*lrw1;
   }
@@ -458,8 +458,8 @@ static int CVBandPrecFree(CVodeMem cv_mem)
   SUNLinSolFree(pdata->LS);
   SUNMatDestroy(pdata->savedP);
   SUNMatDestroy(pdata->savedJ);
-  N_VDestroy(pdata->tmp1);
-  N_VDestroy(pdata->tmp2);
+  NVECTOR_MACRO(Destroy)(pdata->tmp1);
+  NVECTOR_MACRO(Destroy)(pdata->tmp2);
 
   free(pdata);
   pdata = NULL;
@@ -495,20 +495,20 @@ static int CVBandPDQJac(CVBandPrecData pdata, realtype t, N_Vector y,
   cv_mem = (CVodeMem) pdata->cvode_mem;
 
   /* Obtain pointers to the data for various vectors */
-  ewt_data   = N_VGetArrayPointer(cv_mem->cv_ewt);
-  fy_data    = N_VGetArrayPointer(fy);
-  ftemp_data = N_VGetArrayPointer(ftemp);
-  y_data     = N_VGetArrayPointer(y);
-  ytemp_data = N_VGetArrayPointer(ytemp);
+  ewt_data   = NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_ewt);
+  fy_data    = NVECTOR_MACRO(GetArrayPointer)(fy);
+  ftemp_data = NVECTOR_MACRO(GetArrayPointer)(ftemp);
+  y_data     = NVECTOR_MACRO(GetArrayPointer)(y);
+  ytemp_data = NVECTOR_MACRO(GetArrayPointer)(ytemp);
   if (cv_mem->cv_constraintsSet)
-    cns_data = N_VGetArrayPointer(cv_mem->cv_constraints);
+    cns_data = NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_constraints);
 
   /* Load ytemp with y = predicted y vector. */
-  N_VScale(ONE, y, ytemp);
+  NVECTOR_MACRO(Scale)(ONE, y, ytemp);
 
   /* Set minimum increment based on uround and norm of f. */
   srur = SUNRsqrt(cv_mem->cv_uround);
-  fnorm = N_VWrmsNorm(fy, cv_mem->cv_ewt);
+  fnorm = NVECTOR_MACRO(WrmsNorm)(fy, cv_mem->cv_ewt);
   minInc = (fnorm != ZERO) ?
     (MIN_INC_MULT * SUNRabs(cv_mem->cv_h) * cv_mem->cv_uround * pdata->N * fnorm) : ONE;
 
diff --git a/src/cvode/cvode_bbdpre.c b/src/cvode/cvode_bbdpre.c
index f376fe6e1..967fc0d51 100644
--- a/src/cvode/cvode_bbdpre.c
+++ b/src/cvode/cvode_bbdpre.c
@@ -144,7 +144,7 @@ int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal,
   pdata->rlocal = NULL;
   pdata->rlocal = N_VNewEmpty_Serial(Nlocal, cv_mem->cv_sunctx);
   if (pdata->rlocal == NULL) {
-    N_VDestroy(pdata->zlocal);
+    NVECTOR_MACRO(Destroy)(pdata->zlocal);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
     free(pdata); pdata = NULL;
@@ -153,10 +153,10 @@ int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal,
     return(CVLS_MEM_FAIL);
   }
   pdata->tmp1 = NULL;
-  pdata->tmp1 = N_VClone(cv_mem->cv_tempv);
+  pdata->tmp1 = NVECTOR_MACRO(Clone)(cv_mem->cv_tempv);
   if (pdata->tmp1 == NULL) {
-    N_VDestroy(pdata->zlocal);
-    N_VDestroy(pdata->rlocal);
+    NVECTOR_MACRO(Destroy)(pdata->zlocal);
+    NVECTOR_MACRO(Destroy)(pdata->rlocal);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
     free(pdata); pdata = NULL;
@@ -165,11 +165,11 @@ int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal,
     return(CVLS_MEM_FAIL);
   }
   pdata->tmp2 = NULL;
-  pdata->tmp2 = N_VClone(cv_mem->cv_tempv);
+  pdata->tmp2 = NVECTOR_MACRO(Clone)(cv_mem->cv_tempv);
   if (pdata->tmp2 == NULL) {
-    N_VDestroy(pdata->tmp1);
-    N_VDestroy(pdata->zlocal);
-    N_VDestroy(pdata->rlocal);
+    NVECTOR_MACRO(Destroy)(pdata->tmp1);
+    NVECTOR_MACRO(Destroy)(pdata->zlocal);
+    NVECTOR_MACRO(Destroy)(pdata->rlocal);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
     free(pdata); pdata = NULL;
@@ -178,12 +178,12 @@ int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal,
     return(CVLS_MEM_FAIL);
   }
   pdata->tmp3 = NULL;
-  pdata->tmp3 = N_VClone(cv_mem->cv_tempv);
+  pdata->tmp3 = NVECTOR_MACRO(Clone)(cv_mem->cv_tempv);
   if (pdata->tmp3 == NULL) {
-    N_VDestroy(pdata->tmp1);
-    N_VDestroy(pdata->tmp2);
-    N_VDestroy(pdata->zlocal);
-    N_VDestroy(pdata->rlocal);
+    NVECTOR_MACRO(Destroy)(pdata->tmp1);
+    NVECTOR_MACRO(Destroy)(pdata->tmp2);
+    NVECTOR_MACRO(Destroy)(pdata->zlocal);
+    NVECTOR_MACRO(Destroy)(pdata->rlocal);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
     free(pdata); pdata = NULL;
@@ -196,11 +196,11 @@ int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal,
   pdata->LS = NULL;
   pdata->LS = SUNLinSol_Band(pdata->rlocal, pdata->savedP, cv_mem->cv_sunctx);
   if (pdata->LS == NULL) {
-    N_VDestroy(pdata->tmp1);
-    N_VDestroy(pdata->tmp2);
-    N_VDestroy(pdata->tmp3);
-    N_VDestroy(pdata->zlocal);
-    N_VDestroy(pdata->rlocal);
+    NVECTOR_MACRO(Destroy)(pdata->tmp1);
+    NVECTOR_MACRO(Destroy)(pdata->tmp2);
+    NVECTOR_MACRO(Destroy)(pdata->tmp3);
+    NVECTOR_MACRO(Destroy)(pdata->zlocal);
+    NVECTOR_MACRO(Destroy)(pdata->rlocal);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
     free(pdata); pdata = NULL;
@@ -212,11 +212,11 @@ int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal,
   /* initialize band linear solver object */
   flag = SUNLinSolInitialize(pdata->LS);
   if (flag != SUNLS_SUCCESS) {
-    N_VDestroy(pdata->tmp1);
-    N_VDestroy(pdata->tmp2);
-    N_VDestroy(pdata->tmp3);
-    N_VDestroy(pdata->zlocal);
-    N_VDestroy(pdata->rlocal);
+    NVECTOR_MACRO(Destroy)(pdata->tmp1);
+    NVECTOR_MACRO(Destroy)(pdata->tmp2);
+    NVECTOR_MACRO(Destroy)(pdata->tmp3);
+    NVECTOR_MACRO(Destroy)(pdata->zlocal);
+    NVECTOR_MACRO(Destroy)(pdata->rlocal);
     SUNMatDestroy(pdata->savedP);
     SUNMatDestroy(pdata->savedJ);
     SUNLinSolFree(pdata->LS);
@@ -237,12 +237,12 @@ int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal,
   pdata->rpwsize = 0;
   pdata->ipwsize = 0;
   if (cv_mem->cv_tempv->ops->nvspace) {
-    N_VSpace(cv_mem->cv_tempv, &lrw1, &liw1);
+    NVECTOR_MACRO(Space)(cv_mem->cv_tempv, &lrw1, &liw1);
     pdata->rpwsize += 3*lrw1;
     pdata->ipwsize += 3*liw1;
   }
   if (pdata->rlocal->ops->nvspace) {
-    N_VSpace(pdata->rlocal, &lrw1, &liw1);
+    NVECTOR_MACRO(Space)(pdata->rlocal, &lrw1, &liw1);
     pdata->rpwsize += 2*lrw1;
     pdata->ipwsize += 2*liw1;
   }
@@ -549,16 +549,16 @@ static int CVBBDPrecSolve(realtype t, N_Vector y, N_Vector fy,
   pdata = (CVBBDPrecData) bbd_data;
 
   /* Attach local data arrays for r and z to rlocal and zlocal */
-  N_VSetArrayPointer(N_VGetArrayPointer(r), pdata->rlocal);
-  N_VSetArrayPointer(N_VGetArrayPointer(z), pdata->zlocal);
+  NVECTOR_MACRO(SetArrayPointer)(NVECTOR_MACRO(GetArrayPointer)(r), pdata->rlocal);
+  NVECTOR_MACRO(SetArrayPointer)(NVECTOR_MACRO(GetArrayPointer)(z), pdata->zlocal);
 
   /* Call banded solver object to do the work */
   retval = SUNLinSolSolve(pdata->LS, pdata->savedP, pdata->zlocal,
                           pdata->rlocal, ZERO);
 
   /* Detach local data arrays from rlocal and zlocal */
-  N_VSetArrayPointer(NULL, pdata->rlocal);
-  N_VSetArrayPointer(NULL, pdata->zlocal);
+  NVECTOR_MACRO(SetArrayPointer)(NULL, pdata->rlocal);
+  NVECTOR_MACRO(SetArrayPointer)(NULL, pdata->zlocal);
 
   return(retval);
 }
@@ -576,11 +576,11 @@ static int CVBBDPrecFree(CVodeMem cv_mem)
   pdata = (CVBBDPrecData) cvls_mem->P_data;
 
   SUNLinSolFree(pdata->LS);
-  N_VDestroy(pdata->tmp1);
-  N_VDestroy(pdata->tmp2);
-  N_VDestroy(pdata->tmp3);
-  N_VDestroy(pdata->zlocal);
-  N_VDestroy(pdata->rlocal);
+  NVECTOR_MACRO(Destroy)(pdata->tmp1);
+  NVECTOR_MACRO(Destroy)(pdata->tmp2);
+  NVECTOR_MACRO(Destroy)(pdata->tmp3);
+  NVECTOR_MACRO(Destroy)(pdata->zlocal);
+  NVECTOR_MACRO(Destroy)(pdata->rlocal);
   SUNMatDestroy(pdata->savedP);
   SUNMatDestroy(pdata->savedJ);
 
@@ -621,7 +621,7 @@ static int CVBBDDQJac(CVBBDPrecData pdata, realtype t, N_Vector y,
   cv_mem = (CVodeMem) pdata->cvode_mem;
 
   /* Load ytemp with y = predicted solution vector */
-  N_VScale(ONE, y, ytemp);
+  NVECTOR_MACRO(Scale)(ONE, y, ytemp);
 
   /* Call cfn and gloc to get base value of g(t,y) */
   if (pdata->cfn != NULL) {
@@ -635,16 +635,16 @@ static int CVBBDDQJac(CVBBDPrecData pdata, realtype t, N_Vector y,
   if (retval != 0) return(retval);
 
   /* Obtain pointers to the data for various vectors */
-  y_data     =  N_VGetArrayPointer(y);
-  gy_data    =  N_VGetArrayPointer(gy);
-  ewt_data   =  N_VGetArrayPointer(cv_mem->cv_ewt);
-  ytemp_data =  N_VGetArrayPointer(ytemp);
-  gtemp_data =  N_VGetArrayPointer(gtemp);
+  y_data     =  NVECTOR_MACRO(GetArrayPointer)(y);
+  gy_data    =  NVECTOR_MACRO(GetArrayPointer)(gy);
+  ewt_data   =  NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_ewt);
+  ytemp_data =  NVECTOR_MACRO(GetArrayPointer)(ytemp);
+  gtemp_data =  NVECTOR_MACRO(GetArrayPointer)(gtemp);
   if (cv_mem->cv_constraintsSet)
-    cns_data =  N_VGetArrayPointer(cv_mem->cv_constraints);
+    cns_data =  NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_constraints);
 
   /* Set minimum increment based on uround and norm of g */
-  gnorm = N_VWrmsNorm(gy, cv_mem->cv_ewt);
+  gnorm = NVECTOR_MACRO(WrmsNorm)(gy, cv_mem->cv_ewt);
   minInc = (gnorm != ZERO) ?
     (MIN_INC_MULT * SUNRabs(cv_mem->cv_h) *
      cv_mem->cv_uround * pdata->n_local * gnorm) : ONE;
diff --git a/src/cvode/cvode_diag.c b/src/cvode/cvode_diag.c
index 8f9c124a8..f75eb8391 100644
--- a/src/cvode/cvode_diag.c
+++ b/src/cvode/cvode_diag.c
@@ -152,26 +152,26 @@ int CVDiag(void *cvode_mem)
 
   /* Allocate memory for M, bit, and bitcomp */
 
-  M = N_VClone(vec_tmpl);
+  M = NVECTOR_MACRO(Clone)(vec_tmpl);
   if (M == NULL) {
     cvProcessError(cv_mem, CVDIAG_MEM_FAIL, "CVDIAG", "CVDiag", MSGDG_MEM_FAIL);
     free(cvdiag_mem); cvdiag_mem = NULL;
     return(CVDIAG_MEM_FAIL);
   }
 
-  bit = N_VClone(vec_tmpl);
+  bit = NVECTOR_MACRO(Clone)(vec_tmpl);
   if (bit == NULL) {
     cvProcessError(cv_mem, CVDIAG_MEM_FAIL, "CVDIAG", "CVDiag", MSGDG_MEM_FAIL);
-    N_VDestroy(M);
+    NVECTOR_MACRO(Destroy)(M);
     free(cvdiag_mem); cvdiag_mem = NULL;
     return(CVDIAG_MEM_FAIL);
   }
 
-  bitcomp = N_VClone(vec_tmpl);
+  bitcomp = NVECTOR_MACRO(Clone)(vec_tmpl);
   if (bitcomp == NULL) {
     cvProcessError(cv_mem, CVDIAG_MEM_FAIL, "CVDIAG", "CVDiag", MSGDG_MEM_FAIL);
-    N_VDestroy(M);
-    N_VDestroy(bit);
+    NVECTOR_MACRO(Destroy)(M);
+    NVECTOR_MACRO(Destroy)(bit);
     free(cvdiag_mem); cvdiag_mem = NULL;
     return(CVDIAG_MEM_FAIL);
   }
@@ -364,8 +364,8 @@ static int CVDiagSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
   else
 #endif
   {
-    N_VLinearSum(h, fpred, -ONE, zn[1], ftemp);
-    N_VLinearSum(r, ftemp, ONE, ypred, y);
+    NVECTOR_MACRO(LinearSum)(h, fpred, -ONE, zn[1], ftemp);
+    NVECTOR_MACRO(LinearSum)(r, ftemp, ONE, ypred, y);
   }
 
   /* Evaluate f at perturbed y */
@@ -390,21 +390,21 @@ static int CVDiagSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
   else
 #endif
   {
-    N_VLinearSum(ONE, M, -ONE, fpred, M);
-    N_VLinearSum(FRACT, ftemp, -h, M, M);
-    N_VProd(ftemp, ewt, y);
+    NVECTOR_MACRO(LinearSum)(ONE, M, -ONE, fpred, M);
+    NVECTOR_MACRO(LinearSum)(FRACT, ftemp, -h, M, M);
+    NVECTOR_MACRO(Prod)(ftemp, ewt, y);
     /* Protect against deltay_i being at roundoff level */
-    N_VCompare(uround, y, bit);
-    N_VAddConst(bit, -ONE, bitcomp);
-    N_VProd(ftemp, bit, y);
-    N_VLinearSum(FRACT, y, -ONE, bitcomp, y);
-    N_VDiv(M, y, M);
-    N_VProd(M, bit, M);
-    N_VLinearSum(ONE, M, -ONE, bitcomp, M);
+    NVECTOR_MACRO(Compare)(uround, y, bit);
+    NVECTOR_MACRO(AddConst)(bit, -ONE, bitcomp);
+    NVECTOR_MACRO(Prod)(ftemp, bit, y);
+    NVECTOR_MACRO(LinearSum)(FRACT, y, -ONE, bitcomp, y);
+    NVECTOR_MACRO(Div)(M, y, M);
+    NVECTOR_MACRO(Prod)(M, bit, M);
+    NVECTOR_MACRO(LinearSum)(ONE, M, -ONE, bitcomp, M);
   }
 
   /* Invert M with test for zero components */
-  invOK = N_VInvTest(M, M);
+  invOK = NVECTOR_MACRO(InvTest)(M, M);
   if (!invOK) {
     last_flag = CVDIAG_INV_FAIL;
     return(1);
@@ -447,12 +447,12 @@ static int CVDiagSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
     else
 #endif
     {
-      N_VInv(M, M);
-      N_VAddConst(M, -ONE, M);
-      N_VScale(r, M, M);
-      N_VAddConst(M, ONE, M);
+      NVECTOR_MACRO(Inv)(M, M);
+      NVECTOR_MACRO(AddConst)(M, -ONE, M);
+      NVECTOR_MACRO(Scale)(r, M, M);
+      NVECTOR_MACRO(AddConst)(M, ONE, M);
     }
-    invOK = N_VInvTest(M, M);
+    invOK = NVECTOR_MACRO(InvTest)(M, M);
     if (!invOK) {
       last_flag = CVDIAG_INV_FAIL;
       return (1);
@@ -461,7 +461,7 @@ static int CVDiagSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
   }
 
   /* Apply M-inverse to b */
-  N_VProd(b, M, b);
+  NVECTOR_MACRO(Prod)(b, M, b);
 
   last_flag = CVDIAG_SUCCESS;
   return(0);
@@ -481,9 +481,9 @@ static int CVDiagFree(CVodeMem cv_mem)
 
   cvdiag_mem = (CVDiagMem) lmem;
 
-  N_VDestroy(M);
-  N_VDestroy(bit);
-  N_VDestroy(bitcomp);
+  NVECTOR_MACRO(Destroy)(M);
+  NVECTOR_MACRO(Destroy)(bit);
+  NVECTOR_MACRO(Destroy)(bitcomp);
   free(cvdiag_mem);
   cv_mem->cv_lmem = NULL;
 
diff --git a/src/cvode/cvode_fused_stubs.c b/src/cvode/cvode_fused_stubs.c
index 4df3c5a12..6c057d1c6 100644
--- a/src/cvode/cvode_fused_stubs.c
+++ b/src/cvode/cvode_fused_stubs.c
@@ -39,13 +39,13 @@ int cvEwtSetSS_fused(const booleantype atolmin0,
                      N_Vector tempv,
                      N_Vector weight)
 {
-  N_VAbs(ycur, tempv);
-  N_VScale(reltol, tempv, tempv);
-  N_VAddConst(tempv, Sabstol, tempv);
+  NVECTOR_MACRO(Abs)(ycur, tempv);
+  NVECTOR_MACRO(Scale)(reltol, tempv, tempv);
+  NVECTOR_MACRO(AddConst)(tempv, Sabstol, tempv);
   if (atolmin0) {
-    if (N_VMin(tempv) <= ZERO) return(-1);
+    if (NVECTOR_MACRO(Min)(tempv) <= ZERO) return(-1);
   }
-  N_VInv(tempv, weight);
+  NVECTOR_MACRO(Inv)(tempv, weight);
   return 0;
 }
 
@@ -63,13 +63,13 @@ int cvEwtSetSV_fused(const booleantype atolmin0,
                      N_Vector tempv,
                      N_Vector weight)
 {
-  N_VAbs(ycur, tempv);
-  N_VLinearSum(reltol, tempv, ONE,
+  NVECTOR_MACRO(Abs)(ycur, tempv);
+  NVECTOR_MACRO(LinearSum)(reltol, tempv, ONE,
                Vabstol, tempv);
   if (atolmin0) {
-    if (N_VMin(tempv) <= ZERO) return(-1);
+    if (NVECTOR_MACRO(Min)(tempv) <= ZERO) return(-1);
   }
-  N_VInv(tempv, weight);
+  NVECTOR_MACRO(Inv)(tempv, weight);
   return 0;
 }
 
@@ -88,11 +88,11 @@ int cvCheckConstraints_fused(const N_Vector c,
                              const N_Vector mm,
                              N_Vector tmp)
 {
-  N_VCompare(ONEPT5, c, tmp);           /* a[i]=1 when |c[i]|=2  */
-  N_VProd(tmp, c, tmp);                 /* a * c                 */
-  N_VDiv(tmp, ewt, tmp);                /* a * c * wt            */
-  N_VLinearSum(ONE, y, -PT1, tmp, tmp); /* y - 0.1 * a * c * wt  */
-  N_VProd(tmp, mm, tmp);                /* v = mm*(y-0.1*a*c*wt) */
+  NVECTOR_MACRO(Compare)(ONEPT5, c, tmp);           /* a[i]=1 when |c[i]|=2  */
+  NVECTOR_MACRO(Prod)(tmp, c, tmp);                 /* a * c                 */
+  NVECTOR_MACRO(Div)(tmp, ewt, tmp);                /* a * c * wt            */
+  NVECTOR_MACRO(LinearSum)(ONE, y, -PT1, tmp, tmp); /* y - 0.1 * a * c * wt  */
+  NVECTOR_MACRO(Prod)(tmp, mm, tmp);                /* v = mm*(y-0.1*a*c*wt) */
   return 0;
 }
 
@@ -111,8 +111,8 @@ int cvNlsResid_fused(const realtype rl1,
                      const N_Vector ftemp,
                      N_Vector res)
 {
-  N_VLinearSum(rl1, zn1, ONE, ycor, res);
-  N_VLinearSum(ngamma, ftemp, ONE, res, res);
+  NVECTOR_MACRO(LinearSum)(rl1, zn1, ONE, ycor, res);
+  NVECTOR_MACRO(LinearSum)(ngamma, ftemp, ONE, res, res);
   return 0;
 }
 
@@ -130,8 +130,8 @@ int cvDiagSetup_formY(const realtype h,
                       N_Vector ftemp,
                       N_Vector y)
 {
-  N_VLinearSum(h, fpred, -ONE, zn1, ftemp);
-  N_VLinearSum(r, ftemp, ONE, ypred, y);
+  NVECTOR_MACRO(LinearSum)(h, fpred, -ONE, zn1, ftemp);
+  NVECTOR_MACRO(LinearSum)(r, ftemp, ONE, ypred, y);
   return 0;
 }
 
@@ -153,17 +153,17 @@ int cvDiagSetup_buildM(const realtype fract,
                        N_Vector y,
                        N_Vector M)
 {
-  N_VLinearSum(ONE, M, -ONE, fpred, M);
-  N_VLinearSum(FRACT, ftemp, -h, M, M);
-  N_VProd(ftemp, ewt, y);
+  NVECTOR_MACRO(LinearSum)(ONE, M, -ONE, fpred, M);
+  NVECTOR_MACRO(LinearSum)(FRACT, ftemp, -h, M, M);
+  NVECTOR_MACRO(Prod)(ftemp, ewt, y);
   /* Protect against deltay_i being at roundoff level */
-  N_VCompare(uround, y, bit);
-  N_VAddConst(bit, -ONE, bitcomp);
-  N_VProd(ftemp, bit, y);
-  N_VLinearSum(FRACT, y, -ONE, bitcomp, y);
-  N_VDiv(M, y, M);
-  N_VProd(M, bit, M);
-  N_VLinearSum(ONE, M, -ONE, bitcomp, M);
+  NVECTOR_MACRO(Compare)(uround, y, bit);
+  NVECTOR_MACRO(AddConst)(bit, -ONE, bitcomp);
+  NVECTOR_MACRO(Prod)(ftemp, bit, y);
+  NVECTOR_MACRO(LinearSum)(FRACT, y, -ONE, bitcomp, y);
+  NVECTOR_MACRO(Div)(M, y, M);
+  NVECTOR_MACRO(Prod)(M, bit, M);
+  NVECTOR_MACRO(LinearSum)(ONE, M, -ONE, bitcomp, M);
   return 0;
 }
 
@@ -176,9 +176,9 @@ int cvDiagSetup_buildM(const realtype fract,
 
 int cvDiagSolve_updateM(const realtype r, N_Vector M)
 {
-  N_VInv(M, M);
-  N_VAddConst(M, -ONE, M);
-  N_VScale(r, M, M);
-  N_VAddConst(M, ONE, M);
+  NVECTOR_MACRO(Inv)(M, M);
+  NVECTOR_MACRO(AddConst)(M, -ONE, M);
+  NVECTOR_MACRO(Scale)(r, M, M);
+  NVECTOR_MACRO(AddConst)(M, ONE, M);
   return 0;
 }
\ No newline at end of file
diff --git a/src/cvode/cvode_impl.h b/src/cvode/cvode_impl.h
index 9996a03d8..360edaea0 100644
--- a/src/cvode/cvode_impl.h
+++ b/src/cvode/cvode_impl.h
@@ -26,6 +26,7 @@
 #include "sundials_context_impl.h"
 #include "sundials_logger_impl.h"
 #include "sundials/sundials_math.h"
+#include <nvector/nvector_serial.h>
 
 #ifdef __cplusplus  /* wrapper to enable C++ usage */
 extern "C" {
diff --git a/src/cvode/cvode_io.c b/src/cvode/cvode_io.c
index 1c2ac7593..493b79eb3 100644
--- a/src/cvode/cvode_io.c
+++ b/src/cvode/cvode_io.c
@@ -915,7 +915,7 @@ int CVodeSetConstraints(void *cvode_mem, N_Vector constraints)
   /* If there are no constraints, destroy data structures */
   if (constraints == NULL) {
     if (cv_mem->cv_constraintsMallocDone) {
-      N_VDestroy(cv_mem->cv_constraints);
+      NVECTOR_MACRO(Destroy)(cv_mem->cv_constraints);
       cv_mem->cv_lrw -= cv_mem->cv_lrw1;
       cv_mem->cv_liw -= cv_mem->cv_liw1;
     }
@@ -936,21 +936,21 @@ int CVodeSetConstraints(void *cvode_mem, N_Vector constraints)
   }
 
   /* Check the constraints vector */
-  temptest = N_VMaxNorm(constraints);
+  temptest = NVECTOR_MACRO(MaxNorm)(constraints);
   if ((temptest > TWOPT5) || (temptest < HALF)) {
     cvProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeSetConstraints", MSGCV_BAD_CONSTR);
     return(CV_ILL_INPUT);
   }
 
   if ( !(cv_mem->cv_constraintsMallocDone) ) {
-    cv_mem->cv_constraints = N_VClone(constraints);
+    cv_mem->cv_constraints = NVECTOR_MACRO(Clone)(constraints);
     cv_mem->cv_lrw += cv_mem->cv_lrw1;
     cv_mem->cv_liw += cv_mem->cv_liw1;
     cv_mem->cv_constraintsMallocDone = SUNTRUE;
   }
 
   /* Load the constraints vector */
-  N_VScale(ONE, constraints, cv_mem->cv_constraints);
+  NVECTOR_MACRO(Scale)(ONE, constraints, cv_mem->cv_constraints);
 
   cv_mem->cv_constraintsSet = SUNTRUE;
 
@@ -979,7 +979,7 @@ int CVodeSetUseIntegratorFusedKernels(void *cvode_mem, booleantype onoff)
   cv_mem = (CVodeMem) cvode_mem;
 
 #ifdef SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS
-  id = N_VGetVectorID(cv_mem->cv_ewt);
+  id = NVECTOR_MACRO(GetVectorID)(cv_mem->cv_ewt);
   if (!cv_mem->cv_MallocDone ||
       (id != SUNDIALS_NVEC_CUDA && id != SUNDIALS_NVEC_HIP)) {
     cvProcessError(cv_mem, CV_ILL_INPUT, "CVODE",
@@ -1332,7 +1332,7 @@ int CVodeGetErrWeights(void *cvode_mem, N_Vector eweight)
 
   cv_mem = (CVodeMem) cvode_mem;
 
-  N_VScale(ONE, cv_mem->cv_ewt, eweight);
+  NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_ewt, eweight);
 
   return(CV_SUCCESS);
 }
@@ -1354,7 +1354,7 @@ int CVodeGetEstLocalErrors(void *cvode_mem, N_Vector ele)
 
   cv_mem = (CVodeMem) cvode_mem;
 
-  N_VScale(ONE, cv_mem->cv_acor, ele);
+  NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_acor, ele);
 
   return(CV_SUCCESS);
 }
diff --git a/src/cvode/cvode_ls.c b/src/cvode/cvode_ls.c
index f9bf138ec..6cebf7c17 100644
--- a/src/cvode/cvode_ls.c
+++ b/src/cvode/cvode_ls.c
@@ -228,7 +228,7 @@ int CVodeSetLinearSolver(void *cvode_mem, SUNLinearSolver LS,
   }
 
   /* Allocate memory for ytemp and x */
-  cvls_mem->ytemp = N_VClone(cv_mem->cv_tempv);
+  cvls_mem->ytemp = NVECTOR_MACRO(Clone)(cv_mem->cv_tempv);
   if (cvls_mem->ytemp == NULL) {
     cvProcessError(cv_mem, CVLS_MEM_FAIL, "CVLS",
                     "CVodeSetLinearSolver", MSG_LS_MEM_FAIL);
@@ -236,18 +236,18 @@ int CVodeSetLinearSolver(void *cvode_mem, SUNLinearSolver LS,
     return(CVLS_MEM_FAIL);
   }
 
-  cvls_mem->x = N_VClone(cv_mem->cv_tempv);
+  cvls_mem->x = NVECTOR_MACRO(Clone)(cv_mem->cv_tempv);
   if (cvls_mem->x == NULL) {
     cvProcessError(cv_mem, CVLS_MEM_FAIL, "CVLS",
                     "CVodeSetLinearSolver", MSG_LS_MEM_FAIL);
-    N_VDestroy(cvls_mem->ytemp);
+    NVECTOR_MACRO(Destroy)(cvls_mem->ytemp);
     free(cvls_mem); cvls_mem = NULL;
     return(CVLS_MEM_FAIL);
   }
 
   /* For iterative LS, compute default norm conversion factor */
   if (iterative)
-    cvls_mem->nrmfac = SUNRsqrt( N_VGetLength(cvls_mem->ytemp) );
+    cvls_mem->nrmfac = SUNRsqrt( NVECTOR_MACRO(GetLength)(cvls_mem->ytemp) );
 
   /* Check if solution scaling should be enabled */
   if (matrixbased && cv_mem->cv_lmm == CV_BDF)
@@ -373,11 +373,11 @@ int CVodeSetLSNormFactor(void *cvode_mem, realtype nrmfac)
     cvls_mem->nrmfac = nrmfac;
   } else if (nrmfac < ZERO) {
     /* compute factor for WRMS norm with dot product */
-    N_VConst(ONE, cvls_mem->ytemp);
-    cvls_mem->nrmfac = SUNRsqrt(N_VDotProd(cvls_mem->ytemp, cvls_mem->ytemp));
+    NVECTOR_MACRO(Const)(ONE, cvls_mem->ytemp);
+    cvls_mem->nrmfac = SUNRsqrt(NVECTOR_MACRO(DotProd)(cvls_mem->ytemp, cvls_mem->ytemp));
   } else {
     /* compute default factor for WRMS norm from vector legnth */
-    cvls_mem->nrmfac = SUNRsqrt(N_VGetLength(cvls_mem->ytemp));
+    cvls_mem->nrmfac = SUNRsqrt(NVECTOR_MACRO(GetLength)(cvls_mem->ytemp));
   }
 
   return(CVLS_SUCCESS);
@@ -605,7 +605,7 @@ int CVodeGetLinWorkSpace(void *cvode_mem, long int *lenrwLS,
 
   /* add NVector sizes */
   if (cv_mem->cv_tempv->ops->nvspace) {
-    N_VSpace(cv_mem->cv_tempv, &lrw1, &liw1);
+    NVECTOR_MACRO(Space)(cv_mem->cv_tempv, &lrw1, &liw1);
     *lenrwLS += 2*lrw1;
     *leniwLS += 2*liw1;
   }
@@ -893,7 +893,7 @@ int cvLsATimes(void *cvode_mem, N_Vector v, N_Vector z)
   if (retval != 0) return(retval);
 
   /* add contribution from identity matrix */
-  N_VLinearSum(ONE, v, -cv_mem->cv_gamma, z, z);
+  NVECTOR_MACRO(LinearSum)(ONE, v, -cv_mem->cv_gamma, z, z);
 
   return(0);
 }
@@ -1054,24 +1054,24 @@ int cvLsDenseDQJac(realtype t, N_Vector y, N_Vector fy,
   ftemp = tmp1;
 
   /* Create an empty vector for matrix column calculations */
-  jthCol = N_VCloneEmpty(tmp1);
+  jthCol = NVECTOR_MACRO(CloneEmpty)(tmp1);
 
   /* Obtain pointers to the data for ewt, y */
-  ewt_data = N_VGetArrayPointer(cv_mem->cv_ewt);
-  y_data   = N_VGetArrayPointer(y);
+  ewt_data = NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_ewt);
+  y_data   = NVECTOR_MACRO(GetArrayPointer)(y);
   if (cv_mem->cv_constraintsSet)
-    cns_data = N_VGetArrayPointer(cv_mem->cv_constraints);
+    cns_data = NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_constraints);
 
   /* Set minimum increment based on uround and norm of f */
   srur = SUNRsqrt(cv_mem->cv_uround);
-  fnorm = N_VWrmsNorm(fy, cv_mem->cv_ewt);
+  fnorm = NVECTOR_MACRO(WrmsNorm)(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);
+    NVECTOR_MACRO(SetArrayPointer)(SUNDenseMatrix_Column(Jac,j), jthCol);
 
     yjsaved = y_data[j];
     inc = SUNMAX(srur*SUNRabs(yjsaved), minInc/ewt_data[j]);
@@ -1092,13 +1092,13 @@ int cvLsDenseDQJac(realtype t, N_Vector y, N_Vector fy,
     y_data[j] = yjsaved;
 
     inc_inv = ONE/inc;
-    N_VLinearSum(inc_inv, ftemp, -inc_inv, fy, jthCol);
+    NVECTOR_MACRO(LinearSum)(inc_inv, ftemp, -inc_inv, fy, jthCol);
 
   }
 
   /* Destroy jthCol vector */
-  N_VSetArrayPointer(NULL, jthCol);  /* SHOULDN'T BE NEEDED */
-  N_VDestroy(jthCol);
+  NVECTOR_MACRO(SetArrayPointer)(NULL, jthCol);  /* SHOULDN'T BE NEEDED */
+  NVECTOR_MACRO(Destroy)(jthCol);
 
   return(retval);
 }
@@ -1143,20 +1143,20 @@ int cvLsBandDQJac(realtype t, N_Vector y, N_Vector fy, SUNMatrix Jac,
   ytemp = tmp2;
 
   /* Obtain pointers to the data for ewt, fy, ftemp, y, ytemp */
-  ewt_data   = N_VGetArrayPointer(cv_mem->cv_ewt);
-  fy_data    = N_VGetArrayPointer(fy);
-  ftemp_data = N_VGetArrayPointer(ftemp);
-  y_data     = N_VGetArrayPointer(y);
-  ytemp_data = N_VGetArrayPointer(ytemp);
+  ewt_data   = NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_ewt);
+  fy_data    = NVECTOR_MACRO(GetArrayPointer)(fy);
+  ftemp_data = NVECTOR_MACRO(GetArrayPointer)(ftemp);
+  y_data     = NVECTOR_MACRO(GetArrayPointer)(y);
+  ytemp_data = NVECTOR_MACRO(GetArrayPointer)(ytemp);
   if (cv_mem->cv_constraintsSet)
-    cns_data = N_VGetArrayPointer(cv_mem->cv_constraints);
+    cns_data = NVECTOR_MACRO(GetArrayPointer)(cv_mem->cv_constraints);
 
   /* Load ytemp with y = predicted y vector */
-  N_VScale(ONE, y, ytemp);
+  NVECTOR_MACRO(Scale)(ONE, y, ytemp);
 
   /* Set minimum increment based on uround and norm of f */
   srur = SUNRsqrt(cv_mem->cv_uround);
-  fnorm = N_VWrmsNorm(fy, cv_mem->cv_ewt);
+  fnorm = NVECTOR_MACRO(WrmsNorm)(fy, cv_mem->cv_ewt);
   minInc = (fnorm != ZERO) ?
     (MIN_INC_MULT * SUNRabs(cv_mem->cv_h) * cv_mem->cv_uround * N * fnorm) : ONE;
 
@@ -1234,12 +1234,12 @@ int cvLsDQJtimes(N_Vector v, N_Vector Jv, realtype t,
   if (retval != CVLS_SUCCESS)  return(retval);
 
   /* Initialize perturbation to 1/||v|| */
-  sig = ONE/N_VWrmsNorm(v, cv_mem->cv_ewt);
+  sig = ONE/NVECTOR_MACRO(WrmsNorm)(v, cv_mem->cv_ewt);
 
   for (iter=0; iter<MAX_DQITERS; iter++) {
 
     /* Set work = y + sig*v */
-    N_VLinearSum(sig, v, ONE, y, work);
+    NVECTOR_MACRO(LinearSum)(sig, v, ONE, y, work);
 
     /* Set Jv = f(tn, y+sig*v) */
     retval = cvls_mem->jt_f(t, work, Jv, cv_mem->cv_user_data);
@@ -1256,7 +1256,7 @@ int cvLsDQJtimes(N_Vector v, N_Vector Jv, realtype t,
 
   /* Replace Jv by (Jv - fy)/sig */
   siginv = ONE/sig;
-  N_VLinearSum(siginv, Jv, -siginv, fy, Jv);
+  NVECTOR_MACRO(LinearSum)(siginv, Jv, -siginv, fy, Jv);
 
   return(0);
 }
@@ -1615,9 +1615,9 @@ int cvLsSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
      set linear solver tolerance (in left/right scaled 2-norm) */
   if (cvls_mem->iterative) {
     deltar = cvls_mem->eplifac * cv_mem->cv_tq[4];
-    bnorm = N_VWrmsNorm(b, weight);
+    bnorm = NVECTOR_MACRO(WrmsNorm)(b, weight);
     if (bnorm <= deltar) {
-      if (curiter > 0) N_VConst(ZERO, b);
+      if (curiter > 0) NVECTOR_MACRO(Const)(ZERO, b);
       cvls_mem->last_flag = CVLS_SUCCESS;
       return(cvls_mem->last_flag);
     }
@@ -1660,14 +1660,14 @@ int cvLsSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
      So we compute w_mean = ||w||_RMS = ||w||_2 and scale the desired tolerance accordingly. */
   } else if (cvls_mem->iterative) {
 
-    N_VConst(ONE, cvls_mem->x);
-    w_mean = N_VWrmsNorm(weight, cvls_mem->x);
+    NVECTOR_MACRO(Const)(ONE, cvls_mem->x);
+    w_mean = NVECTOR_MACRO(WrmsNorm)(weight, cvls_mem->x);
     delta /= w_mean;
 
   }
 
   /* Set initial guess x = 0 to LS */
-  N_VConst(ZERO, cvls_mem->x);
+  NVECTOR_MACRO(Const)(ZERO, cvls_mem->x);
 
   /* Set zero initial guess flag */
   retval = SUNLinSolSetZeroGuess(cvls_mem->LS, SUNTRUE);
@@ -1691,12 +1691,12 @@ int cvLsSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
 
   /* Call solver, and copy x to b */
   retval = SUNLinSolSolve(cvls_mem->LS, cvls_mem->A, cvls_mem->x, b, delta);
-  N_VScale(ONE, cvls_mem->x, b);
+  NVECTOR_MACRO(Scale)(ONE, cvls_mem->x, b);
 
   /* If using a direct or matrix-iterative solver, BDF method, and gamma has changed,
      scale the correction to account for change in gamma */
   if (cvls_mem->scalesol && cv_mem->cv_gamrat != ONE)
-    N_VScale(TWO/(ONE + cv_mem->cv_gamrat), b, b);
+    NVECTOR_MACRO(Scale)(TWO/(ONE + cv_mem->cv_gamrat), b, b);
 
   /* Retrieve statistics from iterative linear solvers */
 #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO
@@ -1792,11 +1792,11 @@ int cvLsFree(CVodeMem cv_mem)
 
   /* Free N_Vector memory */
   if (cvls_mem->ytemp) {
-    N_VDestroy(cvls_mem->ytemp);
+    NVECTOR_MACRO(Destroy)(cvls_mem->ytemp);
     cvls_mem->ytemp = NULL;
   }
   if (cvls_mem->x) {
-    N_VDestroy(cvls_mem->x);
+    NVECTOR_MACRO(Destroy)(cvls_mem->x);
     cvls_mem->x = NULL;
   }
 
diff --git a/src/cvode/cvode_nls.c b/src/cvode/cvode_nls.c
index 6f5ea0436..f5e1b10b8 100644
--- a/src/cvode/cvode_nls.c
+++ b/src/cvode/cvode_nls.c
@@ -319,7 +319,7 @@ static int cvNlsConvTest(SUNNonlinearSolver NLS, N_Vector ycor, N_Vector delta,
   cv_mem = (CVodeMem) cvode_mem;
 
   /* compute the norm of the correction */
-  del = N_VWrmsNorm(delta, ewt);
+  del = NVECTOR_MACRO(WrmsNorm)(delta, ewt);
 
   /* get the current nonlinear solver iteration count */
   retval = SUNNonlinSolGetCurIter(NLS, &m);
@@ -333,7 +333,7 @@ static int cvNlsConvTest(SUNNonlinearSolver NLS, N_Vector ycor, N_Vector delta,
   dcon = del * SUNMIN(ONE, cv_mem->cv_crate) / tol;
 
   if (dcon <= ONE) {
-    cv_mem->cv_acnrm = (m==0) ? del : N_VWrmsNorm(ycor, ewt);
+    cv_mem->cv_acnrm = (m==0) ? del : NVECTOR_MACRO(WrmsNorm)(ycor, ewt);
     cv_mem->cv_acnrmcur = SUNTRUE;
     return(CV_SUCCESS); /* Nonlinear system was solved successfully */
   }
@@ -361,7 +361,7 @@ static int cvNlsResidual(N_Vector ycor, N_Vector res, void* cvode_mem)
   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);
+  NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[0], ONE, ycor, cv_mem->cv_y);
 
   /* evaluate the rhs function */
   retval = cv_mem->nls_f(cv_mem->cv_tn, cv_mem->cv_y, cv_mem->cv_ftemp,
@@ -379,8 +379,8 @@ static int cvNlsResidual(N_Vector ycor, N_Vector res, void* cvode_mem)
   else
 #endif
   {
-    N_VLinearSum(cv_mem->cv_rl1, cv_mem->cv_zn[1], ONE, ycor, res);
-    N_VLinearSum(-cv_mem->cv_gamma, cv_mem->cv_ftemp, ONE, res, res);
+    NVECTOR_MACRO(LinearSum)(cv_mem->cv_rl1, cv_mem->cv_zn[1], ONE, ycor, res);
+    NVECTOR_MACRO(LinearSum)(-cv_mem->cv_gamma, cv_mem->cv_ftemp, ONE, res, res);
   }
 
   return(CV_SUCCESS);
@@ -399,7 +399,7 @@ static int cvNlsFPFunction(N_Vector ycor, N_Vector res, void* cvode_mem)
   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);
+  NVECTOR_MACRO(LinearSum)(ONE, cv_mem->cv_zn[0], ONE, ycor, cv_mem->cv_y);
 
   /* evaluate the rhs function */
   retval = cv_mem->nls_f(cv_mem->cv_tn, cv_mem->cv_y, res,
@@ -408,8 +408,8 @@ static int cvNlsFPFunction(N_Vector ycor, N_Vector res, void* cvode_mem)
   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);
+  NVECTOR_MACRO(LinearSum)(cv_mem->cv_h, res, -ONE, cv_mem->cv_zn[1], res);
+  NVECTOR_MACRO(Scale)(cv_mem->cv_rl1, res, res);
 
   return(CV_SUCCESS);
 }
diff --git a/src/cvode/cvode_proj.c b/src/cvode/cvode_proj.c
index 11ed32a61..ca19f2e29 100644
--- a/src/cvode/cvode_proj.c
+++ b/src/cvode/cvode_proj.c
@@ -321,7 +321,7 @@ int cvDoProjection(CVodeMem cv_mem, int *nflagPtr, realtype saved_t,
     errP = NULL;
 
   /* Copy acor into errP (if projecting the error) */
-  if (proj_mem->err_proj) N_VScale(ONE, cv_mem->cv_acor, errP);
+  if (proj_mem->err_proj) NVECTOR_MACRO(Scale)(ONE, cv_mem->cv_acor, errP);
 
   /* Call the user projection function */
   retval = proj_mem->pfun(cv_mem->cv_tn, cv_mem->cv_y, acorP,
@@ -336,7 +336,7 @@ int cvDoProjection(CVodeMem cv_mem, int *nflagPtr, realtype saved_t,
   {
     /* Recompute acnrm to be used in error test (if projecting the error) */
     if (proj_mem->err_proj)
-      cv_mem->cv_acnrm = N_VWrmsNorm(errP, cv_mem->cv_ewt);
+      cv_mem->cv_acnrm = NVECTOR_MACRO(WrmsNorm)(errP, cv_mem->cv_ewt);
 
     /* The projection was successful, return now */
     cv_mem->proj_applied = SUNTRUE;
