Porting continues.
[mipav.git] / mipav / src / gov / nih / mipav / model / algorithms / CVODES.java
index 9a25348..0061a79 100644 (file)
@@ -207,18 +207,23 @@ public abstract class CVODES {
        final int MXNEF = 7;\r
        final double CORTES = 0.1;\r
        final double ZERO = 0.0;\r
+       final double TINY = 1.0E-10;\r
        final double PT1 = 0.1;\r
        final double POINT2 = 0.2;\r
+       final double FOURTH = 0.25;\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 THREE = 3.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 ETAMXF = 0.2;\r
-       final double ETAMX1 = 10000.0; \r
+       final double ETAMX1 = 10000.0;\r
+       final double ETAMX2 = 10.0;\r
+       final double ETAMX3 = 10.0;\r
        final double HUB_FACTOR = 0.1;\r
        final double HLB_FACTOR = 100.0;\r
        final double FUZZ_FACTOR = 100.0;\r
@@ -234,6 +239,10 @@ public abstract class CVODES {
     final int MSBP = 20;\r
     // ONEPSM (1+epsilon) used in testing if the step size is below its bound\r
     final double ONEPSM = 1.000001;\r
+    // THRESH      if eta < THRESH reject a change in step size or order\r
+    final double THRESH = 1.5;\r
+    // SMALL_NST   nst > SMALL_NST => use ETAMX3 \r
+    final int SMALL_NST = 10;\r
     final double ETAMIN = 0.1;\r
     // MXNEF1      max no. of error test failures before forcing a reduction of order\r
     // SMALL_NEF   if an error failure occurs and SMALL_NEF <= nef <= MXNEF1, then\r
@@ -636,11 +645,13 @@ public abstract class CVODES {
                abstol.setData(abstolr);\r
                CVodeMemRec cvode_mem;\r
                int flag;\r
+               int flagr;\r
                double A[][];\r
                SUNLinearSolver LS;\r
                double tout;\r
                int iout;\r
                double t[] = new double[1];\r
+               int rootsfound[] = new int[2];\r
                \r
                // Call CVodeCreate to create the solver memory and specify the\r
                // Backward Differentiation Formula and the use of a Newton\r
@@ -708,7 +719,39 @@ public abstract class CVODES {
                \r
                while (true) {\r
                        flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
+                       System.out.println("At t = " + t + " y[0] = " + y.data[0] + " y[1] = " + y.data[1] + " y[2] = " + y.data[2]);\r
+                       \r
+                       if (flag == CV_ROOT_RETURN) {\r
+                           flagr = CVodeGetRootInfo(cvode_mem, rootsfound)     ;\r
+                           if (flagr == CV_MEM_NULL) {\r
+                               return;\r
+                           }\r
+                           System.out.println("Roots found are " + rootsfound[0] + " and " + rootsfound[1]);\r
+                       }\r
+                       \r
+                       if (flag < 0) {\r
+                               System.out.println("CVode failed with flag = " + flag);\r
+                               break;\r
+                       }\r
+                       \r
+                       if (flag == CV_SUCCESS) {\r
+                               iout++;\r
+                               tout *= TMULT;\r
+                       }\r
+                       \r
+                       if (iout == NOUT) {\r
+                               break;\r
+                       }\r
                } // while (true)\r
+               \r
+               // Print some final statistics\r
+               PrintFinalStats(cvode_mem);\r
+       }\r
+       \r
+       private void PrintFinalStats(CVodeMemRec cv_mem) {\r
+           System.out.println("Number of integrations steps = " + cv_mem.cv_nst);\r
+           System.out.println("Number of calls to f = " + cv_mem.cv_nfe);\r
+           System.out.println("Number of calls to the linear solver setup routine = " + cv_mem.cv_nsetups);\r
        }\r
        \r
        private void N_VNew_Serial(NVector y, int length) {\r
@@ -1049,8 +1092,8 @@ public abstract class CVODES {
 \r
          long cv_netf[] = new long[1];        /* number of error test failures                   */\r
          long cv_netfQ[] = new long[1];       /* number of quadr. error test failures            */\r
-         long cv_netfS;       /* number of sensi. error test failures            */\r
-         long cv_netfQS;      /* number of quadr. sensi. error test failures     */\r
+         long cv_netfS[] = new long[1];       /* number of sensi. error test failures            */\r
+         long cv_netfQS[] = new long[1];      /* number of quadr. sensi. error test failures     */\r
 \r
          long cv_nsetups;     /* number of setup calls                           */\r
          long cv_nsetupsS;    /* number of setup calls due to sensitivities      */\r
@@ -4836,13 +4879,15 @@ else                return(snrm);
    double dsm[] = new double[1];\r
    double saved_t;\r
    double dsmQ[] = new double[1];\r
-   double dsmS, dsmQS;\r
+   double dsmS[] = new double[1];\r
+   double dsmQS[] = new double[1];\r
    boolean do_sensi_stg, do_sensi_stg1;\r
    int ncf[] = new int[1]; \r
    int ncfS[] = new int[1];\r
    int nef[] = new int[1];\r
    int nefQ[] = new int[1]; \r
-   int nefS, nefQS;\r
+   int nefS[] = new int[1];\r
+   int nefQS[] = new int[1];\r
    int nflag[] = new int [1];\r
    int kflag, eflag;\r
    int retval, is;\r
@@ -4855,8 +4900,8 @@ else                return(snrm);
    /* Initialize local counters for convergence and error test failures */\r
 \r
    ncf[0]  = nef[0]  = 0;\r
-   nefQ[0] = nefQS = 0;\r
-   ncfS[0] = nefS = 0;\r
+   nefQ[0] = nefQS[0] = 0;\r
+   ncfS[0] = nefS[0] = 0;\r
    if (do_sensi_stg1) {\r
         cv_mem.cv_ncfS1 = new int[cv_mem.cv_Ns][1];\r
      for (is=0; is<cv_mem.cv_Ns; is++)\r
@@ -4930,101 +4975,101 @@ else                return(snrm);
 \r
      /* ------ Correct the sensitivity variables (STAGGERED or STAGGERED1) ------- */\r
 \r
-     /*if (do_sensi_stg || do_sensi_stg1) { \r
+     if (do_sensi_stg || do_sensi_stg1) { \r
 \r
-       ncf[0] = nef[0] = 0;   */     /* reset counters for states     */\r
-       /*if (cv_mem.cv_quadr) nefQ = 0; */ /* reset counter for quadratures */\r
+       ncf[0] = nef[0] = 0;       /* reset counters for states     */\r
+       if (cv_mem.cv_quadr) nefQ[0] = 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 = f(cv_mem.cv_tn, cv_mem.cv_y,\r
+      retval = 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
+         nflag[0] = PREV_CONV_FAIL;\r
          continue;\r
        }\r
 \r
-       if (do_sensi_stg) { */\r
+       if (do_sensi_stg) {\r
          /* Nonlinear solve for sensitivities (all-at-once) */\r
-        /* nflag[0] = cvStgrNls(cv_mem);\r
+         nflag[0] = cvStgrNls(cv_mem);\r
          kflag = cvHandleNFlag(cv_mem, nflag, saved_t, ncfS,\r
                                cv_mem.cv_ncfnS);\r
-       } else { */\r
+       } else { \r
          /* Nonlinear solve for sensitivities (one-by-one) */\r
-         /*for (is=0; is<cv_mem.cv_Ns; is++) { \r
+         for (is=0; is<cv_mem.cv_Ns; is++) { \r
            nflag[0] = 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
+                                 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
+       if (kflag != DO_ERROR_TEST) return(kflag); \r
 \r
        /* Error test on sensitivities */\r
-     /*  if (cv_mem->cv_errconS) {\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
+           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
+         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
+         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
+         if (dsmS[0] > dsm[0]) dsm[0] = dsmS[0];\r
 \r
        }\r
 \r
-     } */\r
+     } \r
 \r
      /* ------ Correct the quadrature sensitivity variables ------ */\r
 \r
-   /*  if (cv_mem->cv_quadr_sensi) { */\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
+       ncf[0] = nef[0] = 0;\r
+       if (cv_mem.cv_quadr) nefQ[0] = 0;\r
+       if (do_sensi_stg) ncfS[0] = nefS[0] = 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
+         for (is=0; is<cv_mem.cv_Ns; is++)\r
+           cv_mem.cv_ncfS1[is][0] = 0;\r
+         nefS[0] = 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
+       nflag[0] = 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
+       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
+       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
+         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
+         if (dsmQS[0] > dsm[0]) dsm[0] = dsmQS[0];\r
        }\r
 \r
 \r
-     }*/\r
+     }\r
 \r
 \r
      /* Everything went fine; exit loop */ \r
@@ -5035,32 +5080,32 @@ else                return(snrm);
    /* 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
+   cvCompleteStep(cv_mem); \r
 \r
-   cvPrepareNextStep(cv_mem, dsm);*/ \r
+   cvPrepareNextStep(cv_mem, dsm[0]);\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
+   if (cv_mem.cv_sldeton) cvBDFStab(cv_mem);\r
 \r
-   cv_mem->cv_etamax = (cv_mem->cv_nst <= SMALL_NST) ? ETAMX2 : ETAMX3;*/\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
+   N_VScale_Serial(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
+   if (cv_mem.cv_quadr)\r
+     N_VScale_Serial(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
+   if (cv_mem.cv_sensi)\r
+     for (is=0; is<cv_mem.cv_Ns; is++)\r
+       N_VScale_Serial(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
+   if (cv_mem.cv_quadr_sensi)\r
+     for (is=0; is<cv_mem.cv_Ns; is++)\r
+       N_VScale_Serial(cv_mem.cv_tq[2], cv_mem.cv_acorQS[is], cv_mem.cv_acorQS[is]);\r
 \r
    return(CV_SUCCESS);\r
        \r
@@ -7842,6 +7887,762 @@ else                return(snrm);
           \r
           return(0);\r
         }\r
+        \r
+        /*\r
+         * cvQuadSensNls\r
+         * \r
+         * This routine solves for the quadrature sensitivity variables\r
+         * at the new step. It does not solve a nonlinear system, but \r
+         * rather updates the quadrature variables. The name for this\r
+         * function is just for uniformity purposes.\r
+         *\r
+         * Possible return values (interpreted by cvHandleNFlag)\r
+         *\r
+         *   CV_SUCCESS        -> continue with error test\r
+         *   CV_QSRHSFUNC_FAIL -> halt the integration \r
+         *   QSRHSFUNC_RECVR   -> predict again or stop if too many\r
+         *   \r
+         */\r
+\r
+        private int cvQuadSensNls(CVodeMemRec cv_mem)\r
+        {\r
+          int is, retval;\r
+\r
+          /* Save quadrature correction in acorQ */\r
+          retval = cvQuadSensRhsInternalDQ(cv_mem.cv_Ns, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                  cv_mem.cv_yS, cv_mem.cv_ftempQ,\r
+                                  cv_mem.cv_acorQS, cv_mem.cv_user_data,\r
+                                  cv_mem.cv_tempv, cv_mem.cv_tempvQ);\r
+          cv_mem.cv_nfQSe++;\r
+          if (retval < 0) return(CV_QSRHSFUNC_FAIL);\r
+          if (retval > 0) return(QSRHSFUNC_RECVR);\r
+\r
+\r
+          for (is=0; is<cv_mem.cv_Ns; is++) {\r
+            N_VLinearSum_Serial(cv_mem.cv_h, cv_mem.cv_acorQS[is], -ONE,\r
+                         cv_mem.cv_znQS[1][is], cv_mem.cv_acorQS[is]);\r
+            N_VScale_Serial(cv_mem.cv_rl1, cv_mem.cv_acorQS[is], cv_mem.cv_acorQS[is]);\r
+            /* Apply correction to quadrature sensitivity variables */\r
+            N_VLinearSum_Serial(ONE, cv_mem.cv_znQS[0][is], ONE,\r
+                         cv_mem.cv_acorQS[is], cv_mem.cv_yQS[is]);\r
+          }\r
+\r
+          return(CV_SUCCESS);\r
+        }\r
+\r
+        /* \r
+         * -----------------------------------------------------------------\r
+         * Functions called after a successful step\r
+         * -----------------------------------------------------------------\r
+         */\r
+\r
+        /*\r
+         * cvCompleteStep\r
+         *\r
+         * This routine performs various update operations when the solution\r
+         * to the nonlinear system has passed the local error test. \r
+         * We increment the step counter nst, record the values hu and qu,\r
+         * update the tau array, and apply the corrections to the zn array.\r
+         * The tau[i] are the last q values of h, with tau[1] the most recent.\r
+         * The counter qwait is decremented, and if qwait == 1 (and q < qmax)\r
+         * we save acor and tq[5] for a possible order increase.\r
+         */\r
+\r
+        private void cvCompleteStep(CVodeMemRec cv_mem)\r
+        {\r
+          int i, j;\r
+          int is;\r
+          \r
+          cv_mem.cv_nst++;\r
+          cv_mem.cv_nscon++;\r
+          cv_mem.cv_hu = cv_mem.cv_h;\r
+          cv_mem.cv_qu = cv_mem.cv_q;\r
+\r
+          for (i=cv_mem.cv_q; i >= 2; i--)\r
+            cv_mem.cv_tau[i] = cv_mem.cv_tau[i-1];\r
+          if ((cv_mem.cv_q==1) && (cv_mem.cv_nst > 1))\r
+            cv_mem.cv_tau[2] = cv_mem.cv_tau[1];\r
+          cv_mem.cv_tau[1] = cv_mem.cv_h;\r
+\r
+          /* Apply correction to column j of zn: l_j * Delta_n */\r
+          for (j=0; j <= cv_mem.cv_q; j++) \r
+            N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_acor, ONE,\r
+                         cv_mem.cv_zn[j], cv_mem.cv_zn[j]);\r
+\r
+          if (cv_mem.cv_quadr) {\r
+            for (j=0; j <= cv_mem.cv_q; j++) \r
+              N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_acorQ, ONE,\r
+                           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=0; j <= cv_mem.cv_q; j++) \r
+                N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_acorS[is], ONE,\r
+                             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=0; j <= cv_mem.cv_q; j++) \r
+                N_VLinearSum_Serial(cv_mem.cv_l[j], cv_mem.cv_acorQS[is], ONE,\r
+                             cv_mem.cv_znQS[j][is], cv_mem.cv_znQS[j][is]);\r
+          }\r
+\r
+\r
+          /* If necessary, store Delta_n in zn[qmax] to be used in order increase.\r
+           * This actually will be Delta_{n-1} in the ELTE at q+1 since it happens at\r
+           * the next to last step of order q before a possible one at order q+1\r
+           */\r
+\r
+          cv_mem.cv_qwait--;\r
+          if ((cv_mem.cv_qwait == 1) && (cv_mem.cv_q != cv_mem.cv_qmax)) {\r
+            \r
+            N_VScale_Serial(ONE, cv_mem.cv_acor, cv_mem.cv_zn[cv_mem.cv_qmax]);\r
+            \r
+            if (cv_mem.cv_quadr)\r
+              N_VScale_Serial(ONE, cv_mem.cv_acorQ, cv_mem.cv_znQ[cv_mem.cv_qmax]);\r
+\r
+            if (cv_mem.cv_sensi)\r
+              for (is=0; is<cv_mem.cv_Ns; is++)\r
+                N_VScale_Serial(ONE, cv_mem.cv_acorS[is], cv_mem.cv_znS[cv_mem.cv_qmax][is]);\r
+            \r
+            if (cv_mem.cv_quadr_sensi)\r
+              for (is=0; is<cv_mem.cv_Ns; is++)\r
+                N_VScale_Serial(ONE, cv_mem.cv_acorQS[is], cv_mem.cv_znQS[cv_mem.cv_qmax][is]);\r
+\r
+            cv_mem.cv_saved_tq5 = cv_mem.cv_tq[5];\r
+            cv_mem.cv_indx_acor = cv_mem.cv_qmax;\r
+          }\r
+\r
+        }\r
+\r
+        /*\r
+         * cvPrepareNextStep\r
+         *\r
+         * This routine handles the setting of stepsize and order for the\r
+         * next step -- hprime and qprime.  Along with hprime, it sets the\r
+         * ratio eta = hprime/h.  It also updates other state variables \r
+         * related to a change of step size or order. \r
+         */\r
+\r
+        private void cvPrepareNextStep(CVodeMemRec cv_mem, double dsm)\r
+        {\r
+          /* If etamax = 1, defer step size or order changes */\r
+          if (cv_mem.cv_etamax == ONE) {\r
+            cv_mem.cv_qwait = Math.max(cv_mem.cv_qwait, 2);\r
+            cv_mem.cv_qprime = cv_mem.cv_q;\r
+            cv_mem.cv_hprime = cv_mem.cv_h;\r
+            cv_mem.cv_eta = ONE;\r
+            return;\r
+          }\r
+\r
+          /* etaq is the ratio of new to old h at the current order */  \r
+          cv_mem.cv_etaq = ONE /(Math.pow(BIAS2*dsm,ONE/cv_mem.cv_L) + ADDON);\r
+          \r
+          /* If no order change, adjust eta and acor in cvSetEta and return */\r
+          if (cv_mem.cv_qwait != 0) {\r
+            cv_mem.cv_eta = cv_mem.cv_etaq;\r
+            cv_mem.cv_qprime = cv_mem.cv_q;\r
+            cvSetEta(cv_mem);\r
+            return;\r
+          }\r
+          \r
+          /* If qwait = 0, consider an order change.   etaqm1 and etaqp1 are \r
+             the ratios of new to old h at orders q-1 and q+1, respectively.\r
+             cvChooseEta selects the largest; cvSetEta adjusts eta and acor */\r
+          cv_mem.cv_qwait = 2;\r
+          cv_mem.cv_etaqm1 = cvComputeEtaqm1(cv_mem);\r
+          cv_mem.cv_etaqp1 = cvComputeEtaqp1(cv_mem);  \r
+          cvChooseEta(cv_mem); \r
+          cvSetEta(cv_mem);\r
+        }\r
+\r
+        /*\r
+         * cvSetEta\r
+         *\r
+         * This routine adjusts the value of eta according to the various\r
+         * heuristic limits and the optional input hmax.\r
+         */\r
+\r
+        private void cvSetEta(CVodeMemRec cv_mem)\r
+        {\r
+\r
+          /* If eta below the threshhold THRESH, reject a change of step size */\r
+          if (cv_mem.cv_eta < THRESH) {\r
+            cv_mem.cv_eta = ONE;\r
+            cv_mem.cv_hprime = cv_mem.cv_h;\r
+          } else {\r
+            /* Limit eta by etamax and hmax, then set hprime */\r
+            cv_mem.cv_eta = Math.min(cv_mem.cv_eta, cv_mem.cv_etamax);\r
+            cv_mem.cv_eta /= Math.max(ONE, Math.abs(cv_mem.cv_h) *\r
+                                     cv_mem.cv_hmax_inv*cv_mem.cv_eta);\r
+            cv_mem.cv_hprime = cv_mem.cv_h * cv_mem.cv_eta;\r
+            if (cv_mem.cv_qprime < cv_mem.cv_q) cv_mem.cv_nscon = 0;\r
+          }\r
+        }\r
+        \r
+        /*\r
+         * cvComputeEtaqm1\r
+         *\r
+         * This routine computes and returns the value of etaqm1 for a\r
+         * possible decrease in order by 1.\r
+         */\r
+\r
+        private double cvComputeEtaqm1(CVodeMemRec cv_mem)\r
+        {\r
+          double ddn;\r
+          \r
+          cv_mem.cv_etaqm1 = ZERO;\r
+\r
+          if (cv_mem.cv_q > 1) {\r
+\r
+            ddn = N_VWrmsNorm_Serial(cv_mem.cv_zn[cv_mem.cv_q], cv_mem.cv_ewt);\r
+\r
+            if ( cv_mem.cv_quadr && cv_mem.cv_errconQ )\r
+              ddn = cvQuadUpdateNorm(cv_mem, ddn, cv_mem.cv_znQ[cv_mem.cv_q],\r
+                                     cv_mem.cv_ewtQ);\r
+\r
+            if ( cv_mem.cv_sensi && cv_mem.cv_errconS )\r
+              ddn = cvSensUpdateNorm(cv_mem, ddn, cv_mem.cv_znS[cv_mem.cv_q],\r
+                                     cv_mem.cv_ewtS);\r
+\r
+            if ( cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS )\r
+              ddn = cvQuadSensUpdateNorm(cv_mem, ddn, cv_mem.cv_znQS[cv_mem.cv_q],\r
+                                         cv_mem.cv_ewtQS);\r
+\r
+            ddn = ddn * cv_mem.cv_tq[1];\r
+            cv_mem.cv_etaqm1 = ONE/(Math.pow(BIAS1*ddn, ONE/cv_mem.cv_q) + ADDON);\r
+          }\r
+\r
+          return(cv_mem.cv_etaqm1);\r
+        }\r
+\r
+        /*\r
+         * cvComputeEtaqp1\r
+         *\r
+         * This routine computes and returns the value of etaqp1 for a\r
+         * possible increase in order by 1.\r
+         */\r
+\r
+        private double cvComputeEtaqp1(CVodeMemRec cv_mem)\r
+        {\r
+          double dup, cquot;\r
+          int is;\r
+          \r
+          cv_mem.cv_etaqp1 = ZERO;\r
+\r
+          if (cv_mem.cv_q != cv_mem.cv_qmax) {\r
+\r
+            if (cv_mem.cv_saved_tq5 == ZERO) return(cv_mem.cv_etaqp1);\r
+\r
+            cquot = (cv_mem.cv_tq[5] / cv_mem.cv_saved_tq5) *\r
+              Math.pow(cv_mem.cv_h/cv_mem.cv_tau[2], cv_mem.cv_L);\r
+            N_VLinearSum_Serial(-cquot, cv_mem.cv_zn[cv_mem.cv_qmax], ONE,\r
+                         cv_mem.cv_acor, cv_mem.cv_tempv);\r
+            dup = 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(-cquot, cv_mem.cv_znQ[cv_mem.cv_qmax], ONE,\r
+                           cv_mem.cv_acorQ, cv_mem.cv_tempvQ);\r
+              dup = cvQuadUpdateNorm(cv_mem, dup, cv_mem.cv_tempvQ, 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(-cquot, cv_mem.cv_znS[cv_mem.cv_qmax][is], ONE,\r
+                             cv_mem.cv_acorS[is], cv_mem.cv_tempvS[is]);\r
+              dup = cvSensUpdateNorm(cv_mem, dup, cv_mem.cv_tempvS, 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(-cquot, cv_mem.cv_znQS[cv_mem.cv_qmax][is], ONE,\r
+                             cv_mem.cv_acorQS[is], cv_mem.cv_tempvQS[is]);\r
+              dup = cvSensUpdateNorm(cv_mem, dup, cv_mem.cv_tempvQS, cv_mem.cv_ewtQS);\r
+            }\r
+\r
+            dup = dup * cv_mem.cv_tq[3];\r
+            cv_mem.cv_etaqp1 = ONE / (Math.pow(BIAS3*dup, ONE/(cv_mem.cv_L+1)) + ADDON);\r
+          }\r
+\r
+          return(cv_mem.cv_etaqp1);\r
+        }\r
+        \r
+        /*\r
+         * cvChooseEta\r
+         * Given etaqm1, etaq, etaqp1 (the values of eta for qprime =\r
+         * q - 1, q, or q + 1, respectively), this routine chooses the \r
+         * maximum eta value, sets eta to that value, and sets qprime to the\r
+         * corresponding value of q.  If there is a tie, the preference\r
+         * order is to (1) keep the same order, then (2) decrease the order,\r
+         * and finally (3) increase the order.  If the maximum eta value\r
+         * is below the threshhold THRESH, the order is kept unchanged and\r
+         * eta is set to 1.\r
+         */\r
+\r
+        private void cvChooseEta(CVodeMemRec cv_mem)\r
+        {\r
+          double etam;\r
+          int is;\r
+          \r
+          etam = Math.max(cv_mem.cv_etaqm1, Math.max(cv_mem.cv_etaq, cv_mem.cv_etaqp1));\r
+          \r
+          if (etam < THRESH) {\r
+            cv_mem.cv_eta = ONE;\r
+            cv_mem.cv_qprime = cv_mem.cv_q;\r
+            return;\r
+          }\r
+\r
+          if (etam == cv_mem.cv_etaq) {\r
+\r
+            cv_mem.cv_eta = cv_mem.cv_etaq;\r
+            cv_mem.cv_qprime = cv_mem.cv_q;\r
+\r
+          } else if (etam == cv_mem.cv_etaqm1) {\r
+\r
+            cv_mem.cv_eta = cv_mem.cv_etaqm1;\r
+            cv_mem.cv_qprime = cv_mem.cv_q - 1;\r
+\r
+          } else {\r
+\r
+            cv_mem.cv_eta = cv_mem.cv_etaqp1;\r
+            cv_mem.cv_qprime = cv_mem.cv_q + 1;\r
+\r
+            if (cv_mem.cv_lmm == CV_BDF) {\r
+\r
+              /* \r
+               * Store Delta_n in zn[qmax] to be used in order increase \r
+               *\r
+               * This happens at the last step of order q before an increase\r
+               * to order q+1, so it represents Delta_n in the ELTE at q+1\r
+               */\r
+              \r
+              N_VScale_Serial(ONE, cv_mem.cv_acor, cv_mem.cv_zn[cv_mem.cv_qmax]);\r
+              \r
+              if (cv_mem.cv_quadr && cv_mem.cv_errconQ)\r
+                N_VScale_Serial(ONE, cv_mem.cv_acorQ, cv_mem.cv_znQ[cv_mem.cv_qmax]);\r
+              \r
+              if (cv_mem.cv_sensi && cv_mem.cv_errconS)\r
+                for (is=0; is<cv_mem.cv_Ns; is++)\r
+                  N_VScale_Serial(ONE, cv_mem.cv_acorS[is], cv_mem.cv_znS[cv_mem.cv_qmax][is]);\r
+\r
+              if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS)\r
+                for (is=0; is<cv_mem.cv_Ns; is++)\r
+                  N_VScale_Serial(ONE, cv_mem.cv_acorQS[is], cv_mem.cv_znQS[cv_mem.cv_qmax][is]);\r
+\r
+            }\r
+          }\r
+        }\r
+\r
+        /*\r
+         * cvBDFStab\r
+         *\r
+         * This routine handles the BDF Stability Limit Detection Algorithm\r
+         * STALD.  It is called if lmm = CV_BDF and the SLDET option is on.\r
+         * If the order is 3 or more, the required norm data is saved.\r
+         * If a decision to reduce order has not already been made, and\r
+         * enough data has been saved, cvSLdet is called.  If it signals\r
+         * a stability limit violation, the order is reduced, and the step\r
+         * size is reset accordingly.\r
+         */\r
+\r
+        private void cvBDFStab(CVodeMemRec cv_mem)\r
+        {\r
+          int i,k, ldflag, factorial;\r
+          double sq, sqm1, sqm2;\r
+              \r
+          /* If order is 3 or greater, then save scaled derivative data,\r
+             push old data down in i, then add current values to top.    */\r
+\r
+          if (cv_mem.cv_q >= 3) {\r
+            for (k = 1; k <= 3; k++)\r
+              for (i = 5; i >= 2; i--)\r
+                cv_mem.cv_ssdat[i][k] = cv_mem.cv_ssdat[i-1][k];\r
+            factorial = 1;\r
+            for (i = 1; i <= cv_mem.cv_q-1; i++) factorial *= i;\r
+            sq = factorial * cv_mem.cv_q * (cv_mem.cv_q+1) *\r
+              cv_mem.cv_acnrm / Math.max(cv_mem.cv_tq[5],TINY);\r
+            sqm1 = factorial * cv_mem.cv_q *\r
+              N_VWrmsNorm_Serial(cv_mem.cv_zn[cv_mem.cv_q], cv_mem.cv_ewt);\r
+            sqm2 = factorial *\r
+              N_VWrmsNorm_Serial(cv_mem.cv_zn[cv_mem.cv_q-1], cv_mem.cv_ewt);\r
+            cv_mem.cv_ssdat[1][1] = sqm2*sqm2;\r
+            cv_mem.cv_ssdat[1][2] = sqm1*sqm1;\r
+            cv_mem.cv_ssdat[1][3] = sq*sq;\r
+          }  \r
+\r
+          if (cv_mem.cv_qprime >= cv_mem.cv_q) {\r
+\r
+            /* If order is 3 or greater, and enough ssdat has been saved,\r
+               nscon >= q+5, then call stability limit detection routine.  */\r
+\r
+            if ( (cv_mem.cv_q >= 3) && (cv_mem.cv_nscon >= cv_mem.cv_q+5) ) {\r
+              ldflag = cvSLdet(cv_mem);\r
+              if (ldflag > 3) {\r
+                /* A stability limit violation is indicated by\r
+                   a return flag of 4, 5, or 6.\r
+                   Reduce new order.                     */\r
+                cv_mem.cv_qprime = cv_mem.cv_q-1;\r
+                cv_mem.cv_eta = cv_mem.cv_etaqm1; \r
+                cv_mem.cv_eta = Math.min(cv_mem.cv_eta,cv_mem.cv_etamax);\r
+                cv_mem.cv_eta = cv_mem.cv_eta /\r
+                  Math.max(ONE,Math.abs(cv_mem.cv_h)*cv_mem.cv_hmax_inv*cv_mem.cv_eta);\r
+                cv_mem.cv_hprime = cv_mem.cv_h * cv_mem.cv_eta;\r
+                cv_mem.cv_nor = cv_mem.cv_nor + 1;\r
+              }\r
+            }\r
+          }\r
+          else {\r
+            /* Otherwise, let order increase happen, and \r
+               reset stability limit counter, nscon.     */\r
+            cv_mem.cv_nscon = 0;\r
+          }\r
+        }\r
+        \r
+        /*\r
+         * cvSLdet\r
+         *\r
+         * This routine detects stability limitation using stored scaled \r
+         * derivatives data. cvSLdet returns the magnitude of the\r
+         * dominate characteristic root, rr. The presents of a stability\r
+         * limit is indicated by rr > "something a little less then 1.0",  \r
+         * and a positive kflag. This routine should only be called if\r
+         * order is greater than or equal to 3, and data has been collected\r
+         * for 5 time steps. \r
+         * \r
+         * Returned values:\r
+         *    kflag = 1 -> Found stable characteristic root, normal matrix case\r
+         *    kflag = 2 -> Found stable characteristic root, quartic solution\r
+         *    kflag = 3 -> Found stable characteristic root, quartic solution,\r
+         *                 with Newton correction\r
+         *    kflag = 4 -> Found stability violation, normal matrix case\r
+         *    kflag = 5 -> Found stability violation, quartic solution\r
+         *    kflag = 6 -> Found stability violation, quartic solution,\r
+         *                 with Newton correction\r
+         *\r
+         *    kflag < 0 -> No stability limitation, \r
+         *                 or could not compute limitation.\r
+         *\r
+         *    kflag = -1 -> Min/max ratio of ssdat too small.\r
+         *    kflag = -2 -> For normal matrix case, vmax > vrrt2*vrrt2\r
+         *    kflag = -3 -> For normal matrix case, The three ratios\r
+         *                  are inconsistent.\r
+         *    kflag = -4 -> Small coefficient prevents elimination of quartics.  \r
+         *    kflag = -5 -> R value from quartics not consistent.\r
+         *    kflag = -6 -> No corrected root passes test on qk values\r
+         *    kflag = -7 -> Trouble solving for sigsq.\r
+         *    kflag = -8 -> Trouble solving for B, or R via B.\r
+         *    kflag = -9 -> R via sigsq[k] disagrees with R from data.\r
+         */\r
+\r
+        private int cvSLdet(CVodeMemRec cv_mem)\r
+        {\r
+          int i, k, j, it, kmin = 0, kflag = 0;\r
+          double rat[][] = new double[5][4];\r
+          double rav[] = new double[4];\r
+          double qkr[] = new double[4];\r
+          double sigsq[] = new double[4];\r
+          double smax[] = new double[4];\r
+          double ssmax[] = new double[4];\r
+          double drr[] = new double[4];\r
+          double rrc[] = new double[4];\r
+          double sqmx[] = new double[4];\r
+          double qjk[][] = new double[4][4];\r
+          double vrat[] = new double[5];\r
+          double qc[][] = new double[6][4];\r
+          double qco[][] = new double[6][4];\r
+          double rr, rrcut, vrrtol, vrrt2, sqtol, rrtol;\r
+          double smink, smaxk, sumrat, sumrsq, vmin, vmax, drrmax, adrr;\r
+          double tem, sqmax, saqk, qp, s, sqmaxk, saqj;\r
+          double sqmin = 0.0;\r
+          double rsa, rsb, rsc, rsd, rd1a, rd1b, rd1c;\r
+          double rd2a, rd2b, rd3a, cest1, corr1; \r
+          double ratp, ratm, qfac1, qfac2, bb, rrb;\r
+\r
+          /* The following are cutoffs and tolerances used by this routine */\r
+\r
+          rrcut  = 0.98;\r
+          vrrtol = 1.0e-4;\r
+          vrrt2  = 5.0e-4;\r
+          sqtol  = 1.0e-3;\r
+          rrtol  = 1.0e-2;\r
+\r
+          rr = ZERO;\r
+\r
+          /*  Index k corresponds to the degree of the interpolating polynomial. */\r
+          /*      k = 1 -> q-1          */\r
+          /*      k = 2 -> q            */\r
+          /*      k = 3 -> q+1          */\r
+\r
+          /*  Index i is a backward-in-time index, i = 1 -> current time, */\r
+          /*      i = 2 -> previous step, etc    */\r
+\r
+          /* get maxima, minima, and variances, and form quartic coefficients  */\r
+\r
+          for (k=1; k<=3; k++) {\r
+            smink = cv_mem.cv_ssdat[1][k];\r
+            smaxk = ZERO;\r
+            \r
+            for (i=1; i<=5; i++) {\r
+              smink = Math.min(smink,cv_mem.cv_ssdat[i][k]);\r
+              smaxk = Math.max(smaxk,cv_mem.cv_ssdat[i][k]);\r
+            }\r
+            \r
+            if (smink < TINY*smaxk) {\r
+              kflag = -1;  \r
+              return(kflag);\r
+            }\r
+            smax[k] = smaxk;\r
+            ssmax[k] = smaxk*smaxk;\r
+\r
+            sumrat = ZERO;\r
+            sumrsq = ZERO;\r
+            for (i=1; i<=4; i++) {\r
+              rat[i][k] = cv_mem.cv_ssdat[i][k] / cv_mem.cv_ssdat[i+1][k];\r
+              sumrat = sumrat + rat[i][k];\r
+              sumrsq = sumrsq + rat[i][k]*rat[i][k];\r
+            } \r
+            rav[k] = FOURTH*sumrat;\r
+            vrat[k] = Math.abs(FOURTH*sumrsq - rav[k]*rav[k]);\r
+\r
+            qc[5][k] = cv_mem.cv_ssdat[1][k] * cv_mem.cv_ssdat[3][k] -\r
+              cv_mem.cv_ssdat[2][k] * cv_mem.cv_ssdat[2][k];\r
+            qc[4][k] = cv_mem.cv_ssdat[2][k] * cv_mem.cv_ssdat[3][k] -\r
+              cv_mem.cv_ssdat[1][k] * cv_mem.cv_ssdat[4][k];\r
+            qc[3][k] = ZERO;\r
+            qc[2][k] = cv_mem.cv_ssdat[2][k] * cv_mem.cv_ssdat[5][k] -\r
+              cv_mem.cv_ssdat[3][k] * cv_mem.cv_ssdat[4][k];\r
+            qc[1][k] = cv_mem.cv_ssdat[4][k] * cv_mem.cv_ssdat[4][k] -\r
+              cv_mem.cv_ssdat[3][k] * cv_mem.cv_ssdat[5][k];\r
+\r
+            for (i=1; i<=5; i++) {\r
+              qco[i][k] = qc[i][k];\r
+            }\r
+          }                            /* End of k loop */\r
+          \r
+          /* Isolate normal or nearly-normal matrix case. Three quartic will\r
+             have common or nearly-common roots in this case. \r
+             Return a kflag = 1 if this procedure works. If three root \r
+             differ more than vrrt2, return error kflag = -3.    */\r
+          \r
+          vmin = Math.min(vrat[1],Math.min(vrat[2],vrat[3]));\r
+          vmax = Math.max(vrat[1],Math.max(vrat[2],vrat[3]));\r
+          \r
+          if (vmin < vrrtol*vrrtol) {\r
+\r
+            if (vmax > vrrt2*vrrt2) {\r
+              kflag = -2;  \r
+              return(kflag);\r
+            } else {\r
+              rr = (rav[1] + rav[2] + rav[3])/THREE;\r
+              drrmax = ZERO;\r
+              for (k = 1;k<=3;k++) {\r
+                adrr = Math.abs(rav[k] - rr);\r
+                drrmax = Math.max(drrmax, adrr);\r
+              }\r
+              if (drrmax > vrrt2) kflag = -3;    \r
+              kflag = 1;\r
+\r
+              /*  can compute charactistic root, drop to next section   */\r
+            }\r
+\r
+          } else {\r
+\r
+            /* use the quartics to get rr. */\r
+\r
+            if (Math.abs(qco[1][1]) < TINY*ssmax[1]) {\r
+              kflag = -4;    \r
+              return(kflag);\r
+            }\r
+\r
+            tem = qco[1][2]/qco[1][1];\r
+            for (i=2; i<=5; i++) {\r
+              qco[i][2] = qco[i][2] - tem*qco[i][1];\r
+            }\r
+\r
+            qco[1][2] = ZERO;\r
+            tem = qco[1][3]/qco[1][1];\r
+            for (i=2; i<=5; i++) {\r
+              qco[i][3] = qco[i][3] - tem*qco[i][1];\r
+            }\r
+            qco[1][3] = ZERO;\r
+\r
+            if (Math.abs(qco[2][2]) < TINY*ssmax[2]) {\r
+              kflag = -4;    \r
+              return(kflag);\r
+            }\r
+\r
+            tem = qco[2][3]/qco[2][2];\r
+            for (i=3; i<=5; i++) {\r
+              qco[i][3] = qco[i][3] - tem*qco[i][2];\r
+            }\r
+\r
+            if (Math.abs(qco[4][3]) < TINY*ssmax[3]) {\r
+              kflag = -4;    \r
+              return(kflag);\r
+            }\r
+\r
+            rr = -qco[5][3]/qco[4][3];\r
+\r
+            if (rr < TINY || rr > HUNDRED) {\r
+              kflag = -5;   \r
+              return(kflag);\r
+            }\r
+\r
+            for (k=1; k<=3; k++)\r
+              qkr[k] = qc[5][k] + rr*(qc[4][k] + rr*rr*(qc[2][k] + rr*qc[1][k]));\r
+\r
+            sqmax = ZERO;\r
+            for (k=1; k<=3; k++) {\r
+              saqk = Math.abs(qkr[k])/ssmax[k];\r
+              if (saqk > sqmax) sqmax = saqk;\r
+            } \r
+\r
+            if (sqmax < sqtol) {\r
+              kflag = 2;\r
+\r
+              /*  can compute charactistic root, drop to "given rr,etc"   */\r
+\r
+            } else {\r
+\r
+              /* do Newton corrections to improve rr.  */\r
+\r
+              for (it=1; it<=3; it++) {\r
+                for (k=1; k<=3; k++) {\r
+                  qp = qc[4][k] + rr*rr*(THREE*qc[2][k] + rr*FOUR*qc[1][k]);\r
+                  drr[k] = ZERO;\r
+                  if (Math.abs(qp) > TINY*ssmax[k]) drr[k] = -qkr[k]/qp;\r
+                  rrc[k] = rr + drr[k];\r
+                } \r
+\r
+                for (k=1; k<=3; k++) {\r
+                  s = rrc[k];\r
+                  sqmaxk = ZERO;\r
+                  for (j=1; j<=3; j++) {\r
+                    qjk[j][k] = qc[5][j] + s*(qc[4][j] + s*s*(qc[2][j] + s*qc[1][j]));\r
+                    saqj = Math.abs(qjk[j][k])/ssmax[j];\r
+                    if (saqj > sqmaxk) sqmaxk = saqj;\r
+                  } \r
+                  sqmx[k] = sqmaxk;\r
+                }\r
+\r
+                sqmin = sqmx[1] + ONE;\r
+                for (k=1; k<=3; k++) {\r
+                  if (sqmx[k] < sqmin) {\r
+                    kmin = k;\r
+                    sqmin = sqmx[k];\r
+                  }\r
+                } \r
+                rr = rrc[kmin];\r
+\r
+                if (sqmin < sqtol) {\r
+                  kflag = 3;\r
+                  /*  can compute charactistic root   */\r
+                  /*  break out of Newton correction loop and drop to "given rr,etc" */ \r
+                  break;\r
+                } else {\r
+                  for (j=1; j<=3; j++) {\r
+                    qkr[j] = qjk[j][kmin];\r
+                  }\r
+                }     \r
+              } /*  end of Newton correction loop  */ \r
+\r
+              if (sqmin > sqtol) {\r
+                kflag = -6;\r
+                return(kflag);\r
+              }\r
+            } /*  end of if (sqmax < sqtol) else   */\r
+          } /*  end of if (vmin < vrrtol*vrrtol) else, quartics to get rr. */\r
+\r
+          /* given rr, find sigsq[k] and verify rr.  */\r
+          /* All positive kflag drop to this section  */\r
+          \r
+          for (k=1; k<=3; k++) {\r
+            rsa = cv_mem.cv_ssdat[1][k];\r
+            rsb = cv_mem.cv_ssdat[2][k]*rr;\r
+            rsc = cv_mem.cv_ssdat[3][k]*rr*rr;\r
+            rsd = cv_mem.cv_ssdat[4][k]*rr*rr*rr;\r
+            rd1a = rsa - rsb;\r
+            rd1b = rsb - rsc;\r
+            rd1c = rsc - rsd;\r
+            rd2a = rd1a - rd1b;\r
+            rd2b = rd1b - rd1c;\r
+            rd3a = rd2a - rd2b;\r
+            \r
+            if (Math.abs(rd1b) < TINY*smax[k]) {\r
+              kflag = -7;\r
+              return(kflag);\r
+            }\r
+            \r
+            cest1 = -rd3a/rd1b;\r
+            if (cest1 < TINY || cest1 > FOUR) {\r
+              kflag = -7;\r
+              return(kflag);\r
+            }\r
+            corr1 = (rd2b/cest1)/(rr*rr);\r
+            sigsq[k] = cv_mem.cv_ssdat[3][k] + corr1;\r
+          }\r
+          \r
+          if (sigsq[2] < TINY) {\r
+            kflag = -8;\r
+            return(kflag);\r
+          }\r
+          \r
+          ratp = sigsq[3]/sigsq[2];\r
+          ratm = sigsq[1]/sigsq[2];\r
+          qfac1 = FOURTH*(cv_mem.cv_q*cv_mem.cv_q - ONE);\r
+          qfac2 = TWO/(cv_mem.cv_q - ONE);\r
+          bb = ratp*ratm - ONE - qfac1*ratp;\r
+          tem = ONE - qfac2*bb;\r
+          \r
+          if (Math.abs(tem) < TINY) {\r
+            kflag = -8;\r
+            return(kflag);\r
+          }\r
+          \r
+          rrb = ONE/tem;\r
+          \r
+          if (Math.abs(rrb - rr) > rrtol) {\r
+            kflag = -9;\r
+            return(kflag);\r
+          }\r
+          \r
+          /* Check to see if rr is above cutoff rrcut  */\r
+          if (rr > rrcut) {\r
+            if (kflag == 1) kflag = 4;\r
+            if (kflag == 2) kflag = 5;\r
+            if (kflag == 3) kflag = 6;\r
+          }\r
+          \r
+          /* All positive kflag returned at this point  */\r
+          \r
+          return(kflag);\r
+          \r
+        }\r
+        \r
+        /* \r
+         * CVodeGetRootInfo\r
+         *\r
+         * Returns pointer to array rootsfound showing roots found\r
+         */\r
+\r
+        private int CVodeGetRootInfo(CVodeMemRec cv_mem, int rootsfound[])\r
+        {\r
+          int i, nrt;\r
+\r
+          if (cv_mem==null) {\r
+            cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetRootInfo", MSGCV_NO_MEM);\r
+            return(CV_MEM_NULL);\r
+          }\r
+\r
+          nrt = cv_mem.cv_nrtfn;\r
+\r
+          for (i=0; i<nrt; i++) rootsfound[i] = cv_mem.cv_iroots[i];\r
+\r
+          return(CV_SUCCESS);\r
+        }\r
+\r
 \r
 \r
 }
\ No newline at end of file