Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 28 Feb 2018 23:20:05 +0000 (23:20 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 28 Feb 2018 23:20:05 +0000 (23:20 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15393 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/CVODES.java

index a899409..b7b5af2 100644 (file)
@@ -104,6 +104,19 @@ public abstract class CVODES {
        //        and return the solution at the point reached by that step.\r
        final int CV_NORMAL = 1;\r
        final int CV_ONE_STEP = 2;\r
+       \r
+       //ism:   This parameter specifies the sensitivity corrector type\r
+       //        to be used. In the CV_SIMULTANEOUS case, the nonlinear\r
+       //        systems for states and all sensitivities are solved\r
+       //        simultaneously. In the CV_STAGGERED case, the nonlinear\r
+       //        system for states is solved first and then, the\r
+       //        nonlinear systems for all sensitivities are solved\r
+       //        at the same time. Finally, in the CV_STAGGERED1 approach\r
+       //        all nonlinear systems are solved in a sequence.\r
+       final int CV_SIMULTANEOUS = 1;\r
+       final int CV_STAGGERED = 2;\r
+       final int CV_STAGGERED1 = 3;\r
+\r
 \r
        \r
        /*\r
@@ -194,8 +207,25 @@ public abstract class CVODES {
        final int MXNEF = 7;\r
        final double CORTES = 0.1;\r
        final double ZERO = 0.0;\r
+       final double PT1 = 0.1;\r
+       final double POINT2 = 0.2;\r
+       final double HALF = 0.5;\r
+       final double H_BIAS = HALF;\r
        final double ONE = 1.0;\r
+       final double TWO = 2.0;\r
+       final double FOUR = 4.0;\r
+       final double FIVE = 5.0;\r
+       final double TWELVE = 12.0;\r
+       final double HUNDRED = 100.0;\r
        final double ETAMX1 = 10000.0; \r
+       final double HUB_FACTOR = 0.1;\r
+       final double HLB_FACTOR = 100.0;\r
+       final double FUZZ_FACTOR = 100.0;\r
+       final int MAX_ITERS = 4;\r
+       \r
+       final int RTFOUND = +1;\r
+       final int CLOSERT = +3;\r
+\r
 \r
        \r
        final String MSG_TIME = "t = %lg";\r
@@ -372,6 +402,23 @@ public abstract class CVODES {
        final String MSGD_LMEMB_NULL = "Linear solver memory is NULL for the backward integration.";\r
        final String MSGD_BAD_TINTERP = "Bad t for interpolation.";\r
        \r
+       final int DO_ERROR_TEST = +2;\r
+       final int PREDICT_AGAIN = +3;\r
+\r
+       final int CONV_FAIL = +4; \r
+       final int TRY_AGAIN = +5;\r
+\r
+       final int FIRST_CALL = +6;\r
+       final int PREV_CONV_FAIL = +7;\r
+       final int PREV_ERR_FAIL = +8;\r
+\r
+       final int RHSFUNC_RECVR = +9;\r
+\r
+       final int QRHSFUNC_RECVR = +11;\r
+       final int SRHSFUNC_RECVR = +12;\r
+       final int QSRHSFUNC_RECVR = +13;\r
+\r
+       \r
        /*-----------------------------------------------------------------\r
         * IV. SUNLinearSolver error codes\r
         * ---------------------------------------------------------------\r
@@ -594,6 +641,7 @@ public abstract class CVODES {
         * @param t\r
         * @param yv\r
         * @param ydotv\r
+        * @param user_data\r
         * @return\r
         */\r
        private int f(double t, NVector yv, NVector ydotv, CVodeMemRec user_data) {\r
@@ -616,9 +664,10 @@ public abstract class CVODES {
         * @param t\r
         * @param yv\r
         * @param gout\r
+        * @param user_data\r
         * @return\r
         */\r
-       private int g(double t, NVector yv, double gout[]) {\r
+       private int g(double t, NVector yv, double gout[], CVodeMemRec user_data) {\r
                double y[] = yv.getData();\r
                if (problem == 1) {\r
                    gout[0] = y[0] - 0.0001;\r
@@ -654,6 +703,23 @@ public abstract class CVODES {
            return 0;\r
        }\r
        \r
+       private int fS(int cv_Ns, double time, NVector ycur, NVector fcur, NVector yScur[], \r
+            NVector fScur[], CVodeMemRec cv_fS_data, NVector temp1, NVector temp2) {\r
+               return 0;\r
+       }\r
+       \r
+       private int fS1(int cv_Ns, double time, NVector ycur, NVector fcur, int is, NVector yScur, \r
+            NVector fScur, CVodeMemRec cv_fS_data, NVector temp1, NVector temp2) {\r
+               return 0;\r
+       }\r
+       \r
+       private int fQS(int cv_Ns, double cv_tn, NVector cv_zn,\r
+                                NVector cv_znS[], NVector cv_znQ,\r
+                                NVector cv_znQS[], CVodeMemRec cv_fQS_data,\r
+                                NVector cv_tempv, NVector cv_tempvQ) {\r
+               return 0;\r
+       }\r
+       \r
        // Types: struct CVodeMemRec, CVodeMem\r
        // -----------------------------------------------------------------\r
        // The type CVodeMem is type pointer to struct CVodeMemRec.\r
@@ -726,6 +792,7 @@ public abstract class CVODES {
          //CVSensRhsFn cv_fS;          /* fS = (df/dy)*yS + (df/dp)                    */\r
          //CVSensRhs1Fn cv_fS1;        /* fS1 = (df/dy)*yS_i + (df/dp)                 */\r
          //void *cv_fS_data;           /* data pointer passed to fS                    */\r
+         CVodeMemRec cv_fS_data;\r
          boolean cv_fSDQ;        /* SUNTRUE if using internal DQ functions       */\r
          int cv_ifS;                 /* ifS = ALLSENS or ONESENS                     */\r
 \r
@@ -750,6 +817,7 @@ public abstract class CVODES {
 \r
          //CVQuadSensRhsFn cv_fQS;     /* fQS = (dfQ/dy)*yS + (dfQ/dp)                 */\r
          //void *cv_fQS_data;          /* data pointer passed to fQS                   */\r
+         CVodeMemRec cv_fQS_data;\r
          boolean cv_fQSDQ;       /* SUNTRUE if using internal DQ functions       */\r
 \r
          boolean cv_errconQS;    /* SUNTRUE if yQS are considered in err. con.   */\r
@@ -1236,7 +1304,7 @@ public abstract class CVODES {
 \r
      void cvProcessError(CVodeMemRec cv_mem, \r
                          int error_code,  String module, String fname, \r
-                         String msgfmt)\r
+                         String msgfmt, double... numbers)\r
      {\r
         MipavUtil.displayError("In module " + module + " in function " + fname + " msgfmt");\r
        //va_list ap;\r
@@ -2043,14 +2111,14 @@ public abstract class CVODES {
      * ----------------------------------------\r
      */\r
 \r
-    /*if (cv_mem.cv_nst == 0) {\r
+    if (cv_mem.cv_nst == 0) {\r
 \r
-      cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;*/\r
+      cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
 \r
       /* Check inputs for corectness */\r
 \r
-      /*ier = cvInitialSetup(cv_mem);\r
-      if (ier!= CV_SUCCESS) return(ier);*/\r
+      ier = cvInitialSetup(cv_mem);\r
+      if (ier!= CV_SUCCESS) return(ier);\r
 \r
       /* \r
        * Call f at (t0,y0), set zn[1] = y'(t0). \r
@@ -2059,13 +2127,12 @@ public abstract class CVODES {
        * If computing quadr. sensi., call fQS at (t0,y0,yS0), set znQS[1][is] = yQS'(t0), is=1,...,Ns.\r
        */\r
 \r
-      /*retval = f(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+      retval = f(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
                             cv_mem.cv_zn[1], cv_mem.cv_user_data); \r
       cv_mem.cv_nfe++;\r
       if (retval < 0) {\r
         cvProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODES", "CVode",\r
-                       //MSGCV_RHSFUNC_FAILED, cv_mem.cv_tn);\r
-                         MSGCV_RHSFUNC_FAILED);\r
+                       MSGCV_RHSFUNC_FAILED, cv_mem.cv_tn);\r
         return(CV_RHSFUNC_FAIL);\r
       }\r
       if (retval > 0) {\r
@@ -2080,7 +2147,7 @@ public abstract class CVODES {
         cv_mem.cv_nfQe++;\r
         if (retval < 0) {\r
           cvProcessError(cv_mem, CV_QRHSFUNC_FAIL, "CVODES", "CVode",\r
-                         MSGCV_QRHSFUNC_FAILED, cv_mem->cv_tn);\r
+                         MSGCV_QRHSFUNC_FAILED, cv_mem.cv_tn);\r
           return(CV_QRHSFUNC_FAIL);\r
         }\r
         if (retval > 0) {\r
@@ -2091,13 +2158,13 @@ public abstract class CVODES {
       }\r
 \r
       if (cv_mem.cv_sensi) {\r
-        retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, cv_mem->cv_zn[0],\r
-                                  cv_mem->cv_zn[1], cv_mem->cv_znS[0],\r
-                                  cv_mem->cv_znS[1], cv_mem->cv_tempv,\r
-                                  cv_mem->cv_ftemp);\r
+        retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                                  cv_mem.cv_zn[1], cv_mem.cv_znS[0],\r
+                                  cv_mem.cv_znS[1], cv_mem.cv_tempv,\r
+                                  cv_mem.cv_ftemp);\r
         if (retval < 0) {\r
           cvProcessError(cv_mem, CV_SRHSFUNC_FAIL, "CVODES", "CVode",\r
-                         MSGCV_SRHSFUNC_FAILED, cv_mem->cv_tn);\r
+                         MSGCV_SRHSFUNC_FAILED, cv_mem.cv_tn);\r
           return(CV_SRHSFUNC_FAIL);\r
         } \r
         if (retval > 0) {\r
@@ -2107,15 +2174,15 @@ public abstract class CVODES {
         }\r
       }\r
 \r
-      if (cv_mem->cv_quadr_sensi) {\r
-        retval = cv_mem->cv_fQS(cv_mem->cv_Ns, cv_mem->cv_tn, cv_mem->cv_zn[0],\r
-                                cv_mem->cv_znS[0], cv_mem->cv_znQ[1],\r
-                                cv_mem->cv_znQS[1], cv_mem->cv_fQS_data,\r
-                                cv_mem->cv_tempv, cv_mem->cv_tempvQ); \r
-        cv_mem->cv_nfQSe++;\r
+      if (cv_mem.cv_quadr_sensi) {\r
+        retval = fQS(cv_mem.cv_Ns, cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                                cv_mem.cv_znS[0], cv_mem.cv_znQ[1],\r
+                                cv_mem.cv_znQS[1], cv_mem.cv_fQS_data,\r
+                                cv_mem.cv_tempv, cv_mem.cv_tempvQ); \r
+        cv_mem.cv_nfQSe++;\r
         if (retval < 0) {\r
           cvProcessError(cv_mem, CV_QSRHSFUNC_FAIL, "CVODES", "CVode",\r
-                         MSGCV_QSRHSFUNC_FAILED, cv_mem->cv_tn);\r
+                         MSGCV_QSRHSFUNC_FAILED, cv_mem.cv_tn);\r
           return(CV_QSRHSFUNC_FAIL);\r
         } \r
         if (retval > 0) {\r
@@ -2123,47 +2190,47 @@ public abstract class CVODES {
                          "CVode", MSGCV_QSRHSFUNC_FIRST);\r
           return(CV_FIRST_QSRHSFUNC_ERR);\r
         }\r
-      } */\r
+      } \r
 \r
       /* Test input tstop for legality. */\r
 \r
-      /*if (cv_mem->cv_tstopset) {\r
-        if ( (cv_mem->cv_tstop - cv_mem->cv_tn)*(tout - cv_mem->cv_tn) <= ZERO ) {\r
+      if (cv_mem.cv_tstopset) {\r
+        if ( (cv_mem.cv_tstop - cv_mem.cv_tn)*(tout - cv_mem.cv_tn) <= ZERO ) {\r
           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                         MSGCV_BAD_TSTOP, cv_mem->cv_tstop, cv_mem->cv_tn);\r
+                         MSGCV_BAD_TSTOP, cv_mem.cv_tstop, cv_mem.cv_tn);\r
           return(CV_ILL_INPUT);\r
         }\r
-      }*/\r
+      }\r
 \r
       /* Set initial h (from H0 or cvHin). */\r
       \r
-      /*cv_mem->cv_h = cv_mem->cv_hin;\r
-      if ( (cv_mem->cv_h != ZERO) && ((tout-cv_mem->cv_tn)*cv_mem->cv_h < ZERO) ) {\r
+      cv_mem.cv_h = cv_mem.cv_hin;\r
+      if ( (cv_mem.cv_h != ZERO) && ((tout-cv_mem.cv_tn)*cv_mem.cv_h < ZERO) ) {\r
         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode", MSGCV_BAD_H0);\r
         return(CV_ILL_INPUT);\r
       }\r
-      if (cv_mem->cv_h == ZERO) {\r
+      if (cv_mem.cv_h == ZERO) {\r
         tout_hin = tout;\r
-        if ( cv_mem->cv_tstopset &&\r
-             (tout-cv_mem->cv_tn)*(tout-cv_mem->cv_tstop) > ZERO )\r
-          tout_hin = cv_mem->cv_tstop; \r
+        if ( cv_mem.cv_tstopset &&\r
+             (tout-cv_mem.cv_tn)*(tout-cv_mem.cv_tstop) > ZERO )\r
+          tout_hin = cv_mem.cv_tstop; \r
         hflag = cvHin(cv_mem, tout_hin);\r
         if (hflag != CV_SUCCESS) {\r
           istate = cvHandleFailure(cv_mem, hflag);\r
           return(istate);\r
         }\r
       }\r
-      rh = SUNRabs(cv_mem->cv_h)*cv_mem->cv_hmax_inv;\r
-      if (rh > ONE) cv_mem->cv_h /= rh;\r
-      if (SUNRabs(cv_mem->cv_h) < cv_mem->cv_hmin)\r
-        cv_mem->cv_h *= cv_mem->cv_hmin/SUNRabs(cv_mem->cv_h);*/\r
+      rh = Math.abs(cv_mem.cv_h)*cv_mem.cv_hmax_inv;\r
+      if (rh > ONE) cv_mem.cv_h /= rh;\r
+      if (Math.abs(cv_mem.cv_h) < cv_mem.cv_hmin)\r
+        cv_mem.cv_h *= cv_mem.cv_hmin/Math.abs(cv_mem.cv_h);\r
 \r
       /* Check for approach to tstop */\r
 \r
-      /*if (cv_mem->cv_tstopset) {\r
-        if ( (cv_mem->cv_tn + cv_mem->cv_h - cv_mem->cv_tstop)*cv_mem->cv_h > ZERO ) \r
-          cv_mem->cv_h = (cv_mem->cv_tstop - cv_mem->cv_tn)*(ONE-FOUR*cv_mem->cv_uround);\r
-      }*/\r
+      if (cv_mem.cv_tstopset) {\r
+        if ( (cv_mem.cv_tn + cv_mem.cv_h - cv_mem.cv_tstop)*cv_mem.cv_h > ZERO ) \r
+          cv_mem.cv_h = (cv_mem.cv_tstop - cv_mem.cv_tn)*(ONE-FOUR*cv_mem.cv_uround);\r
+      }\r
 \r
       /* \r
        * Scale zn[1] by h.\r
@@ -2172,38 +2239,38 @@ public abstract class CVODES {
        * If computing quadrature sensitivities,  scale znQS[1][is] by h. \r
        */\r
 \r
-      /*cv_mem->cv_hscale = cv_mem->cv_h;\r
-      cv_mem->cv_h0u    = cv_mem->cv_h;\r
-      cv_mem->cv_hprime = cv_mem->cv_h;\r
+      cv_mem.cv_hscale = cv_mem.cv_h;\r
+      cv_mem.cv_h0u    = cv_mem.cv_h;\r
+      cv_mem.cv_hprime = cv_mem.cv_h;\r
 \r
-      N_VScale(cv_mem->cv_h, cv_mem->cv_zn[1], cv_mem->cv_zn[1]);\r
+      N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_zn[1], cv_mem.cv_zn[1]);\r
       \r
-      if (cv_mem->cv_quadr)\r
-        N_VScale(cv_mem->cv_h, cv_mem->cv_znQ[1], cv_mem->cv_znQ[1]);\r
+      if (cv_mem.cv_quadr)\r
+        N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_znQ[1], cv_mem.cv_znQ[1]);\r
 \r
-      if (cv_mem->cv_sensi)\r
-        for (is=0; is<cv_mem->cv_Ns; is++) \r
-          N_VScale(cv_mem->cv_h, cv_mem->cv_znS[1][is], cv_mem->cv_znS[1][is]);\r
+      if (cv_mem.cv_sensi)\r
+        for (is=0; is<cv_mem.cv_Ns; is++) \r
+          N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_znS[1][is], cv_mem.cv_znS[1][is]);\r
 \r
-      if (cv_mem->cv_quadr_sensi)\r
-        for (is=0; is<cv_mem->cv_Ns; is++) \r
-          N_VScale(cv_mem->cv_h, cv_mem->cv_znQS[1][is], cv_mem->cv_znQS[1][is]);*/\r
+      if (cv_mem.cv_quadr_sensi)\r
+        for (is=0; is<cv_mem.cv_Ns; is++) \r
+          N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_znQS[1][is], cv_mem.cv_znQS[1][is]);\r
       \r
       /* Check for zeros of root function g at and near t0. */\r
 \r
-      /*if (cv_mem->cv_nrtfn > 0) {\r
+      if (cv_mem.cv_nrtfn > 0) {\r
 \r
         retval = cvRcheck1(cv_mem);\r
 \r
         if (retval == CV_RTFUNC_FAIL) {\r
           cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck1",\r
-                         MSGCV_RTFUNC_FAILED, cv_mem->cv_tn);\r
+                         MSGCV_RTFUNC_FAILED, cv_mem.cv_tn);\r
           return(CV_RTFUNC_FAIL);\r
         }\r
 \r
       }\r
 \r
-    }*/ /* end first call block */ // if (cv_mem.cv_nst == 0)\r
+    } /* end first call block */ // if (cv_mem.cv_nst == 0)\r
 \r
     /*\r
      * ------------------------------------------------------\r
@@ -2217,66 +2284,66 @@ public abstract class CVODES {
      * -------------------------------------------------------\r
      */\r
 \r
-   /* if (cv_mem->cv_nst > 0) {*/\r
+    if (cv_mem.cv_nst > 0) {\r
 \r
       /* Estimate an infinitesimal time interval to be used as\r
          a roundoff for time quantities (based on current time \r
          and step size) */\r
-      /*troundoff = FUZZ_FACTOR * cv_mem->cv_uround *\r
-        (SUNRabs(cv_mem->cv_tn) + SUNRabs(cv_mem->cv_h));*/\r
+      troundoff = FUZZ_FACTOR * cv_mem.cv_uround *\r
+        (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h));\r
 \r
       /* First check for a root in the last step taken, other than the\r
          last root found, if any.  If itask = CV_ONE_STEP and y(tn) was not\r
          returned because of an intervening root, return y(tn) now.     */\r
-      /*if (cv_mem->cv_nrtfn > 0) {\r
+      if (cv_mem.cv_nrtfn > 0) {\r
         \r
-        irfndp = cv_mem->cv_irfnd;\r
+        irfndp = cv_mem.cv_irfnd;\r
         \r
         retval = cvRcheck2(cv_mem);\r
 \r
         if (retval == CLOSERT) {\r
           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvRcheck2",\r
-                         MSGCV_CLOSE_ROOTS, cv_mem->cv_tlo);\r
+                         MSGCV_CLOSE_ROOTS, cv_mem.cv_tlo);\r
           return(CV_ILL_INPUT);\r
         } else if (retval == CV_RTFUNC_FAIL) {\r
           cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck2",\r
-                         MSGCV_RTFUNC_FAILED, cv_mem->cv_tlo);\r
+                         MSGCV_RTFUNC_FAILED, cv_mem.cv_tlo);\r
           return(CV_RTFUNC_FAIL);\r
         } else if (retval == RTFOUND) {\r
-          cv_mem->cv_tretlast = *tret = cv_mem->cv_tlo;\r
+          cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tlo;\r
           return(CV_ROOT_RETURN);\r
-        }*/\r
+        }\r
         \r
         /* If tn is distinct from tretlast (within roundoff),\r
            check remaining interval for roots */\r
-        /*if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tretlast) > troundoff ) {\r
+        if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tretlast) > troundoff ) {\r
 \r
           retval = cvRcheck3(cv_mem);\r
 \r
-          if (retval == CV_SUCCESS) {*/     /* no root found */\r
-            /*cv_mem->cv_irfnd = 0;\r
+          if (retval == CV_SUCCESS) {     /* no root found */\r
+            cv_mem.cv_irfnd = 0;\r
             if ((irfndp == 1) && (itask == CV_ONE_STEP)) {\r
-              cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-              N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+              cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+              N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
               return(CV_SUCCESS);\r
             }\r
-          } else if (retval == RTFOUND) {*/  /* a new root was found */\r
-            /*cv_mem->cv_irfnd = 1;\r
-            cv_mem->cv_tretlast = *tret = cv_mem->cv_tlo;\r
+          } else if (retval == RTFOUND) {  /* a new root was found */\r
+            cv_mem.cv_irfnd = 1;\r
+            cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tlo;\r
             return(CV_ROOT_RETURN);\r
-          } else if (retval == CV_RTFUNC_FAIL) { */ /* g failed */\r
-            /*cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3", \r
-                           MSGCV_RTFUNC_FAILED, cv_mem->cv_tlo);\r
+          } else if (retval == CV_RTFUNC_FAIL) {  /* g failed */\r
+            cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3", \r
+                           MSGCV_RTFUNC_FAILED, cv_mem.cv_tlo);\r
             return(CV_RTFUNC_FAIL);\r
           }\r
 \r
         }\r
         \r
-      } *//* end of root stop check */\r
+      } /* end of root stop check */\r
       \r
       /* In CV_NORMAL mode, test if tout was reached */\r
-     /* if ( (itask == CV_NORMAL) && ((cv_mem->cv_tn-tout)*cv_mem->cv_h >= ZERO) ) {\r
-        cv_mem->cv_tretlast = *tret = tout;\r
+      if ( (itask == CV_NORMAL) && ((cv_mem.cv_tn-tout)*cv_mem.cv_h >= ZERO) ) {\r
+        cv_mem.cv_tretlast = tret[0] = tout;\r
         ier =  CVodeGetDky(cv_mem, tout, 0, yout);\r
         if (ier != CV_SUCCESS) {\r
           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
@@ -2284,40 +2351,40 @@ public abstract class CVODES {
           return(CV_ILL_INPUT);\r
         }\r
         return(CV_SUCCESS);\r
-      }*/\r
+      }\r
       \r
       /* In CV_ONE_STEP mode, test if tn was returned */\r
-     /* if ( itask == CV_ONE_STEP &&\r
-           SUNRabs(cv_mem->cv_tn - cv_mem->cv_tretlast) > troundoff ) {\r
-        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+      if ( itask == CV_ONE_STEP &&\r
+           Math.abs(cv_mem.cv_tn - cv_mem.cv_tretlast) > troundoff ) {\r
+        cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+        N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
         return(CV_SUCCESS);\r
       }\r
       \r
       /* Test for tn at tstop or near tstop */\r
-     /* if ( cv_mem->cv_tstopset ) {\r
+      if ( cv_mem.cv_tstopset ) {\r
         \r
-        if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tstop) <= troundoff ) {\r
-          ier =  CVodeGetDky(cv_mem, cv_mem->cv_tstop, 0, yout);\r
+        if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tstop) <= troundoff ) {\r
+          ier =  CVodeGetDky(cv_mem, cv_mem.cv_tstop, 0, yout);\r
           if (ier != CV_SUCCESS) {\r
             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_BAD_TSTOP, cv_mem->cv_tstop, cv_mem->cv_tn);\r
+                           MSGCV_BAD_TSTOP, cv_mem.cv_tstop, cv_mem.cv_tn);\r
             return(CV_ILL_INPUT);\r
           }\r
-          cv_mem->cv_tretlast = *tret = cv_mem->cv_tstop;\r
-          cv_mem->cv_tstopset = SUNFALSE;\r
+          cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tstop;\r
+          cv_mem.cv_tstopset = false;\r
           return(CV_TSTOP_RETURN);\r
-        }*/\r
+        }\r
         \r
         /* If next step would overtake tstop, adjust stepsize */\r
-        /*if ( (cv_mem->cv_tn + cv_mem->cv_hprime - cv_mem->cv_tstop)*cv_mem->cv_h > ZERO ) {\r
-          cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn)*(ONE-FOUR*cv_mem->cv_uround);\r
-          cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h;\r
+        if ( (cv_mem.cv_tn + cv_mem.cv_hprime - cv_mem.cv_tstop)*cv_mem.cv_h > ZERO ) {\r
+          cv_mem.cv_hprime = (cv_mem.cv_tstop - cv_mem.cv_tn)*(ONE-FOUR*cv_mem.cv_uround);\r
+          cv_mem.cv_eta = cv_mem.cv_hprime / cv_mem.cv_h;\r
         }\r
         \r
       }\r
       \r
-    }*/ /* end stopping tests block at nst>0 */  \r
+    } /* end stopping tests block at nst>0 */  \r
     \r
     /*\r
      * --------------------------------------------------\r
@@ -2335,224 +2402,223 @@ public abstract class CVODES {
      * --------------------------------------------------\r
      */  \r
     \r
-    /*nstloc = 0;\r
+    nstloc = 0;\r
     for(;;) {\r
      \r
-      cv_mem->cv_next_h = cv_mem->cv_h;\r
-      cv_mem->cv_next_q = cv_mem->cv_q;*/\r
+      cv_mem.cv_next_h = cv_mem.cv_h;\r
+      cv_mem.cv_next_q = cv_mem.cv_q;\r
       \r
       /* Reset and check ewt, ewtQ, ewtS */   \r
-      /*if (cv_mem->cv_nst > 0) {\r
+      if (cv_mem.cv_nst > 0) {\r
 \r
-        ier = cv_mem->cv_efun(cv_mem->cv_zn[0], cv_mem->cv_ewt, cv_mem->cv_e_data);\r
+        ier = cv_efun(cv_mem.cv_zn[0], cv_mem.cv_ewt, cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
         if(ier != 0) {\r
-          if (cv_mem->cv_itol == CV_WF)\r
+          if (cv_mem.cv_itol == CV_WF)\r
             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_EWT_NOW_FAIL, cv_mem->cv_tn);\r
+                           MSGCV_EWT_NOW_FAIL, cv_mem.cv_tn);\r
           else\r
             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_EWT_NOW_BAD, cv_mem->cv_tn);\r
+                           MSGCV_EWT_NOW_BAD, cv_mem.cv_tn);\r
           istate = CV_ILL_INPUT;\r
-          cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-          N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+          cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+          N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
           break;\r
         }\r
 \r
-        if (cv_mem->cv_quadr && cv_mem->cv_errconQ) {\r
-          ier = cvQuadEwtSet(cv_mem, cv_mem->cv_znQ[0], cv_mem->cv_ewtQ);\r
+        if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
+          ier = cvQuadEwtSet(cv_mem, cv_mem.cv_znQ[0], cv_mem.cv_ewtQ);\r
           if(ier != 0) {\r
             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_EWTQ_NOW_BAD, cv_mem->cv_tn);\r
+                           MSGCV_EWTQ_NOW_BAD, cv_mem.cv_tn);\r
             istate = CV_ILL_INPUT;\r
-            cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-            N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+            cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+            N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
             break;\r
           }\r
         }\r
 \r
-        if (cv_mem->cv_sensi) {\r
-          ier = cvSensEwtSet(cv_mem, cv_mem->cv_znS[0], cv_mem->cv_ewtS);\r
+        if (cv_mem.cv_sensi) {\r
+          ier = cvSensEwtSet(cv_mem, cv_mem.cv_znS[0], cv_mem.cv_ewtS);\r
           if (ier != 0) {\r
             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_EWTS_NOW_BAD, cv_mem->cv_tn);\r
+                           MSGCV_EWTS_NOW_BAD, cv_mem.cv_tn);\r
             istate = CV_ILL_INPUT;\r
-            cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-            N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+            cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+            N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
             break;\r
           }\r
         }\r
 \r
-        if (cv_mem->cv_quadr_sensi && cv_mem->cv_errconQS) {\r
-          ier = cvQuadSensEwtSet(cv_mem, cv_mem->cv_znQS[0], cv_mem->cv_ewtQS);\r
+        if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS) {\r
+          ier = cvQuadSensEwtSet(cv_mem, cv_mem.cv_znQS[0], cv_mem.cv_ewtQS);\r
           if (ier != 0) {\r
             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_EWTQS_NOW_BAD, cv_mem->cv_tn);\r
+                           MSGCV_EWTQS_NOW_BAD, cv_mem.cv_tn);\r
             istate = CV_ILL_INPUT;\r
-            cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-            N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+            cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+            N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
             break;\r
           }\r
         }\r
 \r
-      }*/\r
+      }\r
 \r
       /* Check for too many steps */\r
-      /*if ( (cv_mem->cv_mxstep>0) && (nstloc >= cv_mem->cv_mxstep) ) {\r
+      if ( (cv_mem.cv_mxstep>0) && (nstloc >= cv_mem.cv_mxstep) ) {\r
         cvProcessError(cv_mem, CV_TOO_MUCH_WORK, "CVODES", "CVode",\r
-                       MSGCV_MAX_STEPS, cv_mem->cv_tn);\r
+                       MSGCV_MAX_STEPS, cv_mem.cv_tn);\r
         istate = CV_TOO_MUCH_WORK;\r
-        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+        cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+        N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
         break;\r
-      }*/\r
+      }\r
 \r
       /* Check for too much accuracy requested */\r
-      /*nrm = N_VWrmsNorm(cv_mem->cv_zn[0], cv_mem->cv_ewt);\r
-      if (cv_mem->cv_quadr && cv_mem->cv_errconQ) {\r
-        nrm = cvQuadUpdateNorm(cv_mem, nrm, cv_mem->cv_znQ[0], cv_mem->cv_ewtQ); \r
+      nrm = N_VWrmsNorm_Serial(cv_mem.cv_zn[0], cv_mem.cv_ewt);\r
+      if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
+        nrm = cvQuadUpdateNorm(cv_mem, nrm, cv_mem.cv_znQ[0], cv_mem.cv_ewtQ); \r
       }\r
-      if (cv_mem->cv_sensi && cv_mem->cv_errconS) {\r
-        nrm = cvSensUpdateNorm(cv_mem, nrm, cv_mem->cv_znS[0], cv_mem->cv_ewtS);\r
+      if (cv_mem.cv_sensi && cv_mem.cv_errconS) {\r
+        nrm = cvSensUpdateNorm(cv_mem, nrm, cv_mem.cv_znS[0], cv_mem.cv_ewtS);\r
       }\r
-      if (cv_mem->cv_quadr_sensi && cv_mem->cv_errconQS) {\r
-        nrm = cvQuadSensUpdateNorm(cv_mem, nrm, cv_mem->cv_znQS[0], cv_mem->cv_ewtQS);\r
+      if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS) {\r
+        nrm = cvQuadSensUpdateNorm(cv_mem, nrm, cv_mem.cv_znQS[0], cv_mem.cv_ewtQS);\r
       }\r
-      cv_mem->cv_tolsf = cv_mem->cv_uround * nrm;\r
-      if (cv_mem->cv_tolsf > ONE) {\r
+      cv_mem.cv_tolsf = cv_mem.cv_uround * nrm;\r
+      if (cv_mem.cv_tolsf > ONE) {\r
         cvProcessError(cv_mem, CV_TOO_MUCH_ACC, "CVODES", "CVode",\r
-                       MSGCV_TOO_MUCH_ACC, cv_mem->cv_tn);\r
+                       MSGCV_TOO_MUCH_ACC, cv_mem.cv_tn);\r
         istate = CV_TOO_MUCH_ACC;\r
-        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
-        cv_mem->cv_tolsf *= TWO;\r
+        cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+        N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
+        cv_mem.cv_tolsf *= TWO;\r
         break;\r
       } else {\r
-        cv_mem->cv_tolsf = ONE;\r
-      }*/\r
+        cv_mem.cv_tolsf = ONE;\r
+      }\r
       \r
       /* Check for h below roundoff level in tn */\r
-     /* if (cv_mem->cv_tn + cv_mem->cv_h == cv_mem->cv_tn) {\r
-        cv_mem->cv_nhnil++;\r
-        if (cv_mem->cv_nhnil <= cv_mem->cv_mxhnil) \r
+      if (cv_mem.cv_tn + cv_mem.cv_h == cv_mem.cv_tn) {\r
+        cv_mem.cv_nhnil++;\r
+        if (cv_mem.cv_nhnil <= cv_mem.cv_mxhnil) \r
           cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode", MSGCV_HNIL,\r
-                         cv_mem->cv_tn, cv_mem->cv_h);\r
-        if (cv_mem->cv_nhnil == cv_mem->cv_mxhnil) \r
+                         cv_mem.cv_tn, cv_mem.cv_h);\r
+        if (cv_mem.cv_nhnil == cv_mem.cv_mxhnil) \r
           cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode", MSGCV_HNIL_DONE);\r
-      }*/\r
+      }\r
 \r
       /* Call cvStep to take a step */\r
-      /*kflag = cvStep(cv_mem); */\r
+      kflag = cvStep(cv_mem);\r
 \r
       /* Process failed step cases, and exit loop */\r
-     /* if (kflag != CV_SUCCESS) {\r
+      if (kflag != CV_SUCCESS) {\r
         istate = cvHandleFailure(cv_mem, kflag);\r
-        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
+        cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+        N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
         break;\r
       }\r
       \r
-      nstloc++;*/\r
+      nstloc++;\r
 \r
       /* If tstop is set and was reached, reset tn = tstop */\r
-      /*if ( cv_mem->cv_tstopset ) {\r
-        troundoff = FUZZ_FACTOR * cv_mem->cv_uround *\r
-          (SUNRabs(cv_mem->cv_tn) + SUNRabs(cv_mem->cv_h));\r
-        if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tstop) <= troundoff)\r
-          cv_mem->cv_tn = cv_mem->cv_tstop;\r
-      }*/\r
+      if ( cv_mem.cv_tstopset ) {\r
+        troundoff = FUZZ_FACTOR * cv_mem.cv_uround *\r
+          (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h));\r
+        if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tstop) <= troundoff)\r
+          cv_mem.cv_tn = cv_mem.cv_tstop;\r
+      }\r
 \r
       /* Check for root in last step taken. */    \r
-      /*if (cv_mem->cv_nrtfn > 0) {\r
+      if (cv_mem.cv_nrtfn > 0) {\r
         \r
         retval = cvRcheck3(cv_mem);\r
         \r
-        if (retval == RTFOUND) { */ /* A new root was found */\r
-         /* cv_mem->cv_irfnd = 1;\r
+        if (retval == RTFOUND) {  /* A new root was found */\r
+          cv_mem.cv_irfnd = 1;\r
           istate = CV_ROOT_RETURN;\r
-          cv_mem->cv_tretlast = *tret = cv_mem->cv_tlo;\r
+          cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tlo;\r
           break;\r
-        } else if (retval == CV_RTFUNC_FAIL) {*/ /* g failed */\r
-         /* cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3",\r
-                         MSGCV_RTFUNC_FAILED, cv_mem->cv_tlo);\r
+        } else if (retval == CV_RTFUNC_FAIL) { /* g failed */\r
+          cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3",\r
+                         MSGCV_RTFUNC_FAILED, cv_mem.cv_tlo);\r
           istate = CV_RTFUNC_FAIL;\r
           break;\r
-        }*/\r
+        }\r
 \r
         /* If we are at the end of the first step and we still have\r
          * some event functions that are inactive, issue a warning\r
          * as this may indicate a user error in the implementation\r
          * of the root function. */\r
 \r
-       /* if (cv_mem->cv_nst==1) {\r
-          inactive_roots = SUNFALSE;\r
-          for (ir=0; ir<cv_mem->cv_nrtfn; ir++) { \r
-            if (!cv_mem->cv_gactive[ir]) {\r
-              inactive_roots = SUNTRUE;\r
+        if (cv_mem.cv_nst==1) {\r
+          inactive_roots = false;\r
+          for (ir=0; ir<cv_mem.cv_nrtfn; ir++) { \r
+            if (!cv_mem.cv_gactive[ir]) {\r
+              inactive_roots = true;\r
               break;\r
             }\r
           }\r
-          if ((cv_mem->cv_mxgnull > 0) && inactive_roots) {\r
+          if ((cv_mem.cv_mxgnull > 0) && inactive_roots) {\r
             cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode",\r
                            MSGCV_INACTIVE_ROOTS);\r
           }\r
         }\r
 \r
-      }*/\r
+      }\r
 \r
       /* In NORMAL mode, check if tout reached */\r
-     /* if ( (itask == CV_NORMAL) &&  (cv_mem->cv_tn-tout)*cv_mem->cv_h >= ZERO ) {\r
+      if ( (itask == CV_NORMAL) &&  (cv_mem.cv_tn-tout)*cv_mem.cv_h >= ZERO ) {\r
         istate = CV_SUCCESS;\r
-        cv_mem->cv_tretlast = *tret = tout;\r
-        (void) CVodeGetDky(cv_mem, tout, 0, yout);\r
-        cv_mem->cv_next_q = cv_mem->cv_qprime;\r
-        cv_mem->cv_next_h = cv_mem->cv_hprime;\r
+        cv_mem.cv_tretlast = tret[0] = tout;\r
+        CVodeGetDky(cv_mem, tout, 0, yout);\r
+        cv_mem.cv_next_q = cv_mem.cv_qprime;\r
+        cv_mem.cv_next_h = cv_mem.cv_hprime;\r
         break;\r
-      }*/\r
+      }\r
 \r
       /* Check if tn is at tstop, or about to pass tstop */\r
-     /* if ( cv_mem->cv_tstopset ) {\r
-\r
-        troundoff = FUZZ_FACTOR * cv_mem->cv_uround *\r
-          (SUNRabs(cv_mem->cv_tn) + SUNRabs(cv_mem->cv_h));\r
-        if ( SUNRabs(cv_mem->cv_tn - cv_mem->cv_tstop) <= troundoff) {\r
-          (void) CVodeGetDky(cv_mem, cv_mem->cv_tstop, 0, yout);\r
-          cv_mem->cv_tretlast = *tret = cv_mem->cv_tstop;\r
-          cv_mem->cv_tstopset = SUNFALSE;\r
+      if ( cv_mem.cv_tstopset ) {\r
+\r
+        troundoff = FUZZ_FACTOR * cv_mem.cv_uround *\r
+          (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h));\r
+        if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tstop) <= troundoff) {\r
+          CVodeGetDky(cv_mem, cv_mem.cv_tstop, 0, yout);\r
+          cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tstop;\r
+          cv_mem.cv_tstopset = false;\r
           istate = CV_TSTOP_RETURN;\r
           break;\r
         }\r
 \r
-        if ( (cv_mem->cv_tn + cv_mem->cv_hprime - cv_mem->cv_tstop)*cv_mem->cv_h > ZERO ) {\r
-          cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn)*(ONE-FOUR*cv_mem->cv_uround);\r
-          cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h;\r
+        if ( (cv_mem.cv_tn + cv_mem.cv_hprime - cv_mem.cv_tstop)*cv_mem.cv_h > ZERO ) {\r
+          cv_mem.cv_hprime = (cv_mem.cv_tstop - cv_mem.cv_tn)*(ONE-FOUR*cv_mem.cv_uround);\r
+          cv_mem.cv_eta = cv_mem.cv_hprime / cv_mem.cv_h;\r
         }\r
 \r
-      }*/\r
+      }\r
 \r
       /* In ONE_STEP mode, copy y and exit loop */\r
-      /*if (itask == CV_ONE_STEP) {\r
+      if (itask == CV_ONE_STEP) {\r
         istate = CV_SUCCESS;\r
-        cv_mem->cv_tretlast = *tret = cv_mem->cv_tn;\r
-        N_VScale(ONE, cv_mem->cv_zn[0], yout);\r
-        cv_mem->cv_next_q = cv_mem->cv_qprime;\r
-        cv_mem->cv_next_h = cv_mem->cv_hprime;\r
+        cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+        N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
+        cv_mem.cv_next_q = cv_mem.cv_qprime;\r
+        cv_mem.cv_next_h = cv_mem.cv_hprime;\r
         break;\r
       }\r
 \r
-    } *//* end looping for internal steps */\r
+    } /* end looping for internal steps */\r
     \r
     /* Load optional output */\r
-   /* if (cv_mem->cv_sensi && (cv_mem->cv_ism==CV_STAGGERED1)) { \r
-      cv_mem->cv_nniS  = 0;\r
-      cv_mem->cv_ncfnS = 0;\r
-      for (is=0; is<cv_mem->cv_Ns; is++) {\r
-        cv_mem->cv_nniS  += cv_mem->cv_nniS1[is];\r
-        cv_mem->cv_ncfnS += cv_mem->cv_ncfnS1[is];\r
+    if (cv_mem.cv_sensi && (cv_mem.cv_ism==CV_STAGGERED1)) { \r
+      cv_mem.cv_nniS  = 0;\r
+      cv_mem.cv_ncfnS = 0;\r
+      for (is=0; is<cv_mem.cv_Ns; is++) {\r
+        cv_mem.cv_nniS  += cv_mem.cv_nniS1[is];\r
+        cv_mem.cv_ncfnS += cv_mem.cv_ncfnS1[is];\r
       }\r
-    }*/\r
+    }\r
     \r
-    /*return(istate);*/\r
-    return -100;\r
+    return(istate);\r
 \r
   }\r
   \r
@@ -3143,6 +3209,81 @@ public abstract class CVODES {
 \r
     return;\r
   }\r
+  \r
+  private void N_VDiv_Serial(NVector x, NVector y, NVector z)\r
+  {\r
+         int i, N;\r
+         double xd[], yd[], zd[];\r
+\r
+         xd = yd = zd = null;\r
+\r
+         xd = x.data;\r
+         yd = y.data;\r
+         zd = z.data;\r
+         N = xd.length;\r
+\r
+    for (i = 0; i < N; i++)\r
+      zd[i] = xd[i]/yd[i];\r
+\r
+    return;\r
+  }\r
+  \r
+  private double N_VMaxNorm_Serial(NVector x)\r
+  {\r
+    int i, N;\r
+    double max, xd[];\r
+\r
+    max = ZERO;\r
+    xd = null;\r
+\r
+    xd = x.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++) {\r
+      if (Math.abs(xd[i]) > max) max = Math.abs(xd[i]);\r
+    }\r
+\r
+    return(max);\r
+  }\r
+\r
+  \r
+  private double N_VWrmsNorm_Serial(NVector x, NVector w)\r
+  {\r
+    int i, N;\r
+    double sum, prodi, xd[],wd[];\r
+\r
+    sum = ZERO;\r
+    xd = wd = null;\r
+\r
+    xd = x.data;\r
+    wd = w.data;\r
+    N = xd.length;\r
+\r
+    for (i = 0; i < N; i++) {\r
+      prodi = xd[i]*wd[i];\r
+      sum += (prodi*prodi);\r
+    }\r
+\r
+    return(Math.sqrt(sum/N));\r
+  }\r
+  \r
+  private void N_VConst_Serial(double c, NVector z)\r
+  {\r
+    int i, N;\r
+    double zd[];\r
+\r
+    zd = null;\r
+\r
+    zd = z.data;\r
+    N = zd.length;\r
+\r
+    for (i = 0; i < N; i++) zd[i] = c;\r
+\r
+    return;\r
+  }\r
+\r
+\r
+\r
 \r
   private int cv_efun(NVector ycur, NVector weight, CVodeMemRec cv_mem, int cv_efun_select) {\r
          int flag = 0;\r
@@ -3486,5 +3627,2289 @@ public abstract class CVODES {
    return(CVDLS_SUCCESS);\r
  }\r
 \r
+ /*\r
+  * cvSensRhsWrapper\r
+  *\r
+  * CVSensRhs is a high level routine that returns right hand side \r
+  * of sensitivity equations. Depending on the 'ifS' flag, it either \r
+  * calls directly the fS routine (ifS=CV_ALLSENS) or (if ifS=CV_ONESENS) \r
+  * calls the fS1 routine in a loop over all sensitivities.\r
+  *\r
+  * CVSensRhs is called:\r
+  *  (*) by CVode at the first step\r
+  *  (*) by cvYddNorm if errcon=SUNTRUE\r
+  *  (*) by cvNlsFunctional, cvNlsNewton, and cvNewtonIteration\r
+  *      if ism=CV_SIMULTANEOUS\r
+  *  (*) by cvDoErrorTest when restarting from scratch\r
+  *  (*) in the corrector loop if ism=CV_STAGGERED\r
+  *  (*) by cvStgrDoErrorTest when restarting from scratch \r
+  *\r
+  * The return value is that of the sensitivity RHS function fS,\r
+  *\r
+  */\r
+\r
+ int cvSensRhsWrapper(CVodeMemRec cv_mem, double time, \r
+                      NVector ycur, NVector fcur, \r
+                      NVector yScur[], NVector fScur[],\r
+                      NVector temp1, NVector temp2)\r
+ {\r
+   int retval=0, is;\r
+\r
+   if (cv_mem.cv_ifS==CV_ALLSENS) {\r
+     retval = fS(cv_mem.cv_Ns, time, ycur, fcur, yScur, \r
+                            fScur, cv_mem.cv_fS_data, temp1, temp2);\r
+     cv_mem.cv_nfSe++;\r
+   } else {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+       retval = fS1(cv_mem.cv_Ns, time, ycur, fcur, is, yScur[is], \r
+                               fScur[is], cv_mem.cv_fS_data, temp1, temp2);\r
+       cv_mem.cv_nfSe++;\r
+       if (retval != 0) break;\r
+     }\r
+   }\r
+\r
+   return(retval);\r
+ }\r
+\r
\r
+ /*\r
+  * cvHin\r
+  *\r
+  * This routine computes a tentative initial step size h0. \r
+  * If tout is too close to tn (= t0), then cvHin returns CV_TOO_CLOSE\r
+  * and h remains uninitialized. Note that here tout is either the value\r
+  * passed to CVode at the first call or the value of tstop (if tstop is \r
+  * enabled and it is closer to t0=tn than tout).\r
+  * If any RHS function fails unrecoverably, cvHin returns CV_*RHSFUNC_FAIL.\r
+  * If any RHS function fails recoverably too many times and recovery is\r
+  * not possible, cvHin returns CV_REPTD_*RHSFUNC_ERR.\r
+  * Otherwise, cvHin sets h to the chosen value h0 and returns CV_SUCCESS.\r
+  *\r
+  * The algorithm used seeks to find h0 as a solution of\r
+  *       (WRMS norm of (h0^2 ydd / 2)) = 1, \r
+  * where ydd = estimated second derivative of y. Here, y includes\r
+  * all variables considered in the error test.\r
+  *\r
+  * We start with an initial estimate equal to the geometric mean of the\r
+  * lower and upper bounds on the step size.\r
+  *\r
+  * Loop up to MAX_ITERS times to find h0.\r
+  * Stop if new and previous values differ by a factor < 2.\r
+  * Stop if hnew/hg > 2 after one iteration, as this probably means\r
+  * that the ydd value is bad because of cancellation error.        \r
+  *  \r
+  * For each new proposed hg, we allow MAX_ITERS attempts to\r
+  * resolve a possible recoverable failure from f() by reducing\r
+  * the proposed stepsize by a factor of 0.2. If a legal stepsize\r
+  * still cannot be found, fall back on a previous value if possible,\r
+  * or else return CV_REPTD_RHSFUNC_ERR.\r
+  *\r
+  * Finally, we apply a bias (0.5) and verify that h0 is within bounds.\r
+  */\r
+\r
+ private int cvHin(CVodeMemRec cv_mem, double tout)\r
+ {\r
+   int retval = 0;\r
+   int sign, count1, count2;\r
+   double tdiff, tdist, tround, hlb, hub;\r
+   double hnew = 0.0;\r
+   double hg, hgs, hs, hrat, h0;\r
+   double yddnrm[] = new double[1];\r
+   boolean hgOK, hnewOK;\r
+\r
+   /* If tout is too close to tn, give up */\r
+   \r
+   if ((tdiff = tout-cv_mem.cv_tn) == ZERO) return(CV_TOO_CLOSE);\r
+   \r
+   sign = (tdiff > ZERO) ? 1 : -1;\r
+   tdist = Math.abs(tdiff);\r
+   tround = cv_mem.cv_uround * Math.max(Math.abs(cv_mem.cv_tn), Math.abs(tout));\r
+\r
+   if (tdist < TWO*tround) return(CV_TOO_CLOSE);\r
+   \r
+   /* \r
+      Set lower and upper bounds on h0, and take geometric mean \r
+      as first trial value.\r
+      Exit with this value if the bounds cross each other.\r
+   */\r
+\r
+   hlb = HLB_FACTOR * tround;\r
+   hub = cvUpperBoundH0(cv_mem, tdist);\r
+\r
+   hg  = Math.sqrt(hlb*hub);\r
+\r
+   if (hub < hlb) {\r
+     if (sign == -1) cv_mem.cv_h = -hg;\r
+     else            cv_mem.cv_h =  hg;\r
+     return(CV_SUCCESS);\r
+   }\r
+   \r
+   /* Outer loop */\r
+\r
+   hnewOK = false;\r
+   hs = hg;         /* safeguard against 'uninitialized variable' warning */\r
+\r
+   for(count1 = 1; count1 <= MAX_ITERS; count1++) {\r
+\r
+     /* Attempts to estimate ydd */\r
+\r
+     hgOK = false;\r
+\r
+     for (count2 = 1; count2 <= MAX_ITERS; count2++) {\r
+       hgs = hg*sign;\r
+       retval = cvYddNorm(cv_mem, hgs, yddnrm);\r
+       /* If a RHS function failed unrecoverably, give up */\r
+       if (retval < 0) return(retval);\r
+       /* If successful, we can use ydd */\r
+       if (retval == CV_SUCCESS) {hgOK = true; break;}\r
+       /* A RHS function failed recoverably; cut step size and test it again */\r
+       hg *= POINT2;\r
+     }\r
+\r
+     /* If a RHS function failed recoverably MAX_ITERS times */\r
+\r
+     if (!hgOK) {\r
+       /* Exit if this is the first or second pass. No recovery possible */\r
+       if (count1 <= 2) \r
+         if (retval == RHSFUNC_RECVR)  return(CV_REPTD_RHSFUNC_ERR);\r
+         if (retval == QRHSFUNC_RECVR) return(CV_REPTD_QRHSFUNC_ERR);\r
+         if (retval == SRHSFUNC_RECVR) return(CV_REPTD_SRHSFUNC_ERR);\r
+       /* We have a fall-back option. The value hs is a previous hnew which\r
+          passed through f(). Use it and break */\r
+       hnew = hs;\r
+       break;\r
+     }\r
+\r
+     /* The proposed step size is feasible. Save it. */\r
+     hs = hg;\r
+\r
+     /* If the stopping criteria was met, or if this is the last pass, stop */\r
+     if ( (hnewOK) || (count1 == MAX_ITERS))  {hnew = hg; break;}\r
+\r
+     /* Propose new step size */\r
+     hnew = (yddnrm[0]*hub*hub > TWO) ? Math.sqrt(TWO/yddnrm[0]) : Math.sqrt(hg*hub);\r
+     hrat = hnew/hg;\r
+     \r
+     /* Accept hnew if it does not differ from hg by more than a factor of 2 */\r
+     if ((hrat > HALF) && (hrat < TWO)) {\r
+       hnewOK = true;\r
+     }\r
+\r
+     /* After one pass, if ydd seems to be bad, use fall-back value. */\r
+     if ((count1 > 1) && (hrat > TWO)) {\r
+       hnew = hg;\r
+       hnewOK = true;\r
+     }\r
+\r
+     /* Send this value back through f() */\r
+     hg = hnew;\r
+\r
+   }\r
+\r
+   /* Apply bounds, bias factor, and attach sign */\r
+\r
+   h0 = H_BIAS*hnew;\r
+   if (h0 < hlb) h0 = hlb;\r
+   if (h0 > hub) h0 = hub;\r
+   if (sign == -1) h0 = -h0;\r
+   cv_mem.cv_h = h0;\r
+\r
+   return(CV_SUCCESS);\r
+ }\r
\r
+ /*\r
+  * cvUpperBoundH0\r
+  *\r
+  * This routine sets an upper bound on abs(h0) based on\r
+  * tdist = tn - t0 and the values of y[i]/y'[i].\r
+  */\r
+\r
+ private double cvUpperBoundH0(CVodeMemRec cv_mem, double tdist)\r
+ {\r
+   double hub_inv, hubQ_inv, hubS_inv, hubQS_inv, hub;\r
+   NVector temp1, temp2;\r
+   NVector tempQ1, tempQ2;\r
+   NVector tempS1[];\r
+   NVector tempQS1[];\r
+   int is;\r
+\r
+   /* \r
+    * Bound based on |y|/|y'| -- allow at most an increase of\r
+    * HUB_FACTOR in y0 (based on a forward Euler step). The weight \r
+    * factor is used as a safeguard against zero components in y0. \r
+    */\r
+\r
+   temp1 = cv_mem.cv_tempv;\r
+   temp2 = cv_mem.cv_acor;\r
+\r
+   N_VAbs_Serial(cv_mem.cv_zn[0], temp2);\r
+   cv_efun(cv_mem.cv_zn[0], temp1, cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
+   N_VInv_Serial(temp1, temp1);\r
+   N_VLinearSum_Serial(HUB_FACTOR, temp2, ONE, temp1, temp1);\r
+\r
+   N_VAbs_Serial(cv_mem.cv_zn[1], temp2);\r
+\r
+   N_VDiv_Serial(temp2, temp1, temp1);\r
+   hub_inv = N_VMaxNorm_Serial(temp1);\r
+\r
+   /* Bound based on |yQ|/|yQ'| */\r
+   \r
+   if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
+\r
+     tempQ1 = cv_mem.cv_tempvQ;\r
+     tempQ2 = cv_mem.cv_acorQ;\r
+\r
+     N_VAbs_Serial(cv_mem.cv_znQ[0], tempQ2);\r
+     cvQuadEwtSet(cv_mem, cv_mem.cv_znQ[0], tempQ1);\r
+     N_VInv_Serial(tempQ1, tempQ1);\r
+     N_VLinearSum_Serial(HUB_FACTOR, tempQ2, ONE, tempQ1, tempQ1);\r
+     \r
+     N_VAbs_Serial(cv_mem.cv_znQ[1], tempQ2);\r
+     \r
+     N_VDiv_Serial(tempQ2, tempQ1, tempQ1);\r
+     hubQ_inv = N_VMaxNorm_Serial(tempQ1);\r
+\r
+     if (hubQ_inv > hub_inv) hub_inv = hubQ_inv;\r
+\r
+   }\r
+\r
+   /* Bound based on |yS|/|yS'| */\r
+\r
+   if (cv_mem.cv_sensi && cv_mem.cv_errconS) {\r
+\r
+     tempS1 = cv_mem.cv_acorS;\r
+     cvSensEwtSet(cv_mem, cv_mem.cv_znS[0], tempS1);\r
+\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+\r
+       N_VAbs_Serial(cv_mem.cv_znS[0][is], temp2);\r
+       N_VInv_Serial(tempS1[is], temp1);\r
+       N_VLinearSum_Serial(HUB_FACTOR, temp2, ONE, temp1, temp1);\r
+       \r
+       N_VAbs_Serial(cv_mem.cv_znS[1][is], temp2);\r
+       \r
+       N_VDiv_Serial(temp2, temp1, temp1);\r
+       hubS_inv = N_VMaxNorm_Serial(temp1);\r
+\r
+       if (hubS_inv > hub_inv) hub_inv = hubS_inv;\r
+\r
+     }\r
+\r
+   }\r
+\r
+   /* Bound based on |yQS|/|yQS'| */\r
+\r
+   if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS) {\r
+\r
+     tempQ1 = cv_mem.cv_tempvQ;\r
+     tempQ2 = cv_mem.cv_acorQ;\r
+\r
+     tempQS1 = cv_mem.cv_acorQS;\r
+     cvQuadSensEwtSet(cv_mem, cv_mem.cv_znQS[0], tempQS1);\r
+\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+\r
+       N_VAbs_Serial(cv_mem.cv_znQS[0][is], tempQ2);\r
+       N_VInv_Serial(tempQS1[is], tempQ1);\r
+       N_VLinearSum_Serial(HUB_FACTOR, tempQ2, ONE, tempQ1, tempQ1);\r
+       \r
+       N_VAbs_Serial(cv_mem.cv_znQS[1][is], tempQ2);\r
+       \r
+       N_VDiv_Serial(tempQ2, tempQ1, tempQ1);\r
+       hubQS_inv = N_VMaxNorm_Serial(tempQ1);\r
+\r
+       if (hubQS_inv > hub_inv) hub_inv = hubQS_inv;\r
+\r
+     }\r
+\r
+   }\r
+\r
+\r
+   /*\r
+    * bound based on tdist -- allow at most a step of magnitude\r
+    * HUB_FACTOR * tdist\r
+    */\r
+   \r
+   hub = HUB_FACTOR*tdist;\r
+\r
+   /* Use the smaler of the two */\r
+\r
+   if (hub*hub_inv > ONE) hub = ONE/hub_inv;\r
+\r
+   return(hub);\r
+ }\r
\r
+ /*\r
+  * cvYddNorm\r
+  *\r
+  * This routine computes an estimate of the second derivative of Y\r
+  * using a difference quotient, and returns its WRMS norm.\r
+  *\r
+  * Y contains all variables included in the error test. \r
+  */\r
+\r
+ private int cvYddNorm(CVodeMemRec cv_mem, double hg, double yddnrm[])\r
+ {\r
+   int retval, is;\r
+   NVector wrk1, wrk2;\r
+   \r
+   /* y <- h*y'(t) + y(t) */\r
+   \r
+   N_VLinearSum_Serial(hg, cv_mem.cv_zn[1], ONE, cv_mem.cv_zn[0], cv_mem.cv_y);\r
+   \r
+   if (cv_mem.cv_sensi && cv_mem.cv_errconS) \r
+     for (is=0; is<cv_mem.cv_Ns; is++)\r
+       N_VLinearSum_Serial(hg, cv_mem.cv_znS[1][is], ONE,\r
+                    cv_mem.cv_znS[0][is], cv_mem.cv_yS[is]);\r
+   \r
+   /* tempv <- f(t+h, h*y'(t)+y(t)) */\r
+\r
+   retval = f(cv_mem.cv_tn+hg, cv_mem.cv_y,\r
+                         cv_mem.cv_tempv, cv_mem.cv_user_data);\r
+   cv_mem.cv_nfe++;\r
+   if (retval < 0) return(CV_RHSFUNC_FAIL);\r
+   if (retval > 0) return(RHSFUNC_RECVR);\r
+\r
+   if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
+     retval = fQ(cv_mem.cv_tn+hg, cv_mem.cv_y,\r
+                            cv_mem.cv_tempvQ, cv_mem.cv_user_data);\r
+     cv_mem.cv_nfQe++;\r
+     if (retval < 0) return(CV_QRHSFUNC_FAIL);\r
+     if (retval > 0) return(QRHSFUNC_RECVR);\r
+   }\r
+\r
+   if (cv_mem.cv_sensi && cv_mem.cv_errconS) {\r
+     wrk1 = cv_mem.cv_ftemp;\r
+     wrk2 = cv_mem.cv_acor;\r
+     retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn+hg, cv_mem.cv_y,\r
+                               cv_mem.cv_tempv, cv_mem.cv_yS,\r
+                               cv_mem.cv_tempvS, wrk1, wrk2);\r
+     if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+     if (retval > 0) return(SRHSFUNC_RECVR);\r
+   }  \r
+\r
+   if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS) {\r
+     wrk1 = cv_mem.cv_ftemp;\r
+     wrk2 = cv_mem.cv_acorQ;\r
+     retval = fQS(cv_mem.cv_Ns, cv_mem.cv_tn+hg,\r
+                             cv_mem.cv_y, cv_mem.cv_yS,\r
+                             cv_mem.cv_tempvQ, cv_mem.cv_tempvQS,\r
+                             cv_mem.cv_fQS_data, wrk1, wrk2);\r
+\r
+     cv_mem.cv_nfQSe++;\r
+     if (retval < 0) return(CV_QSRHSFUNC_FAIL);\r
+     if (retval > 0) return(QSRHSFUNC_RECVR);\r
+   } \r
+\r
+   /* Load estimate of ||y''|| into tempv:\r
+    * tempv <-  (1/h) * f(t+h, h*y'(t)+y(t)) - y'(t) */\r
+   \r
+   N_VLinearSum_Serial(ONE, cv_mem.cv_tempv, -ONE, cv_mem.cv_zn[1], cv_mem.cv_tempv);\r
+   N_VScale_Serial(ONE/hg, cv_mem.cv_tempv, cv_mem.cv_tempv);\r
+   yddnrm[0] = N_VWrmsNorm_Serial(cv_mem.cv_tempv, cv_mem.cv_ewt);\r
+\r
+   if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
+     N_VLinearSum_Serial(ONE, cv_mem.cv_tempvQ, -ONE, \r
+                  cv_mem.cv_znQ[1], cv_mem.cv_tempvQ);\r
+     N_VScale_Serial(ONE/hg, cv_mem.cv_tempvQ, cv_mem.cv_tempvQ);\r
+     yddnrm[0] = cvQuadUpdateNorm(cv_mem, yddnrm[0], cv_mem.cv_tempvQ,\r
+                                cv_mem.cv_ewtQ);\r
+   }\r
+\r
+   if (cv_mem.cv_sensi && cv_mem.cv_errconS) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+       N_VLinearSum_Serial(ONE, cv_mem.cv_tempvS[is], -ONE,\r
+                    cv_mem.cv_znS[1][is], cv_mem.cv_tempvS[is]);\r
+       N_VScale_Serial(ONE/hg, cv_mem.cv_tempvS[is], cv_mem.cv_tempvS[is]);\r
+     }\r
+     yddnrm[0] = cvSensUpdateNorm(cv_mem, yddnrm[0], cv_mem.cv_tempvS,\r
+                                cv_mem.cv_ewtS);\r
+   }\r
+\r
+   if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+       N_VLinearSum_Serial(ONE, cv_mem.cv_tempvQS[is], -ONE,\r
+                    cv_mem.cv_znQS[1][is], cv_mem.cv_tempvQS[is]);\r
+       N_VScale_Serial(ONE/hg, cv_mem.cv_tempvQS[is], cv_mem.cv_tempvQS[is]);\r
+     }\r
+     yddnrm[0] = cvQuadSensUpdateNorm(cv_mem, yddnrm[0], cv_mem.cv_tempvQS,\r
+                                    cv_mem.cv_ewtQS);\r
+   }\r
+\r
+   return(CV_SUCCESS);\r
+ }\r
+\r
+ /*\r
+  * cvQuadUpdateNorm\r
+  *\r
+  * Updates the norm old_nrm to account for all quadratures.\r
+  */\r
+\r
+ private double cvQuadUpdateNorm(CVodeMemRec cv_mem, double old_nrm,\r
+                                  NVector xQ, NVector wQ)\r
+ {\r
+   double qnrm;\r
+\r
+   qnrm = N_VWrmsNorm_Serial(xQ, wQ);\r
+   if (old_nrm > qnrm) return(old_nrm);\r
+   else                return(qnrm);\r
+ }\r
+\r
+ /*\r
+  * cvSensUpdateNorm\r
+  *\r
+  * Updates the norm old_nrm to account for all sensitivities.\r
+  */\r
+\r
+ private double cvSensUpdateNorm(CVodeMemRec cv_mem, double old_nrm,\r
+                                  NVector xS[], NVector wS[])\r
+ {\r
+   double snrm;\r
+   \r
+   snrm = cvSensNorm(cv_mem, xS, wS);\r
+   if (old_nrm > snrm) return(old_nrm);\r
+   else                return(snrm);\r
+ }\r
\r
+ private double cvSensNorm(CVodeMemRec cv_mem, NVector xS[], NVector wS[])\r
+ {\r
+   int is;\r
+   double nrm, snrm;\r
+\r
+   nrm = N_VWrmsNorm_Serial(xS[0],wS[0]);\r
+   for (is=1; is<cv_mem.cv_Ns; is++) {\r
+     snrm = N_VWrmsNorm_Serial(xS[is],wS[is]);\r
+     if ( snrm > nrm ) nrm = snrm;\r
+   }\r
+\r
+   return(nrm);\r
+ }\r
+\r
+ private double cvQuadSensUpdateNorm(CVodeMemRec cv_mem, double old_nrm,\r
+         NVector xQS[], NVector wQS[])\r
+{\r
+double snrm;\r
+\r
+snrm = cvQuadSensNorm(cv_mem, xQS, wQS);\r
+if (old_nrm > snrm) return(old_nrm);\r
+else                return(snrm);\r
+}\r
\r
+ /*\r
+  * cvQuadSensNorm\r
+  *\r
+  * This routine returns the maximum over the weighted root mean \r
+  * square norm of xQS with weight vectors wQS:\r
+  *\r
+  *  max { wrms(xQS[0],wS[0]) ... wrms(xQS[Ns-1],wS[Ns-1]) }    \r
+  *\r
+  * Called by cvQuadSensUpdateNorm.\r
+  */\r
+\r
+ private double cvQuadSensNorm(CVodeMemRec cv_mem, NVector xQS[], NVector wQS[])\r
+ {\r
+   int is;\r
+   double nrm, snrm;\r
+\r
+   nrm = N_VWrmsNorm_Serial(xQS[0],wQS[0]);\r
+   for (is=1; is<cv_mem.cv_Ns; is++) {\r
+     snrm = N_VWrmsNorm_Serial(xQS[is],wQS[is]);\r
+     if ( snrm > nrm ) nrm = snrm;\r
+   }\r
+\r
+   return(nrm);\r
+ }\r
+\r
+ /* \r
+  * -----------------------------------------------------------------\r
+  * Function to handle failures\r
+  * -----------------------------------------------------------------\r
+  */\r
+\r
+ /*\r
+  * cvHandleFailure\r
+  *\r
+  * This routine prints error messages for all cases of failure by\r
+  * cvHin or cvStep. \r
+  * It returns to CVode the value that CVode is to return to the user.\r
+  */\r
+\r
+ private int cvHandleFailure(CVodeMemRec cv_mem, int flag)\r
+ {\r
+\r
+   /* Set vector of  absolute weighted local errors */\r
+   /*\r
+   N_VProd(acor, ewt, tempv);\r
+   N_VAbs(tempv, tempv);\r
+   */\r
+\r
+   /* Depending on flag, print error message and return error flag */\r
+   switch (flag) {\r
+   case CV_ERR_FAILURE:\r
+     cvProcessError(cv_mem, CV_ERR_FAILURE, "CVODES", "CVode",\r
+                    MSGCV_ERR_FAILS, cv_mem.cv_tn, cv_mem.cv_h);\r
+     break;\r
+   case CV_CONV_FAILURE:\r
+     cvProcessError(cv_mem, CV_CONV_FAILURE, "CVODES", "CVode",\r
+                    MSGCV_CONV_FAILS, cv_mem.cv_tn, cv_mem.cv_h);\r
+     break;\r
+   case CV_LSETUP_FAIL:\r
+     cvProcessError(cv_mem, CV_LSETUP_FAIL, "CVODES", "CVode",\r
+                    MSGCV_SETUP_FAILED, cv_mem.cv_tn);\r
+     break;\r
+   case CV_LSOLVE_FAIL:\r
+     cvProcessError(cv_mem, CV_LSOLVE_FAIL, "CVODES", "CVode",\r
+                    MSGCV_SOLVE_FAILED, cv_mem.cv_tn);\r
+     break;\r
+   case CV_RHSFUNC_FAIL:\r
+     cvProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODES", "CVode",\r
+                    MSGCV_RHSFUNC_FAILED, cv_mem.cv_tn);\r
+     break;\r
+   case CV_UNREC_RHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_UNREC_RHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_RHSFUNC_UNREC, cv_mem.cv_tn);\r
+     break;\r
+   case CV_REPTD_RHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_REPTD_RHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_RHSFUNC_REPTD, cv_mem.cv_tn);\r
+     break;\r
+   case CV_RTFUNC_FAIL:\r
+     cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "CVode",\r
+                    MSGCV_RTFUNC_FAILED, cv_mem.cv_tn);\r
+     break;\r
+   case CV_QRHSFUNC_FAIL:\r
+     cvProcessError(cv_mem, CV_QRHSFUNC_FAIL, "CVODES", "CVode",\r
+                    MSGCV_QRHSFUNC_FAILED, cv_mem.cv_tn);\r
+     break;\r
+   case CV_UNREC_QRHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_UNREC_QRHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_QRHSFUNC_UNREC, cv_mem.cv_tn);\r
+     break;\r
+   case CV_REPTD_QRHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_REPTD_QRHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_QRHSFUNC_REPTD, cv_mem.cv_tn);\r
+     break;\r
+   case CV_SRHSFUNC_FAIL:\r
+     cvProcessError(cv_mem, CV_SRHSFUNC_FAIL, "CVODES", "CVode",\r
+                    MSGCV_SRHSFUNC_FAILED, cv_mem.cv_tn);\r
+     break;\r
+   case CV_UNREC_SRHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_UNREC_SRHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_SRHSFUNC_UNREC, cv_mem.cv_tn);\r
+     break;\r
+   case CV_REPTD_SRHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_REPTD_SRHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_SRHSFUNC_REPTD, cv_mem.cv_tn);\r
+     break;\r
+   case CV_QSRHSFUNC_FAIL:\r
+     cvProcessError(cv_mem, CV_QSRHSFUNC_FAIL, "CVODES", "CVode",\r
+                    MSGCV_QSRHSFUNC_FAILED, cv_mem.cv_tn);\r
+     break;\r
+   case CV_UNREC_QSRHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_UNREC_QSRHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_QSRHSFUNC_UNREC, cv_mem.cv_tn);\r
+     break;\r
+   case CV_REPTD_QSRHSFUNC_ERR:\r
+     cvProcessError(cv_mem, CV_REPTD_QSRHSFUNC_ERR, "CVODES", "CVode",\r
+                    MSGCV_QSRHSFUNC_REPTD, cv_mem.cv_tn);\r
+     break;\r
+   case CV_TOO_CLOSE:\r
+     cvProcessError(cv_mem, CV_TOO_CLOSE, "CVODES", "CVode",\r
+                    MSGCV_TOO_CLOSE);\r
+     break;\r
+   default:\r
+     return(CV_SUCCESS);\r
+   }\r
+\r
+   return(flag);\r
+ }\r
\r
+ /* \r
+  * cvRcheck1\r
+  *\r
+  * This routine completes the initialization of rootfinding memory\r
+  * information, and checks whether g has a zero both at and very near\r
+  * the initial point of the IVP.\r
+  *\r
+  * This routine returns an int equal to:\r
+  *  CV_RTFUNC_FAIL < 0 if the g function failed, or\r
+  *  CV_SUCCESS     = 0 otherwise.\r
+  */\r
+\r
+ private int cvRcheck1(CVodeMemRec cv_mem)\r
+ {\r
+   int i, retval;\r
+   double smallh, hratio, tplus;\r
+   boolean zroot;\r
+\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++)\r
+     cv_mem.cv_iroots[i] = 0;\r
+   cv_mem.cv_tlo = cv_mem.cv_tn;\r
+   cv_mem.cv_ttol = (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h)) *\r
+     cv_mem.cv_uround*HUNDRED;\r
+\r
+   /* Evaluate g at initial t and check for zero values. */\r
+   retval = g(cv_mem.cv_tlo, cv_mem.cv_zn[0],\r
+                            cv_mem.cv_glo, cv_mem.cv_user_data);\r
+   cv_mem.cv_nge = 1;\r
+   if (retval != 0) return(CV_RTFUNC_FAIL);\r
+\r
+   zroot = false;\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++) {\r
+     if (Math.abs(cv_mem.cv_glo[i]) == ZERO) {\r
+       zroot = true;\r
+       cv_mem.cv_gactive[i] = false;\r
+     }\r
+   }\r
+   if (!zroot) return(CV_SUCCESS);\r
+\r
+   /* Some g_i is zero at t0; look at g at t0+(small increment). */\r
+   hratio = Math.max(cv_mem.cv_ttol/Math.abs(cv_mem.cv_h), PT1);\r
+   smallh = hratio*cv_mem.cv_h;\r
+   tplus = cv_mem.cv_tlo + smallh;\r
+   N_VLinearSum_Serial(ONE, cv_mem.cv_zn[0], hratio, cv_mem.cv_zn[1], cv_mem.cv_y);\r
+   retval = g(tplus, cv_mem.cv_y,\r
+                            cv_mem.cv_ghi, cv_mem.cv_user_data);\r
+   cv_mem.cv_nge++;\r
+   if (retval != 0) return(CV_RTFUNC_FAIL);\r
+\r
+   /* We check now only the components of g which were exactly 0.0 at t0\r
+    * to see if we can 'activate' them. */\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++) {\r
+     if (!cv_mem.cv_gactive[i] && Math.abs(cv_mem.cv_ghi[i]) != ZERO) {\r
+       cv_mem.cv_gactive[i] = true;\r
+       cv_mem.cv_glo[i] = cv_mem.cv_ghi[i];\r
+     }\r
+   }\r
+   return(CV_SUCCESS);\r
+ }\r
+\r
+ /*\r
+  * cvRcheck2\r
+  *\r
+  * This routine checks for exact zeros of g at the last root found,\r
+  * if the last return was a root.  It then checks for a close pair of\r
+  * zeros (an error condition), and for a new root at a nearby point.\r
+  * The array glo = g(tlo) at the left endpoint of the search interval\r
+  * is adjusted if necessary to assure that all g_i are nonzero\r
+  * there, before returning to do a root search in the interval.\r
+  *\r
+  * On entry, tlo = tretlast is the last value of tret returned by\r
+  * CVode.  This may be the previous tn, the previous tout value,\r
+  * or the last root location.\r
+  *\r
+  * This routine returns an int equal to:\r
+  *     CV_RTFUNC_FAIL  < 0 if the g function failed, or\r
+  *     CLOSERT         = 3 if a close pair of zeros was found, or\r
+  *     RTFOUND         = 1 if a new zero of g was found near tlo, or\r
+  *     CV_SUCCESS      = 0 otherwise.\r
+  */\r
+\r
+ private int cvRcheck2(CVodeMemRec cv_mem)\r
+ {\r
+   int i, retval;\r
+   double smallh, hratio, tplus;\r
+   boolean zroot;\r
+\r
+   if (cv_mem.cv_irfnd == 0) return(CV_SUCCESS);\r
+\r
+   CVodeGetDky(cv_mem, cv_mem.cv_tlo, 0, cv_mem.cv_y);\r
+   retval = g(cv_mem.cv_tlo, cv_mem.cv_y,\r
+                            cv_mem.cv_glo, cv_mem.cv_user_data);\r
+   cv_mem.cv_nge++;\r
+   if (retval != 0) return(CV_RTFUNC_FAIL);\r
+\r
+   zroot = false;\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++)\r
+     cv_mem.cv_iroots[i] = 0;\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++) {\r
+     if (!cv_mem.cv_gactive[i]) continue;\r
+     if (Math.abs(cv_mem.cv_glo[i]) == ZERO) {\r
+       zroot = true;\r
+       cv_mem.cv_iroots[i] = 1;\r
+     }\r
+   }\r
+   if (!zroot) return(CV_SUCCESS);\r
+\r
+   /* One or more g_i has a zero at tlo.  Check g at tlo+smallh. */\r
+   cv_mem.cv_ttol = (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h)) *\r
+     cv_mem.cv_uround*HUNDRED;\r
+   smallh = (cv_mem.cv_h > ZERO) ? cv_mem.cv_ttol : -cv_mem.cv_ttol;\r
+   tplus = cv_mem.cv_tlo + smallh;\r
+   if ( (tplus - cv_mem.cv_tn)*cv_mem.cv_h >= ZERO) {\r
+     hratio = smallh/cv_mem.cv_h;\r
+     N_VLinearSum_Serial(ONE, cv_mem.cv_y, hratio, cv_mem.cv_zn[1], cv_mem.cv_y);\r
+   } else {\r
+     CVodeGetDky(cv_mem, tplus, 0, cv_mem.cv_y);\r
+   }\r
+   retval = g(tplus, cv_mem.cv_y,\r
+                            cv_mem.cv_ghi, cv_mem.cv_user_data);\r
+   cv_mem.cv_nge++;\r
+   if (retval != 0) return(CV_RTFUNC_FAIL);\r
+\r
+   /* Check for close roots (error return), for a new zero at tlo+smallh,\r
+   and for a g_i that changed from zero to nonzero. */\r
+   zroot = false;\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++) {\r
+     if (!cv_mem.cv_gactive[i]) continue;\r
+     if (Math.abs(cv_mem.cv_ghi[i]) == ZERO) {\r
+       if (cv_mem.cv_iroots[i] == 1) return(CLOSERT);\r
+       zroot = true;\r
+       cv_mem.cv_iroots[i] = 1;\r
+     } else {\r
+       if (cv_mem.cv_iroots[i] == 1)\r
+         cv_mem.cv_glo[i] = cv_mem.cv_ghi[i];\r
+     }\r
+   }\r
+   if (zroot) return(RTFOUND);\r
+   return(CV_SUCCESS);\r
+ }\r
\r
+ /*\r
+  * CVodeGetDky\r
+  *\r
+  * This routine computes the k-th derivative of the interpolating\r
+  * polynomial at the time t and stores the result in the vector dky.\r
+  * The formula is:\r
+  *         q \r
+  *  dky = SUM c(j,k) * (t - tn)^(j-k) * h^(-j) * zn[j] , \r
+  *        j=k \r
+  * where c(j,k) = j*(j-1)*...*(j-k+1), q is the current order, and\r
+  * zn[j] is the j-th column of the Nordsieck history array.\r
+  *\r
+  * This function is called by CVode with k = 0 and t = tout, but\r
+  * may also be called directly by the user.\r
+  */\r
+\r
+ private int CVodeGetDky(CVodeMemRec cv_mem, double t, int k, NVector dky)\r
+ {\r
+   double s, c, r;\r
+   double tfuzz, tp, tn1;\r
+   int i, j;\r
+   \r
+   /* Check all inputs for legality */\r
+  \r
+   if (cv_mem == null) {\r
+     cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetDky", MSGCV_NO_MEM);\r
+     return(CV_MEM_NULL);\r
+   }\r
+\r
+   if (dky == null) {\r
+     cvProcessError(cv_mem, CV_BAD_DKY, "CVODES", "CVodeGetDky", MSGCV_NULL_DKY);\r
+     return(CV_BAD_DKY);\r
+   }\r
+\r
+   if ((k < 0) || (k > cv_mem.cv_q)) {\r
+     cvProcessError(cv_mem, CV_BAD_K, "CVODES", "CVodeGetDky", MSGCV_BAD_K);\r
+     return(CV_BAD_K);\r
+   }\r
+   \r
+   /* Allow for some slack */\r
+   tfuzz = FUZZ_FACTOR * cv_mem.cv_uround *\r
+     (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_hu));\r
+   if (cv_mem.cv_hu < ZERO) tfuzz = -tfuzz;\r
+   tp = cv_mem.cv_tn - cv_mem.cv_hu - tfuzz;\r
+   tn1 = cv_mem.cv_tn + tfuzz;\r
+   if ((t-tp)*(t-tn1) > ZERO) {\r
+     cvProcessError(cv_mem, CV_BAD_T, "CVODES", "CVodeGetDky", MSGCV_BAD_T,\r
+                    t, cv_mem.cv_tn-cv_mem.cv_hu, cv_mem.cv_tn);\r
+     return(CV_BAD_T);\r
+   }\r
+\r
+   /* Sum the differentiated interpolating polynomial */\r
+\r
+   s = (t - cv_mem.cv_tn) / cv_mem.cv_h;\r
+   for (j=cv_mem.cv_q; j >= k; j--) {\r
+     c = ONE;\r
+     for (i=j; i >= j-k+1; i--) c *= i;\r
+     if (j == cv_mem.cv_q) {\r
+       N_VScale_Serial(c, cv_mem.cv_zn[cv_mem.cv_q], dky);\r
+     } else {\r
+       N_VLinearSum_Serial(c, cv_mem.cv_zn[j], s, dky, dky);\r
+     }\r
+   }\r
+   if (k == 0) return(CV_SUCCESS);\r
+   r = Math.pow(cv_mem.cv_h, -k);\r
+   N_VScale_Serial(r, dky, dky);\r
+   return(CV_SUCCESS);\r
+ }\r
+\r
+ /*\r
+  * cvRcheck3\r
+  *\r
+  * This routine interfaces to cvRootfind to look for a root of g\r
+  * between tlo and either tn or tout, whichever comes first.\r
+  * Only roots beyond tlo in the direction of integration are sought.\r
+  *\r
+  * This routine returns an int equal to:\r
+  *     CV_RTFUNC_FAIL  < 0 if the g function failed, or\r
+  *     RTFOUND         = 1 if a root of g was found, or\r
+  *     CV_SUCCESS      = 0 otherwise.\r
+  */\r
+\r
+ private int cvRcheck3(CVodeMemRec cv_mem)\r
+ {\r
+   int i, ier, retval;\r
+\r
+   /* Set thi = tn or tout, whichever comes first; set y = y(thi). */\r
+   if (cv_mem.cv_taskc == CV_ONE_STEP) {\r
+     cv_mem.cv_thi = cv_mem.cv_tn;\r
+     N_VScale_Serial(ONE, cv_mem.cv_zn[0], cv_mem.cv_y);\r
+   }\r
+   if (cv_mem.cv_taskc == CV_NORMAL) {\r
+     if ( (cv_mem.cv_toutc - cv_mem.cv_tn)*cv_mem.cv_h >= ZERO) {\r
+       cv_mem.cv_thi = cv_mem.cv_tn; \r
+       N_VScale_Serial(ONE, cv_mem.cv_zn[0], cv_mem.cv_y);\r
+     } else {\r
+       cv_mem.cv_thi = cv_mem.cv_toutc;\r
+       CVodeGetDky(cv_mem, cv_mem.cv_thi, 0, cv_mem.cv_y);\r
+     }\r
+   }\r
+\r
+   /* Set ghi = g(thi) and call cvRootfind to search (tlo,thi) for roots. */\r
+   retval = g(cv_mem.cv_thi, cv_mem.cv_y,\r
+                            cv_mem.cv_ghi, cv_mem.cv_user_data);\r
+   cv_mem.cv_nge++;\r
+   if (retval != 0) return(CV_RTFUNC_FAIL);\r
+\r
+   cv_mem.cv_ttol = (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h)) *\r
+     cv_mem.cv_uround*HUNDRED;\r
+   ier = cvRootfind(cv_mem);\r
+   if (ier == CV_RTFUNC_FAIL) return(CV_RTFUNC_FAIL);\r
+   for(i=0; i<cv_mem.cv_nrtfn; i++) {\r
+     if(!cv_mem.cv_gactive[i] && cv_mem.cv_grout[i] != ZERO)\r
+       cv_mem.cv_gactive[i] = true;\r
+   }\r
+   cv_mem.cv_tlo = cv_mem.cv_trout;\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++)\r
+     cv_mem.cv_glo[i] = cv_mem.cv_grout[i];\r
+\r
+   /* If no root found, return CV_SUCCESS. */  \r
+   if (ier == CV_SUCCESS) return(CV_SUCCESS);\r
+\r
+   /* If a root was found, interpolate to get y(trout) and return.  */\r
+   CVodeGetDky(cv_mem, cv_mem.cv_trout, 0, cv_mem.cv_y);\r
+   return(RTFOUND);\r
+ }\r
+\r
+ /*\r
+  * cvRootfind\r
+  *\r
+  * This routine solves for a root of g(t) between tlo and thi, if\r
+  * one exists.  Only roots of odd multiplicity (i.e. with a change\r
+  * of sign in one of the g_i), or exact zeros, are found.\r
+  * Here the sign of tlo - thi is arbitrary, but if multiple roots\r
+  * are found, the one closest to tlo is returned.\r
+  *\r
+  * The method used is the Illinois algorithm, a modified secant method.\r
+  * Reference: Kathie L. Hiebert and Lawrence F. Shampine, Implicitly\r
+  * Defined Output Points for Solutions of ODEs, Sandia National\r
+  * Laboratory Report SAND80-0180, February 1980.\r
+  *\r
+  * This routine uses the following parameters for communication:\r
+  *\r
+  * nrtfn    = number of functions g_i, or number of components of\r
+  *            the vector-valued function g(t).  Input only.\r
+  *\r
+  * gfun     = user-defined function for g(t).  Its form is\r
+  *            (void) gfun(t, y, gt, user_data)\r
+  *\r
+  * rootdir  = in array specifying the direction of zero-crossings.\r
+  *            If rootdir[i] > 0, search for roots of g_i only if\r
+  *            g_i is increasing; if rootdir[i] < 0, search for\r
+  *            roots of g_i only if g_i is decreasing; otherwise\r
+  *            always search for roots of g_i.\r
+  *\r
+  * gactive  = array specifying whether a component of g should\r
+  *            or should not be monitored. gactive[i] is initially\r
+  *            set to SUNTRUE for all i=0,...,nrtfn-1, but it may be\r
+  *            reset to SUNFALSE if at the first step g[i] is 0.0\r
+  *            both at the I.C. and at a small perturbation of them.\r
+  *            gactive[i] is then set back on SUNTRUE only after the \r
+  *            corresponding g function moves away from 0.0.\r
+  *\r
+  * nge      = cumulative counter for gfun calls.\r
+  *\r
+  * ttol     = a convergence tolerance for trout.  Input only.\r
+  *            When a root at trout is found, it is located only to\r
+  *            within a tolerance of ttol.  Typically, ttol should\r
+  *            be set to a value on the order of\r
+  *               100 * UROUND * max (SUNRabs(tlo), SUNRabs(thi))\r
+  *            where UROUND is the unit roundoff of the machine.\r
+  *\r
+  * tlo, thi = endpoints of the interval in which roots are sought.\r
+  *            On input, these must be distinct, but tlo - thi may\r
+  *            be of either sign.  The direction of integration is\r
+  *            assumed to be from tlo to thi.  On return, tlo and thi\r
+  *            are the endpoints of the final relevant interval.\r
+  *\r
+  * glo, ghi = arrays of length nrtfn containing the vectors g(tlo)\r
+  *            and g(thi) respectively.  Input and output.  On input,\r
+  *            none of the glo[i] should be zero.\r
+  *\r
+  * trout    = root location, if a root was found, or thi if not.\r
+  *            Output only.  If a root was found other than an exact\r
+  *            zero of g, trout is the endpoint thi of the final\r
+  *            interval bracketing the root, with size at most ttol.\r
+  *\r
+  * grout    = array of length nrtfn containing g(trout) on return.\r
+  *\r
+  * iroots   = int array of length nrtfn with root information.\r
+  *            Output only.  If a root was found, iroots indicates\r
+  *            which components g_i have a root at trout.  For\r
+  *            i = 0, ..., nrtfn-1, iroots[i] = 1 if g_i has a root\r
+  *            and g_i is increasing, iroots[i] = -1 if g_i has a\r
+  *            root and g_i is decreasing, and iroots[i] = 0 if g_i\r
+  *            has no roots or g_i varies in the direction opposite\r
+  *            to that indicated by rootdir[i].\r
+  *\r
+  * This routine returns an int equal to:\r
+  *      CV_RTFUNC_FAIL  < 0 if the g function failed, or\r
+  *      RTFOUND         = 1 if a root of g was found, or\r
+  *      CV_SUCCESS      = 0 otherwise.\r
+  */\r
+\r
+ private int cvRootfind(CVodeMemRec cv_mem)\r
+ {\r
+   double alph, tmid, gfrac, maxfrac, fracint, fracsub;\r
+   int i, retval, imax, side, sideprev;\r
+   boolean zroot, sgnchg;\r
+\r
+   imax = 0;\r
+\r
+   /* First check for change in sign in ghi or for a zero in ghi. */\r
+   maxfrac = ZERO;\r
+   zroot = false;\r
+   sgnchg = false;\r
+   for (i = 0;  i < cv_mem.cv_nrtfn; i++) {\r
+     if(!cv_mem.cv_gactive[i]) continue;\r
+     if (Math.abs(cv_mem.cv_ghi[i]) == ZERO) {\r
+       if(cv_mem.cv_rootdir[i]*cv_mem.cv_glo[i] <= ZERO) {\r
+         zroot = true;\r
+       }\r
+     } else {\r
+       if ( (cv_mem.cv_glo[i]*cv_mem.cv_ghi[i] < ZERO) &&\r
+            (cv_mem.cv_rootdir[i]*cv_mem.cv_glo[i] <= ZERO) ) {\r
+         gfrac = Math.abs(cv_mem.cv_ghi[i]/(cv_mem.cv_ghi[i] - cv_mem.cv_glo[i]));\r
+         if (gfrac > maxfrac) {\r
+           sgnchg = true;\r
+           maxfrac = gfrac;\r
+           imax = i;\r
+         }\r
+       }\r
+     }\r
+   }\r
+\r
+   /* If no sign change was found, reset trout and grout.  Then return\r
+      CV_SUCCESS if no zero was found, or set iroots and return RTFOUND.  */ \r
+   if (!sgnchg) {\r
+     cv_mem.cv_trout = cv_mem.cv_thi;\r
+     for (i = 0; i < cv_mem.cv_nrtfn; i++)\r
+       cv_mem.cv_grout[i] = cv_mem.cv_ghi[i];\r
+     if (!zroot) return(CV_SUCCESS);\r
+     for (i = 0; i < cv_mem.cv_nrtfn; i++) {\r
+       cv_mem.cv_iroots[i] = 0;\r
+       if(!cv_mem.cv_gactive[i]) continue;\r
+       if ( (Math.abs(cv_mem.cv_ghi[i]) == ZERO) &&\r
+            (cv_mem.cv_rootdir[i]*cv_mem.cv_glo[i] <= ZERO) )\r
+         cv_mem.cv_iroots[i] = cv_mem.cv_glo[i] > 0 ? -1:1;\r
+     }\r
+     return(RTFOUND);\r
+   }\r
+\r
+   /* Initialize alph to avoid compiler warning */\r
+   alph = ONE;\r
+\r
+   /* A sign change was found.  Loop to locate nearest root. */\r
+\r
+   side = 0;  sideprev = -1;\r
+   for(;;) {                                    /* Looping point */\r
+\r
+     /* If interval size is already less than tolerance ttol, break. */\r
+       if (Math.abs(cv_mem.cv_thi - cv_mem.cv_tlo) <= cv_mem.cv_ttol) break;\r
+\r
+     /* Set weight alph.\r
+        On the first two passes, set alph = 1.  Thereafter, reset alph\r
+        according to the side (low vs high) of the subinterval in which\r
+        the sign change was found in the previous two passes.\r
+        If the sides were opposite, set alph = 1.\r
+        If the sides were the same, then double alph (if high side),\r
+        or halve alph (if low side).\r
+        The next guess tmid is the secant method value if alph = 1, but\r
+        is closer to cv_mem->cv_tlo if alph < 1, and closer to thi if alph > 1.    */\r
+\r
+     if (sideprev == side) {\r
+       alph = (side == 2) ? alph*TWO : alph*HALF;\r
+     } else {\r
+       alph = ONE;\r
+     }\r
+\r
+     /* Set next root approximation tmid and get g(tmid).\r
+        If tmid is too close to tlo or thi, adjust it inward,\r
+        by a fractional distance that is between 0.1 and 0.5.  */\r
+     tmid = cv_mem.cv_thi - (cv_mem.cv_thi - cv_mem.cv_tlo) *\r
+       cv_mem.cv_ghi[imax] / (cv_mem.cv_ghi[imax] - alph*cv_mem.cv_glo[imax]);\r
+     if (Math.abs(tmid - cv_mem.cv_tlo) < HALF*cv_mem.cv_ttol) {\r
+       fracint = Math.abs(cv_mem.cv_thi - cv_mem.cv_tlo)/cv_mem.cv_ttol;\r
+       fracsub = (fracint > FIVE) ? PT1 : HALF/fracint;\r
+       tmid = cv_mem.cv_tlo + fracsub*(cv_mem.cv_thi - cv_mem.cv_tlo);\r
+     }\r
+     if (Math.abs(cv_mem.cv_thi - tmid) < HALF*cv_mem.cv_ttol) {\r
+       fracint = Math.abs(cv_mem.cv_thi - cv_mem.cv_tlo)/cv_mem.cv_ttol;\r
+       fracsub = (fracint > FIVE) ? PT1 : HALF/fracint;\r
+       tmid = cv_mem.cv_thi - fracsub*(cv_mem.cv_thi - cv_mem.cv_tlo);\r
+     }\r
+\r
+     CVodeGetDky(cv_mem, tmid, 0, cv_mem.cv_y);\r
+     retval = g(tmid, cv_mem.cv_y, cv_mem.cv_grout,\r
+                              cv_mem.cv_user_data);\r
+     cv_mem.cv_nge++;\r
+     if (retval != 0) return(CV_RTFUNC_FAIL);\r
+\r
+     /* Check to see in which subinterval g changes sign, and reset imax.\r
+        Set side = 1 if sign change is on low side, or 2 if on high side.  */  \r
+     maxfrac = ZERO;\r
+     zroot = false;\r
+     sgnchg = false;\r
+     sideprev = side;\r
+     for (i = 0;  i < cv_mem.cv_nrtfn; i++) {\r
+       if(!cv_mem.cv_gactive[i]) continue;\r
+       if (Math.abs(cv_mem.cv_grout[i]) == ZERO) {\r
+         if(cv_mem.cv_rootdir[i]*cv_mem.cv_glo[i] <= ZERO) zroot = true;\r
+       } else {\r
+         if ( (cv_mem.cv_glo[i]*cv_mem.cv_grout[i] < ZERO) &&\r
+              (cv_mem.cv_rootdir[i]*cv_mem.cv_glo[i] <= ZERO) ) {\r
+           gfrac = Math.abs(cv_mem.cv_grout[i] /\r
+                           (cv_mem.cv_grout[i] - cv_mem.cv_glo[i]));\r
+           if (gfrac > maxfrac) {\r
+             sgnchg = true;\r
+             maxfrac = gfrac;\r
+             imax = i;\r
+           }\r
+         }\r
+       }\r
+     }\r
+     if (sgnchg) {\r
+       /* Sign change found in (tlo,tmid); replace thi with tmid. */\r
+       cv_mem.cv_thi = tmid;\r
+       for (i = 0; i < cv_mem.cv_nrtfn; i++)\r
+         cv_mem.cv_ghi[i] = cv_mem.cv_grout[i];\r
+       side = 1;\r
+       /* Stop at root thi if converged; otherwise loop. */\r
+       if (Math.abs(cv_mem.cv_thi - cv_mem.cv_tlo) <= cv_mem.cv_ttol) break;\r
+       continue;  /* Return to looping point. */\r
+     }\r
+\r
+     if (zroot) {\r
+       /* No sign change in (tlo,tmid), but g = 0 at tmid; return root tmid. */\r
+       cv_mem.cv_thi = tmid;\r
+       for (i = 0; i < cv_mem.cv_nrtfn; i++)\r
+         cv_mem.cv_ghi[i] = cv_mem.cv_grout[i];\r
+       break;\r
+     }\r
+\r
+     /* No sign change in (tlo,tmid), and no zero at tmid.\r
+        Sign change must be in (tmid,thi).  Replace tlo with tmid. */\r
+     cv_mem.cv_tlo = tmid;\r
+     for (i = 0; i < cv_mem.cv_nrtfn; i++)\r
+       cv_mem.cv_glo[i] = cv_mem.cv_grout[i];\r
+     side = 2;\r
+     /* Stop at root thi if converged; otherwise loop back. */\r
+     if (Math.abs(cv_mem.cv_thi - cv_mem.cv_tlo) <= cv_mem.cv_ttol) break;\r
+\r
+   } /* End of root-search loop */\r
+\r
+   /* Reset trout and grout, set iroots, and return RTFOUND. */\r
+   cv_mem.cv_trout = cv_mem.cv_thi;\r
+   for (i = 0; i < cv_mem.cv_nrtfn; i++) {\r
+     cv_mem.cv_grout[i] = cv_mem.cv_ghi[i];\r
+     cv_mem.cv_iroots[i] = 0;\r
+     if(!cv_mem.cv_gactive[i]) continue;\r
+     if ( (Math.abs(cv_mem.cv_ghi[i]) == ZERO) &&\r
+          (cv_mem.cv_rootdir[i]*cv_mem.cv_glo[i] <= ZERO) )\r
+       cv_mem.cv_iroots[i] = cv_mem.cv_glo[i] > 0 ? -1:1;\r
+     if ( (cv_mem.cv_glo[i]*cv_mem.cv_ghi[i] < ZERO) &&\r
+          (cv_mem.cv_rootdir[i]*cv_mem.cv_glo[i] <= ZERO) ) \r
+       cv_mem.cv_iroots[i] = cv_mem.cv_glo[i] > 0 ? -1:1;\r
+   }\r
+   return(RTFOUND);\r
+ }\r
+\r
+ /* \r
+  * -----------------------------------------------------------------\r
+  * Main cvStep function\r
+  * -----------------------------------------------------------------\r
+  */\r
+\r
+ /* \r
+  * cvStep\r
+  *\r
+  * This routine performs one internal cvode step, from tn to tn + h.\r
+  * It calls other routines to do all the work.\r
+  *\r
+  * The main operations done here are as follows:\r
+  * - preliminary adjustments if a new step size was chosen;\r
+  * - prediction of the Nordsieck history array zn at tn + h;\r
+  * - setting of multistep method coefficients and test quantities;\r
+  * - solution of the nonlinear system;\r
+  * - testing the local error;\r
+  * - updating zn and other state data if successful;\r
+  * - resetting stepsize and order for the next step.\r
+  * - if SLDET is on, check for stability, reduce order if necessary.\r
+  * On a failure in the nonlinear system solution or error test, the\r
+  * step may be reattempted, depending on the nature of the failure.\r
+  */\r
+\r
+ private int cvStep(CVodeMemRec cv_mem)\r
+ {\r
+   double saved_t, dsm, dsmQ, dsmS, dsmQS;\r
+   boolean do_sensi_stg, do_sensi_stg1;\r
+   int ncf, ncfS;\r
+   int nef, nefQ, nefS, nefQS;\r
+   int nflag, kflag, eflag;\r
+   int retval, is;\r
+\r
+   /* Are we computing sensitivities with a staggered approach? */\r
+\r
+   do_sensi_stg  = (cv_mem.cv_sensi && (cv_mem.cv_ism==CV_STAGGERED));\r
+   do_sensi_stg1 = (cv_mem.cv_sensi && (cv_mem.cv_ism==CV_STAGGERED1));\r
+\r
+   /* Initialize local counters for convergence and error test failures */\r
+\r
+   ncf  = nef  = 0;\r
+   nefQ = nefQS = 0;\r
+   ncfS = nefS = 0;\r
+   if (do_sensi_stg1) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++)\r
+       cv_mem.cv_ncfS1[is] = 0;\r
+   }\r
+\r
+   /* If needed, adjust method parameters */\r
+\r
+   if ((cv_mem.cv_nst > 0) && (cv_mem.cv_hprime != cv_mem.cv_h))\r
+     cvAdjustParams(cv_mem);\r
+\r
+   /* Looping point for attempts to take a step */\r
+\r
+   saved_t = cv_mem.cv_tn;\r
+   nflag = FIRST_CALL;\r
+\r
+   for(;;) {  \r
+\r
+     cvPredict(cv_mem);  \r
+     cvSet(cv_mem);\r
+\r
+     /* ------ Correct state variables ------ */\r
+     \r
+     nflag = cvNls(cv_mem, nflag);\r
+   /*  kflag = cvHandleNFlag(cv_mem, &nflag, saved_t, &ncf, &(cv_mem->cv_ncfn));*/\r
+\r
+     /* Go back in loop if we need to predict again (nflag=PREV_CONV_FAIL) */\r
+    /* if (kflag == PREDICT_AGAIN) continue;*/\r
+\r
+     /* Return if nonlinear solve failed and recovery not possible. */\r
+    /* if (kflag != DO_ERROR_TEST) return(kflag);*/\r
+\r
+     /* Perform error test (nflag=CV_SUCCESS) */\r
+   /*  eflag = cvDoErrorTest(cv_mem, &nflag, saved_t, cv_mem->cv_acnrm,\r
+                           &nef, &(cv_mem->cv_netf), &dsm);*/\r
+\r
+     /* Go back in loop if we need to predict again (nflag=PREV_ERR_FAIL) */\r
+    /* if (eflag == TRY_AGAIN) continue;*/\r
+\r
+     /* Return if error test failed and recovery not possible. */\r
+    /* if (eflag != CV_SUCCESS) return(eflag); */\r
+\r
+     /* Error test passed (eflag=CV_SUCCESS, nflag=CV_SUCCESS), go on */\r
+\r
+     /* ------ Correct the quadrature variables ------ */\r
+\r
+    /* if (cv_mem->cv_quadr) {\r
+\r
+       ncf = nef = 0;*/ /* reset counters for states */\r
+\r
+     /*  nflag = cvQuadNls(cv_mem);\r
+       kflag = cvHandleNFlag(cv_mem, &nflag, saved_t, &ncf, &(cv_mem->cv_ncfn));\r
+\r
+       if (kflag == PREDICT_AGAIN) continue;\r
+       if (kflag != DO_ERROR_TEST) return(kflag);*/\r
+\r
+       /* Error test on quadratures */\r
+      /* if (cv_mem->cv_errconQ) {\r
+         cv_mem->cv_acnrmQ = N_VWrmsNorm(cv_mem->cv_acorQ, cv_mem->cv_ewtQ);\r
+         eflag = cvDoErrorTest(cv_mem, &nflag, saved_t, cv_mem->cv_acnrmQ,\r
+                               &nefQ, &(cv_mem->cv_netfQ), &dsmQ);\r
+\r
+         if (eflag == TRY_AGAIN) continue;\r
+         if (eflag != CV_SUCCESS) return(eflag);*/\r
+\r
+         /* Set dsm = max(dsm, dsmQ) to be used in cvPrepareNextStep */\r
+        /* if (dsmQ > dsm) dsm = dsmQ;\r
+       }\r
+\r
+     }*/\r
+\r
+     /* ------ Correct the sensitivity variables (STAGGERED or STAGGERED1) ------- */\r
+\r
+    /* if (do_sensi_stg || do_sensi_stg1) {\r
+\r
+       ncf = nef = 0;    */    /* reset counters for states     */\r
+      /* if (cv_mem->cv_quadr) nefQ = 0;  *//* reset counter for quadratures */\r
+\r
+       /* Evaluate f at converged y, needed for future evaluations of sens. RHS \r
+        * If f() fails recoverably, treat it as a convergence failure and \r
+        * attempt the step again */\r
+\r
+      /* retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_y,\r
+                             cv_mem->cv_ftemp, cv_mem->cv_user_data);\r
+       cv_mem->cv_nfe++;\r
+       if (retval < 0) return(CV_RHSFUNC_FAIL);\r
+       if (retval > 0) {\r
+         nflag = PREV_CONV_FAIL;\r
+         continue;\r
+       }\r
+\r
+       if (do_sensi_stg) {*/\r
+         /* Nonlinear solve for sensitivities (all-at-once) */\r
+        /* nflag = cvStgrNls(cv_mem);\r
+         kflag = cvHandleNFlag(cv_mem, &nflag, saved_t, &ncfS,\r
+                               &(cv_mem->cv_ncfnS));\r
+       } else {*/\r
+         /* Nonlinear solve for sensitivities (one-by-one) */\r
+        /* for (is=0; is<cv_mem->cv_Ns; is++) { \r
+           nflag = cvStgr1Nls(cv_mem, is); \r
+           kflag = cvHandleNFlag(cv_mem, &nflag, saved_t,\r
+                                 &(cv_mem->cv_ncfS1[is]),\r
+                                 &(cv_mem->cv_ncfnS1[is]));\r
+           if (kflag != DO_ERROR_TEST) break; \r
+         }\r
+       }\r
+\r
+       if (kflag == PREDICT_AGAIN) continue;\r
+       if (kflag != DO_ERROR_TEST) return(kflag);*/\r
+\r
+       /* Error test on sensitivities */\r
+     /*  if (cv_mem->cv_errconS) {\r
+\r
+         if (do_sensi_stg1)\r
+           cv_mem->cv_acnrmS = cvSensNorm(cv_mem, cv_mem->cv_acorS, cv_mem->cv_ewtS);\r
+\r
+         eflag = cvDoErrorTest(cv_mem, &nflag, saved_t, cv_mem->cv_acnrmS,\r
+                               &nefS, &(cv_mem->cv_netfS), &dsmS);\r
+\r
+         if (eflag == TRY_AGAIN)  continue;\r
+         if (eflag != CV_SUCCESS) return(eflag);*/\r
+\r
+         /* Set dsm = max(dsm, dsmS) to be used in cvPrepareNextStep */\r
+      /*   if (dsmS > dsm) dsm = dsmS;\r
+\r
+       }\r
+\r
+     } */\r
+\r
+     /* ------ Correct the quadrature sensitivity variables ------ */\r
+\r
+   /*  if (cv_mem->cv_quadr_sensi) { */\r
+\r
+       /* Reset local convergence and error test failure counters */\r
+    /*   ncf = nef = 0;\r
+       if (cv_mem->cv_quadr) nefQ = 0;\r
+       if (do_sensi_stg) ncfS = nefS = 0;\r
+       if (do_sensi_stg1) {\r
+         for (is=0; is<cv_mem->cv_Ns; is++)\r
+           cv_mem->cv_ncfS1[is] = 0;\r
+         nefS = 0;\r
+       }*/\r
+\r
+       /* Note that ftempQ contains yQdot evaluated at the converged y\r
+        * (stored in cvQuadNls) and can be used in evaluating fQS */\r
+\r
+     /*  nflag = cvQuadSensNls(cv_mem);\r
+       kflag = cvHandleNFlag(cv_mem, &nflag, saved_t, &ncf, &(cv_mem->cv_ncfn));\r
+\r
+       if (kflag == PREDICT_AGAIN) continue;\r
+       if (kflag != DO_ERROR_TEST) return(kflag);*/\r
+\r
+       /* Error test on quadrature sensitivities */\r
+     /*  if (cv_mem->cv_errconQS) {\r
+         cv_mem->cv_acnrmQS = cvQuadSensNorm(cv_mem, cv_mem->cv_acorQS,\r
+                                             cv_mem->cv_ewtQS);\r
+         eflag = cvDoErrorTest(cv_mem, &nflag, saved_t, cv_mem->cv_acnrmQS,\r
+                               &nefQS, &(cv_mem->cv_netfQS), &dsmQS);\r
+\r
+         if (eflag == TRY_AGAIN) continue;\r
+         if (eflag != CV_SUCCESS) return(eflag);*/\r
+\r
+         /* Set dsm = max(dsm, dsmQS) to be used in cvPrepareNextStep */\r
+       /*  if (dsmQS > dsm) dsm = dsmQS;\r
+       }\r
+\r
+\r
+     }*/\r
+\r
+\r
+     /* Everything went fine; exit loop */ \r
+     break;\r
+\r
+   }\r
+\r
+   /* Nonlinear system solve and error test were both successful.\r
+      Update data, and consider change of step and/or order.       */\r
+\r
+   /*cvCompleteStep(cv_mem); \r
+\r
+   cvPrepareNextStep(cv_mem, dsm);*/ \r
+\r
+   /* If Stablilty Limit Detection is turned on, call stability limit\r
+      detection routine for possible order reduction. */\r
+\r
+  /* if (cv_mem->cv_sldeton) cvBDFStab(cv_mem);\r
+\r
+   cv_mem->cv_etamax = (cv_mem->cv_nst <= SMALL_NST) ? ETAMX2 : ETAMX3;*/\r
+\r
+   /*  Finally, we rescale the acor array to be the \r
+       estimated local error vector. */\r
+\r
+ /*  N_VScale(cv_mem->cv_tq[2], cv_mem->cv_acor, cv_mem->cv_acor);\r
+\r
+   if (cv_mem->cv_quadr)\r
+     N_VScale(cv_mem->cv_tq[2], cv_mem->cv_acorQ, cv_mem->cv_acorQ);\r
+\r
+   if (cv_mem->cv_sensi)\r
+     for (is=0; is<cv_mem->cv_Ns; is++)\r
+       N_VScale(cv_mem->cv_tq[2], cv_mem->cv_acorS[is], cv_mem->cv_acorS[is]);\r
+\r
+   if (cv_mem->cv_quadr_sensi)\r
+     for (is=0; is<cv_mem->cv_Ns; is++)\r
+       N_VScale(cv_mem->cv_tq[2], cv_mem->cv_acorQS[is], cv_mem->cv_acorQS[is]);*/\r
+\r
+   return(CV_SUCCESS);\r
+       \r
+ }\r
\r
+ /*\r
+  * cvAdjustParams\r
+  *\r
+  * This routine is called when a change in step size was decided upon,\r
+  * and it handles the required adjustments to the history array zn.\r
+  * If there is to be a change in order, we call cvAdjustOrder and reset\r
+  * q, L = q+1, and qwait.  Then in any case, we call cvRescale, which\r
+  * resets h and rescales the Nordsieck array.\r
+  */\r
+\r
+ private void cvAdjustParams(CVodeMemRec cv_mem)\r
+ {\r
+   if (cv_mem.cv_qprime != cv_mem.cv_q) {\r
+     cvAdjustOrder(cv_mem, cv_mem.cv_qprime-cv_mem.cv_q);\r
+     cv_mem.cv_q = cv_mem.cv_qprime;\r
+     cv_mem.cv_L = cv_mem.cv_q+1;\r
+     cv_mem.cv_qwait = cv_mem.cv_L;\r
+   }\r
+   cvRescale(cv_mem);\r
+ }\r
\r
+ /*\r
+  * cvAdjustOrder\r
+  *\r
+  * This routine is a high level routine which handles an order\r
+  * change by an amount deltaq (= +1 or -1). If a decrease in order\r
+  * is requested and q==2, then the routine returns immediately.\r
+  * Otherwise cvAdjustAdams or cvAdjustBDF is called to handle the\r
+  * order change (depending on the value of lmm).\r
+  */\r
+\r
+ private void cvAdjustOrder(CVodeMemRec cv_mem, int deltaq)\r
+ {\r
+   if ((cv_mem.cv_q==2) && (deltaq != 1)) return;\r
+   \r
+   switch(cv_mem.cv_lmm){\r
+   case CV_ADAMS:\r
+     cvAdjustAdams(cv_mem, deltaq);\r
+     break;\r
+   case CV_BDF:\r
+     cvAdjustBDF(cv_mem, deltaq);\r
+     break;\r
+   }\r
+ }\r
+\r
+ /*\r
+  * cvAdjustAdams\r
+  *\r
+  * This routine adjusts the history array on a change of order q by\r
+  * deltaq, in the case that lmm == CV_ADAMS.\r
+  */\r
+\r
+ private void cvAdjustAdams(CVodeMemRec cv_mem, int deltaq)\r
+ {\r
+   int i, j;\r
+   int is;\r
+   double xi, hsum;\r
+\r
+   /* On an order increase, set new column of zn to zero and return */\r
+   \r
+   if (deltaq==1) {\r
+     N_VConst_Serial(ZERO, cv_mem.cv_zn[cv_mem.cv_L]);\r
+     if (cv_mem.cv_quadr)\r
+       N_VConst_Serial(ZERO, cv_mem.cv_znQ[cv_mem.cv_L]);\r
+     if (cv_mem.cv_sensi)\r
+       for (is=0; is<cv_mem.cv_Ns; is++)\r
+         N_VConst_Serial(ZERO, cv_mem.cv_znS[cv_mem.cv_L][is]);\r
+     return;\r
+   }\r
+\r
+   /*\r
+    * On an order decrease, each zn[j] is adjusted by a multiple of zn[q].\r
+    * The coeffs. in the adjustment are the coeffs. of the polynomial:\r
+    *        x\r
+    * q * INT { u * ( u + xi_1 ) * ... * ( u + xi_{q-2} ) } du \r
+    *        0\r
+    * where xi_j = [t_n - t_(n-j)]/h => xi_0 = 0\r
+    */\r
+\r
+   for (i=0; i <= cv_mem.cv_qmax; i++) cv_mem.cv_l[i] = ZERO;\r
+   cv_mem.cv_l[1] = ONE;\r
+   hsum = ZERO;\r
+   for (j=1; j <= cv_mem.cv_q-2; j++) {\r
+     hsum += cv_mem.cv_tau[j];\r
+     xi = hsum / cv_mem.cv_hscale;\r
+     for (i=j+1; i >= 1; i--)\r
+       cv_mem.cv_l[i] = cv_mem.cv_l[i]*xi + cv_mem.cv_l[i-1];\r
+   }\r
+   \r
+   for (j=1; j <= cv_mem.cv_q-2; j++)\r
+     cv_mem.cv_l[j+1] = cv_mem.cv_q * (cv_mem.cv_l[j] / (j+1));\r
+   \r
+   for (j=2; j < cv_mem.cv_q; j++)\r
+     N_VLinearSum_Serial(-cv_mem.cv_l[j], cv_mem.cv_zn[cv_mem.cv_q], ONE,\r
+                  cv_mem.cv_zn[j], cv_mem.cv_zn[j]);\r
+\r
+   if (cv_mem.cv_quadr)\r
+     for (j=2; j < cv_mem.cv_q; j++)\r
+       N_VLinearSum_Serial(-cv_mem.cv_l[j], cv_mem.cv_znQ[cv_mem.cv_q], ONE,\r
+                    cv_mem.cv_znQ[j], cv_mem.cv_znQ[j]);\r
+\r
+   if (cv_mem.cv_sensi)\r
+     for (is=0; is<cv_mem.cv_Ns; is++)\r
+       for (j=2; j < cv_mem.cv_q; j++)\r
+         N_VLinearSum_Serial(-cv_mem.cv_l[j], cv_mem.cv_znS[cv_mem.cv_q][is],\r
+                      ONE, cv_mem.cv_znS[j][is], cv_mem.cv_znS[j][is]);\r
+\r
+ }\r
+\r
+ /*\r
+  * cvAdjustBDF\r
+  *\r
+  * This is a high level routine which handles adjustments to the\r
+  * history array on a change of order by deltaq in the case that \r
+  * lmm == CV_BDF.  cvAdjustBDF calls cvIncreaseBDF if deltaq = +1 and \r
+  * cvDecreaseBDF if deltaq = -1 to do the actual work.\r
+  */\r
+\r
+ private void cvAdjustBDF(CVodeMemRec cv_mem, int deltaq)\r
+ {\r
+   switch(deltaq) {\r
+   case 1: \r
+     cvIncreaseBDF(cv_mem);\r
+     return;\r
+   case -1: \r
+     cvDecreaseBDF(cv_mem);\r
+     return;\r
+   }\r
+ }\r
+\r
+ /*\r
+  * cvIncreaseBDF\r
+  *\r
+  * This routine adjusts the history array on an increase in the \r
+  * order q in the case that lmm == CV_BDF.  \r
+  * A new column zn[q+1] is set equal to a multiple of the saved \r
+  * vector (= acor) in zn[indx_acor].  Then each zn[j] is adjusted by\r
+  * a multiple of zn[q+1].  The coefficients in the adjustment are the \r
+  * coefficients of the polynomial x*x*(x+xi_1)*...*(x+xi_j),\r
+  * where xi_j = [t_n - t_(n-j)]/h.\r
+  */\r
+\r
+ private void cvIncreaseBDF(CVodeMemRec cv_mem)\r
+ {\r
+   double alpha0, alpha1, prod, xi, xiold, hsum, A1;\r
+   int i, j;\r
+   int is;\r
+\r
+   for (i=0; i <= cv_mem.cv_qmax; i++)\r
+     cv_mem.cv_l[i] = ZERO;\r
+   cv_mem.cv_l[2] = alpha1 = prod = xiold = ONE;\r
+   alpha0 = -ONE;\r
+   hsum = cv_mem.cv_hscale;\r
+   if (cv_mem.cv_q > 1) {\r
+     for (j=1; j < cv_mem.cv_q; j++) {\r
+       hsum += cv_mem.cv_tau[j+1];\r
+       xi = hsum / cv_mem.cv_hscale;\r
+       prod *= xi;\r
+       alpha0 -= ONE / (j+1);\r
+       alpha1 += ONE / xi;\r
+       for (i=j+2; i >= 2; i--)\r
+         cv_mem.cv_l[i] = cv_mem.cv_l[i]*xiold + cv_mem.cv_l[i-1];\r
+       xiold = xi;\r
+     }\r
+   }\r
+   A1 = (-alpha0 - alpha1) / prod;\r
+\r
+   /* \r
+      zn[indx_acor] contains the value Delta_n = y_n - y_n(0) \r
+      This value was stored there at the previous successful\r
+      step (in cvCompleteStep) \r
+      \r
+      A1 contains dbar = (1/xi* - 1/xi_q)/prod(xi_j)\r
+   */\r
+   \r
+   N_VScale_Serial(A1, cv_mem.cv_zn[cv_mem.cv_indx_acor], cv_mem.cv_zn[cv_mem.cv_L]);\r
+   for (j=2; j <= cv_mem.cv_q; j++)\r
+     N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_zn[cv_mem.cv_L],\r
+                  ONE, cv_mem.cv_zn[j], cv_mem.cv_zn[j]);\r
+\r
+   if (cv_mem.cv_quadr) {\r
+     N_VScale_Serial(A1, cv_mem.cv_znQ[cv_mem.cv_indx_acor], cv_mem.cv_znQ[cv_mem.cv_L]);\r
+     for (j=2; j <= cv_mem.cv_q; j++)\r
+       N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_znQ[cv_mem.cv_L],\r
+                    ONE, cv_mem.cv_znQ[j], cv_mem.cv_znQ[j]);\r
+   }\r
+\r
+   if (cv_mem.cv_sensi) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+       N_VScale_Serial(A1, cv_mem.cv_znS[cv_mem.cv_indx_acor][is],\r
+                cv_mem.cv_znS[cv_mem.cv_L][is]);\r
+       for (j=2; j <= cv_mem.cv_q; j++)\r
+         N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_znS[cv_mem.cv_L][is],\r
+                      ONE, cv_mem.cv_znS[j][is], cv_mem.cv_znS[j][is]);\r
+     }\r
+   }\r
+\r
+   if (cv_mem.cv_quadr_sensi) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+       N_VScale_Serial(A1, cv_mem.cv_znQS[cv_mem.cv_indx_acor][is],\r
+                cv_mem.cv_znQS[cv_mem.cv_L][is]);\r
+       for (j=2; j <= cv_mem.cv_q; j++)\r
+         N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_znQS[cv_mem.cv_L][is],\r
+                      ONE, cv_mem.cv_znQS[j][is], cv_mem.cv_znQS[j][is]);\r
+     }\r
+   }\r
+\r
+ }\r
+\r
+ /*\r
+  * cvDecreaseBDF\r
+  *\r
+  * This routine adjusts the history array on a decrease in the \r
+  * order q in the case that lmm == CV_BDF.  \r
+  * Each zn[j] is adjusted by a multiple of zn[q].  The coefficients\r
+  * in the adjustment are the coefficients of the polynomial\r
+  *   x*x*(x+xi_1)*...*(x+xi_j), where xi_j = [t_n - t_(n-j)]/h.\r
+  */\r
+\r
+ private void cvDecreaseBDF(CVodeMemRec cv_mem)\r
+ {\r
+   double hsum, xi;\r
+   int i, j;\r
+   int is;\r
+   \r
+   for (i=0; i <= cv_mem.cv_qmax; i++)\r
+     cv_mem.cv_l[i] = ZERO;\r
+   cv_mem.cv_l[2] = ONE;\r
+   hsum = ZERO;\r
+   for (j=1; j <= cv_mem.cv_q-2; j++) {\r
+     hsum += cv_mem.cv_tau[j];\r
+     xi = hsum / cv_mem.cv_hscale;\r
+     for (i=j+2; i >= 2; i--)\r
+       cv_mem.cv_l[i] = cv_mem.cv_l[i]*xi + cv_mem.cv_l[i-1];\r
+   }\r
+   \r
+   for (j=2; j < cv_mem.cv_q; j++)\r
+     N_VLinearSum_Serial(-cv_mem.cv_l[j], cv_mem.cv_zn[cv_mem.cv_q],\r
+                  ONE, cv_mem.cv_zn[j], cv_mem.cv_zn[j]);\r
+\r
+   if (cv_mem.cv_quadr) {\r
+     for (j=2; j < cv_mem.cv_q; j++)\r
+       N_VLinearSum_Serial(-cv_mem.cv_l[j], cv_mem.cv_znQ[cv_mem.cv_q],\r
+                    ONE, cv_mem.cv_znQ[j], cv_mem.cv_znQ[j]);\r
+   }\r
+\r
+   if (cv_mem.cv_sensi) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) \r
+       for (j=2; j < cv_mem.cv_q; j++)\r
+         N_VLinearSum_Serial(-cv_mem.cv_l[j], cv_mem.cv_znS[cv_mem.cv_q][is],\r
+                      ONE, cv_mem.cv_znS[j][is], cv_mem.cv_znS[j][is]);\r
+   }\r
+\r
+   if (cv_mem.cv_quadr_sensi) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) \r
+       for (j=2; j < cv_mem.cv_q; j++)\r
+         N_VLinearSum_Serial(-cv_mem.cv_l[j], cv_mem.cv_znQS[cv_mem.cv_q][is],\r
+                      ONE, cv_mem.cv_znQS[j][is], cv_mem.cv_znQS[j][is]);\r
+   }\r
+ }\r
+\r
+ /*\r
+  * cvRescale\r
+  *\r
+  * This routine rescales the Nordsieck array by multiplying the\r
+  * jth column zn[j] by eta^j, j = 1, ..., q.  Then the value of\r
+  * h is rescaled by eta, and hscale is reset to h.\r
+  */\r
+\r
+ private void cvRescale(CVodeMemRec cv_mem)\r
+ {\r
+   int j;\r
+   int is;\r
+   double factor;\r
+\r
+   factor = cv_mem.cv_eta;\r
+   for (j=1; j <= cv_mem.cv_q; j++) {\r
+\r
+     N_VScale_Serial(factor, cv_mem.cv_zn[j], cv_mem.cv_zn[j]);\r
+\r
+     if (cv_mem.cv_quadr)\r
+       N_VScale_Serial(factor, cv_mem.cv_znQ[j], cv_mem.cv_znQ[j]);\r
+\r
+     if (cv_mem.cv_sensi)\r
+       for (is=0; is<cv_mem.cv_Ns; is++)\r
+         N_VScale_Serial(factor, cv_mem.cv_znS[j][is], cv_mem.cv_znS[j][is]);\r
+\r
+     if (cv_mem.cv_quadr_sensi)\r
+       for (is=0; is<cv_mem.cv_Ns; is++)\r
+         N_VScale_Serial(factor, cv_mem.cv_znQS[j][is], cv_mem.cv_znQS[j][is]);\r
+\r
+     factor *= cv_mem.cv_eta;\r
+\r
+   }\r
+\r
+   cv_mem.cv_h = cv_mem.cv_hscale * cv_mem.cv_eta;\r
+   cv_mem.cv_next_h = cv_mem.cv_h;\r
+   cv_mem.cv_hscale = cv_mem.cv_h;\r
+   cv_mem.cv_nscon = 0;\r
+\r
+ }\r
\r
+ /*\r
+  * cvPredict\r
+  *\r
+  * This routine advances tn by the tentative step size h, and computes\r
+  * the predicted array z_n(0), which is overwritten on zn.  The\r
+  * prediction of zn is done by repeated additions.\r
+  * If tstop is enabled, it is possible for tn + h to be past tstop by roundoff,\r
+  * and in that case, we reset tn (after incrementing by h) to tstop.\r
+  */\r
+\r
+ private void cvPredict(CVodeMemRec cv_mem)\r
+ {\r
+   int j, k;\r
+   int is;\r
+\r
+   cv_mem.cv_tn += cv_mem.cv_h;\r
+   if (cv_mem.cv_tstopset) {\r
+     if ((cv_mem.cv_tn - cv_mem.cv_tstop)*cv_mem.cv_h > ZERO)\r
+       cv_mem.cv_tn = cv_mem.cv_tstop;\r
+   }\r
+\r
+   for (k = 1; k <= cv_mem.cv_q; k++)\r
+     for (j = cv_mem.cv_q; j >= k; j--) \r
+       N_VLinearSum_Serial(ONE, cv_mem.cv_zn[j-1], ONE,\r
+                    cv_mem.cv_zn[j], cv_mem.cv_zn[j-1]); \r
+\r
+   if (cv_mem.cv_quadr) {\r
+     for (k = 1; k <= cv_mem.cv_q; k++)\r
+       for (j = cv_mem.cv_q; j >= k; j--) \r
+         N_VLinearSum_Serial(ONE, cv_mem.cv_znQ[j-1], ONE,\r
+                      cv_mem.cv_znQ[j], cv_mem.cv_znQ[j-1]);\r
+   }\r
+\r
+   if (cv_mem.cv_sensi) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+       for (k = 1; k <= cv_mem.cv_q; k++)\r
+         for (j = cv_mem.cv_q; j >= k; j--) \r
+           N_VLinearSum_Serial(ONE, cv_mem.cv_znS[j-1][is], ONE,\r
+                        cv_mem.cv_znS[j][is], cv_mem.cv_znS[j-1][is]);\r
+     }\r
+   }\r
+\r
+   if (cv_mem.cv_quadr_sensi) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++) {\r
+       for (k = 1; k <= cv_mem.cv_q; k++)\r
+         for (j = cv_mem.cv_q; j >= k; j--) \r
+           N_VLinearSum_Serial(ONE, cv_mem.cv_znQS[j-1][is], ONE,\r
+                        cv_mem.cv_znQS[j][is], cv_mem.cv_znQS[j-1][is]);\r
+     }\r
+   }\r
+\r
+ }\r
+\r
+ /*\r
+  * cvSet\r
+  *\r
+  * This routine is a high level routine which calls cvSetAdams or\r
+  * cvSetBDF to set the polynomial l, the test quantity array tq, \r
+  * and the related variables  rl1, gamma, and gamrat.\r
+  *\r
+  * The array tq is loaded with constants used in the control of estimated\r
+  * local errors and in the nonlinear convergence test.  Specifically, while\r
+  * running at order q, the components of tq are as follows:\r
+  *   tq[1] = a coefficient used to get the est. local error at order q-1\r
+  *   tq[2] = a coefficient used to get the est. local error at order q\r
+  *   tq[3] = a coefficient used to get the est. local error at order q+1\r
+  *   tq[4] = constant used in nonlinear iteration convergence test\r
+  *   tq[5] = coefficient used to get the order q+2 derivative vector used in\r
+  *           the est. local error at order q+1\r
+  */\r
+\r
+ private void cvSet(CVodeMemRec cv_mem)\r
+ {\r
+   switch(cv_mem.cv_lmm) {\r
+   case CV_ADAMS:\r
+     cvSetAdams(cv_mem);\r
+     break;\r
+   case CV_BDF:\r
+     cvSetBDF(cv_mem);\r
+     break;\r
+   }\r
+   cv_mem.cv_rl1 = ONE / cv_mem.cv_l[1];\r
+   cv_mem.cv_gamma = cv_mem.cv_h * cv_mem.cv_rl1;\r
+   if (cv_mem.cv_nst == 0) cv_mem.cv_gammap = cv_mem.cv_gamma;\r
+   cv_mem.cv_gamrat = (cv_mem.cv_nst > 0) ?\r
+     cv_mem.cv_gamma / cv_mem.cv_gammap : ONE;  /* protect x / x != 1.0 */\r
+ }\r
\r
+ /*\r
+  * cvSetAdams\r
+  *\r
+  * This routine handles the computation of l and tq for the\r
+  * case lmm == CV_ADAMS.\r
+  *\r
+  * The components of the array l are the coefficients of a\r
+  * polynomial Lambda(x) = l_0 + l_1 x + ... + l_q x^q, given by\r
+  *                          q-1\r
+  * (d/dx) Lambda(x) = c * PRODUCT (1 + x / xi_i) , where\r
+  *                          i=1\r
+  *  Lambda(-1) = 0, Lambda(0) = 1, and c is a normalization factor.\r
+  * Here xi_i = [t_n - t_(n-i)] / h.\r
+  *\r
+  * The array tq is set to test quantities used in the convergence\r
+  * test, the error test, and the selection of h at a new order.\r
+  */\r
+\r
+ private void cvSetAdams(CVodeMemRec cv_mem)\r
+ {\r
+   double m[] = new double[L_MAX];\r
+   double M[] = new double[3];\r
+   double hsum;\r
+   \r
+   if (cv_mem.cv_q == 1) {\r
+     cv_mem.cv_l[0] = cv_mem.cv_l[1] = cv_mem.cv_tq[1] = cv_mem.cv_tq[5] = ONE;\r
+     cv_mem.cv_tq[2] = HALF;\r
+     cv_mem.cv_tq[3] = ONE/TWELVE;\r
+     cv_mem.cv_tq[4] = cv_mem.cv_nlscoef / cv_mem.cv_tq[2];       /* = 0.1 / tq[2] */\r
+     return;\r
+   }\r
+   \r
+   hsum = cvAdamsStart(cv_mem, m);\r
+   \r
+   M[0] = cvAltSum(cv_mem.cv_q-1, m, 1);\r
+   M[1] = cvAltSum(cv_mem.cv_q-1, m, 2);\r
+   \r
+   cvAdamsFinish(cv_mem, m, M, hsum);\r
+ }\r
+\r
+ /*\r
+  * cvAdamsStart\r
+  *\r
+  * This routine generates in m[] the coefficients of the product\r
+  * polynomial needed for the Adams l and tq coefficients for q > 1.\r
+  */\r
+\r
+ private double cvAdamsStart(CVodeMemRec cv_mem, double m[])\r
+ {\r
+   double hsum, xi_inv, sum;\r
+   int i, j;\r
+   \r
+   hsum = cv_mem.cv_h;\r
+   m[0] = ONE;\r
+   for (i=1; i <= cv_mem.cv_q; i++) m[i] = ZERO;\r
+   for (j=1; j < cv_mem.cv_q; j++) {\r
+     if ((j==cv_mem.cv_q-1) && (cv_mem.cv_qwait == 1)) {\r
+       sum = cvAltSum(cv_mem.cv_q-2, m, 2);\r
+       cv_mem.cv_tq[1] = cv_mem.cv_q * sum / m[cv_mem.cv_q-2];\r
+     }\r
+     xi_inv = cv_mem.cv_h / hsum;\r
+     for (i=j; i >= 1; i--)\r
+       m[i] += m[i-1] * xi_inv;\r
+     hsum += cv_mem.cv_tau[j];\r
+     /* The m[i] are coefficients of product(1 to j) (1 + x/xi_i) */\r
+   }\r
+   return(hsum);\r
+ }\r
+\r
+ /*\r
+  * cvAdamsFinish\r
+  *\r
+  * This routine completes the calculation of the Adams l and tq.\r
+  */\r
+\r
+ private void cvAdamsFinish(CVodeMemRec cv_mem, double m[], double M[], double hsum)\r
+ {\r
+   int i;\r
+   double M0_inv, xi, xi_inv;\r
+   \r
+   M0_inv = ONE / M[0];\r
+   \r
+   cv_mem.cv_l[0] = ONE;\r
+   for (i=1; i <= cv_mem.cv_q; i++)\r
+     cv_mem.cv_l[i] = M0_inv * (m[i-1] / i);\r
+   xi = hsum / cv_mem.cv_h;\r
+   xi_inv = ONE / xi;\r
+   \r
+   cv_mem.cv_tq[2] = M[1] * M0_inv / xi;\r
+   cv_mem.cv_tq[5] = xi / cv_mem.cv_l[cv_mem.cv_q];\r
+\r
+   if (cv_mem.cv_qwait == 1) {\r
+     for (i=cv_mem.cv_q; i >= 1; i--)\r
+       m[i] += m[i-1] * xi_inv;\r
+     M[2] = cvAltSum(cv_mem.cv_q, m, 2);\r
+     cv_mem.cv_tq[3] = M[2] * M0_inv / cv_mem.cv_L;\r
+   }\r
+\r
+   cv_mem.cv_tq[4] = cv_mem.cv_nlscoef / cv_mem.cv_tq[2];\r
+ }\r
+\r
+ /*  \r
+  * cvAltSum\r
+  *\r
+  * cvAltSum returns the value of the alternating sum\r
+  *   sum (i= 0 ... iend) [ (-1)^i * (a[i] / (i + k)) ].\r
+  * If iend < 0 then cvAltSum returns 0.\r
+  * This operation is needed to compute the integral, from -1 to 0,\r
+  * of a polynomial x^(k-1) M(x) given the coefficients of M(x).\r
+  */\r
+\r
+ private double cvAltSum(int iend, double a[], int k)\r
+ {\r
+   int i, sign;\r
+   double sum;\r
+   \r
+   if (iend < 0) return(ZERO);\r
+   \r
+   sum = ZERO;\r
+   sign = 1;\r
+   for (i=0; i <= iend; i++) {\r
+     sum += sign * (a[i] / (i+k));\r
+     sign = -sign;\r
+   }\r
+   return(sum);\r
+ }\r
+\r
+ /*\r
+  * cvSetBDF\r
+  *\r
+  * This routine computes the coefficients l and tq in the case\r
+  * lmm == CV_BDF.  cvSetBDF calls cvSetTqBDF to set the test\r
+  * quantity array tq. \r
+  * \r
+  * The components of the array l are the coefficients of a\r
+  * polynomial Lambda(x) = l_0 + l_1 x + ... + l_q x^q, given by\r
+  *                                 q-1\r
+  * Lambda(x) = (1 + x / xi*_q) * PRODUCT (1 + x / xi_i) , where\r
+  *                                 i=1\r
+  *  xi_i = [t_n - t_(n-i)] / h.\r
+  *\r
+  * The array tq is set to test quantities used in the convergence\r
+  * test, the error test, and the selection of h at a new order.\r
+  */\r
+\r
+ private void cvSetBDF(CVodeMemRec cv_mem)\r
+ {\r
+   double alpha0, alpha0_hat, xi_inv, xistar_inv, hsum;\r
+   int i,j;\r
+   \r
+   cv_mem.cv_l[0] = cv_mem.cv_l[1] = xi_inv = xistar_inv = ONE;\r
+   for (i=2; i <= cv_mem.cv_q; i++) cv_mem.cv_l[i] = ZERO;\r
+   alpha0 = alpha0_hat = -ONE;\r
+   hsum = cv_mem.cv_h;\r
+   if (cv_mem.cv_q > 1) {\r
+     for (j=2; j < cv_mem.cv_q; j++) {\r
+       hsum += cv_mem.cv_tau[j-1];\r
+       xi_inv = cv_mem.cv_h / hsum;\r
+       alpha0 -= ONE / j;\r
+       for (i=j; i >= 1; i--)\r
+         cv_mem.cv_l[i] += cv_mem.cv_l[i-1]*xi_inv;\r
+       /* The l[i] are coefficients of product(1 to j) (1 + x/xi_i) */\r
+     }\r
+     \r
+     /* j = q */\r
+     alpha0 -= ONE / cv_mem.cv_q;\r
+     xistar_inv = -cv_mem.cv_l[1] - alpha0;\r
+     hsum += cv_mem.cv_tau[cv_mem.cv_q-1];\r
+     xi_inv = cv_mem.cv_h / hsum;\r
+     alpha0_hat = -cv_mem.cv_l[1] - xi_inv;\r
+     for (i=cv_mem.cv_q; i >= 1; i--)\r
+       cv_mem.cv_l[i] += cv_mem.cv_l[i-1]*xistar_inv;\r
+   }\r
+\r
+   cvSetTqBDF(cv_mem, hsum, alpha0, alpha0_hat, xi_inv, xistar_inv);\r
+ }\r
+\r
+ /*\r
+  * cvSetTqBDF\r
+  *\r
+  * This routine sets the test quantity array tq in the case\r
+  * lmm == CV_BDF.\r
+  */\r
+\r
+ private void cvSetTqBDF(CVodeMemRec cv_mem, double hsum, double alpha0,\r
+                        double alpha0_hat, double xi_inv, double xistar_inv)\r
+ {\r
+   double A1, A2, A3, A4, A5, A6;\r
+   double C, Cpinv, Cppinv;\r
+   \r
+   A1 = ONE - alpha0_hat + alpha0;\r
+   A2 = ONE + cv_mem.cv_q * A1;\r
+   cv_mem.cv_tq[2] = Math.abs(A1 / (alpha0 * A2));\r
+   cv_mem.cv_tq[5] = Math.abs(A2 * xistar_inv / (cv_mem.cv_l[cv_mem.cv_q] * xi_inv));\r
+   if (cv_mem.cv_qwait == 1) {\r
+     if (cv_mem.cv_q > 1) {\r
+       C = xistar_inv / cv_mem.cv_l[cv_mem.cv_q];\r
+       A3 = alpha0 + ONE / cv_mem.cv_q;\r
+       A4 = alpha0_hat + xi_inv;\r
+       Cpinv = (ONE - A4 + A3) / A3;\r
+       cv_mem.cv_tq[1] = Math.abs(C * Cpinv);\r
+     }\r
+     else cv_mem.cv_tq[1] = ONE;\r
+     hsum += cv_mem.cv_tau[cv_mem.cv_q];\r
+     xi_inv = cv_mem.cv_h / hsum;\r
+     A5 = alpha0 - (ONE / (cv_mem.cv_q+1));\r
+     A6 = alpha0_hat - xi_inv;\r
+     Cppinv = (ONE - A6 + A5) / A2;\r
+     cv_mem.cv_tq[3] = Math.abs(Cppinv / (xi_inv * (cv_mem.cv_q+2) * A5));\r
+   }\r
+   cv_mem.cv_tq[4] = cv_mem.cv_nlscoef / cv_mem.cv_tq[2];\r
+ }\r
+\r
+ /*\r
+  * cvNls\r
+  *\r
+  * This routine attempts to solve the nonlinear system associated\r
+  * with a single implicit step of the linear multistep method.\r
+  * Depending on iter, it calls cvNlsFunctional or cvNlsNewton\r
+  * to do the work.\r
+  */\r
+\r
+ private int cvNls(CVodeMemRec cv_mem, int nflag)\r
+ {\r
+   int flag = CV_SUCCESS;\r
+\r
+   switch(cv_mem.cv_iter) {\r
+   case CV_FUNCTIONAL:\r
+     flag = cvNlsFunctional(cv_mem);\r
+     break;\r
+   case CV_NEWTON:\r
+     flag = cvNlsNewton(cv_mem, nflag);\r
+     break;\r
+   }\r
+   \r
+   return(flag);\r
+\r
+ }\r
+\r
+ /*\r
+  * cvNlsFunctional\r
+  *\r
+  * This routine attempts to solve the nonlinear system using \r
+  * functional iteration (no matrices involved).\r
+  * \r
+  * This routine also handles the functional iteration of the\r
+  * combined system (states + sensitivities) when sensitivities are \r
+  * computed using the CV_SIMULTANEOUS approach.\r
+  *\r
+  * Possible return values are:\r
+  *\r
+  *   CV_SUCCESS       ---> continue with error test\r
+  *\r
+  *   CV_RHSFUNC_FAIL  -+   \r
+  *   CV_SRHSFUNC_FAIL -+-> halt the integration \r
+  *\r
+  *   CONV_FAIL        -+\r
+  *   RHSFUNC_RECVR     |-> predict again or stop if too many\r
+  *   SRHSFUNC_RECVR   -+\r
+  *\r
+  */\r
+\r
+ private int cvNlsFunctional(CVodeMemRec cv_mem)\r
+ {\r
+   int m;\r
+   double del, delS, Del, Delp, dcon;\r
+   int retval, is;\r
+   boolean do_sensi_sim;\r
+   NVector wrk1, wrk2;\r
+\r
+   /* Are we computing sensitivities with the CV_SIMULTANEOUS approach? */\r
+   do_sensi_sim = (cv_mem.cv_sensi && (cv_mem.cv_ism==CV_SIMULTANEOUS));\r
+\r
+   /* Initialize counter and evaluate f at predicted y */\r
+   cv_mem.cv_crate = ONE;\r
+   m = 0;\r
+\r
+   /* Initialize delS and Delp to avoid compiler warning message */\r
+   delS = Delp = ZERO;\r
+\r
+   retval = f(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                         cv_mem.cv_tempv, cv_mem.cv_user_data);\r
+   cv_mem.cv_nfe++;\r
+   if (retval < 0) return(CV_RHSFUNC_FAIL);\r
+   if (retval > 0) return(RHSFUNC_RECVR);\r
+\r
+   if (do_sensi_sim) {\r
+     wrk1 = cv_mem.cv_ftemp;\r
+     wrk2 = cv_mem.cv_ftempS[0];\r
+     retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                               cv_mem.cv_tempv, cv_mem.cv_znS[0],\r
+                               cv_mem.cv_tempvS, wrk1, wrk2);\r
+     if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+     if (retval > 0) return(SRHSFUNC_RECVR);\r
+   }\r
+\r
+   /* Initialize correction to zero */\r
+\r
+   N_VConst_Serial(ZERO, cv_mem.cv_acor);\r
+   if (do_sensi_sim) {\r
+     for (is=0; is<cv_mem.cv_Ns; is++)\r
+       N_VConst_Serial(ZERO, cv_mem.cv_acorS[is]);\r
+   }\r
+\r
+   /* Loop until convergence; accumulate corrections in acor */\r
+\r
+   for(;;) {\r
+\r
+     cv_mem.cv_nni++;\r
+\r
+     /* Correct y directly from the last f value */\r
+\r
+     N_VLinearSum_Serial(cv_mem.cv_h, cv_mem.cv_tempv, -ONE,\r
+                  cv_mem.cv_zn[1], cv_mem.cv_tempv);\r
+     N_VScale_Serial(cv_mem.cv_rl1, cv_mem.cv_tempv, cv_mem.cv_tempv);\r
+     N_VLinearSum_Serial(ONE, cv_mem.cv_zn[0], ONE, cv_mem.cv_tempv, cv_mem.cv_y);\r
+\r
+     if (do_sensi_sim)\r
+       for (is=0; is<cv_mem.cv_Ns; is++) {\r
+         N_VLinearSum_Serial(cv_mem.cv_h, cv_mem.cv_tempvS[is], -ONE,\r
+                      cv_mem.cv_znS[1][is], cv_mem.cv_tempvS[is]);\r
+         N_VScale_Serial(cv_mem.cv_rl1, cv_mem.cv_tempvS[is], cv_mem.cv_tempvS[is]);\r
+         N_VLinearSum_Serial(ONE, cv_mem.cv_znS[0][is], ONE,\r
+                      cv_mem.cv_tempvS[is], cv_mem.cv_yS[is]);\r
+       }\r
+     \r
+     /* Get WRMS norm of current correction to use in convergence test */\r
+\r
+     N_VLinearSum_Serial(ONE, cv_mem.cv_tempv, -ONE, cv_mem.cv_acor, cv_mem.cv_acor);\r
+     if (do_sensi_sim)\r
+       for (is=0; is<cv_mem.cv_Ns; is++)\r
+         N_VLinearSum_Serial(ONE, cv_mem.cv_tempvS[is], -ONE,\r
+                      cv_mem.cv_acorS[is], cv_mem.cv_acorS[is]);\r
+\r
+     del = N_VWrmsNorm_Serial(cv_mem.cv_acor, cv_mem.cv_ewt);\r
+     if (do_sensi_sim)\r
+       delS = cvSensUpdateNorm(cv_mem, del, cv_mem.cv_acorS, cv_mem.cv_ewtS);\r
+\r
+     /*N_VScale(ONE, cv_mem->cv_tempv, cv_mem->cv_acor);\r
+     if (do_sensi_sim) \r
+       for (is=0; is<cv_mem->cv_Ns; is++)\r
+         N_VScale(ONE, cv_mem->cv_tempvS[is], cv_mem->cv_acorS[is]);*/\r
+     \r
+     /* Test for convergence.  If m > 0, an estimate of the convergence\r
+        rate constant is stored in crate, and used in the test. \r
+\r
+        Recall that, even when errconS=SUNFALSE, all variables are used in the\r
+        convergence test. Hence, we use Del (and not del). However, acnrm\r
+        is used in the error test and thus it has different forms\r
+        depending on errconS (and this explains why we have to carry around\r
+        del and delS)\r
+     */\r
+     \r
+     /*Del = (do_sensi_sim) ? delS : del;\r
+     if (m > 0) cv_mem->cv_crate = SUNMAX(CRDOWN * cv_mem->cv_crate, Del / Delp);\r
+     dcon = Del * SUNMIN(ONE, cv_mem->cv_crate) / cv_mem->cv_tq[4];\r
+\r
+     if (dcon <= ONE) {\r
+       if (m == 0)\r
+         if (do_sensi_sim && cv_mem->cv_errconS)\r
+           cv_mem->cv_acnrm = delS;\r
+         else\r
+           cv_mem->cv_acnrm = del;\r
+       else {\r
+         cv_mem->cv_acnrm = N_VWrmsNorm(cv_mem->cv_acor, cv_mem->cv_ewt);\r
+         if (do_sensi_sim && cv_mem->cv_errconS)\r
+           cv_mem->cv_acnrm = cvSensUpdateNorm(cv_mem, cv_mem->cv_acnrm,\r
+                                               cv_mem->cv_acorS, cv_mem->cv_ewtS);\r
+       }\r
+       return(CV_SUCCESS);*/  /* Convergence achieved */\r
+     /*}*/\r
+\r
+     /* Stop at maxcor iterations or if iter. seems to be diverging */\r
+\r
+     /*m++;\r
+     if ((m==cv_mem->cv_maxcor) || ((m >= 2) && (Del > RDIV * Delp)))\r
+       return(CONV_FAIL);*/\r
+\r
+     /* Save norm of correction, evaluate f, and loop again */\r
+\r
+    /* Delp = Del;\r
+\r
+     retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_y,\r
+                           cv_mem->cv_tempv, cv_mem->cv_user_data);\r
+     cv_mem->cv_nfe++;\r
+     if (retval < 0) return(CV_RHSFUNC_FAIL);\r
+     if (retval > 0) return(RHSFUNC_RECVR);\r
+\r
+     if (do_sensi_sim) {\r
+       wrk1 = cv_mem->cv_ftemp;\r
+       wrk2 = cv_mem->cv_ftempS[0];\r
+       retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, cv_mem->cv_y,\r
+                                 cv_mem->cv_tempv, cv_mem->cv_yS,\r
+                                 cv_mem->cv_tempvS, wrk1, wrk2);\r
+       if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+       if (retval > 0) return(SRHSFUNC_RECVR);\r
+     }  */    \r
+\r
+   } /* end loop */\r
+\r
+ }\r
+\r
+ /*\r
+  * cvNlsNewton\r
+  *\r
+  * This routine handles the Newton iteration. It calls lsetup if \r
+  * indicated, calls cvNewtonIteration to perform the iteration, and \r
+  * retries a failed attempt at Newton iteration if that is indicated.\r
+  * See return values at top of this file.\r
+  *\r
+  * This routine also handles the Newton iteration of the combined \r
+  * system when sensitivities are computed using the CV_SIMULTANEOUS\r
+  * approach. Since in that case we use a quasi-Newton on the \r
+  * combined system (by approximating the Jacobian matrix by its\r
+  * block diagonal) and thus only solve linear systems with \r
+  * multiple right hand sides (all sharing the same coefficient\r
+  * matrix - whatever iteration matrix we decide on) we set-up\r
+  * the linear solver to handle N equations at a time.\r
+  *\r
+  * Possible return values:\r
+  *\r
+  *   CV_SUCCESS       ---> continue with error test\r
+  *\r
+  *   CV_RHSFUNC_FAIL  -+  \r
+  *   CV_LSETUP_FAIL    |\r
+  *   CV_LSOLVE_FAIL    |-> halt the integration \r
+  *   CV_SRHSFUNC_FAIL -+\r
+  *\r
+  *   CONV_FAIL        -+\r
+  *   RHSFUNC_RECVR     |-> predict again or stop if too many\r
+  *   SRHSFUNC_RECVR   -+\r
+  *\r
+  */\r
+\r
+ private int cvNlsNewton(CVodeMemRec cv_mem, int nflag)\r
+ {\r
+  /* N_Vector vtemp1, vtemp2, vtemp3, wrk1, wrk2;\r
+   int convfail, ier;\r
+   booleantype callSetup, do_sensi_sim;\r
+   int retval, is;*/\r
+   \r
+   /* Are we computing sensitivities with the CV_SIMULTANEOUS approach? */\r
+  /* do_sensi_sim = (cv_mem->cv_sensi && (cv_mem->cv_ism==CV_SIMULTANEOUS));\r
+\r
+   vtemp1 = cv_mem->cv_acor; */ /* rename acor as vtemp1 for readability  */\r
+  /* vtemp2 = cv_mem->cv_y;  */   /* rename y as vtemp2 for readability     */\r
+  /* vtemp3 = cv_mem->cv_tempv; *//* rename tempv as vtemp3 for readability */\r
+   \r
+   /* Set flag convfail, input to lsetup for its evaluation decision */\r
+  /* convfail = ((nflag == FIRST_CALL) || (nflag == PREV_ERR_FAIL)) ?\r
+     CV_NO_FAILURES : CV_FAIL_OTHER;*/\r
+\r
+   /* Decide whether or not to call setup routine (if one exists) */\r
+   /*if (cv_mem->cv_lsetup) {      \r
+     callSetup = (nflag == PREV_CONV_FAIL) || (nflag == PREV_ERR_FAIL) ||\r
+       (cv_mem->cv_nst == 0) ||\r
+       (cv_mem->cv_nst >= cv_mem->cv_nstlp + MSBP) ||\r
+       (SUNRabs(cv_mem->cv_gamrat-ONE) > DGMAX);*/\r
+\r
+     /* Decide whether to force a call to setup */\r
+    /* if (cv_mem->cv_forceSetup) {\r
+       callSetup = SUNTRUE;\r
+       convfail = CV_FAIL_OTHER;\r
+     }\r
+\r
+   } else {  \r
+     cv_mem->cv_crate = ONE;\r
+     cv_mem->cv_crateS = ONE; */ /* if NO lsetup all conv. rates are set to ONE */\r
+   /*  callSetup = SUNFALSE;\r
+   }*/\r
+   \r
+   /* Looping point for the solution of the nonlinear system.\r
+      Evaluate f at the predicted y, call lsetup if indicated, and\r
+      call cvNewtonIteration for the Newton iteration itself.      */\r
+\r
+   /*for(;;) {\r
+\r
+     retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_zn[0],\r
+                           cv_mem->cv_ftemp, cv_mem->cv_user_data);\r
+     cv_mem->cv_nfe++; \r
+     if (retval < 0) return(CV_RHSFUNC_FAIL);\r
+     if (retval > 0) return(RHSFUNC_RECVR);\r
+\r
+     if (do_sensi_sim) {\r
+       wrk1 = cv_mem->cv_tempv;\r
+       wrk2 = cv_mem->cv_tempvS[0];\r
+       retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, cv_mem->cv_zn[0],\r
+                                 cv_mem->cv_ftemp, cv_mem->cv_znS[0],\r
+                                 cv_mem->cv_ftempS, wrk1, wrk2);\r
+       if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+       if (retval > 0) return(SRHSFUNC_RECVR);\r
+     }\r
+\r
+     if (callSetup) {\r
+       ier = cv_mem->cv_lsetup(cv_mem, convfail, cv_mem->cv_zn[0],\r
+                               cv_mem->cv_ftemp, &(cv_mem->cv_jcur), \r
+                               vtemp1, vtemp2, vtemp3);\r
+       cv_mem->cv_nsetups++;\r
+       callSetup = SUNFALSE;\r
+       cv_mem->cv_forceSetup = SUNFALSE;\r
+       cv_mem->cv_gamrat = ONE; \r
+       cv_mem->cv_gammap = cv_mem->cv_gamma;\r
+       cv_mem->cv_crate = ONE;\r
+       cv_mem->cv_crateS = ONE;*/ /* after lsetup all conv. rates are reset to ONE */\r
+       /*cv_mem->cv_nstlp = cv_mem->cv_nst;*/\r
+       /* Return if lsetup failed */\r
+       /*if (ier < 0) return(CV_LSETUP_FAIL);\r
+       if (ier > 0) return(CONV_FAIL);\r
+     }*/\r
+\r
+     /* Set acor to zero and load prediction into y vector */\r
+    /* N_VConst(ZERO, cv_mem->cv_acor);\r
+     N_VScale(ONE, cv_mem->cv_zn[0], cv_mem->cv_y);\r
+\r
+     if (do_sensi_sim)\r
+       for (is=0; is<cv_mem->cv_Ns; is++) {\r
+         N_VConst(ZERO, cv_mem->cv_acorS[is]);\r
+         N_VScale(ONE, cv_mem->cv_znS[0][is], cv_mem->cv_yS[is]);\r
+       }*/\r
+\r
+     /* Do the Newton iteration */\r
+    /* ier = cvNewtonIteration(cv_mem);*/\r
+\r
+     /* If there is a convergence failure and the Jacobian-related \r
+        data appears not to be current, loop again with a call to lsetup\r
+        in which convfail=CV_FAIL_BAD_J.  Otherwise return.                 */\r
+    /* if (ier != TRY_AGAIN) return(ier);\r
+     \r
+     callSetup = SUNTRUE;\r
+     convfail = CV_FAIL_BAD_J;\r
+   }*/\r
+        return -100;\r
+ }\r
+\r
 \r
 }
\ No newline at end of file